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 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.backproject   = 'yes' or 'no',  determines when reducerank is applied whether the
%                       lower rank leadfield is projected back onto the original linear
%                       subspace, or not (default = 'yes')
%
% 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

% 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, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729
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', 'all');   % for comp input
cfg.frequency       = ft_getopt(cfg, 'frequency');          % for freq input
cfg.latency         = ft_getopt(cfg, 'latency', 'all');     % for timelock 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

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'});

% determine data type
iscomp = ft_datatype(data, 'comp');           % it can also be raw+comp, timelock+comp or freq+comp
isfreq = ft_datatype(data, 'freq');           % it might also be freq+comp, in that case it should be treated as component data
istimelock = ft_datatype(data, 'timelock');   % it might also be timelock+comp, in that case it should be treated as component data

% the default for this depends on the data type
if ~isfield(cfg, 'model')
  if iscomp
    % each component is fitted independently
    cfg.model = 'moving';
  elseif isfreq
    % fit the data with a dipole at one location
    cfg.model = 'regional';
  elseif istimelock
    % 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_INVERSE_DIPOLEFIT 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_INVERSE_DIPOLEFIT 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_INVERSE_DIPOLEFIT 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 iscomp
  % transform the data into a representation on which the timelocked dipole fit can perform its trick
  data = comp2timelock(cfg, data);

  % default component selection is all components
  if ischar(cfg.component) && strcmp(cfg.component, 'all')
    cfg.component = (1:size(data.avg, 2));
  end
elseif isfreq
  % transform the data into a representation on which the timelocked dipole fit can perform its trick
  data = freq2timelock(cfg, data);
elseif istimelock
  % no transformation is needed
end

% collect and preprocess the electrodes/gradiometer and head model
% this will also update cfg.channel to match the electrodes/gradiometers
[headmodel, sens, cfg] = prepare_headmodel(cfg, data);

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD and FT_INVERSE_DIPOLEFIT
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(cfg, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(cfg, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(cfg, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(cfg, 'weight'));

% construct the low-level options for the dipole fitting as key-value pairs, these are passed to FT_INVERSE_DIPOLEFIT
dipfitopt = ft_cfg2keyval(cfg.dipfit);

% select the desired channels, ordered according to the sensor structure or configuration
[selcfg, seldata] = match_str(cfg.channel, 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(:, :);
elseif istimelock
  % 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);
ft_info('selected %d channels\n', nchans);
ft_info('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

  if isfield(cfg.sourcemodel, 'leadfield')
    ft_notice('using precomputed leadfields for the gridsearch');

    sourcemodel = keepfields(cfg.sourcemodel, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label'});

    % select the channels corresponding to the data and the user configuration
    tmpcfg = keepfields(cfg, 'channel');
    sourcemodel = ft_selectdata(tmpcfg, sourcemodel);

    % sort the channels to be consistent with the data
    [dum, chansel] = match_str(data.label, sourcemodel.label);
    sourcemodel.label = sourcemodel.label(chansel);
    for i=1:numel(sourcemodel.leadfield)
      if ~isempty(sourcemodel.leadfield{i})
        sourcemodel.leadfield{i} = sourcemodel.leadfield{i}(chansel, :);
      end
    end

    % ensure that the channels are consistent with the data
    assert(isequal(sourcemodel.label, cfg.channel), 'cannot match the channels in the sourcemodel to those in the data')

  else
    ft_notice('computing the leadfields for the gridsearch on the fly');

    % construct the dipole positions on which the source reconstruction will be done
    tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid' 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
    tmpcfg.headmodel = headmodel;
    if ft_senstype(sens, 'eeg')
      tmpcfg.elec = sens;
    elseif ft_senstype(sens, 'meg')
      tmpcfg.grad = sens;
    end
    sourcemodel = ft_prepare_sourcemodel(tmpcfg);

  end % if precomputed leadfield or not

  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
        ft_info('found minimum after scanning on grid point [%g %g %g]\n', dip.pos(1), dip.pos(2), dip.pos(3));
      elseif cfg.numdipoles==2
        ft_info('found minimum after scanning on grid point [%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

    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
          ft_info('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
          ft_info('found minimum after scanning for topography %d on grid point [%g %g %g; %g %g %g]\n', t, dip(t).pos(1,1), dip(t).pos(1,2), dip(t).pos(1,3), dip(t).pos(2,1), dip(t).pos(2,2), dip(t).pos(2,3));
        end
      end

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

elseif strcmp(cfg.gridsearch, 'no')
  % there is no grid needed for dipole scanning
  sourcemodel = [];
  % use the initial guess supplied in the configuration for the remainder
  switch cfg.model
    case 'regional'
      dip = cfg.dip;
    case 'moving'
      for t=1:ntime
        dip(t) = cfg.dip;
      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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 = ft_inverse_dipolefit(dip, sens, headmodel, Vdata, dipfitopt{:}, leadfieldopt{:});
        success = 1;
        if cfg.numdipoles==1
          ft_info('found minimum after non-linear optimization on [%g %g %g]\n', dip.pos(1), dip.pos(2), dip.pos(3));
        elseif cfg.numdipoles==2
          ft_info('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) = ft_inverse_dipolefit(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) = ft_inverse_dipolefit(dipin(t), sens, headmodel, Vdata(:,t), dipfitopt{:}, leadfieldopt{:});
          success(t) = 1;
          if cfg.numdipoles==1
            ft_info('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
            ft_info('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)
  if isfield(headmodel, 'unit')
    source.dip(i).unit = headmodel.unit;
  elseif isfield(sourcemodel, 'unit')
    source.dip(i).unit = sourcemodel.unit;
  end
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
elseif istimelock
  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 previous   data
ft_postamble provenance source
ft_postamble history    source
ft_postamble savevar    source


text:
function sim_2prob(value2trans) % autogenerated function wrapper
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun(mfilename('fullpath'));

    % needed variables
    % BigCatalog        big catalog from which to take the eqs randomly, produced by translating
    %                   consists of 100000 eqs
    % ni                number of earthquakes in a bin, i.e. sample size
    % NuBins            number of bins
    % BinLength         1/length(xt), length of shortest possible interval
    % winlen_days               length of interval in times shortest
    % NuRep             number of repetitions

    delta=winlen_days/NuBins;

    for nto=1:NuRep
        disp(nto);

        which=ceil(100000*(rand(ni)));
        for i=1:ni
            rancata(i)=BigCatalog(which(i));
        end
        clear i which;
        rancata=ceil(rancata*NuBins);

        for i=1:NuBins
            l=sum(rancata==i); Bins(i,1)=sum(l); clear l;
        end
        clear rancata i;

        FirstBin=ceil(rand(1)*(NuBins-winlen_days+1));

        zin=Bins(FirstBin:FirstBin+winlen_days-1); zout=[Bins(1:FirstBin-1,1); Bins(FirstBin+winlen_days:NuBins,1)];
        ToBeFitted(nto,1)=nto;
        % calculating beta
        ToBeFitted(nto,2)=(sum(zin)-ni*delta)/(sqrt(ni*delta*(1-delta)));
        % calculating z
        ToBeFitted(nto,3)=(mean(zout)-mean(zin))/(sqrt(var(zin)/sum(zin)+var(zout)/sum(zout)));
        clear Bins FirstBin zin zout;
    end
    clear BigCatalog nto delta;

    [meanval, std] =normfit(ToBeFitted(:,2)); IsFitted(1,1)=meanval; IsFitted(1,2)=std;
    [meanval, std] =normfit(ToBeFitted(:,3)); IsFitted(2,1)=meanval; IsFitted(2,2)=std;
    clear meanval std;
    clear ToBeFitted;

    switch value2trans
        case 'beta'
            Pbeta = normcdf(BetaValues,IsFitted(1,1),IsFitted(1,2));
            l = Pbeta == 0; Pbeta(l) = nan;
        case 'z'
            Pbeta = normcdf(BetaValues,IsFitted(2,1),IsFitted(2,2));
            l = Pbeta == 0; Pbeta(l) = nan;
    end

    % plot the resuts
    figure
    pq = -log10(1-Pbeta); l = isinf(pq);pq(l) = 18 ;
    pl1 = plot(xt,pq,'color',[0.0 0.5 0.9]);
    hold on
    l = pq < 1.3; pq(l) = nan;
    pl3 = plot(xt,pq,'b','Linewidth',2);

    pq = -log10(Pbeta);l = isinf(pq);pq(l) = 18 ;
    pl2 = plot(xt,pq,'color',[0.8 0.6 0.8]);
    l = pq < 1.3; pq(l) = nan;
    pl4 = plot(xt,pq,'r','Linewidth',2);

    maxd = [get(pl1,'Ydata') get(pl2,'ydata') ]; maxd(isinf(maxd)) = []; maxd = max(maxd);
    if maxd < 5 ; maxd = 5; end
    if isnan(maxd) == 1 ; maxd = 10; end

    legend([pl3 pl4],'Rate increases','Rate decreases');
    set(gca,'Ylim',[0 maxd+1])
    set(gca,'YTick',[1.3 2 3 4 5])
    set(gca,'YTickLabel',[ '    5%' ; '    1%' ;  '  0.1%' ;  ' 0.01%' ; '0.001%'])
    set(gca,'TickDir','out','Ticklength',[0.02 0.02],'pos',[0.2 0.2 0.7 0.7]);
    xlabel('Time [years]')
    ylabel('Significance level');
    set(gcf,'color','w')
    grid

    uicontrol('Units','normal',...
        'Position',[.8 .0 .1 .05],'String','Explain ... ',...
        'callback',@(~,~)showweb('explproba'));

    delete(probut)

end

author: CelsoReyes


text:
function [ C ] = gsp_ngwft(G,f,g, param)
%GSP_NGWFT Normalized graph windowed Fourier transform
%   Usage:  G = gsp_ngwft(G,f,g, param);
%           G = gsp_ngwft(G,f,g);
%
%   Input parameters:
%         G     : Graph
%         f     : Graph signal
%         g     : window
%         param : Structure of optional parameter
%   Output parameters:
%         C     : Coefficient
%
%   This function compute the normalized graph windowed Fourier transform
%   of a signal *f* with the window *g*. The function returns a matrix of
%   size N^2*N.
%
%   *param* a Matlab structure containing the following fields:
%
%     * *param.verbose* : 0 no log, 1 print main steps, 2 print all steps.
%       By default, it is 1.
%     * *param.lowmemory* : use less memory. By default, it is 1.
%

% Author : Nathanael Perraudin
% Testing: test_ngwft

% Optional input arguments
if nargin<4, param=struct; end
if ~isfield(param, 'verbose'), param.verbose=1 ; end
if ~isfield(param, 'lowmemory'), param.lowmemory=1 ; end

if ~isfield(G,'U')
   error(['GSP_NGWFT: You need first to compute the Fourier basis\n',...
       'You can do it with the function gsp_compute_fourier_basis']);
end

if ~param.lowmemory
    % Compute the Frame into a big matrix
    Frame=gsp_ngwft_frame_matrix(G,g,param);

    C=Frame'*f;

    C=reshape(C,G.N,G.N);
else
    % Compute the translate of g
    ghat=G.U'*g;
    Ftrans=sqrt(G.N)*G.U*(repmat(ghat,1,G.N).*G.U');
    C=zeros(G.N);
    for ii=1:G.N
         atoms = (repmat(1./G.U(:,1),1,G.N) .* ...
             G.U.*repmat(Ftrans(:,ii),1,G.N))';
         % normalization
         atoms = atoms./ repmat(sqrt(sum(abs(atoms).^2,2)),1,G.N);
         C(:,ii)=atoms*f;
    end

end

repo: gspbox


text:
function [inputImages, inputSignals, choices] = signalSorter(inputImages,inputSignals,varargin)
	% Displays a GUI for sorting images (e.g. cells) and their associated signals (e.g. fluorescence activity traces). Also does preliminary sorting based on image/signal properties if requested by user.
	% See following URL for details of GUI and tips on manual sorting: https://github.com/bahanonu/calciumImagingAnalysis/wiki/Manual-cell-sorting-of-cell-extraction-outputs.
	% Biafra Ahanonu
	% started: 2013.10.08
	% Dependent code
		% getOptions.m, createObjMap.m, removeSmallICs.m, identifySpikes.m, etc., see repository
	% inputs
		% inputImages - [x y N] matrix where N = number of images, x/y are dimensions. Use permute(inputImages,[2 3 1]) if you use [N x y] for matrix indexing.
			% Alternatively, make inputImages = give path to NWB file and inputSignals = [] for signalSorter to automatically load NWB files.
		% inputSignals - [N time] matrix where N = number of signals (traces) and time = frames.
		% inputID - obsolete, kept for compatibility, just input empty []
		% nSignals - obsolete, kept for compatibility
	% outputs
		% inputImages - [N x y] matrix where N = number of images, x/y are dimensions with only manual choices kept.
		% inputSignals
		% choices

	[inputImages, inputSignals, choices] = ciapkg.classification.signalSorter(inputImages,inputSignals,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function filename = savecplexlp(varargin)
%SAVCPLEXLP Saves a problem definition in CPLEX-LP format
%
%    SAVCPLEXLP(F,h,'filename')    Saves the problem min(h(x)), F(x)>0 to the file filename
%    SAVCPLEXLP(F,h)               A "Save As"- box will be opened
%

F = varargin{1};
h = varargin{2};

[aux1,aux2,aux3,model] = export(F,h);

% Check so that it really is an LP
% if any(any(model.Q)) | any(model.variabletype)
% Only check the variabletype
if any(model.variabletype)
    error('This is not an LP or QP');
end

c = model.c;
Q = 2*full(model.Q);
b = model.F_struc(:,1);
A = -model.F_struc(:,2:end);
if model.K.f>0
    Aeq = A(1:model.K.f,:);
    beq = b(1:model.K.f,:);
    A(1:model.K.f,:) = [];
    b(1:model.K.f) = [];
else
    Aeq = [];
    beq = [];
end

lb = model.lb;
ub = model.ub;

[lb,ub,A,b] = remove_bounds_from_Ab(A,b,lb,ub);
[lb,ub,Aeq,beq] = remove_bounds_from_Aeqbeq(Aeq,beq,lb,ub);

% Is a filename supplied
if nargin<3
    [filename, pathname] = uiputfile('*.lp', 'Save LP format file');
    if isa(filename,'double')
        return % User cancelled
    else
        % Did the user change the extension
        if isempty(strfind(filename,'.'))
            filename = [pathname filename '.lp'];
        else
            filename = [pathname filename];
        end
    end
else
    filename = varargin{3};
end

fid = fopen(filename,'w');

obj = strrep(lptext(c(:)'),'+ -','-');

fprintf(fid,['Minimize\r\n obj:  ' obj(1:end-2) '']);
if any(any(Q))
    obj = qptext(Q);
    fprintf(fid,[' + [' obj(1:end-2) '] / 2']);
end
fprintf(fid,'\r\n');

fprintf(fid,['\r\n']);
fprintf(fid,['Subject To\r\n']);

for i = 1:length(b)
    rowtext = lptext(-A(i,:));
    rhs = sprintf('%0.20g',full(-b(i)));
    rowtext = [rowtext(1:end-2) ' >= ' rhs];
    fprintf(fid,[' c%i: ' strrep(rowtext,'+ -','-') ''],i);
    fprintf(fid,'\r\n');
end
for i = 1:length(beq)
    rowtext = lptext(-Aeq(i,:));
    rowtext = [rowtext(1:end-2) '== ' sprintf('%0.20g',full(-beq(i)))];
    fprintf(fid,[' eq%i: ' strrep(rowtext,'+ -','-') ''],i);
    fprintf(fid,'\r\n');
end

if length(c)>length(model.binary_variables)
    fprintf(fid,['\r\nBounds\r\n']);
    for i = 1:length(c)
        %        if ~ismember(i,model.binary_variables)
        if isinf(lb(i)) & isinf(ub(i))
            fprintf(fid,[' x%i free\n\r'],i);
        elseif lb(i)==0 & isinf(ub(i))
            % Standard non-negative variable
        elseif isinf(ub(i))
            s = strrep(sprintf(['%0.20g <= x%i \r\n'],[lb(i) i ]),'Inf','inf');
            fprintf(fid,s);
        else
            s = strrep(sprintf(['%0.20g <= x%i <= %0.20g \r\n'],[lb(i) i ub(i)]),'Inf','inf');
            fprintf(fid,s);
        end
        %        end
    end
end

if length(model.binary_variables)>0
    fprintf(fid,['\r\n']);
    fprintf(fid,['Binary\r\n']);
    for i = 1:length(model.binary_variables)
        fprintf(fid,[' x%i\r\n'],model.binary_variables(i));
    end
end

if length(model.integer_variables)>0
    fprintf(fid,['\r\n']);
    fprintf(fid,['Integer\r\n']);
    for i = 1:length(model.integer_variables)
        fprintf(fid,[' x%i\r\n'],model.integer_variables(i));
    end
end

fprintf(fid,['\r\nEnd']);
fclose(fid);

function rowtext = lptext(a)
[aux,poss,vals] = find(a);
rowtext = sprintf('%0.20g x%d + ',reshape([vals(:) poss(:)]',[],1));
%rowtext = strrep(rowtext,'+ -','- ');
%rowtext(isspace(rowtext))=[];
%rowtext = strrep(rowtext,'+-','-');
%rowtext = strrep(rowtext,'-1x','-x');
%rowtext = strrep(rowtext,'+1x','+x');

function rowtext = qptext(Q)
n=size(Q,2);
q = diag(Q);
if any(q)
    i = find(q);
    rowtext = sprintf('%0.20g x%d ^2 + ',reshape([q(i) i]',[],1));
end
for i=1:n
    for j=i+1:n
        if ~(Q(i,j)+Q(j,i)==0)
        rowtext = [rowtext sprintf('%0.20g x%d * x%d + ',Q(i,j)+Q(j,i),i,j)];
        end
    end
end
rowtext = strrep(rowtext,'+ -','- ');

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
disp('Starting tests...')
test_digamma
test_gammaln
test_trigamma
test_int_hist
%test_inv_posdef
test_java
test_logmulexp
test_logsumexp
test_mutable
%test_ndsum
%test_normcdf
test_normpdf
test_randgamma
test_randbeta
test_randbinom
test_randwishart
test_repmat
test_row_sum
test_sameobject
test_scale
test_solve_tri
test_sorted
test_sparse
test_sqdist
disp('All tests completed')

author: andrewssobral

repo: lrslibrary


text:
function [enums,Matrices] = mpt_enumerate_binary(Matrices)

binary_var_index = Matrices.binary_var_index;
notbinary_var_index = setdiff(1:Matrices.nu,binary_var_index);

% Detect and extract pure binary equalities. These are used
% to detect SOS constraints, and for pruning
nbin = length(binary_var_index);
only_binary = ~any(Matrices.Aeq(:,notbinary_var_index),2);
Aeq_bin = Matrices.Aeq(find(only_binary),binary_var_index);
%Beq_bin = Beq(find(only_binary),binary_var_index);
beq_bin = Matrices.beq(find(only_binary),:);

% Keep the rest
Matrices.Aeq = Matrices.Aeq(find(~only_binary),:);
Matrices.Beq = Matrices.Beq(find(~only_binary),:);
Matrices.beq = Matrices.beq(find(~only_binary),:);

% Dummy pre-solve to avoid problems
if length(beq_bin)>0
    [ii,jj,kk] = unique([Aeq_bin beq_bin],'rows');
    Aeq_bin = Aeq_bin(jj,:);
    beq_bin = beq_bin(jj,:);
end

% Normalize...
neg_b = find(beq_bin < 0);
Aeq_bin(neg_b,:) = -Aeq_bin(neg_b,:);
beq_bin(neg_b,:) = -beq_bin(neg_b,:);

% Detect and extract simple binary LP constraints
only_binary_lp = find(~any([Matrices.G(:,notbinary_var_index) Matrices.E],2));
Abin = Matrices.G(only_binary_lp,binary_var_index);
bbin = Matrices.W(only_binary_lp);

% Remove pure binary constraints
Matrices.G(only_binary_lp,:) = [];
Matrices.W(only_binary_lp,:) = [];
Matrices.E(only_binary_lp,:) = [];

% if ~isempty(Matrices.binary_var_index)
%     fixed_up = find(Matrices.lb(Matrices.binary_var_index) == 1);
%     fixed_down = find(Matrices.ub(Matrices.binary_var_index) == 0);
%     fixed = [fixed_up fixed_down];
%     if ~isempty(fixed)
%         % Ok, some of the binary are fixed.
%         % Do not enumerate over these
%         % Fix them and correct the constraints that we use
%         % for enumeration and pruning
%         bbin = bbin - Abin(:,fixed)*Matrices.lb(Matrices.binary_var_index(fixed));
%         Abin(:,fixed) = [];
%         beq_bin = beq_bin - Aeq_bin(:,fixed)*Matrices.lb(Matrices.binary_var_index(fixed));
%         Aeq_bin(:,fixed) = [];
%         Matrices.binary_var_index(fixed) = [];
%         binary_var_index = Matrices.binary_var_index;
%         notbinary_var_index = setdiff(1:Matrices.nu,binary_var_index);
%     end
% end

% Detect groups with constraints sum(d_i) == d_j
SOS = [];
variables_in_sos = [];
for i = 1:size(Aeq_bin,1)
    if beq_bin(i) == 0
        [ix,jx,sx] = find(Aeq_bin(i,:));
        if all(abs(sx) == 1)
            j = find(sx == -1);
            if length(j) == 1 & length(jx)>2
                % Aha, we have sum binary(i) == binary(j)
                j = jx(j);
                jx = setdiff(jx,j);
                this_sos = sparse(1:length(jx),jx,1,length(jx),length(binary_var_index));
                this_sos(1:end,j)= 1;
                this_sos(end+1,1)=0;
                if isempty(SOS)
                    SOS = this_sos;
                else
%                    new_sos  = [];
                    SOS = kron(ones(size(SOS,1),1),this_sos) | kron(SOS,ones(size(this_sos,1),1));
%                     for k = 1:size(SOS,1)
%                         for r = 1:size(this_sos,1)
%                             new_sos = [new_sos;this_sos(r,:) | SOS(k,:)];
%                         end
%                     end
%                     if ~isequal(new_sos,new_sos2)
%                         error
%                     end
%                    SOS = new_sos;
                end
                variables_in_sos = [variables_in_sos jx j];
            end
        end
    end
end

% Detect groups with constraints sum(d_i) == 1
for i = 1:size(Aeq_bin,1)
    if beq_bin(i) == 1
        [ix,jx,sx] = find(Aeq_bin(i,:));
        if all(sx == 1)
            % Aha, we have sum binary(jx) == 1
            this_sos = sparse(1:length(jx),jx,1,length(jx),length(binary_var_index));
            if isempty(SOS)
                SOS = this_sos;
            else
               % new_sos  = [];
                SOS = kron(ones(size(SOS,1),1),this_sos) | kron(SOS,ones(size(this_sos,1),1));
%                 for k = 1:size(SOS,1)
%                     for r = 1:size(this_sos,1)
%                         new_sos = [new_sos;this_sos(r,:) | SOS(k,:)];
%                     end
%                 end
%                 SOS = new_sos;
            end
            variables_in_sos = [variables_in_sos jx];
        end
    end
end

variables_not_in_sos = setdiff(1:length(binary_var_index),variables_in_sos);
if ~isempty(variables_not_in_sos)
    % we need to add some more binaries
    n_left = length(variables_not_in_sos);
    all_perms = dec2decbin(0:2^n_left-1,n_left);
    [ix,jx,sx] = find(all_perms);jx = variables_not_in_sos(jx);
    this_sos = sparse(ix,jx,sx,size(all_perms,1),length(binary_var_index));
    if isempty(SOS)
        SOS = this_sos;
    else
        new_sos  = [];
        for k = 1:size(SOS,1)
            for r = 1:size(this_sos,1)
                new_sos = [new_sos;this_sos(r,:) | SOS(k,:)];
            end
        end
        SOS = new_sos;
    end
end

enums = unique(SOS,'rows')';

% FIX :FULL, F**N Linux 6.5 repmat bug
feasible = find(~any(Aeq_bin*enums-repmat(full(beq_bin),1,size(enums,2)),1));
enums = enums(:,feasible);

% Remove binary vertices violating these
if ~isempty(Abin)
    remove = find(any(Abin*enums - repmat(bbin,1,size(enums,2)) >0,1));
%    keep = setdiff(1:size(enums,2),remove);
%    mapper = find(ismember(1:size(enums,2),keep));
%    keep(remove)=0;keep2 = 1:size(enums,2);keep2(find(keep)) = find(keep);
%     map = find(keep);
%     remove = find(ismember(j,remove));
%     i(remove) = [];
%     j(remove) = [];
%     k(remove) = [];
%     enums = sparse(i,
%
%     [i,j,k] = find(enums);
    enums(:,remove)=[];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function shim = cvx_sedumi( shim )

% CVX_SOLVER_SHIM	SeDuMi interface for CVX.
%   This procedure returns a 'shim': a structure containing the necessary
%   information CVX needs to use this solver in its modeling framework.

global cvx___
if ~isempty( shim.solve ),
    return
end
if isempty( shim.name ),
    fname = 'sedumi.m';
    fs = cvx___.fs;
    ps = cvx___.ps;
    int_path = [ cvx___.where, fs ];
    int_plen = length( int_path );
    shim.name = 'SeDuMi';
    shim.dualize = true;
    flen = length(fname);
    fpaths = { [ int_path, 'sedumi', fs, fname ] };
    fpaths = [ fpaths ; which( fname, '-all' ) ];
    old_dir = pwd;
    oshim = shim;
    shim = [];
    for k = 1 : length(fpaths),
        fpath = fpaths{k};
        if ~exist( fpath, 'file' ) || any( strcmp( fpath, fpaths(1:k-1) ) ),
            continue
        end
        new_dir = fpath(1:end-flen-1);
        cd( new_dir );
        tshim = oshim;
        tshim.fullpath = fpath;
        tshim.version = 'unknown';
        is_internal = strncmp( new_dir, int_path, int_plen );
        if is_internal,
            tshim.location = [ '{cvx}', new_dir(int_plen:end) ];
        else
            tshim.location = new_dir;
        end
        try
            fid = fopen(fname);
            otp = fread(fid,Inf,'uint8=>char')';
            fclose(fid);
        catch errmsg
            tshim.error = sprintf( 'Unexpected error:\n%s\n', errmsg.message );
        end
        if isempty( tshim.error ),
            otp = regexp( otp, 'SeDuMi \d\S+', 'match' );
            if ~isempty(otp), tshim.version = otp{end}(8:end); end
            vnum = str2double( tshim.version );
            tshim.check = @check;
            tshim.solve = @solve;
            tshim.eargs = { vnum >= 1.3 && vnum < 1.32 };
            if k ~= 2,
                tshim.path = [ new_dir, ps ];
                if ~isempty(cvx___.msub) && exist([new_dir,fs,cvx___.msub],'dir'),
                   tshim.path = [ new_dir, fs, cvx___.msub, ps, tshim.path ];
                end
            end
        end
        shim = [ shim, tshim ]; %#ok
    end
    cd( old_dir );
    if isempty( shim ),
        shim = oshim;
        shim.error = 'Could not find a SeDuMi installation.';
    end
else
    shim.check = @check;
    shim.solve = @solve;
    vnum = str2double( shim.version );
    shim.eargs = { vnum >= 1.3 && vnum < 1.32 };
end

function found_bad = check( nonls ) %#ok
found_bad = false;

function [ x, status, tol, iters, y, z ] = solve( At, b, c, nonls, quiet, prec, settings, nocplx )

n = length( c );
m = length( b );
K = struct( 'f', 0, 'l', 0, 'q', [], 'r', [], 's', [], 'scomplex', [], 'ycomplex', [] );
reord = struct( 'n', 0, 'r', [], 'c', [], 'v', [] );
reord = struct( 'f', reord, 'l', reord, 'a', reord, 'q', reord, 'r', reord, 's', reord, 'h', reord );
reord.f.n = n;
zinv = [];
for k = 1 : length( nonls ),
    temp = nonls( k ).indices;
    nn = size( temp, 1 );
    nv = size( temp, 2 );
    nnv = nn * nv;
    tt = nonls( k ).type;
    reord.f.n = reord.f.n - nnv;
    if strncmp( tt, 'i_', 2 ),
        error( 'SeDuMi does not support integer variables.' );
    elseif nn == 1 || isequal( tt, 'nonnegative' ),
        reord.l.r = [ reord.l.r ; temp(:) ];
        reord.l.c = [ reord.l.c ; reord.l.n + ( 1 : nnv )' ];
        reord.l.v = [ reord.l.v ; ones( nnv, 1 ) ];
        reord.l.n = reord.l.n + nnv;
    elseif isequal( tt, 'lorentz' ),
        if nn == 2,
            rr = [ temp ; temp ];
            cc = reshape( floor( 1 : 0.5 : 2 * nv + 0.5 ), 4, nv );
            vv = [1;1;-1;1]; vv = vv(:,ones(1,nv));
            reord.a.r = [ reord.a.r ; rr(:) ];
            reord.a.c = [ reord.a.c ; cc(:) + reord.a.n ];
            reord.a.v = [ reord.a.v ; vv(:) ];
            reord.a.n = reord.a.n + nnv;
            zinv = [ zinv ; temp(:) ]; %#ok
        else
            temp = temp( [ end, 1 : end - 1 ], : );
            reord.q.r = [ reord.q.r ; temp(:) ];
            reord.q.c = [ reord.q.c ; reord.q.n + ( 1 : nnv )' ];
            reord.q.v = [ reord.q.v ; ones(nnv,1) ];
            reord.q.n = reord.q.n + nnv;
            K.q = [ K.q, nn * ones( 1, nv ) ];
        end
    elseif isequal( tt, 'semidefinite' ),
        if nn == 3,
            temp = temp( [1,3,2], : );
            tempv = [sqrt(2);sqrt(2);1] * ones(1,nv);
            reord.r.r = [ reord.r.r ; temp(:) ];
            reord.r.c = [ reord.r.c ; reord.r.n + ( 1 : nnv )' ];
            reord.r.v = [ reord.r.v ; tempv(:) ];
            reord.r.n = reord.r.n + nnv;
            K.r = [ K.r, 3 * ones( 1, nv ) ];
            temp = temp(1:2,:);
            zinv = [ zinv ; temp(:) ]; %#ok
        else
            nn = 0.5 * ( sqrt( 8 * nn + 1 ) - 1 );
            str = cvx_create_structure( [ nn, nn, nv ], 'symmetric' );
            K.s = [ K.s, nn * ones( 1, nv ) ];
            [ cc, rr, vv ] = find( cvx_invert_structure( str, 'compact' ) );
            rr = temp( rr );
            reord.s.r = [ reord.s.r; rr( : ) ];
            reord.s.c = [ reord.s.c; cc( : ) + reord.s.n ];
            reord.s.v = [ reord.s.v; vv( : ) ];
            reord.s.n = reord.s.n + nn * nn * nv;
            reord.s.z = reord.s.v;
        end
    elseif isequal( tt, 'hermitian-semidefinite' ),
        if nn == 4,
            temp = temp( [1,4,2,3], : );
            tempv = [sqrt(2);sqrt(2);1;1] * ones(1,nv);
            reord.r.r = [ reord.r.r ; temp(:) ];
            reord.r.c = [ reord.r.c ; reord.r.n + ( 1 : nnv )' ];
            reord.r.v = [ reord.r.v ; tempv(:) ];
            reord.r.n = reord.r.n + nnv;
            reord.r.z = reord.r.v;
            K.r = [ K.r, 4 * ones( 1, nv ) ];
            temp = temp(1:2,:);
            zinv = [ zinv ; temp(:) ]; %#ok
        elseif nocplx,
            % SeDuMi's complex SDP support was broken with the 1.3 update. So
            % we must use the following complex-to-real SDP conversion to work
            % around it, at a modest cost of problem size.
            %   X >= 0 <==> exists [ Y1, Y2^T ; Y2, Y3 ] >= 0 s.t.
            %               Y1 + Y3 == real(X), Y2 - Y2^T == imag(X)
            nsq = nn; nn = sqrt( nn );
            str = cvx_create_structure( [ nn, nn, nv ], 'hermitian' );
            [ cc, rr, vv ] = find( cvx_invert_structure( str, 'compact' ) );
            cc = cc - 1;
            mm = floor( cc / nsq );
            cc = cc - mm * nsq;
            jj = floor( cc / nn );
            ii = cc - jj * nn + 1;
            jj = jj + 1;
            mm = mm + 1;
            vr = real( vv );
            vi = imag( vv );
            ii = [ ii + nn * ~vr ; ii + nn * ~vi ];
            jj = [ jj ; jj + nn ]; %#ok
            vv = sqrt( 0.5 ) * [ vr + vi ; vr - vi ];
            rr = [ rr ; rr ]; %#ok
            mm = [ mm ; mm ]; %#ok
            [ jj, ii ] = deal( min( ii, jj ), max( ii, jj ) );
            cc = ii + ( jj - 1 ) * ( 2 * nn ) + ( mm - 1 ) * ( 4 * nsq );
            K.s = [ K.s, 2 * nn * ones( 1, nv ) ];
            rr = temp( rr );
            reord.s.r = [ reord.s.r; rr( : ) ];
            reord.s.c = [ reord.s.c; cc( : ) + reord.s.n ];
            reord.s.v = [ reord.s.v; vv( : ) ];
            reord.s.n = reord.s.n + 4 * nsq * nv;
            reord.s.z = reord.s.v;
        else
            % SeDuMi's complex SDP support was restored in v1.33.
            K.scomplex = [ K.scomplex, length( K.s ) + ( 1 : nv ) ];
            nn = sqrt( nn );
            str = cvx_create_structure( [ nn, nn, nv ], 'hermitian' );
            K.s = [ K.s, nn * ones( 1, nv ) ];
            stri = cvx_invert_structure( str, 'compact' )';
            [ rr, cc, vv ] = find( stri );
            rr = temp( rr );
            reord.s.r = [ reord.s.r; rr( : ) ];
            reord.s.c = [ reord.s.c; cc( : ) + reord.s.n ];
            reord.s.v = [ reord.s.v; vv( : ) ];
            reord.s.n = reord.s.n + size( stri, 2 );
            reord.s.z = reord.s.v;
        end
    else
        error( 'Unsupported nonlinearity: %s', tt );
    end
end
if reord.f.n > 0,
    reord.f.r = ( 1 : n )';
    reord.f.r( [ reord.l.r ; reord.a.r ; reord.q.r ; reord.r.r ; reord.s.r ] ) = [];
    reord.f.c = ( 1 : reord.f.n )';
    reord.f.v = ones(reord.f.n,1);
end
n_d = max( m - n - reord.f.n + 1, isempty( At ) );
if n_d,
    reord.l.n = reord.l.n + n_d;
end
K.f = reord.f.n;
K.l = reord.l.n + reord.a.n;
n_out = reord.f.n;
reord.l.c = reord.l.c + n_out; n_out = n_out + reord.l.n;
reord.a.c = reord.a.c + n_out; n_out = n_out + reord.a.n;
reord.q.c = reord.q.c + n_out; n_out = n_out + reord.q.n;
reord.r.c = reord.r.c + n_out; n_out = n_out + reord.r.n;
reord.s.c = reord.s.c + n_out; n_out = n_out + reord.s.n;
reord = sparse( ...
    [ reord.f.r ; reord.l.r ; reord.a.r ; reord.q.r ; reord.r.r ; reord.s.r ], ...
    [ reord.f.c ; reord.l.c ; reord.a.c ; reord.q.c ; reord.r.c ; reord.s.c ], ...
    [ reord.f.v ; reord.l.v ; reord.a.v ; reord.q.v ; reord.r.v ; reord.s.v ], ...
    n, n_out );

At = reord' * At;
c  = reord' * c;
pars.free = K.f > 1 && nnz( K.q );
pars.eps     = prec(1);
pars.bigeps  = prec(3);
if quiet,
    pars.fid = 0;
end
add_row = isempty( At );
if add_row,
    K.f = K.f + 1;
    At = sparse( 1, 1, 1, n_out + 1, 1 );
    b = 1;
    c = [ 0 ; c ];
end
[ xx, yy, info ] = cvx_run_solver( @sedumi, At, b, c, K, pars, 'xx', 'yy', 'info', settings, 5 );
if add_row,
    xx = xx(2:end);
    yy = zeros(0,1);
    At = zeros(n_out,0);
    % b  = zeros(0,1);
    c  = c(2:end);
end
if ~isfield( info, 'r0' ) && info.pinf,
    info.r0 = 0;
    info.iter = 0;
    info.numerr = 0;
end
tol = info.r0;
iters = info.iter;
xx = full( xx );
yy = full( yy );
status = '';
if info.pinf ~= 0,
    status = 'Infeasible';
    x = NaN * ones( n, 1 );
    y = yy;
    z = - real( reord * ( At * yy ) );
    if add_row, y = zeros( 0, 1 ); end
elseif info.dinf ~= 0
    status = 'Unbounded';
    y = NaN * ones( m, 1 );
    z = NaN * ones( n, 1 );
    x = real( reord * xx );
else
    x = real( reord * xx );
    y = yy;
    z = real( reord * ( c - At * yy ) );
    if add_row, y = zeros( 0, 1 ); end
end
if ~isempty(zinv),
    z(zinv) = z(zinv) * 0.5;
end
if info.numerr == 2,
    status = 'Failed';
    if any( K.q == 2 ),
        warning( 'CVX:SeDuMi', cvx_error_format( 'This solver failure may possibly be due to a known bug in the SeDuMi solver. Try switching to SDPT3 by inserting "cvx_solver sdpt3" into your model.', ...
            [66,75], false, '' ) );
    end
else
    if isempty( status ),
        status = 'Solved';
    end
    if info.numerr == 1 && info.r0 > prec(2),
        status = [ 'Inaccurate/', status ];
    end
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
 function st = nufft_table_init(om, Nd, Jd, Kd, n_shift, Ld, varargin)
%function st = nufft_table_init(om, Nd, Jd, Kd, n_shift, Ld, ...)
%|
%| Initialize structure for d-dimension NUFFT using table-based interpolator,
%| This should be called only by nufft_init with its 'table' argument!
%| (So this m-file belongs in the private directory.)
%|
%| in
%|	om, Nd, Jd, Kd, n_shift		see nufft_init.m
%|	Ld [d]				table over-sampling factor(s)
%|
%| optional arguments (same as nufft_init.m, to specify interp type)
%| out
%|	st.interp_table		function: X = st.interp_table(st, Xk);
%|	st.interp_table_adj	function: Xk = st.interp_table_adj(st, Xk);
%|	st.phase_shift		phase shift due to n_shift
%|	st.*			same as in nufft_init.m
%|
%| Copyright 2004-3-30, Jeff Fessler and Yingying Zhang, University of Michigan

if nargin < 5, help(mfilename), error args, end

% dimensionality of input space (usually 2 or 3)
dd = length(Nd);
if dd > 1 && length(Ld) == 1
	Ld = Ld * ones(1,dd); % allow scalar L to be used for all dimensions
end
if dd ~= length(Jd) || dd ~= length(Kd) || dd ~= length(Ld)
	printm('dd:'), disp(dd)
	printm('Jd:'), disp(Jd)
	printm('Kd:'), disp(Kd)
	printm('Ld:'), disp(Ld)
	error 'inconsistent dim'
end
if dd ~= size(om,2), error('omega needs %d columns', dd), end

% D-dim NUFFT structure with desired scaling factors, using fake omega=0
st = nufft_init(zeros(size(Nd)), Nd, Jd, Kd, 0*n_shift, varargin{:});
st.p = [];
st.om = om;
st.Ld = Ld;
st.M = size(om,1);
st.n_shift = n_shift; % trick: store it in structure for history

% phase shift associated with indexing n_shift
% trick: because we put this phase shift here, we use "0" for n_shift below
% when building the tables.  in fact, n_shift definitely cannot be in table.
if any(n_shift ~= 0)
	st.phase_shift = exp(1i * (om * n_shift(:))); % [M 1]
end

% 1D tables for each dimension
for id = 1:dd
	J = Jd(id);
	L = Ld(id);
	K = Kd(id);

	% This is a slow and inefficient way to get the table
	% since it builds a huge sparse matrix but only uses 1 column!
	t0 = [-J*L/2:J*L/2]' / L; % [J*L+1]
	if 0
		om0 = t0 * 2*pi/K; % gam(id)
		s1 = nufft_init(om0, Nd(id), J, K, ...
			0*n_shift(id), varargin{:}); % user's interpolator
		h0 = full(s1.p(:,1));
	end

	% This efficient way uses only "J" columns of sparse matrix!
	% The trick to this is to use fake small values for N and K,
	% which works for interpolators that depend only on the ratio K/N.
	ktype = varargin{1};
	if 0 && (streq(ktype, 'minmax:kb') || streq(ktype, 'kaiser'))

		if streq(ktype, 'minmax:kb')
			args = {'minmax:user', {st.alpha{id}}, {st.beta{id}}};
		elseif streq(ktype, 'kaiser')
			args = {'kaiser'};
			if length(varargin) > 1, warning 'bug?', end
		else
			error 'bug'
		end

		Nfake = J;
		Kfake = Nfake * K/Nd(id);
		t1 = J/2-1 + [0:(L-1)]' / L;	% [L]
		om1 = t1 * 2*pi/Kfake;		% "gam"
		s1 = nufft_init(om1, Nfake, J, Kfake, 0, args{:});

		h = [];
		for jj=J:-1:1
			h = [h; full(s1.p(:,jj))]; % stack up pieces
		end
		h = [h; 0]; % finally J*L+1 long.
		h = h .* exp(1i * pi * t0 * (1/K - 1/Kfake)); % fix phase

		if 0	% testing
			clf, subplot(211)
			plot(t0, real(h0), 'c-', t0, real(h), 'y.')
			plot(t0, real(h0-h), 'y.')
			subplot(212)
			plot(t0, imag(h0), 'c-', t0, imag(h), 'y.')
			plot(t0, imag(h0-h), 'y.')
			printf('diff=%g', max_percent_diff(h0, h))
			keyboard
		end
		st.h{id} = double(h); % table mex files want double

	% This way is "J times faster" than the slow way, but still not ideal.
	% It works for any user-specified interpolator.
	else
		t1 = J/2-1 + [0:(L-1)]' / L; % [L]
		om1 = t1 * 2*pi/K; % gam(id)

		% trick: handle {'kaiser', [kb_alf], [kb_m]} in multidim case
		if streq(varargin{1}, 'kaiser') && length(varargin) == 3
			kb_alf = varargin{2};
			kb_m = varargin{3};
			if any(kb_alf ~= kb_alf(1)) || any(kb_m ~= kb_m(1))
				error 'table based needs identical interpolator for each dimension'
			end
			varargin{2} = kb_alf(1);
			varargin{3} = kb_m(1);
		end

		s1 = nufft_init(om1, Nd(id), J, K, 0*n_shift(id), varargin{:});

		h = [];
		for jj=J:-1:1
			h = [h; full(s1.p(:,jj))];
		end
		h = [h; 0];

		st.h{id} = double(h); % table mex files want double
		if 0	% testing
			clf, subplot(211)
			plot(t0, real(h0), 'c-', t0, real(h), 'y.')
			subplot(212)
			plot(t0, imag(h0), 'c-', t0, imag(h), 'y.')
			printf('diff=%g', max_percent_diff(h0, h))
			keyboard
		end
	end

	if isreal(st.h{id})
		if J == 1
			st.h{id} = complex(st.h{id});
		else
			whos
			warning 'real kernel?'
			keyboard
		end
	end

end

% interface routines
% st.interp_table = inline('nufft_table_interp(st, Xk)', 'st', 'Xk');
st.interp_table		= @nufft_table_interp;
st.interp_table_adj	= @nufft_table_adj;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function det = r8mat_det_4d ( a )

%*****************************************************************************80
%
%% R8MAT_DET_4D computes the determinant of a 4 by 4 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(4,4), the matrix whose determinant is desired.
%
%    Output, real DET, the determinant of the matrix.
%
  det = ...
      a(1,1) * ( ...
        a(2,2) * ( a(3,3) * a(4,4) - a(3,4) * a(4,3) ) ...
      - a(2,3) * ( a(3,2) * a(4,4) - a(3,4) * a(4,2) ) ...
      + a(2,4) * ( a(3,2) * a(4,3) - a(3,3) * a(4,2) ) ) ...
    - a(1,2) * ( ...
        a(2,1) * ( a(3,3) * a(4,4) - a(3,4) * a(4,3) ) ...
      - a(2,3) * ( a(3,1) * a(4,4) - a(3,4) * a(4,1) ) ...
      + a(2,4) * ( a(3,1) * a(4,3) - a(3,3) * a(4,1) ) ) ...
    + a(1,3) * ( ...
        a(2,1) * ( a(3,2) * a(4,4) - a(3,4) * a(4,2) ) ...
      - a(2,2) * ( a(3,1) * a(4,4) - a(3,4) * a(4,1) ) ...
      + a(2,4) * ( a(3,1) * a(4,2) - a(3,2) * a(4,1) ) ) ...
    - a(1,4) * ( ...
        a(2,1) * ( a(3,2) * a(4,3) - a(3,3) * a(4,2) ) ...
      - a(2,2) * ( a(3,1) * a(4,3) - a(3,3) * a(4,1) ) ...
      + a(2,3) * ( a(3,1) * a(4,2) - a(3,2) * a(4,1) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x=cover(A,B)
x.t=min(A.t,B.t);
x.b=max(A.b,B.b);
x.l=min(A.l,B.l);
x.r=max(A.r,B.r);

author: Zzh-tju

repo: DIoU

save_path: github-repos/MATLAB/Zzh-tju-DIoU


text:
%  Figure 7.33      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 7.33
% fig7_33.m
clf;
f=[0 1;-1 0];
g=[0;1];
h=[1 0];
K=[3 4];
fc=f-g*K;
gc=[4*g [1;0]];
hc=[h;[0 1]; -K/4];
jc=[0 0;0 0; 1 0];
A=[fc [0;0]; -101*h-K -10];
B=[0 1;4 0;4 0];
C=eye(3);
C(3,1) = 10;
D=[0 0;0 0;0 0];
t=0:.1:4;
[y]=impulse(A,B,C,D,2,t);
plot(t,y,'LineWidth',2);
xlabel('Time (sec)');
ylabel('Amplitude');
text(.2,3,'$$\mathsf{\hat x_2}$$','interpreter','latex','FontSize',11);
text(.75,1,'x_1');
text(.25,-1,'x_2');
title('Fig.7.33  Initial condition response with reduced order estimator')
nicegrid

author: Sable

repo: mcbench-benchmarks


text:
function A = SOMPerr(D,X,errorGoal,BETA,LL)

[n,P]=size(X);
[n,K]=size(D);
%%%%%%%%%%%%%%%guiyihua
BETA_matrix1=repmat(BETA,[1 P]);
BETA_matrix2=repmat(BETA,[1 K]);
X_scale=X.*BETA_matrix1;
D_scale=D.*BETA_matrix2;
X=X_scale;
D=D_scale;

%%%%%%%%%%%%%%%%%%%%
E2 = errorGoal^2*n;
% maxNumCoef = K/3;
maxNumCoef = LL;
A = sparse(K,P);
a = [];
residual = X;
indx = [];
currResNorm2 = sum(sum(residual.^2))/P;
j = 0;

while currResNorm2 >E2 && j < maxNumCoef

        j = j + 1;
        proj = mean(D' * residual,2) ;
%         proj = sum(proj,2);
        pos = find(abs(proj) == max(abs(proj)));
        pos = pos(1);
        indx(j) = pos;
        a = pinv(D(:,indx(1:j)))*X;
        residual = X - D(:,indx(1:j))*a;
		currResNorm2 = sum(sum(residual.^2))/P;
end
if (~isempty(indx))
       A(indx,:) = a;
end

return;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function z = cinf_1D(d,len,params)

% Generating sensor signals for a 1D sensor array in a cylindrically
% isotropic noise field [1]
%
%    z = cinf_1D(d,len,params)
%
% Input parameters:
%    d            : sensor distances
%    len          : desired data length
%    params.fs    : sample frequency
%    params.c     : sound velocity in m/s
%    params.N_phi : number of cylindrical angles
%
% Output parameters:
%    z            : output signal
%
% References:
%    [1] E.A.P. Habets and S. Gannot, 'Generating sensor signals
%        in isotropic noise fields', The Journal of the Acoustical
%        Society of America, Vol. 122, Issue 6, pp. 3464-3470, Dec. 2007.
%
% Authors:  E.A.P. Habets and S. Gannot
%
% History:  2007-11-02 - Initial version
%           2010-09-16 - Minor corrections
%           2017-06-29 - Use native waitbar
%
% Copyright (C) 2007-2017 E.A.P. Habets
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M = length(d);                        % Number of sensors
NFFT = 2^ceil(log2(len));             % Number of frequency bins
X = zeros(M,NFFT/2+1);

if ~isfield(params,'fs')
    fs = 8000;                        % Default
else
    fs = params.fs;
end
if ~isfield(params,'c')
    c = 340;                          % Default
else
    c = params.c;
end
if ~isfield(params,'N_phi')
    N_phi = 64;                       % Default
else
    N_phi = params.N_phi;
end

w = 2*pi*fs*(0:NFFT/2)/NFFT;
phi = 2*pi*(0:1/N_phi:1-1/N_phi);

% Calculate relative sensor distances
d_rel = d - d(1);

% Initialize waitbar
h = waitbar(0,'Generating sensor signals...');

% Calculate sensor signals in the frequency domain
for phi_idx = 1:N_phi
    waitbar(phi_idx/N_phi);
    X_prime = randn(1,NFFT/2+1) + 1i*randn(1,NFFT/2+1);
    X(1,:) = X(1,:) + X_prime;
    for m = 2:M
        Delta = d_rel(m)*cos(phi(phi_idx));
        X(m,:) = X(m,:) + X_prime.*exp(-1i*Delta*w/c);
    end
end
X = X/sqrt(N_phi);

% Transform to time domain
X = [sqrt(NFFT)*real(X(:,1)), sqrt(NFFT/2)*X(:,2:NFFT/2),...
    sqrt(NFFT)*real(X(:,NFFT/2+1)), sqrt(NFFT/2)*conj(X(:,NFFT/2:-1:2))];
z = real(ifft(X,NFFT,2));

% Truncate output signals
z = z(:,1:len);

% Close waitbar
close(h);

author: ZitengWang

repo: MASP

save_path: github-repos/MATLAB/ZitengWang-MASP


text:
%kkabsdiff 'Output = Absolute Value of (Input 1 - (Input 2 or Constant))'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kabsdiff.pane file
%
% Parameters:
% InputFile: i1 'Input 1', required: 'First input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% InputFile: i2 'Input 2', optional: 'Second input data object'
%
% Example: o = kkabsdiff({i1, i2}, {'i1','';'o','';'i2',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kabsdiff - Output = Absolute Value of (Input 1 - (Input 2 or Constant))
%
%  DESCRIPTION
% Output = Absolute Value(Input 1 - (Input 2 or a Constant))
%
% The "Absolute Difference" operator computes the absolute value
% of the result of subtracting either \fBInput 2"
% or the \fBConstant" value, which ever is specified by the user,
% from \fBInput 1".
%
% Executing "Absolute Difference" runs the program \fIkarith2\fP
% with the -absdiff flag.
%
%
%  "Data Type - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Data Type - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_2input
%
%  "Map Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "Map Data - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_2input
%
%  "Validity Mask" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_2input
%
%  "Input Objects of Different Sizes" 5
% .cI $DATAMANIP/repos/shared/man/sections/resize_2input
% The values used to pad the data when input files are not the same size are
% (0.0, 0.0).
%
%  "Explicit Location and Time Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Explicit Location and Time Data - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_2input
%
%  "Failure Modes - Single Input"
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
%  "Failure Modes - Two Input Objects"
% .cI $DATAMANIP/repos/shared/man/sections/fail_2input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::karith2
%
%  RESTRICTIONS
% Operations on complex data are not supported at this time.  If the
% input object is complex, only the real component of the data is
% processed, and the output object is real.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkabsdiff(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkabsdiff(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output';'i2', '__input'};
maxval={0,0,1};
minval={0,0,1};
istoggle=[0,0,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','InputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'karith2"  -absdiff'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef PTKReporting < MimReporting
    % PTKReporting. Legacy support function. Replaced by MimReporting
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    methods
        function obj = PTKReporting(progress_dialog, verbose_mode, log_file_name)
            obj = obj@MimReporting(progress_dialog, verbose_mode, log_file_name);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function op = proj_l2group( q,group_indices )
%PROJ_L2GROUP   Projection onto the intersection of scaled 2-norm balls.
%    OP = PROJ_L2GROUP( Q, GROUP_INDICES ) returns an operator implementing the
%    indicator function for the intersection of 2-norm ball of size q_k,
%    i.e. intersection_k B_k
%    where B_k = { X |  norm( X(indK) ) <= q(k) }
%    and indK is indexed by group_indices. Group_indices keeps track
%    of the last index in each index set.
%    The index sets must be non-overlapping.
%    For example, if K=2, and ind1 = 1:10 and ind2 = 11:20,
%    then group_indices = [10,20].
%    Q must be a vector of length K or a vector of length N.
%
% With contributions from Joseph Salmon
%
% See also: proj_l2.m

if isempty(q), q = 1; end
% Make sure q is a column vector
if size(q,1) == 1 && size(q,2) > 1, q = q.'; end
q = expand_q(q,group_indices);
op = @(varargin)proj_l2_q( q,group_indices, varargin{:} );
end

function [ v, x ] = proj_l2_q( q,group_indices, x, t )
v = 0;
nrm=group_norm(x,group_indices);

switch nargin,
	case 3,
		if nargout == 2,
			error( 'This function is not differentiable.' );
		elseif norm( nrm, 'inf') > q,
			v = Inf;
		end
	case 4,
            x = x .* min(1,( q ./ nrm ));
	otherwise,
		error( 'Not enough arguments.' );
end
end

function y=group_norm(x,group_indices)
% the groups are indexed by their end point
    K=length(group_indices);
    y=zeros(size(x));
    j=1; % start of the group
    for k=1:K
        end_point=group_indices(k);
        y(j:end_point)=norm(x(j:end_point));
        j=end_point+1;
    end
end

function Q = expand_q(q,group_indices)
% If we have 2 groups, each of size 10
% (group_indices = [10,20]), then q should
% be a vector of size 20. But it's more convenient
% for the user to pass in a vector of size 2, so this function
% will convert it...
n = group_indices(end);
K = length(group_indices);
if length(q) < n
    Q = ones(n,1);
    j=1; % start of the group
    for k=1:K
        end_point=group_indices(k);
        Q(j:end_point) = q(k);
        j=end_point+1;
    end
else
    Q = q;
end
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
classdef RotatedRect
    %ROTATEDRECT  The class represents rotated (i.e. not up-right) rectangles on a plane
    %
    % Each rectangle is specified by the center point (mass center), length of
    % each side (represented by `[width,height]`) and the rotation angle in
    % degrees.
    %
    % The sample `RotatedRect_demo.m` demonstrates how to use RotatedRect.
    %
    % See also: cv.CamShift, cv.fitEllipse, cv.minAreaRect
    %

    methods (Static)
        function rrect = from3points(pt1, pt2, pt3)
            %FROM3POINTS  Create a rotated rectangle from 3 points
            %
            %     rrect = cv.RotatedRect.from3points(pt1, pt2, pt3)
            %
            % ## Input
            % * __pt1__, __pt2__, __pt3__ Any 3 end points `[x,y]` of the
            %   rotated rectangle. They must be given in order (either
            %   clockwise or anticlockwise). By definition, the two sides
            %   formed by these three points must be perpendicular.
            %
            % ## Output
            % * __rrect__ output rotated rectangle. A structure with the
            %   following fields:
            %   * __center__ The rectangle mass center `[x,y]`.
            %   * __size__ Width and height of the rectangle `[w,h]`.
            %   * __angle__ The rotation angle in a clockwise direction. When
            %     the angle is 0, 90, 180, 270 etc., the rectangle becomes an
            %     up-right rectangle.
            %
            rrect = RotatedRect_('from3points', pt1, pt2, pt3);
        end

        function pts = points(rrect)
            %POINTS  Returns 4 vertices of the rectangle
            %
            %     pts = cv.RotatedRect.points(rrect)
            %
            % ## Input
            % * __rrect__ rotated rectangle. A structure with the following
            %   fields:
            %   * __center__ The rectangle mass center `[x,y]`.
            %   * __size__ Width and height of the rectangle `[w,h]`.
            %   * __angle__ The rotation angle in a clockwise direction. When
            %     the angle is 0, 90, 180, 270 etc., the rectangle becomes an
            %     up-right rectangle.
            %
            % ## Output
            % * __pts__ 4-by-2 points matrix of the rectangle vertices
            %   `[x1 y1; x2 y2; x3 y3; x4 y4]`. The order is bottom-left,
            %   top-left, top-right, bottom-right.
            %
            % See also: cv.boxPoints, bbox2points
            %
            pts = RotatedRect_('points', rrect);
        end

        function rect = boundingRect(rrect)
            %BOUNDINGRECT  Returns the minimal up-right integer rectangle containing the rotated rectangle
            %
            %     rect = cv.RotatedRect.boundingRect(rrect)
            %
            % ## Input
            % * __rrect__ rotated rectangle. A structure with the following
            %   fields:
            %   * __center__ The rectangle mass center `[x,y]`.
            %   * __size__ Width and height of the rectangle `[w,h]`.
            %   * __angle__ The rotation angle in a clockwise direction. When
            %     the angle is 0, 90, 180, 270 etc., the rectangle becomes an
            %     up-right rectangle.
            %
            % ## Output
            % * __rect__ bounding rectangle, a 1-by-4 vector `[x, y, w, h]`
            %
            % See also: cv.RotatedRect.boundingRect2f
            %
            rect = RotatedRect_('boundingRect', rrect);
        end

        function rect = boundingRect2f(rrect)
            %BOUNDINGRECT  returns the minimal (exact) floating point rectangle containing the rotated rectangle
            %
            %     rect = cv.RotatedRect.boundingRect2f(rrect)
            %
            % ## Input
            % * __rrect__ rotated rectangle. A structure with the following
            %   fields:
            %   * __center__ The rectangle mass center `[x,y]`.
            %   * __size__ Width and height of the rectangle `[w,h]`.
            %   * __angle__ The rotation angle in a clockwise direction. When
            %     the angle is 0, 90, 180, 270 etc., the rectangle becomes an
            %     up-right rectangle.
            %
            % ## Output
            % * __rect__ bounding rectangle, a 1-by-4 vector `[x, y, w, h]`
            %
            % Not intended for use with images.
            %
            % See also: cv.RotatedRect.boundingRect
            %
            rect = RotatedRect_('boundingRect2f', rrect);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef BodyCenteredInertialFrame < AbstractReferenceFrame
    %BodyCenteredInertialFrame Summary of this class goes here
    %   Detailed explanation goes here

    properties
        bodyInfo KSPTOT_BodyInfo
        celBodyData

        bodyRotMatFromGlobalInertialToBodyInertial(3,3) double
    end

    properties(Constant)
        typeEnum = ReferenceFrameEnum.BodyCenteredInertial
    end

    methods
        function obj = BodyCenteredInertialFrame(bodyInfo, celBodyData)
            obj.bodyInfo = bodyInfo;
            obj.celBodyData = celBodyData;
            obj.bodyRotMatFromGlobalInertialToBodyInertial = obj.bodyInfo.bodyRotMatFromGlobalInertialToBodyInertial;
        end

        function [posOffsetOrigin, velOffsetOrigin, angVelWrtOrigin, rotMatToInertial] = getOffsetsWrtInertialOrigin(obj, time, ~, ~)
            [rVectB, vVectB] = getPositOfBodyWRTSun(time, [obj.bodyInfo], obj(1).celBodyData);

            posOffsetOrigin = rVectB;
            velOffsetOrigin = vVectB;

            [angVelWrtOrigin, rotMatToInertial] = obj.getAngVelWrtOriginAndRotMatToInertial(time, [], []);
        end

        function rotMatToInertial = getRotMatToInertialAtTime(obj, time, ~, ~)
            if(numel(obj) == 1)
                rotMatToInertial = repmat(obj.bodyRotMatFromGlobalInertialToBodyInertial, [1, 1, length(time)]);

            else
                rotMatToInertial = NaN([3, 3, length(time)]);

                [uniObj,~,ic] = unique(obj,'stable');
                for(i=1:length(uniObj)) %#ok<*NO4LP>
                    subObj = uniObj(i);
                    bool = i == ic;
                    subTimes = time(bool);

                    rotMatToInertial(:,:,bool) = repmat(subObj.bodyRotMatFromGlobalInertialToBodyInertial, [1, 1, numel(subTimes)]);
                end
            end
        end

        function [angVelWrtOrigin, rotMatToInertial] = getAngVelWrtOriginAndRotMatToInertial(obj, time, vehElemSet, bodyInfoInertialOrigin)
            angVelWrtOrigin = repmat([0;0;0], [1, length(time)]);
            rotMatToInertial = obj.getRotMatToInertialAtTime(time, vehElemSet, bodyInfoInertialOrigin);
        end

        function tf = frameOriginIsACelBody(obj)
            tf = true;
        end

        function bodyInfo = getOriginBody(obj)
            bodyInfo = [obj.bodyInfo];
        end

        function setOriginBody(obj, newBodyInfo)
            obj.bodyInfo = newBodyInfo;
        end

        function nameStr = getNameStr(obj)
            nameStr = sprintf('Inertial Frame (Origin: %s)', obj.bodyInfo.name);
        end

        function newFrame = editFrameDialogUI(obj, ~)
            if(isempty(obj.celBodyData))
                obj.celBodyData = obj.bodyInfo.celBodyData;
            end

            [bodiesStr, sortedBodyInfo] = ma_getSortedBodyNames(obj.celBodyData);
            for(i=1:length(sortedBodyInfo))
                sortedBodyInfoArr(i) = sortedBodyInfo{i}; %#ok<AGROW>
            end

            ind = find(obj.bodyInfo == sortedBodyInfoArr,1,'first');

            out = AppDesignerGUIOutput();
            listdlgARH_App('ListString',bodiesStr, 'SelectionMode','single', 'InitialValue',ind, ...
                           'Name','Select Frame Central Body', ...
                           'PromptString', {'Select the central body of the reference frame:'}, ...
                           'ListSize', [250 300], ...
                           'out',out);
            selection = out.output{1};
            ok = out.output{2};

            if(ok == 1)
                newBodyInfo = sortedBodyInfoArr(selection);
                newFrame = BodyCenteredInertialFrame(newBodyInfo, obj.celBodyData);
            else
                newFrame = obj;
            end
        end
    end

    methods
        function bool = eq(A,B)
            try
                if(numel(A) == 1 && numel(B) == 1 && ...
                   A.typeEnum == ReferenceFrameEnum.BodyCenteredInertial && ...
                   B.typeEnum == ReferenceFrameEnum.BodyCenteredInertial)
                    bool = A.bodyInfo == B.bodyInfo;
                else
                    bool = eq@handle(A,B);
                end
            catch ME %#ok<NASGU>
                bool = eq@handle(A,B);
            end
        end

        function bool = ne(A,B)
            bool = not(eq(A,B));
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            obj.bodyRotMatFromGlobalInertialToBodyInertial = obj.bodyInfo.bodyRotMatFromGlobalInertialToBodyInertial;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function r = toep_cholesky_upper ( n, g )

%*****************************************************************************80
%
%% TOEP_CHOLESKY_UPPER: upper Cholesky factor of a compressed Toeplitz matrix.
%
%  Discussion:
%
%    The Toeplitz matrix A is supplied in a compressed form G.
%
%    The Toeplitz matrix must be positive semi-definite.
%
%    After factorization, A = R' * R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2012
%
%  Author:
%
%    Michael Stewart
%
%  Reference:
%
%    Michael Stewart,
%    Cholesky factorization of semi-definite Toeplitz matrices.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real G(2,N), the compressed Toeplitz matrix.
%    G(1,1:N) contains the first row.
%    G(2,2:N) contains the first column.
%
%    Output, real R(N,N), the upper Cholesky factor.
%
  r = zeros ( n, n );
  r(1,1:n) = g(1,1:n);
  g(1,2:n) = g(1,1:n-1);
  g(1,1) = 0.0;
  for i = 2 : n
    rho = - g(2,i) / g(1,i);
    g(:,i:n) = ( [ 1 rho; rho 1 ] * g(:,i:n) ) / sqrt ( ( 1 - rho ) * ( 1 + rho ) );
    r(i,i:n) = g(1,i:n);
    g(1,i+1:n) = g(1,i:n-1);
    g(1,i) = 0.0;
  end

  return
end


text:
function [Ymat Xmat N n m p T k q h]=panel6prelim(data_endo,data_exo,const,lags)

% first compute N, the number of units, as the dimension of the data_endo matrix
N=size(data_endo,3);

% compute p, the number of lags in the model
p=lags;

% then compute n, the number of endogenous variables in the model; it is simply the number of columns in the matrix 'data_endo'
n=size(data_endo,2);

% if the constant has been selected, augment the matrix of exogenous with a column of ones (number of rows equal to number of rows in data_endo)
if const==1
data_exo=[ones(size(data_endo,1),1) data_exo];
% if no constant was included, do nothing
else
end

% compute m, the number of exogenous variables in the model
% if data_exo is empty, set m=0
if isempty(data_exo)==1
m=0;
% if data_exo is not empty, count the number of exogenous variables that will be included in the model
else
m=size(data_exo,2);
% Also, trim a number initial rows equal to the number of lags, as they will be suppressed from the endogenous as well to create initial conditions
data_exo=data_exo(p+1:end,:);
end

% determine k, the number of parameters to estimate in each equation; it is equal to np+m
k=N*n*p+m;

% determine q, the total number of VAR parameters for each unit
q=n*k;

% determine h, the total number of VAR parameters for the whole model
h=N*q;

% obtain Ymat and Xmat
temp=[];
% stack the matrices of endogenous variables to obtain a temporary matrix
for ii=1:N
temp=[temp data_endo(:,:,ii)];
end
% use the lagx function on this matrix
temp=bear.lagx(temp,lags);

% set Ymat as the first Nn columns of the result
Ymat=temp(:,1:N*n);

% to build Xmat, take off the Nn initial columns of temp, and concatenate the exogenous
Xmat=[temp(:,N*n+1:end) data_exo];

% Define T, the number of periods of the model, as the number of rows of X
T=size(Xmat,1);

author: european-central-bank

repo: BEAR-toolbox


text:
function y = shiftl(A,row,shift,type)

%  PURPOSE:
%  --------
%  y = shiftl (A, row, shift) moves #row of matrix A to the left
%                                by #shift positions.
%
%  INPUT ARGUMENTS:
%  ----------------
%  'A' is the input matrix. ('A' can be a vector)
%
%  'row' is the number of the row to be shifted. If 'row' is zero,
%  then all rows in the matrix are shifted.
%
%  'shift' is the number of positions by which the row is shifted
%  to the right.
%
%  'type' is an optional argument.
%
%         The shifted matrix-elements are discarded if this argument
%         is 0 or is omitted,
%         then vacated spaces to the right are filled with zeros.
%
%         The shifted matrix-elements are retained if 'type' is 1
%         or any other non-zero value,
%         then vacated spaces to the right are filled with the shifted
%         row-elements from the left (i.e. "wraparound").
%
%  EXAMPLES:   A = [1 2 3 4 5;
%  ---------        6 7 8 9 0]
%
%              y = shiftl(A,0,2)   --> [3 4 5 0 0;
%                                       8 9 0 0 0]
%              y = shiftl(A,1,2)   --> [3 4 5 0 0;
%                                       6 7 8 9 0]
%              y = shiftl(A,1,2,0) --> [3 4 5 0 0;
%                                       6 7 8 9 0]
%              y = shiftl(A,1,2,1) --> [3 4 5 1 2;
%                                       6 7 8 9 0]
%              y = shiftl(A,0,2,1) --> [3 4 5 1 2;
%                                       8 9 0 6 7]
%              B = [1 2 3 4 5]
%
%              z = shiftl(B,1,3)   --> [4 5 0 0 0]
%
%  SEE ALSO:  shiftr, shiftu, shiftd.

[M,N] = size(A);
if row > M | row < 0, error('Invalid Row'); end
if shift < 0, error('Negative shift value - use "shiftr" instead'); end
if shift > N, error('Shift value exceeds number of columns'); end

if row == 0
   if nargin == 4 & type ~= 0
      A = [A(:,1+shift:N) A(:,1:shift)];
   else
      A = [A(:,1+shift:N) zeros(M,shift)];
   end
else
   if nargin == 4 & type ~= 0
      A(row,:) = [A(row,1+shift:N) A(row,1:shift)];
   else
      A(row,:) = [A(row,1+shift:N) zeros(1,shift)];
   end
end
y = A;


text:
function [result, prediction] = svr_test_linear_shift(test_labels, test_samples, model)

    prediction = test_samples * model.w(1:end-1)' + model.w(end);
%     prediction = predict(test_labels, test_samples, model);

    prediction(~model.success) = 0;

    if(model.cutoff >= 0)
        % perform shifting here per person
        users = unique(model.vid_ids);

        for i=1:numel(users)

            preds_user = prediction(strcmp(model.vid_ids, users(i)));
            sorted = sort(preds_user);

            % alternative, move to histograms and pick the highest one

            shift = sorted(round(end*model.cutoff)+1);

            prediction(strcmp(model.vid_ids, users(i))) = preds_user - shift;

        end
    end

    % Cap the prediction as well
    prediction(prediction<0)=0;
    prediction(prediction>5)=5;

    % using the average of RMS errors
%     result = mean(sqrt(mean((prediction - test_labels).^2)));
    if(~isfield(model, 'eval_ids'))
        result = corr(test_labels, prediction);
        [ ~, ~, ~, ccc, ~, ~ ] = evaluate_regression_results( prediction, test_labels );
        result = ccc;
    else
        eval_ids = unique(model.eval_ids)';
        ccc = 0;
        fprintf('CCC: ');
        for i=eval_ids
            [ ~, ~, ~, ccc_curr, ~, ~ ] = evaluate_regression_results( prediction(model.eval_ids == i), test_labels(model.eval_ids == i) );
            ccc = ccc + ccc_curr;
            fprintf('%.3f ', ccc_curr);
        end
        ccc = ccc / numel(eval_ids);
        fprintf('mean : %.3f\n', ccc);
        result = ccc;
    end

    if(isnan(result))
        result = 0;
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%
for i=1:size(X,1)
    tempMin=inf;
    for j=1:K
        temp =sum((centroids(j,:)-X(i,:)).^2);
        if temp<tempMin
           idx(i)=j;
           tempMin = temp;
        end
    end

end

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

end

author: zzlyw

repo: machine-learning-exercises