text:
function [segmented] = ft_volumesegment(cfg, mri)

% FT_VOLUMESEGMENT segments an anatomical MRI. The behavior depends on the output requested. It can
% return probabilistic tissue maps of gray/white/csf compartments, a skull-stripped anatomy, or
% binary masks representing the brain surface, skull, or scalp surface.
%
% Use as
%   segmented = ft_volumesegment(cfg, mri)
% where the input mri should be a single anatomical volume that was for example read with
% FT_READ_MRI. For the purpose of creating binary masks of the brain or of the skull, you can also
% provide either the anatomical volume or the already segmented volume (with the probabilistic
% tissue maps) as input.
%
% The configuration structure can contain
%   cfg.output         = string or cell-array of strings, see below (default = 'tpm')
%   cfg.spmversion     = string, 'spm2', 'spm8', 'spm12' (default = 'spm12')
%   cfg.spmmethod      = string, 'old', 'new', 'mars', the algorithm used when
%                        cfg.spmversion='spm12' (default = 'old')
%   cfg.opts           = structure  with spm-version specific options. See the
%                        code and/or the SPM-documentation for more detail.
%   cfg.template       = filename of the template anatomical MRI (default =
%                        '/spm2/templates/T1.mnc' or '/spm8/templates/T1.nii')
%   cfg.tpm            = cell-array containing the filenames of the tissue probability maps
%   cfg.name           = string for output filename
%   cfg.write          = 'no' or 'yes' (default = 'no'), writes the probabilistic tissue maps
%                        to SPM compatible analyze (spm2), or nifti (spm8 or spm12) files,
%                        with the following suffix for spm2
%                         _seg1, for the gray matter segmentation
%                         _seg2, for the white matter segmentation
%                         _seg3, for the csf segmentation
%                        or with the following prefix for spm8 and spm12 with spmmethod='old'
%                         c1, for the gray matter segmentation
%                         c2, for the white matter segmentation
%                         c3, for the csf segmentation
%                        and with spm12 with spmmethod='new' there will be 3 additional tissue types
%                         c4, for the bone segmentation
%                         c5, for the soft tissue segmentation
%                         c6, for the air segmentation
%                        When using spm12 with spmmethod='mars', the tpms will be postprocessed
%                        with the mars toolbox, yielding smoother segmentations in general.
%   cfg.brainsmooth    = 'no', or scalar, the FWHM of the gaussian kernel in voxels, (default = 5)
%   cfg.scalpsmooth    = 'no', or scalar, the FWHM of the gaussian kernel in voxels, (default = 5)
%   cfg.skullsmooth    = 'no', or scalar, the FWHM of the gaussian kernel in voxels, (default = 5)
%                        this parameter is only used when the segmentation contains 6 tisuse types,
% %                      including 'bone'
%   cfg.brainthreshold = 'no', or scalar, relative threshold value which is used to threshold the
%                        tpm in order to create a volumetric brainmask (see below), (default = 0.5)
%   cfg.scalpthreshold = 'no', or scalar, relative threshold value which is used to threshold the
%                        anatomical data in order to create a volumetric scalpmask (see below),
%                        (default = 0.1)
%   cfg.skullthreshold = 'no', or scalar, relative threshold value which is used to threshold the
%                        anatomical data in order to create a volumetric scalpmask (see below),
%                        (default = 0.5). this parameter is only used when the segmentation
%                        contains 6 tissue types, including 'bone'
%   cfg.downsample     = integer, amount of downsampling before segmentation (default = 1, which
%                        means no downsampling)
%
% The desired segmentation output is specified with cfg.output as a string or cell-array of strings
% and can contain
%   'tpm'         - tissue probability map for csf, white and gray matter
%   'brain'       - binary representation of the brain (the combination of csf, white and gray matter)
%   'skull'       - binary representation of the skull
%   'scalp'       - binary representation of the scalp
%   'skullstrip'  - anatomy with only the brain
%
% Example use:
%   cfg        = [];
%   segmented  = ft_volumesegment(cfg, mri) will segmented the anatomy and will output the
%                segmentation result as 3 probabilistic masks in gray, white and csf.
%
%   cfg        = [];
%   cfg.output = 'skullstrip';
%   segmented  = ft_volumesegment(cfg, mri) will generate a skull-stripped anatomy based on a
%                brainmask generated from the probabilistic tissue maps. The skull-stripped anatomy
%                is stored in the field segmented.anatomy.
%
%   cfg        = [];
%   cfg.output = {'brain' 'scalp' 'skull'};
%   segmented  = ft_volumesegment(cfg, mri) will produce a volume with 3 binary masks, representing
%                the brain surface, scalp surface, and skull which do not overlap.
%
%   cfg        = [];
%   cfg.output = {'scalp'};
%   segmented  = ft_volumesegment(cfg, mri) will produce a volume with a binary mask (based on the
%                anatomy), representing the border of the scalp surface (i.e., everything inside the
%                surface is also included). Such representation of the scalp is produced faster,
%                because it doesn't require to create the tissue probabilty maps prior to creating
%                the mask.
%
% It is not possible to request tissue-probability maps (tpm) in combination with binary masks
% (brain, scalp or skull) or with a skull-stripped anatomy. The output will return only the probabilistic
% maps in gray, white and csf. However, when a segmentation with the probabilistic gray, white
% and csf representations is available, it is possible to use it as input to create the brain or skull
% binary mask. For example:
%   cfg           = [];
%   cfg.output    = {'tpm'};
%   segment_tpm   = ft_volumesegment(cfg, mri);
%   cfg.output    = {'brain'};
%   segment_brain = ft_volumesegment(cfg, segment_tpm);
%
% For the SPM-based segmentation to work, the coordinate frame of the input MRI needs to be
% approximately coregistered to the templates of the probabilistic tissue maps. The templates
% are defined in SPM/MNI-space. FieldTrip attempts to do an automatic alignment based on the
% coordsys-field in the MRI, and if this is not present, based on the coordsys-field in the cfg.
% If none of them is specified the FT_DETERMINE_COORDSYS function is used to interactively
% assess the coordinate system in which the MRI is expressed.
%
% The template MRI is defined in SPM/MNI-coordinates, see also http://bit.ly/2sw7eC4
%   x-axis pointing to the right ear
%   y-axis along the acpc-line
%   z-axis pointing to the top of the head
%   origin in the anterior commissure.
% Note that the segmentation requires the template MRI to be in SPM coordinates.
%
% 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_READ_MRI, FT_DETERMINE_COORDSYS, FT_PREPARE_HEADMODEL

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

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

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

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

% this is not supported any more as of 26/10/2011
if ischar(mri)
  ft_error('please use cfg.inputfile instead of specifying the input variable as a string');
end

% check if the input data is valid for this function
mri = ft_checkdata(mri, 'datatype', 'volume', 'feedback', 'yes', 'hasunit', 'yes', 'hascoordsys', 'yes');

% ensure that old and unsupported options are not being relied on by the end-user's script
% instead of specifying cfg.coordsys, the user should specify the coordsys in the data
cfg = ft_checkconfig(cfg, 'forbidden', {'units', 'coordsys', 'inputcoord', 'inputcoordsys', 'coordinates'});
cfg = ft_checkconfig(cfg, 'renamedval', {'output', 'skin', 'scalp'});

% set the defaults
cfg.output           = ft_getopt(cfg, 'output',         'tpm');
cfg.downsample       = ft_getopt(cfg, 'downsample',     1);
cfg.spmversion       = ft_getopt(cfg, 'spmversion',     'spm12');
cfg.write            = ft_getopt(cfg, 'write',          'no');
cfg.spmmethod        = ft_getopt(cfg, 'spmmethod',      'old'); % doing old-style in case of spm12

% set default for smooth and threshold
cfg.brainsmooth      = ft_getopt(cfg, 'brainsmooth',    ''); % see also below
cfg.scalpsmooth      = ft_getopt(cfg, 'scalpsmooth',    ''); % see also below
cfg.skullsmooth      = ft_getopt(cfg, 'skullsmooth',    ''); % see also below
cfg.brainthreshold   = ft_getopt(cfg, 'brainthreshold', ''); % see also below
cfg.scalpthreshold   = ft_getopt(cfg, 'scalpthreshold', ''); % see also below

% earlier version of smooth and threshold specification
cfg.smooth           = ft_getopt(cfg, 'smooth',      '');
cfg.threshold        = ft_getopt(cfg, 'threshold',   '');

% chech whether earlier version of smooth and threshold was specified
if ~(isempty(cfg.smooth))
  if isempty(cfg.brainsmooth)
    cfg.brainsmooth = cfg.smooth;
    ft_warning('Smoothing can be specified separately for scalp and brain. User-specified smoothing will be applied for brainmask.')
  end
  if isempty(cfg.scalpsmooth)
    cfg.scalpsmooth = cfg.smooth;
    ft_warning('Smoothing can be specified separately for scalp and brain. User-specified smoothing will be applied for scalpmask.')
  end
end
if ~(isempty(cfg.threshold))
  if isempty(cfg.brainthreshold)
    cfg.brainthreshold = cfg.threshold;
    ft_warning('Threshold can be specified separately for scalp and brain. User-specified threshold will be applied for brainmask.')
  end
  if isempty(cfg.scalpthreshold)
    cfg.scalpthreshold = cfg.threshold;
    ft_warning('Threshold can be specified separately for scalp and brain. User-specified threshold will be applied for scalpmask.')
  end
end
% then set defaults again
cfg.brainsmooth      = ft_getopt(cfg, 'brainsmooth',      5);
cfg.scalpsmooth      = ft_getopt(cfg, 'scalpsmooth',      5);
cfg.skullsmooth      = ft_getopt(cfg, 'skullsmooth',      5);
cfg.brainthreshold   = ft_getopt(cfg, 'brainthreshold',   0.5);
cfg.scalpthreshold   = ft_getopt(cfg, 'scalpthreshold',   0.1);
cfg.skullthreshold   = ft_getopt(cfg, 'skullthreshold',   0.5);

% check that the preferred SPM version is on the path
ft_hastoolbox(cfg.spmversion, 1);

if ~isfield(cfg, 'name')
  if ~strcmp(cfg.write, 'yes')
    tmp = tempname;
    cfg.name = tmp;
  else
    ft_error('you must specify the output filename in cfg.name');
  end
end
[pathstr, name] = fileparts(cfg.name);
cfg.name = fullfile(pathstr, name); % remove any possible file extension, to be added later

if ~iscell(cfg.output)
  % ensure it to be cell, to allow for multiple outputs
  cfg.output = {cfg.output};
end

% check whether SPM is needed to generate tissue probability maps
if numel(cfg.output) == 1 && strcmp('scalp', cfg.output)
  needtpm = 0; % not needed for (cummulative type) scalpmask
else
  needtpm = any(ismember(cfg.output, {'tpm' 'gray' 'white' 'csf' 'brain' 'skull' 'skullstrip'}));
end

tissue  = isfield(mri, {'gray', 'white', 'csf', 'bone', 'softtissue', 'air'});
ntissue = sum(tissue);
if ntissue==6
  % this is new-style segmentation
  hastpm = ~islogical(mri.gray) && ~islogical(mri.white) && ~islogical(mri.csf);  % tpm should be probabilistic and not binary!
elseif all(tissue(1:3)==true)
  hastpm = ~islogical(mri.gray) && ~islogical(mri.white) && ~islogical(mri.csf);  % tpm should be probabilistic and not binary!
elseif any(tissue)
  ft_warning('the input seems to contain tissue probability maps, but is incomplete for this function');
  hastpm = false;
else
  hastpm = false;
end

if needtpm && ~hastpm
  % spm needs to be used for the creation of the tissue probability maps
  dotpm = 1;
else
  dotpm = 0;
end

if cfg.downsample~=1
  % optionally downsample the anatomical and/or functional volumes
  tmpcfg = keepfields(cfg, {'downsample', 'spmversion', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.smooth = 'no'; % smoothing is done in ft_volumesegment itself
  mri = ft_volumedownsample(tmpcfg, mri);
  % restore the provenance information
  [cfg, mri] = rollback_provenance(cfg, mri);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create the tissue probability maps if needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if dotpm

  if isdeployed
    % ensure that these exist in this case, otherwise deal with the defaults below
    cfg = ft_checkconfig(cfg, 'required', {'template' 'tpm'});
  end

  % remember the original transformation matrix coordinate system
  original = [];
  original.transform = mri.transform;
  original.coordsys  = mri.coordsys;
  if isfield(mri, 'unit')
    original.unit = mri.unit;
  else
    mri = ft_determine_units(mri); % guess the unit field if not present
    original.unit = mri.unit;
  end
  mri = ft_convert_units(mri, 'mm');
  if isdeployed
    mri = ft_convert_coordsys(mri, 'acpc', 2, cfg.template);
  else
    mri = ft_convert_coordsys(mri, 'acpc');
  end

  % flip and permute the 3D volume itself, so that the voxel and
  % headcoordinates approximately correspond this improves the convergence
  % of the segmentation algorithm
  [mri, permutevec, flipflags] = align_ijk2xyz(mri);

  switch lower(cfg.spmversion)
    case 'spm2'
      cfg.template = ft_getopt(cfg, 'template', fullfile(spm('Dir'), 'templates', 'T1.mnc'));
      opts          = ft_getopt(cfg, 'opts');
      opts.estimate = ft_getopt(opts, 'estimate');
      opts.write    = ft_getopt(opts, 'write');
      opts.estimate.affreg = ft_getopt(opts.estimate, 'affreg');

      VF = ft_write_mri([cfg.name, '.img'], mri.anatomy, 'transform', mri.transform, 'spmversion', cfg.spmversion, 'dataformat', 'analyze_img');

      % set the spm segmentation defaults (from /opt/spm2/spm_defaults.m script)
      opts.estimate.priors = ft_getopt(opts.estimate, 'priors', ...
        char(fullfile(spm('Dir'), 'apriori', 'gray.mnc'),...
        fullfile(spm('Dir'), 'apriori', 'white.mnc'),...
        fullfile(spm('Dir'), 'apriori', 'csf.mnc')));
      opts.estimate.reg    = ft_getopt(opts.estimate, 'reg',    0.01);
      opts.estimate.cutoff = ft_getopt(opts.estimate, 'cutoff', 30);
      opts.estimate.samp   = ft_getopt(opts.estimate, 'samp',   3);
      opts.estimate.bb     = ft_getopt(opts.estimate, 'bb',     [[-88 88]' [-122 86]' [-60 95]']);
      opts.estimate.affreg.smosrc  = ft_getopt(opts.estimate.affreg, 'smosrc',  8);
      opts.estimate.affreg.regtype = ft_getopt(opts.estimate.affreg, 'regtype', 'mni');
      opts.estimate.affreg.weight  = ft_getopt(opts.estimate.affreg, 'weight', '');
      opts.write.cleanup   = ft_getopt(opts.write, 'cleanup', 1);
      opts.write.wrt_cor   = ft_getopt(opts.write, 'wrt_cor', 1);

      % perform the segmentation
      fprintf('performing the segmentation on the specified volume\n');
      spm_segment(VF, cfg.template, opts);

      % generate the list of filenames that contains the segmented volumes
      filenames = {[cfg.name, '_seg1.img'];...
        [cfg.name, '_seg2.img'];...
        [cfg.name, '_seg3.img']};

    case 'spm8'
      cfg.tpm = ft_getopt(cfg, 'tpm');
      cfg.tpm = char(cfg.tpm(:));
      if isempty(cfg.tpm)
        cfg.tpm = char(fullfile(spm('Dir'),'tpm','grey.nii'),...
          fullfile(spm('Dir'),'tpm','white.nii'),...
          fullfile(spm('Dir'),'tpm','csf.nii'));
      end
      px.tpm = cfg.tpm;

      VF = ft_write_mri([cfg.name, '.img'], mri.anatomy, 'transform', mri.transform, 'spmversion', cfg.spmversion, 'dataformat', 'nifti_spm');

      fprintf('performing the segmentation on the specified volume\n');
      p         = spm_preproc(VF, px);
      [po, dum] = spm_prep2sn(p);

      % this writes a mat file, may be needed for Dartel, not sure yet
      save([cfg.name '_sn.mat'],'-struct','po');

      % These settings were taken from a batch
      opts     = ft_getopt(cfg, 'opts');
      opts.GM  = [0 0 1];
      opts.WM  = [0 0 1];
      opts.CSF = [0 0 1];
      opts.biascor = 1;
      opts.cleanup = 0;

      % write the segmented volumes, -> this can probably be done differently
      spm_preproc_write(po, opts);

      % generate the list of filenames that contains the segmented volumes
      [pathstr, name] = fileparts(cfg.name);
      filenames = {fullfile(pathstr,['c1', name, '.img']);...
        fullfile(pathstr,['c2', name, '.img']);...
        fullfile(pathstr,['c3', name, '.img'])};

    case 'spm12'
      addpath(fullfile(spm('Dir'),'toolbox/OldSeg'));
      if strcmp(cfg.spmmethod, 'old')
        cfg.tpm = ft_getopt(cfg, 'tpm');
        cfg.tpm = char(cfg.tpm(:));
        if isempty(cfg.tpm)
          cfg.tpm = char(fullfile(spm('Dir'),'toolbox/OldSeg','grey.nii'),...
            fullfile(spm('Dir'),'toolbox/OldSeg','white.nii'),...
            fullfile(spm('Dir'),'toolbox/OldSeg','csf.nii'));
        end
        px.tpm = cfg.tpm;

        VF = ft_write_mri([cfg.name, '.nii'], mri.anatomy, 'transform', mri.transform, 'unit', mri.unit, 'spmversion', cfg.spmversion, 'dataformat', 'nifti_spm');

        fprintf('performing the segmentation on the specified volume, using the old-style segmentation\n');
        p         = spm_preproc(VF, px);
        [po, dum] = spm_prep2sn(p);

        % this write a mat file, may be needed for Dartel, not sure yet
        save([cfg.name '_sn.mat'],'-struct','po');

        % These settings are taken from a batch
        opts     = [];
        opts.GM  = [0 0 1];
        opts.WM  = [0 0 1];
        opts.CSF = [0 0 1];
        opts.biascor = 1;
        opts.cleanup = 0;

        % write the segmented volumes, -> this can be done differently
        spm_preproc_write(po, opts);

        % generate the list of filenames that contains the segmented volumes
        [pathstr, name] = fileparts(cfg.name);
        filenames = {fullfile(pathstr,['c1', name, '.nii']);...
          fullfile(pathstr,['c2', name, '.nii']);...
          fullfile(pathstr,['c3', name, '.nii'])};

      elseif strcmp(cfg.spmmethod, 'new') || strcmp(cfg.spmmethod, 'mars')
        cfg.tpm = ft_getopt(cfg, 'tpm');
        cfg.tpm = char(cfg.tpm);
        if ~isfield(cfg, 'tpm') || isempty(cfg.tpm)
          cfg.tpm = fullfile(spm('dir'),'tpm','TPM.nii');
        end

        VF = ft_write_mri([cfg.name, '.nii'], mri.anatomy, 'transform', mri.transform, 'spmversion', cfg.spmversion, 'dataformat', 'nifti_spm');

        fprintf('performing the segmentation on the specified volume, using the new-style segmentation\n');

        % create the structure that is required for spm_preproc8
        opts          = ft_getopt(cfg, 'opts');
        opts.image    = VF;
        opts.tpm      = ft_getopt(opts, 'tpm', spm_load_priors8(cfg.tpm));
        opts.biasreg  = ft_getopt(opts, 'biasreg',  0.0001);
        opts.biasfwhm = ft_getopt(opts, 'biasfwhm', 60);
        opts.lkp      = ft_getopt(opts, 'lkp',      [1 1 2 2 3 3 4 4 4 5 5 5 5 6 6 ]);
        opts.reg      = ft_getopt(opts, 'reg',      [0 0.001 0.5 0.05 0.2]);
        opts.samp     = ft_getopt(opts, 'samp',     3);
        opts.fwhm     = ft_getopt(opts, 'fwhm',     1);

        Affine = spm_maff8(opts.image(1),3,32,opts.tpm,eye(4),'mni');
        Affine = spm_maff8(opts.image(1),3, 1,opts.tpm,Affine,'mni');
        opts.Affine = Affine;

        % run the segmentation
        p = spm_preproc8(opts);

        % this writes the 'native' segmentations
        if strcmp(cfg.spmmethod, 'new')
          spm_preproc_write8(p, [ones(6,2) zeros(6,2)], [0 0], [0 1], 1, 1, nan(2,3), nan);
        elseif strcmp(cfg.spmmethod, 'mars')
          ft_hastoolbox('mars', 1);
          if ~isfield(cfg, 'mars'), cfg.mars = []; end
          beta        = ft_getopt(cfg.mars, 'beta', 0.1);
          convergence = ft_getopt(cfg.mars, 'convergence', 0.1);
          tcm{1}      = fullfile(fileparts(which('spm_mars_mrf')), 'rTCM_BW20_S1.mat');
          p = spm_mars_mrf(p, [ones(6,2) zeros(6,2)], [0 0], [0 1], tcm, beta, convergence, 1);
        end

        % this writes a mat file, may be needed for Dartel, not sure yet
        save([cfg.name '_seg8.mat'],'-struct','p');

        [pathstr, name] = fileparts(cfg.name);
        filenames = {fullfile(pathstr,['c1', name, '.nii']);...
          fullfile(pathstr,['c2', name, '.nii']);...
          fullfile(pathstr,['c3', name, '.nii']);...
          fullfile(pathstr,['c4', name, '.nii']);...
          fullfile(pathstr,['c5', name, '.nii']);...
          fullfile(pathstr,['c6', name, '.nii'])};

      else
        ft_error('cfg.spmmethod should be either ''old'', ''new'' or ''mars''');
      end % if spmmethod is old, new or mars

    otherwise
      ft_error('unsupported SPM version');
  end % switch spmversion

  for k = 1:numel(filenames)
    Vtmp = spm_vol(filenames{k});
    dat  = spm_read_vols(Vtmp);
    Vtmp.dat = dat;
    V(k)     = Vtmp;
  end

  if strcmp(cfg.write, 'no')
    [pathstr, name] = fileparts(cfg.name);
    prefix = {'c1';'c2';'c3';'c3';'c4';'c5';'c6';'m';'y_';''};
    suffix = {'_seg1.hdr';'_seg2.hdr';'_seg3.hdr';'_seg1.img';'_seg2.img';'_seg3.img';'_seg1.mat';'_seg2.mat';'_seg3.mat';'.hdr';'.img';'.nii'};
    for k = 1:numel(prefix)
      for j = 1:numel(suffix)
        warning off; % delete is quite verbose in its warnings
        try, delete(fullfile(pathstr,[prefix{k}, name, suffix{j}])); end
        warning on;
      end
    end
  elseif strcmp(cfg.write, 'yes')
    for k= 1:numel(V)
      % I am not sure whether or why this is needed
      V(k).mat = mri.transform;
      V(k)     = spm_create_vol(V(k));
    end
  end

  % collect the results
  segmented.dim       = reshape(size(V(1).dat), 1, []); % enforce a row vector
  segmented.transform = original.transform; % use the original transform
  segmented.coordsys  = original.coordsys;  % use the original coordsys
  segmented.unit      = original.unit;      % use the original units
  segmented.gray      = V(1).dat;
  if length(V)>1, segmented.white      = V(2).dat; end
  if length(V)>2, segmented.csf        = V(3).dat; end
  if length(V)>3, segmented.bone       = V(4).dat; end
  if length(V)>4, segmented.softtissue = V(5).dat; end
  if length(V)>5, segmented.air        = V(6).dat; end
  segmented.anatomy = mri.anatomy;

  % flip the volumes back according to the changes introduced by align_ijk2xyz
  fn = fieldnames(segmented);
  fn = fn(ismember(fn, {'anatomy', 'gray', 'white', 'csf', 'bone', 'softtissue', 'air'}));
  for k = 1:3
    if flipflags(k)
      for j = 1:numel(fn)
        segmented.(fn{j}) = flipdim(segmented.(fn{j}), k);
      end
    end
  end

  if ~all(permutevec == [1 2 3])
    for j = 1:numel(fn)
      segmented.(fn{j}) = ipermute(segmented.(fn{j}), permutevec);
    end
    segmented.dim  = size(segmented.(fn{1}));
  end

else
  % rename the output data
  segmented = mri;
  clear mri
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now the data contains the tissue probability maps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create the requested output fields

remove = {'anatomy' 'csf' 'gray' 'white' 'bone' 'softtissue' 'air' 'skull' 'skullstrip' 'brain' 'scalp'}; % all possible tissues

% check if smoothing or thresholding is required

dosmooth_brain = ~strcmp(cfg.brainsmooth, 'no');
if dosmooth_brain && ischar(cfg.brainsmooth)
  ft_error('invalid value %s for cfg.brainsmooth', cfg.brainsmooth);
end
dosmooth_scalp = ~strcmp(cfg.scalpsmooth, 'no');
if dosmooth_scalp && ischar(cfg.scalpsmooth)
  ft_error('invalid value %s for cfg.scalpsmooth', cfg.scalpsmooth);
end
dosmooth_skull = ~strcmp(cfg.skullsmooth, 'no');
if dosmooth_skull && ischar(cfg.skullsmooth)
  ft_error('invalid value %s for cfg.skullsmooth', cfg.skullsmooth);
end

dothres_brain = ~strcmp(cfg.brainthreshold, 'no');
if dothres_brain && ischar(cfg.brainthreshold)
  ft_error('invalid value %s for cfg.brainthreshold', cfg.brainthreshold);
end
dothres_scalp = ~strcmp(cfg.scalpthreshold, 'no');
if dothres_scalp && ischar(cfg.scalpthreshold)
  ft_error('invalid value %s for cfg.scalpthreshold', cfg.scalpthreshold);
end
dothres_skull = ~strcmp(cfg.skullthreshold, 'no');
if dothres_skull && ischar(cfg.skullthreshold)
  ft_error('invalid value %s for cfg.skullthreshold', cfg.skullthreshold);
end

outp = cfg.output;

if ~isempty(intersect(outp, 'tpm'))
  % output: probability tissue maps
  remove = intersect(remove, {'anatomy'});
elseif  ~isempty(intersect(outp, {'white' 'gray' 'csf' 'brain' 'skull' 'scalp' 'skullstrip'}))

  createoutputs = true;
  while createoutputs
    % create scalpmask - no tpm or brainmask is required to create it
    if any(strcmp('scalp', outp))

      fprintf('creating scalpmask ... ');
      % let the softtissue mask take precedence
      if isfield(segmented, 'softtissue')
        fprintf('using the softtissue tpm for segmentation\n');
        fname   = 'softtissue';
        anatomy = segmented.softtissue;
      elseif isfield(segmented, 'anatomy')
        fprintf('using the anatomy field for segmentation\n');
        fname   = 'anatomy';
        anatomy = segmented.anatomy;
      else
        ft_error('no appropriate volume is present for the creation of a scalpmask');
      end
      if dosmooth_scalp
        anatomy = volumesmooth(anatomy, cfg.scalpsmooth, fname);
      else
        fprintf('no smoothing applied on %s for scalp segmentation\n', fname);
      end
      if dothres_scalp
        anatomy = volumethreshold(anatomy, cfg.scalpthreshold, fname);
      else
        fprintf('no threshold applied on %s for scalp segmentation\n', fname);
      end

      % fill the slices along each dimension (because using a single one is
      % just arbitrary, and behavior depends on how the voxeldata is in the
      % volume.
      a1 = volumefillholes(anatomy, 1);
      a2 = volumefillholes(anatomy, 2);
      a3 = volumefillholes(anatomy, 3);

      scalpmask = a1 | a2 | a3;
      clear anatomy a1 a2 a3

      % threshold again to remove little parts outside of head
      scalpmask = volumethreshold(scalpmask);

      % output: scalp (cummulative) (if this is the only requested output)
      if numel(outp)==1
        segmented.scalp = scalpmask;
        remove(strcmp(remove, 'scalp')) = []; % keep this
        break
      end
    end % if scalp

    % create the brain from the tpm
    fprintf('creating brainmask ... using the sum of gray, white and csf tpms\n');
    brain = segmented.gray + segmented.white + segmented.csf;
    if dosmooth_brain
      brain = volumesmooth(brain,  cfg.brainsmooth, 'brainmask');
    else
      fprintf('no smoothing applied on brainmask\n')
    end
    if dothres_brain
      brain = volumethreshold(brain, cfg.brainthreshold, 'brainmask');
    else
      fprintf('no threshold applied on brainmask\n')
    end

    % output: skullstrip
    if any(strcmp('skullstrip', outp))
      if ~isfield(segmented, 'anatomy'), ft_error('no anatomy field present'); end
      fprintf('creating skullstripped anatomy ...');
      brain_ss = cast(brain, class(segmented.anatomy));
      segmented.anatomy = segmented.anatomy.*brain_ss;
      clear brain_ss
      remove(strcmp(remove, 'skullstrip')) = []; % keep this
      remove(strcmp(remove, 'anatomy'))    = []; % keep this
      if numel(outp)==1
        break
      end
    end % if skullstrip

    % make binary mask from brain
    brainmask = brain>0;
    clear brain

    if any(strcmp(outp, 'brain'))
      % output: brain
      segmented.brain = brainmask;
      remove(strcmp(remove, 'brain')) = []; % keep this
      if numel(outp)==1
        break
      end

    elseif any(strcmp(outp, 'gray')) || any(strcmp(outp, 'white')) || any(strcmp(outp, 'csf'))
      % output: gray, white, csf
      [dum, tissuetype] = max(cat(4, segmented.csf, segmented.gray, segmented.white), [], 4);
      clear dum
      if any(strcmp(outp, 'white'))
        segmented.white = (tissuetype == 3) & brainmask;
        remove(strcmp(remove, 'white')) = []; % keep this
      end
      if any(strcmp(outp, 'gray'))
        segmented.gray = (tissuetype == 2) & brainmask;
        remove(strcmp(remove, 'gray')) = []; % keep this
      end
      if any(strcmp(outp, 'csf'))
        segmented.csf = (tissuetype == 1) & brainmask;
        remove(strcmp(remove, 'csf')) = []; % keep this
      end

    end % if brain or gray/while/csf

    if any(strcmp('skull', outp)) || any(strcmp('scalp', outp))
      % create the skull from the brain mask
      fprintf('creating skullmask ... ');
      if ~isfield(segmented, 'bone')
        fprintf('using the brainmask\n');
        braindil  = imdilate(brainmask>0, strel_bol(6)); % FIXME this should be user configurable
        skullmask = braindil & ~brainmask;
        clear braindil
      elseif isfield(segmented, 'bone')
        fprintf('using the bone tpm for segmentation\n');
        skull = segmented.bone + segmented.gray + segmented.white + segmented.csf;
        if dosmooth_skull
          skull = volumesmooth(skull, cfg.skullsmooth, 'skull');
        else
          fprintf('no smoothing applied on skull tpm for skull segmentation\n');
        end
        if dothres_skull
          skull = volumethreshold(skull, cfg.skullthreshold, 'skull');
        else
          fprintf('no threshold applied on skull tom for skull segmentation\n')
        end

        a1 = volumefillholes(skull, 1);
        a2 = volumefillholes(skull, 2);
        a3 = volumefillholes(skull, 3);

        skullmask = a1 | a2 | a3;
        skullmask = volumethreshold(skullmask);
        skullmask = skullmask & ~brainmask;

        clear a1 a2 a3
      end % if segmentation without bone

      if any(strcmp(outp, 'skull'))
        segmented.skull = skullmask;
        remove(strcmp(remove, 'skull')) = []; % keep this
        if numel(outp)==1
          break
        end
      end

      % output: scalp (exclusive type)
      if numel(outp) > 1 && any(strcmp('scalp', outp))
        scalpmask(brainmask>0)=0;
        clear brainmask
        scalpmask(skullmask>0)=0;
        clear skullmask
        segmented.scalp = scalpmask;
        remove(strcmp(remove, 'scalp')) = []; % keep this
        clear scalpmask
      end
    end

    createoutputs = false; % exit the while loop
  end % while createoutputs

else
  ft_error('unknown output %s requested\n', cfg.output{:});
end

% remove unnecessary fields
segmented = removefields(segmented, remove);

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   mri
ft_postamble provenance segmented
ft_postamble history    segmented
ft_postamble savevar    segmented

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function structureSet = getStructureSetAssociatedScan(scanSet, planC)
% DK
% LM: APA, 10/22/08: Fixed to output correct associated structure index
%
% copyright (c) 2001-2008, Washington University in St. Louis.
%
% Permission is granted to use or modify only for non-commercial,
% non-treatment-decision applications, and further only if this header is
% not removed from any file. No warranty is expressed or implied for any
% use whatever: use at your own risk.  Users can request use of CERR for
% institutional review board-approved protocols.  Commercial users can
% request a license.  Contact Joe Deasy for more information
% (radonc.wustl.edu@jdeasy, reversed).

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

structureSet = [];
structureSet = sparse(structureSet);

[assocScanUID{1:length(planC{indexS.structureArray})}] = deal(planC{indexS.structureArray}.assocScanUID);

if isempty(assocScanUID)
    structureSet = [];
    return
end

for i = 1:length(scanSet)

    scanUID = planC{indexS.scan}(scanSet(i)).scanUID;
    try
        structureSet(i) = find(strcmpi(scanUID,assocScanUID));
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function r8vec_bracket3_test ( )

%*****************************************************************************80
%
%% R8VEC_BRACKET3_TEST tests R8VEC_BRACKET3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  ntest = 6;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_BRACKET3_TEST\n' );
  fprintf ( 1, '  R8VEC_BRACKET3 finds a pair of entries in a\n' );
  fprintf ( 1, '  sorted real array which bracket a value.\n' );

  xtest(1:6) = [ -10.0, 1.0, 4.5, 5.0, 10.0, 12.0 ];

  x = r8vec_indicator1 ( n );
  x(6) = x(5);

  r8vec_print ( n, x, '  Sorted array:' );

  left = floor ( ( n + 1 ) / 2 );

  for itest = 1 : ntest

    xval = xtest(itest);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Search for XVAL = %f\n', xval );

    fprintf ( 1, '  Starting guess for interval is = %d\n', left );

    left = r8vec_bracket3 ( n, x, xval, left );

    fprintf ( 1, '  Nearest interval:' );
    fprintf ( 1, '    X[%d]= %f\n', left, x(left) );
    fprintf ( 1, '    X[%d]= %f\n', left+1, x(left+1) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function CreateAffinePanel(handles)

    %Affine Panel
    handles.affinePanel = uipanel('Parent',handles.toolPanel,'Units','normalize','FontSize',10,'Title','Registration Parameters',...
                                    'Tag','affinePanel','Clipping','off',...
                                    'Position',[0 0 1 1],...
                                    'Visible','on');
    dy = 160;
    h312 = uicontrol(...
    'Parent',handles.affinePanel,...
    'Units','pixel',...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[9 124+dy 117 19],...
    'String','Translation Scale:',...
    'Style','text',...
    'Tag','text62');

    handles.affine_tscale = uicontrol(...
    'Parent',handles.affinePanel,...
    'Units','pixel',...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[129 122+dy 78 23],...
    'String','0.001',...
    'Style','edit',...
    'Tag','affine_tscale');

    handles.affine_minstep = uicontrol(...
    'Parent',handles.affinePanel,...
    'Units','pixel',...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[129 86+dy 78 24],...
    'String','0.01',...
    'Style','edit',...
    'Tag','affine_minstep');

    handles.affine_maxstep = uicontrol(...
    'Parent',handles.affinePanel,...
    'Units','pixel',...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[129 52+dy 78 23],...
    'String','1.0',...
    'Style','edit',...
    'Tag','affine_maxstep');

    h316 = uicontrol(...
    'Parent',handles.affinePanel,...
    'Units','pixel',...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[8 92+dy 118 17],...
    'String','Minimum Step:',...
    'Style','text',...
    'Tag','text65');

    h317 = uicontrol(...
    'Parent',handles.affinePanel,...
    'Units','pixel',...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[8 58+dy 118 17],...
    'String','Maximum Step:',...
    'Style','text',...
    'Tag','text66');

    h318 = uicontrol(...
    'Parent',handles.affinePanel,...
    'Units','pixel',...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[8 24+dy 118 17],...
    'String','Iteration Number :',...
    'Style','text',...
    'Tag','text67');

    handles.affine_iternum = uicontrol(...
    'Parent',handles.affinePanel,...
    'Units','pixel',...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[129 19+dy 78 23],...
    'String','200',...
    'Style','edit',...
    'Tag','affine_iternum');

    guidata(handles.mainframe, handles);

end


text:
function newcat = excludequarries(orgcatalog,quarries)
    %this function searches with ismember every quarry defined with 'quarries'
    %(typical the variable q) in the orgcatalog (typical a) and removes
    %them in the original catalog.

    foundquars=ismember(orgcatalog,quarries,'rows');

    newcat=orgcatalog(~foundquars,:);

end

author: CelsoReyes


text:
function C = mrdivide(W,B)
   %MRDIVIDE (/) Slash or right matrix divide for WAVEFORMS.

   % VERSION: 1.0 of waveform objects
   % AUTHOR: Celso Reyes (celso@gi.alaska.edu)
   % LASTUPDATE: 3/15/2009

   if ~(isnumeric(B))
      error('Waveform:mrdivide:invalidNumeratorClass',...
         'Cannot divide by a %s, must be numeric', class(B));
   end

   C = set(W, 'data', get(W,'data') / double(B));
end

author: geoscience-community-codes

repo: GISMO


text:
function gray = rgb_to_gray ( rgb, equal )

%*****************************************************************************80
%
%% RGB_TO_GRAY returns a grayscale version of an RGB image.
%
%  Discussion:
%
%    An RGB image is an (M,N,3) array.
%
%    A grayscale image is an (M,N) array.
%
%    A grayscale version of an RGB image can be made by averaging
%    the RGB components of each pixel.
%
%    A more sophisticated approach uses the luminance function:
%
%      GRAY = 0.2126 * R + 0.7152 * G + 0.0722 * B
%
%    which attempts to better model the contributions to brightness
%    of the different colors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, uint8 RGB(M,N,3), the original data.
%
%    Input, logical EQUAL, is TRUE (1) if the R, G and B channels
%    are to be assigned equal weight.  EQUAL is FALSE by default.
%
%    Output, uint8 GRAY(M,N), the grayscale version of the data.
%
  if ( nargin < 1 )
    error ( 'RGB_TO_GRAY - Fatal error! Missing RGB input argument.' )
  end

  if ( nargin < 2 )
    equal = 0;
  end

  if ( equal )
    v = [ 1, 1, 1 ]' / 3;
  else
    v = [ 0.2126, 0.7152, 0.0722 ]';
  end
%
%  This is really a matrix multiply, but the obvious equation
%    gray = rgb * v
%  is illegal according to MATLAB.
%
  gray = uint8 ( double ( rgb(:,:,1) ) * v(1) ...
               + double ( rgb(:,:,2) ) * v(2) ...
               + double ( rgb(:,:,3) ) * v(3) );

  return
end


text:
function[nodes, categories, heights] = parentsToTrees(parents)
% [nodes, categories, heights] = parentsToTrees(parents)
%
% Converts an [s x 2] dimensional cell of strings into a tree where the
% first column indicates nodes and the second column indicates their parent
% nodes.
%
% The output can be plotted using Matlab's treeplot function:
% treeplot(nodes');
% [xs, ys] = treelayout(nodes);
%
% Copyright by Holger Caesar, 2017

% Extract categories and make sure they are unique
categories = parents(:, 1);
assert(numel(categories) == numel(unique(categories)));

% Create pointers to parent nodes
categoryCount = size(categories, 1);
nodes = nan(categoryCount, 1);
heights = nan(categoryCount, 1);
nodes(1) = 0;
heights(1) = 0;
for i = 2 : categoryCount
    childNode = find(strcmp(parents(:, 1), categories{i}));
    if isempty(childNode)
        error('Error: No parent node found for %s!', categories{i});
    end
    parentNode = find(strcmp(categories, parents(childNode, 2)));
    assert(numel(parentNode) == 1, 'Error: Node %s has %d parent labels! Should be 1.', parents{childNode, 2}, numel(parentNode));
    nodes(i) = parentNode;
    heights(i) = heights(parentNode) + 1;
end

author: nightrome


text:
% Get vertex coordinates for the patch type 'T'
% [C, NW, NE, SW, SE] = HealpixGetPatchVertexCoordsT(n, i, j, INFO)
%
% Parameters
% n : grid resolution
% i : ring index
% j : intra-ring index
% INFO : intermediate information (output of HealpixSelectPatchClass())
% C : intra-patch coordinates
% NW : coordinates for north-east vertex
% NE : coordinates for north-west vertex
% SW, SE : coordinates for south vertex (always SW = SE since the patch is triangular shape)

function [C, NW, NE, SW, SE] = HealpixGetPatchVertexCoordsT(n, i, j, INFO)

int_j = fix(j);
if int_j == 0
    int_j = 4 * n;
end
east_int_j = mod(int_j + 1 - 1, 4 * INFO.int_i_n) + 1;

if INFO.is_south_pole == 0
    NW = [n, int_j];
    NE = [n, east_int_j];
    SW = [n + 1, int_j];
    SE = SW;

    decimal_i = INFO.decimal_i_n;
else
    SW = [3 * n, int_j];
    SE = [3 * n, east_int_j];
    NW = [3 * n - 1, int_j];
    NE = NW;

    decimal_i = 1 - INFO.decimal_i_n;
end

decimal_j = INFO.decimal_polar_intra_part;

if 1 > INFO.decimal_i_n
    C = [decimal_i, decimal_j / (1 - INFO.decimal_i_n)] - 0.5;
else
    C = [decimal_i, decimal_j] - 0.5;
end

author: mtex-toolbox

repo: mtex

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


text:
function plotacq(chan)
%PLOTACQ Plots a channel structure generated from loadacq
%  (c) 2003 Mihai Moldovan
%   M.Moldovan@mfi.ku.dk

fig_hndl=figure(...
            'Units','normalized',...
             'Color',get(0,'defaultUicontrolBackgroundColor'),...
            'Name','Acqknowledge import',...
            'NumberTitle','off',...
            'Resize','on',...
            'WindowStyle','normal',...
            'PaperUnits','normalized',...
            'PaperOrientation','Landscape',...
            'PaperPositionMode','manual',...
            'PaperPosition',[0 0 1 1],...
            'PaperType','A4');

% --------------------------------------------------------------------
% plot all channels
% --------------------------------------------------------------------

n=length(chan);

for i=1:n
    ax_hndl=subplot(n,1,i);
    set(ax_hndl,'FontUnits','points','FontName','Arial','FontSize',7,'NextPlot','add','LineWidth',1);

    y=chan(i).data;
    x=(1:length(y));
    x=x.*chan(i).ms/1000;
    plot (x,y,'Color',chan(i).color);

    t=ylabel([chan(i).name ' (' chan(i).units ')']);
    set (t,'Rotation',0,'HorizontalAlignment','right','FontWeight','bold');

    if i~=n
        set(ax_hndl,'XTickLabel',{})
    else
        xlabel ('Time (seconds)')
    end

    mk=chan(i).mdata;
    for j=1:length(mk)
        x=mk(j).*chan(i).ms/1000;
        plot ([x x],ylim,'k--')
    end
end

% --------------------------------------------------------------------
% save the metafile
% --------------------------------------------------------------------

   [ST,I] = dbstack;
   n=length(ST);

   mfile=ST(n).name;

   [file_path mname emxt mversn] = fileparts(mfile);
   meta=[file_path filesep mname '.emf'];

   %print -dwinc for color
   print ( '-dmeta', meta);

%return
%-----------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
%% t_tseriesLoadFromInplane
%
% Illustrates how to load a time series from a functional data set.
%
% Tested 01/05/2011 - MATLAB r2008a, Fedora 12, Current Repos
%
% Stanford VISTA
%

%% Initialize the key variables and data path
% Data directory (where the mrSession file is located)
dataDir = mrtInstallSampleData('functional','vwfaLoc');

%% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

% There can be several data types - name the one you want to plot
dataType = 'MotionComp';

% Which scan number from that data type?
scan = 1;

% Would you like the raw time series?
isRawTSeries = true;

%% Get data structure:
vw = initHiddenInplane(); % Foregoes interface - loads data silently

%% Set data structure properties:
vw = viewSet(vw, 'CurrentDataType', dataType); % Data type

%% Get time series from ROI:
% Format returned is rows x cols x slices x time
tSeries = tSeries4D(vw, scan, [], 'usedefaults', ~isRawTSeries);

%% Show movie of a single slice across the given scan
figure;
colormap autumn;
nSlices = size(tSeries, 3);
nTimePoints = size(tSeries, 4);
rg = [-1 1] * max(abs(tSeries(:)));
for i = 1:nTimePoints
    imagesc(tSeries(:, :, ceil(nSlices/2), i), rg);
    axis image; colormap gray
    title(sprintf('Volume %d', i))
    pause(.1);
end
close;

%% Show movie of slices at a single time point
figure;
colormap autumn;
for i = 1:nSlices
    imagesc(tSeries(:, :, i, ceil(nTimePoints/2)));
    colormap gray
    axis image;
    title(sprintf('Slice %d', i))
    pause(.1);
end
close;

%% Restore original directory
cd(curDir);

%% END

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [data] = ft_denoise_hfc(cfg,data)

% FT_DENOISE_HFC implements harmonic field correction, which models
% interference as a harmonic magnetic field. It is particulaly useful
% for MEG data with low channel numbers (e.g. OPM data). Homogenous field based
% on Tierney et al. (2021) NIMG, 118484. Harmonic expansion based on Tierney
% et al. (2022) NIMG, 119338.
%
% Use as:
%   data = ft_denoise_hfc(cfg,data)
%
% Where cfg is a configuration structure that contains:
%   cfg.channel         = channels for HFC (default = 'all')
% 	cfg.order           = spherical harmonic order:
%                           order = 1 is a homogenous field; order = 2
%                           includes gradients; order = 3 includes
%                           quadratic terms etc. (default = 1)
%   cfg.trials          = which trials do you want to denoise?
%                           (default = 'all')
%   cfg.updatesens      = do you want to update sensor info with projector?
%                           (default = 'yes')
%   cfg.feedback        = do you want feedback (default = 'yes')
%   cfg.residualcheck   = do you want to check channel residuals
%                          (default = 'yes')
%   cfg.residualthresh  = (in pT) what level of residual signal is fine for
%                           quality assurance (default = 50)
% See also FT_DENOISE_SYNTHETIC, FT_DENOISE_PCA, FT_DENOISE_DSSP, FT_DENOISE_TSP

% Copyright (C) 2021-22, Tim Tierney, George O'Neill, Robert Seymour
% Wellcome Centre for Human Neuroimaging, UCL
%
% 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
  % do not continue function execution in case the outputfile is present and the user indicated to keep it
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'trial'}); % prevent accidental typos, see issue 1729

% check the input data
data = ft_checkdata(data, 'datatype', {'raw'});

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729

% Check if the data has a grad structure
if ~isfield(data,'grad')
  error(['Data needs a grad structure']);
end

% set the defaults
cfg.order           = ft_getopt(cfg, 'order', 1);
cfg.channel         = ft_getopt(cfg, 'channel', 'all', 1);
cfg.trials          = ft_getopt(cfg, 'trials', 'all', 1);
cfg.updatesens      = ft_getopt(cfg, 'updatesens', 'yes');
cfg.feedback        = ft_getopt(cfg, 'feedback', 'yes');
cfg.residualcheck   = ft_getopt(cfg, 'residualcheck', 'yes');
cfg.residualthresh  = ft_getopt(cfg, 'residualthresh', 50);

% select data based on cfg.channel and cfg.trials
tmpcfg = keepfields(cfg, {'trials', 'showcallinfo'});

% Select the correct channels using ft_channelselection
tmpcfg.channel = cfg.channel;
data   = ft_selectdata(tmpcfg, data);

% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

% Check match between input data chans and grad chans
[x, dummy] = ismember(data.grad.label,data.label);

% Warn the user if there are chans in grad not in data
num_mismatch = sum(x(:) == 0);

if num_mismatch > 0 && istrue(cfg.feedback)
  ft_warning(['Found ' num2str(num_mismatch) ' channels in grad structure'...
    ' not present in cfg.channel. These channels will NOT be used for '...
    'Harmonic Field Correction']);
end

% Check for Dr. Tim Tierney's OPM toolbox on the path, and add if needed
% See: https://github.com/tierneytim/OPM
ft_hastoolbox('opm', 1);

% generate harmonic basis set
opt = [];
opt.li = cfg.order;
opt.v = data.grad.coilpos(x,:);
opt.o = data.grad.coilori(x,:);
N = spm_opm_vslm(opt);

% Make montage for the next step
montage     = [];
montage.tra = eye(length(N)) - N*pinv(N);
montage.labelold = data.grad.label(x);
montage.labelnew = data.grad.label(x);
montage.chantypeold = data.grad.chantype(x);
montage.chantypenew = data.grad.chantype(x);
montage.chanunitold = data.grad.chanunit(x);
montage.chanunitnew = data.grad.chanunit(x);

labelold = data.label;

% apply montage;
% Tell the user
data = ft_apply_montage(data,montage,'keepunused','yes');
if istrue(cfg.feedback)
  disp('Applied HFC to the data');
end

% Update the tra to account for the g. Essential to correct the lead
% fields going forward.
if istrue(cfg.updatesens)
  data.grad = ft_apply_montage(data.grad, montage, 'keepunused',...
    'yes', 'balancename', 'hfc','warning',false);
  if istrue(cfg.feedback)
    disp('Converted the sensor description to HFC');
  end
end

% reorder the channels to stay close to the original ordering
[dummy, selnew] = match_str(montage.labelold, data.label);
if numel(selnew)==numel(labelold)
  for i=1:numel(data.trial)
    data.trial{i} = data.trial{i}(selnew,:);
  end
  data.label = data.label(selnew);
else
  ft_warning('channel ordering might have changed');
end

% Perform running variance check to identify odd channels
if strcmp(cfg.residualcheck,'yes')
  residual_check(cfg.residualthresh,data,montage.labelold)
end

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

function residual_check(residualthresh,data,oldlabels)
% Script to determine residual variance post HFC

% find corrected channels in the output data
[dummy, selnew2] = match_str(oldlabels, data.label);

trvar = [];
for ii = 1:numel(data.trial)
  tmp = data.trial{ii}(selnew2,:);
  Mk =  tmp(:,1);
  Sk = zeros(size(Mk));
  count = 1;
  for jj = 1:size(tmp,2)
    Xk = tmp(:,jj);
    Mkprev = Mk;
    Mk = Mkprev +(Xk-Mkprev)/count;
    Sk=Sk+(Xk-Mkprev).*(Xk-Mk) ;
    count=count+1;
  end
  trvar(:,ii)=Sk/(count-1);
end

% Identify the most common chanunit
chanunit = ft_chanunit(data);
[s,dummy,j]=unique(chanunit);
chanunit = s{mode(j)};

switch chanunit % some of this are silly, but safety first!
  case 'fT'
    scale = 1e-3;
  case 'pT'
    scale = 1;
  case 'nT'
    scale = 1e3;
  case 'uT'
    scale = 1e6;
  case 'mT'
    scale = 1e9;
  case {'T','T/m'}
    scale = 1e12;
  otherwise
    ft_error('Cannot check residuals due to unknown sensor units!')
end

SD = mean(sqrt(trvar),2)*scale;

fprintf('Checking for unsual channels post-corrections\n')
count = 0;
for ii = 1:length(SD)
  index = selnew2(ii);
  if SD(ii) > residualthresh
    count = count + 1;
    if strcmp(chanunit,'fT')
      fprintf(['Residual on channel ' num2str(index) ', '...
        data.label{index} ': %3.2f pT\n'], SD(ii));
    else
      fprintf(['Residual on channel ' num2str(index) ', '...
        data.label{index} ': %3.2f' chanunit '\n'], SD(ii));
    end
  end
end
if ~count
  fprintf('No unusual channel residuals found!\n')
end


text:
%% check Clebsch Gordan Tensor

%% the reference

% some arbitrary rotation
g = rotation.byEuler(-72*degree,88*degree,134*degree);

% we want to express the product of two wigner D functions
D1 = WignerD(g,'order',1);
D1D1_ref = D1(:) * D1(:).'

%% expansion into Wigner functions of lower order

% zero order component
D0 = WignerD(g,'order',0);
CG0 = ClebschGordanTensor(0);
C0 = D0*EinsteinSum(CG0,[1 3],CG0,[2 4])

% first order component
D1 = WignerD(g,'order',1);
CG1 = ClebschGordanTensor(1);
C1 = EinsteinSum(CG1,[1 3 -1],D1,[-1 -2],CG1,[2 4 -2])

% second order component
D2 = WignerD(g,'order',2);
CG2 = ClebschGordanTensor(2);
C2 = EinsteinSum(CG2,[1 3 -1],D2,[-1 -2],CG2,[2 4 -2])

a = reshape(matrix(C0 + C1 + C2),[9,9]) ./ D;
A = round(a)

imagesc(A)

assert(norm(A.* D - reshape(matrix(C0 + C1 + C2),[9,9]))<=1e-10,'Clebsch Gordan check failed')


text:
%addpath /white/u2/bob/matlab/stats/
volView = getSelectedVolume;

%baseDir = '/home/bob/research/papers/dti/callosal_areaMaps/matlabFigs';
mkdir(tempdir, 'matlabFigs');
baseDir = fullfile(tempdir, 'matlabFigs');
disp(['figures will be saved in ' baseDir]);
[junk,sessionCode] = fileparts(pwd);
fileNameBase = fullfile(baseDir, [mrSESSION.subject '_' sessionCode]);
fileName = [fileNameBase '_' volView.ROIs(viewGet(volView, 'selectedROI')).name];
labelX = 'Left phase (rad)';
labelY = 'Right phase (rad)';
% axisLim = [-2*pi, 2*pi, -2*pi, 2*pi];
% tick = [-3*pi/2, -pi, -pi/2, 0, pi/2, pi, 3*pi/2];
% tickLabel = ['   ';'-pi';'   ';' 0 ';'   ';' pi';'   '];
eccenDeg = [0, pi/2, pi, 3*pi/2, 4*pi/2]/(2*pi) * 16;
axisLim = [-pi/2, 5*pi/2, -pi/2, 5*pi/2];
tick = [0, pi/2, pi, 3*pi/2, 4*pi/2];
tickLabel = [' 0  ';'    ';' pi ';'    ';'2 pi'];
scanNames = {'wedge','ring'};
coThresh = 0.25;
paperFigs = 0;
if(paperFigs)
    fontName = 'Helvetica';
    fontSize = 12;
    markSize = 9;
    markColor = 'k';
    lineWidth = 2;
    lineColor = 'k';
else
    fontName = 'Comic Sans MS';
    fontSize = 16;
    markSize = 12;
    markColor = 'r';
    lineWidth = 3;
    lineColor = 'b';
end
generalNotes = '';
roiCoords = getCurROIcoords(volView);
% We get the last coord, which is the left-right axis (lower = left)
% The following will swap coordinate pairs so that the left-most coord
% is always first.
fiberEndptCoordX = [roiCoords(3,1:2:end); roiCoords(3,2:2:end)];
swapThese = diff(fiberEndptCoordX)<0;
tmp = fiberEndptCoordX(1,swapThese);
fiberEndptCoordX(1,swapThese) = fiberEndptCoordX(2,swapThese);
fiberEndptCoordX(2,swapThese) = tmp;
for(scanNum=1:2)
    co = getCurDataROI(volView,'co',scanNum,roiCoords);
    ph = getCurDataROI(volView,'ph',scanNum,roiCoords);
    fiberEndptPh = [ph(1:2:end); ph(2:2:end)];
    tmp = fiberEndptPh(1,swapThese);
    fiberEndptPh(1,swapThese) = fiberEndptPh(2,swapThese);
    fiberEndptPh(2,swapThese) = tmp;
    goodVals = isfinite(fiberEndptPh(1,:)) & isfinite(fiberEndptPh(2,:)) & co(1:2:end)>coThresh & co(2:2:end)>coThresh;
    fiberEndptPh = fiberEndptPh(:,goodVals);
    % Null hypothesis test
    fiberEndptPh(2,:) = fiberEndptPh(2,randperm(length(fiberEndptPh(2,:))));
    % Unwrap phases
    if(strcmpi(scanNames{scanNum}, 'wedge'))
    else
        labelX = 'Left Eccentricity (deg)';
        labelY = 'Right Eccentricity (deg)';
        tickLabel = num2str(eccenDeg');
        d = diff(fiberEndptPh);
        unwrapThese = d>pi & fiberEndptPh(2,:)>3*pi/2;
        fiberEndptPh(1,unwrapThese) = fiberEndptPh(2,unwrapThese)+2*pi;
        unwrapThese = d<-pi & fiberEndptPh(1,:)>3*pi/2;
        fiberEndptPh(2,unwrapThese) = fiberEndptPh(2,unwrapThese)+2*pi;
    end
    x = squeeze(fiberEndptPh(1,:));
    y = squeeze(fiberEndptPh(2,:));
    %[p, r, df] = statTest(x, y, 'r'); rsq = r^2;
    %fit = polyfit(x, y, 1)
    [b,bint,r,rint,s]=regress(y',[x',ones(size(y'))]);
    fitX = [min(x),max(x)];
    fitY = fitX.*b(1) + b(2);
    p = s(3); rsq = s(1); f = s(2);
    if(paperFigs)
        if(scanNum==1) fh = figure; else figure(fh); end
        set(fh,'Position',[6, 33, 400, 820]);
        subplot(2,1,scanNum);
    else
        fh = figure;
        set(fh,'Position',[7, 50, 400, 430]);
    end
    axes;
    line([axisLim(1),axisLim(2)], [axisLim(3),axisLim(4)], 'Color', 'k', 'LineWidth', 1);
    hold on;
    h = scatter(x, y, markColor);
    hold off;
    set(gca, 'fontName', fontName);
    set(gca, 'fontSize', fontSize);
    xlabel(labelX, 'fontName', fontName, 'fontSize', fontSize, 'fontWeight','normal');
    ylabel(labelY, 'fontName', fontName, 'fontSize', fontSize, 'fontWeight','normal');
    axis(axisLim);
    set(gca, 'XTick', tick);
    set(gca, 'YTick', tick);
    set(gca, 'XTickLabel', tickLabel);
    set(gca, 'YTickLabel', tickLabel);
    set(h, 'MarkerSize', markSize); set(h, 'LineWidth', lineWidth);
    lineH = line(fitX, fitY, 'Color', lineColor, 'LineWidth', lineWidth);
    if(p<0.0001)
        title(sprintf('%s: R^2=%0.4f (p<0.0001)', scanNames{scanNum}, rsq));
    else
        title(sprintf('%s: R^2=%0.4f (p=%0.4f)', scanNames{scanNum}, rsq, p));
    end
    if(~paperFigs)
        pos = get(fh, 'Position');
        set(gca, 'Units', 'pixels');
        %[left bottom width height]
        set(gca, 'Position', [70,70,pos(3)-90,pos(3)-90]);
        set(fh, 'PaperPositionMode', 'auto');
        %print(fh, '-dpng', '-r120', [fileName '_' scanNames{scanNum} '_scatter.png']);
    end

    notes = [sprintf('coThresh=%0.3f; ', coThresh) generalNotes];
    %save([fileName '_' scanNames{scanNum} '_dat.mat'], 'fiberEndptPh', 'co', 'ph', 'roiCoords', 'notes');
end

if(paperFigs)
    set(fh, 'PaperPositionMode', 'auto');
    print(fh, '-depsc', '-tiff', [fileName '_scatter.eps']);
end
refresh;


text:
function [ x, seed ] = beta_binomial_sample ( a, b, c, seed )

%*****************************************************************************80
%
%% BETA_BINOMIAL_SAMPLE samples the Beta Binomial CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, parameters of the PDF.
%    0.0D+00 < A,
%    0.0D+00 < B.
%
%    Input, integer C, a parameter of the PDF.
%    0 <= C.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = beta_binomial_cdf_inv ( cdf, a, b, c );

  return
end


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function run_cases(inputs_all, lambda_vals, test_name, n_proc)
% RUN_CASES: Functions runs the discrete optimization algorithms
% on all test cases specified in inputs_all with the lambda values
% prescribed in lambda_vals vector. The results are saved in the
% results/test_name folder.

% Find number of tests
n_test = length(inputs_all);

%% Parallel Problem setup

parpool(n_proc);
parfor t=1:n_test

    % Set test inputs struct
    inputs_t = inputs_all{t};

    % Declare cells to store results
    rnd  = cell(length(lambda_vals));
    sa   = cell(length(lambda_vals));
    bo   = cell(length(lambda_vals));
    ols  = cell(length(lambda_vals));
    smc  = cell(length(lambda_vals));
    smac = cell(length(lambda_vals));

    bayes = struct;
    bayes.stSA1 = cell(length(lambda_vals));
    bayes.stSA2 = cell(length(lambda_vals));
    bayes.stSA3 = cell(length(lambda_vals));
    bayes.sdp   = cell(length(lambda_vals));

    mle = struct;
    mle.stSA1   = cell(length(lambda_vals));
    mle.stSA2   = cell(length(lambda_vals));
    mle.stSA3   = cell(length(lambda_vals));
    mle.sdp     = cell(length(lambda_vals));

    hs = struct;
    hs.stSA1    = cell(length(lambda_vals));
    hs.stSA2    = cell(length(lambda_vals));
    hs.stSA3    = cell(length(lambda_vals));
    hs.sdp      = cell(length(lambda_vals));

    %% Run optimization
    for l=1:length(lambda_vals)

        % Define objective function with penalty term
        inputs_t.lambda = lambda_vals(l);
        penalty   = @(x) inputs_t.lambda*inputs_t.reg_term(x);
        objective = @(x) inputs_t.model(x) + penalty(x);

        fprintf('--------------------------------------------\n')
        fprintf('Test = %d/%d, Lambda = %f\n\n', t, n_test, lambda_vals(l));

        % Run different ML optimization algorithms
        rnd{l} = random_samp(objective, inputs_t);
        fprintf('Random - Runtime: %f\n', sum(rnd{l}.runTime));

        sa{l}  = simulated_annealing(objective, inputs_t);
        fprintf('SA - Runtime = %f\n', sum(sa{l}.runTime));

        bo{l}  = bayes_opt(objective, inputs_t);
        fprintf('BO - Runtime = %f\n', sum(bo{l}.runTime));

        ols{l} = local_search(objective, inputs_t);
        fprintf('OLS - Runtime = %f\n', sum(ols{l}.runTime));

        smc{l} = binary_smc(objective, inputs_t);
        fprintf('SMC - Runtime = %f\n', sum(smc{l}.runTime));

        smac{l} = run_smac(objective, inputs_t);
        fprintf('SMAC - Runtime = %f\n', sum(smac{l}.runTime));

        % Run BOCS with Bayesian model
        inputs_t.estimator = 'bayes';

        bayes.stSA1{l} = BOCS(inputs_t.model, penalty, inputs_t, 1, 'SA');
        fprintf('Bayes.SA1 - Runtime = %f\n', sum(bayes.stSA1{l}.runTime));

        bayes.stSA2{l} = BOCS(inputs_t.model, penalty, inputs_t, 2, 'SA');
        fprintf('Bayes.SA2 - Runtime = %f\n', sum(bayes.stSA2{l}.runTime));

        bayes.stSA3{l} = BOCS(inputs_t.model, penalty, inputs_t, 3, 'SA');
        fprintf('Bayes.SA3 - Runtime = %f\n', sum(bayes.stSA3{l}.runTime));

        bayes.sdp{l} = BOCS(inputs_t.model, penalty, inputs_t, 2, 'sdp');
        fprintf('Bayes.SDP - Runtime = %f\n', sum(bayes.sdp{l}.runTime));

        % Run BOCS with MLE model
        inputs_t.estimator = 'mle';

        mle.stSA1{l} = BOCS(inputs_t.model, penalty, inputs_t, 1, 'SA');
        fprintf('MLE.SA1 - Runtime = %f\n', sum(mle.stSA1{l}.runTime));

        mle.stSA2{l} = BOCS(inputs_t.model, penalty, inputs_t, 2, 'SA');
        fprintf('MLE.SA2 - Runtime = %f\n', sum(mle.stSA2{l}.runTime));

        mle.stSA3{l} = BOCS(inputs_t.model, penalty, inputs_t, 3, 'SA');
        fprintf('MLE.SA3 - Runtime = %f\n', sum(mle.stSA3{l}.runTime));

        mle.sdp{l} = BOCS(inputs_t.model, penalty, inputs_t, 2, 'sdp');
        fprintf('MLE.SDP - Runtime = %f\n', sum(mle.sdp{l}.runTime));

        % Run BOCS with Horseshoe model
        inputs_t.estimator = 'horseshoe';

        hs.stSA1{l} = BOCS(inputs_t.model, penalty, inputs_t, 1, 'SA');
        fprintf('HS.SA1 - Runtime = %f\n', sum(hs.stSA1{l}.runTime));

        hs.stSA2{l} = BOCS(inputs_t.model, penalty, inputs_t, 2, 'SA');
        fprintf('HS.SA2 - Runtime = %f\n', sum(hs.stSA2{l}.runTime));

        hs.stSA3{l} = BOCS(inputs_t.model, penalty, inputs_t, 3, 'SA');
        fprintf('HS.SA3 - Runtime = %f\n', sum(hs.stSA3{l}.runTime));

        hs.sdp{l} = BOCS(inputs_t.model, penalty, inputs_t, 2, 'sdp');
        fprintf('HS.SDP - Runtime = %f\n', sum(hs.sdp{l}.runTime));

        % Save results
        iSave(sprintf(['../results/' test_name '/test%d.mat'], t), ...
            rnd, sa, bo, ols, smc, smac, bayes, mle, hs, inputs_t);

    end
end

delete(gcp('nocreate'))

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
%Leave One Out Cross Validation with Nearest Neighbor Classifier
%>
%> @param FeatureMatrix: features (dimension iNumFeatures x iNumObservations)
%> @param ClassIdx: vector with class indices (length iNumObservations, starting from 0)
%>
%> @retval Acc overall accuracy after Cross-Validation
% ======================================================================
function [Acc, conf_mat] = ToolLooCrossVal(FeatureMatrix, ClassIdx)

    % initialize
    TP = 0;

    conf_mat = zeros(length(unique(ClassIdx)));

    % loop over observations
    for o = 1:size(FeatureMatrix, 2)
        % remove current observation from 'training set'
        v_train = [FeatureMatrix(:, 1:o-1) FeatureMatrix(:, o+1:end)]';
        C_train = [ClassIdx(1:o-1) ClassIdx(:, o+1:end)]';

        % compute result of Nearest Neighbor Classifier given the traindata
        res = ToolSimpleKnn(FeatureMatrix(:, o)', v_train, C_train, 1);

        conf_mat(ClassIdx(o)+1, res+1) = conf_mat(ClassIdx(o)+1, res+1) + 1;

        % if result is correct increment number of true positives
        if (res == ClassIdx(o))
            TP = TP+1;
        end
    end

    % compute overall (micro) accuracy
    Acc = TP / length(ClassIdx);
end

author: alexanderlerch

repo: ACA-Code


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym piecewise (@var{cond1}, @var{val1}, @var{cond2}, @var{val2}, @dots{})
%% @defmethodx @@sym piecewise (@var{cond1}, @var{val1}, @var{cond2}, @var{val2}, @dots{}, @var{else_val})
%% Construct piecewise function.
%%
%% The returned piecewise function evaluates to @var{val1} if @var{cond1}
%% holds, @var{val2} if @var{cond2} holds, @dots{} etc.  In the case where none
%% of the conditions hold, it evaluates to @var{else_val} if provided.  If
%% @var{else_val} is not provided, it evaluates to @code{nan}.
%%
%% Examples:
%% @example
%% @group
%% syms x real
%% f = piecewise (abs (x) < 1, exp (- 1 / (1 - x^2)), abs (x) >= 1, 0)
%%   @result{} f = (sym)
%%       ⎧  -1
%%       ⎪ ──────
%%       ⎪      2
%%       ⎨ 1 - x
%%       ⎪ℯ        for │x│ < 1
%%       ⎪
%%       ⎩   0      otherwise
%% @end group
%% @end example
%%
%% For this piecewise function, we can omit the redundant condition at the end:
%% @example
%% @group
%% syms x real
%% f = piecewise (abs (x) < 1, exp (- 1 / (1 - x^2)), 0)
%%   @result{} f = (sym)
%%       ⎧  -1
%%       ⎪ ──────
%%       ⎪      2
%%       ⎨ 1 - x
%%       ⎪ℯ        for │x│ < 1
%%       ⎪
%%       ⎩   0      otherwise
%% @end group
%% @end example
%%
%% @seealso{if}
%% @end defmethod

function f = piecewise (varargin)
  if nargin < 1
    print_usage ();
  end

  cmd = {'def pack_args_without_else(args):'
         '    return [(val, cond) for cond, val in unflatten(args, 2)]'
         'def pack_args(args):'
         '    if len(args) % 2 == 0:'
         '        return pack_args_without_else(args)'
         '    else:'
         '        return pack_args_without_else(args[:-1]) + [(args[-1], True)]'
         'return Piecewise(*pack_args(_ins))'
        };

  args = cellfun (@sym, varargin, 'UniformOutput', false);
  f = pycall_sympy__ (cmd, args{:});
end

%!test
%! % basic
%! syms x real
%! f = piecewise (abs (x) < 1, 1);
%! assert (isnan (subs (f, -1)));
%! assert (isequal (subs (f, 0), 1));
%! assert (isnan (subs (f, 1)));

%!test
%! % heaviside
%! syms x real
%! f = rewrite (heaviside (x, 1 / sym (2)), 'Piecewise');
%! g = piecewise (x < 0, 0, x == 0, 1 / sym (2), x > 0, 1);
%! assert (logical (simplify (f == g)));

%% FIXME: expand test suite, add SMT compat tests, ...

repo: octsympy


text:
function [results, loc] = ismember(A,B)
   %ISMEMBER for scnlobjects
   % TRUE for each waveform that matches any scnl in the anythingelse array.
   %
   % see also ISMEMBER
   %
   % WILDCARD is not implemented
   %
   % since we are looking for scnlA in scnlB, then scnlB may have wildcards.
   % the wildcard is a field value of '*' (only).  These haven't been
   % thoroughly tested.

   wildcard = '*';
   if isempty(A) || isempty(B)
      results = [];
      loc = [];
      return
   end

   assert(isa(A,'scnlobject'),class(A))
   tagA = toTag(A);
   tagB = toTag(B);
   [results, loc] = ismember(tagA, tagB);
   return

   results = false(size(A));
   loc = zeros(size(A));
   for i=1:numel(B)
      theseA = A;
      if strcmpi(B(i).station,wildcard)
         theseA = set(theseA,'station',wildcard);
      end
      if strcmpi(B(i).channel,wildcard)
         theseA = set(theseA,'channel',wildcard);
      end
      if strcmpi(B(i).network,wildcard)
         theseA = set(theseA,'network',wildcard);
      end
      if strcmpi(B(i).location,wildcard)
         theseA = set(theseA,'location',wildcard);
      end
      [resultstmp, loctmp] = ismember(getSearchString(theseA),getSearchString(B));
      results = results | resultstmp;
      loc(loc==0) = loctmp(loc==0);
   end
end
function t = toTag(scnls)

   switch class(scnls)
      case 'scnlobject'
         t = [scnls.tag];
      case 'ChannelTag'
         t = [scnls];
      otherwise
         %try to convert it?
         t = ChannelTag(scnls);
   end
end

function ss = getSearchString(scnl)
   sta = get(scnl,'station');
   if ~iscell(sta),
      sta = {sta};
   end;
   sta(strcmpi(sta,'')) = {'*-noStation-*'};
   cha = get(scnl,'channel');
   if ~iscell(cha), cha = {cha}; end;
   cha(strcmpi(cha,'')) = {'*-noChannel-*'};
   net = get(scnl,'network');
   if ~iscell(net), net = {net}; end;
   net(strcmpi(net,'') | strcmpi(net,'--')) = {'*-noNetwork-*'};
   loc = get(scnl,'location');
   if ~iscell(loc), loc = {loc}; end;
   loc(strcmpi(loc,'') | strcmpi(loc,'--')) = {'*-noLocation-*'};

   ss = strcat(sta,'|',cha,'|', net,'|',loc);
end


text:
function [ gradientParams, SigmaInvs, CholDecomps, Sigmas ] = gradientCCRFFull( params, lambda_a, lambda_b, PrecalcBs, x, y, Precalc_yBys, PrecalcBsFlat)
%GRADIENTPRF Summary of this function goes here
%   Detailed explanation goes here

    nExamples = numel(x);

    numBetas = size(PrecalcBsFlat{1},2);
    numAlphas = numel(params) - numBetas;

    alphasInit = params(1:numAlphas);
    betasInit = params(numAlphas+1:end);
    gradientParams = zeros(size(params));

    % These might be use to calculate the LogLikelihood, don't want to
    % recompute them
    SigmaInvs = cell(nExamples, 1);
    CholDecomps = cell(nExamples, 1);
    Sigmas = cell(nExamples, 1);
    gradients = zeros(nExamples, numel(params));
    for q = 1 : nExamples

        yq = y{q};
        xq = x{q};

        PrecalcB = PrecalcBs{q};
        PrecalcB_flat = PrecalcBsFlat{q};

        [ logGradientsAlphas, logGradientsBetas, SigmaInv, CholDecomp, Sigma ] = gradientCCRF_withoutReg(alphasInit, betasInit, PrecalcB, xq, yq, Precalc_yBys(q, :), PrecalcB_flat);
        SigmaInvs{q} = SigmaInv;
        CholDecomps{q} = CholDecomp;
        Sigmas{q} = Sigma;

        gradients(q,:) = [logGradientsAlphas; logGradientsBetas];
    end
    gradientParams = sum(gradients,1)';
    regAlpha = alphasInit * lambda_a;
    regBeta = betasInit * lambda_b;
    gradientParams = gradientParams - [regAlpha; regBeta];
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function marginal = marginal_nodes(engine, query, add_ev)
% MARGINAL_NODES Compute the marginal on the specified query nodes (cond_gauss)
% marginal = marginal_nodes(engine, query, add_ev)
%
% 'query' must be a singleton set
% add_ev is an optional argument; if 1, we will "inflate" the marginal of observed nodes
% to their original size, adding 0s to the positions which contradict the evidence

if nargin < 3, add_ev = 0; end

if length(query) ~= 1
  error('cond_gauss_inf_engine can only handle marginal queries on single nodes')
end
j = query;
bnet = bnet_from_engine(engine);

if myismember(j, bnet.cnodes)
  if ~myismember(j, engine.onodes)
    [m, C] = collapse_mog(engine.mu{j}, engine.Sigma{j}, engine.T);
    marginal.mu = m;
    marginal.Sigma = C;
    marginal.T = 1.0; % single mixture component
  else
    marginal.mu = engine.evidence{j};
    k = bnet.node_sizes(j);
    marginal.Sigma = zeros(k,k);
    marginal.T = 1.0; % since P(E|E)=1
  end
else
  marginal = pot_to_marginal(marginalize_pot(engine.joint_dmarginal, j));
  if add_ev
    marginal = add_ev_to_dmarginal(marginal, engine.evidence, bnet.node_sizes);
  end
end

marginal.domain = query;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function f=chi2pdf(x,k)

% f=chi2pdf(x,k)
%
% Return the pdf of a chi-squared distribution with k degrees of freedom
% evaluated at x.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% 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 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

error(nargchk(2,2,nargin));
error(chk_param(x,'x','numeric','real'));
error(chk_param(k,'k','scalar','real','integer','>=',1));

f=NaN(size(x));
f(x<0)=0;
xgt0=x>=0;
f(xgt0)=x(xgt0).^(k/2-1).*exp(-x(xgt0)/2)/(2^(k/2)*gamma(k/2));

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function [x, f, eflag, output, lambda] = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
%MIQPS_CPLEX  Mixed Integer Quadratic Program Solver based on CPLEX.
%   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
%       MIQPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
%   [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_CPLEX(PROBLEM)
%   A wrapper function providing a standardized interface for using
%   CPLEXQP or CPLEXLP to solve the following QP (quadratic programming)
%   problem:
%
%       min 1/2 X'*H*X + C'*X
%        X
%
%   subject to
%
%       L <= A*X <= U       (linear constraints)
%       XMIN <= X <= XMAX   (variable bounds)
%
%   Inputs (all optional except H, C, A and L):
%       H : matrix (possibly sparse) of quadratic cost coefficients
%       C : vector of linear cost coefficients
%       A, L, U : define the optional linear constraints. Default
%           values for the elements of L and U are -Inf and Inf,
%           respectively.
%       XMIN, XMAX : optional lower and upper bounds on the
%           X variables, defaults are -Inf and Inf, respectively.
%       X0 : optional starting value of optimization vector X
%       VTYPE : character string of length NX (number of elements in X),
%               or 1 (value applies to all variables in x),
%               allowed values are 'C' (continuous), 'B' (binary),
%               'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer).
%       OPT : optional options structure with the following fields,
%           all of which are also optional (default values shown in
%           parentheses)
%           verbose (0) - controls level of progress output displayed
%               0 = no progress output
%               1 = some progress output
%               2 = verbose progress output
%           skip_prices (0) - flag that specifies whether or not to
%               skip the price computation stage, in which the problem
%               is re-solved for only the continuous variables, with all
%               others being constrained to their solved values
%           price_stage_warn_tol (1e-7) - tolerance on the objective fcn
%               value and primal variable relative match required to avoid
%               mis-match warning message
%           cplex_opt - options struct for CPLEX, value in verbose
%                   overrides these options
%       PROBLEM : The inputs can alternatively be supplied in a single
%           PROBLEM struct with fields corresponding to the input arguments
%           described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt
%
%   Outputs:
%       X : solution vector
%       F : final objective function value
%       EXITFLAG : CPLEXQP/CPLEXLP exit flag
%           (see CPLEXQP and CPLEXLP documentation for details)
%       OUTPUT : CPLEXQP/CPLEXLP output struct
%           (see CPLEXQP and CPLEXLP documentation for details)
%       LAMBDA : struct containing the Langrange and Kuhn-Tucker
%           multipliers on the constraints, with fields:
%           mu_l - lower (left-hand) limit on linear constraints
%           mu_u - upper (right-hand) limit on linear constraints
%           lower - lower bound on optimization variables
%           upper - upper bound on optimization variables
%
%   Note the calling syntax is almost identical to that of QUADPROG
%   from MathWorks' Optimization Toolbox. The main difference is that
%   the linear constraints are specified with A, L, U instead of
%   A, B, Aeq, Beq.
%
%   Calling syntax options:
%       [x, f, exitflag, output, lambda] = ...
%           miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
%
%       x = miqps_cplex(H, c, A, l, u)
%       x = miqps_cplex(H, c, A, l, u, xmin, xmax)
%       x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0)
%       x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype)
%       x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
%       x = miqps_cplex(problem), where problem is a struct with fields:
%                       H, c, A, l, u, xmin, xmax, x0, vtype, opt
%                       all fields except 'c', 'A' and 'l' or 'u' are optional
%       x = miqps_cplex(...)
%       [x, f] = miqps_cplex(...)
%       [x, f, exitflag] = miqps_cplex(...)
%       [x, f, exitflag, output] = miqps_cplex(...)
%       [x, f, exitflag, output, lambda] = miqps_cplex(...)
%
%
%   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
%       H = [   1003.1  4.3     6.3     5.9;
%               4.3     2.2     2.1     3.9;
%               6.3     2.1     3.5     4.8;
%               5.9     3.9     4.8     10  ];
%       c = zeros(4,1);
%       A = [   1       1       1       1;
%               0.17    0.11    0.10    0.18    ];
%       l = [1; 0.10];
%       u = [1; Inf];
%       xmin = zeros(4,1);
%       x0 = [1; 0; 0; 1];
%       opt = struct('verbose', 2);
%       [x, f, s, out, lambda] = miqps_cplex(H, c, A, l, u, xmin, [], x0, vtype, opt);
%
%   See also MIQPS_MASTER, CPLEXMIQP, CPLEXMILP, CPLEXQP, CPLEXLP,
%   CPLEX_OPTIONS.

%   MP-Opt-Model
%   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% check for CPLEX
% if ~have_feature('cplexqp')
%     error('miqps_cplex: requires the MATLAB interface for CPLEX');
% end

%%----- input argument handling  -----
%% gather inputs
if nargin == 1 && isstruct(H)       %% problem struct
    p = H;
    if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
    if isfield(p, 'vtype'), vtype = p.vtype;else,   vtype = []; end
    if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
    if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
    if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
    if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
    if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
    if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
    if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
    if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
else                                %% individual args
    if nargin < 10
        opt = [];
        if nargin < 9
            vtype = [];
            if nargin < 8
                x0 = [];
                if nargin < 7
                    xmax = [];
                    if nargin < 6
                        xmin = [];
                    end
                end
            end
        end
    end
end

%% define nx, set default values for missing optional inputs
if isempty(H) || ~any(any(H))
    if isempty(A) && isempty(xmin) && isempty(xmax)
        error('miqps_cplex: LP problem must include constraints or variable bounds');
    else
        if ~isempty(A)
            nx = size(A, 2);
        elseif ~isempty(xmin)
            nx = length(xmin);
        else    % if ~isempty(xmax)
            nx = length(xmax);
        end
    end
else
    nx = size(H, 1);
end
if isempty(c)
    c = zeros(nx, 1);
end
if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
                                 (isempty(u) || all(u == Inf)))
    A = sparse(0,nx);           %% no limits => no linear constraints
end
nA = size(A, 1);                %% number of original linear constraints
if isempty(u)                   %% By default, linear inequalities are ...
    u = Inf(nA, 1);             %% ... unbounded above and ...
end
if isempty(l)
    l = -Inf(nA, 1);            %% ... unbounded below.
end
if isempty(xmin)                %% By default, optimization variables are ...
    xmin = -Inf(nx, 1);         %% ... unbounded below and ...
end
if isempty(xmax)
    xmax = Inf(nx, 1);          %% ... unbounded above.
end
if isempty(x0)
    x0 = zeros(nx, 1);
end

%% default options
if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
    verbose = opt.verbose;
else
    verbose = 0;
end

%% split up linear constraints
ieq = find( abs(u-l) <= eps );          %% equality
igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
Ae = A(ieq, :);
be = u(ieq);
Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];

%% grab some dimensions
nlt = length(ilt);      %% number of upper bounded linear inequalities
ngt = length(igt);      %% number of lower bounded linear inequalities
nbx = length(ibx);      %% number of doubly bounded linear inequalities

%% set up options struct for CPLEX
if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt)
    cplex_opt = cplex_options(opt.cplex_opt);
else
    cplex_opt = cplex_options;
end

vstr = have_feature('cplex', 'vstr');
vnum = have_feature('cplex', 'vnum');
vrb = max([0 verbose-1]);
if vrb && vnum > 12.002 && vnum < 12.007
    cplex_opt.diagnostics   = 'on';
end
if verbose > 2
    cplex_opt.display = 'iter';
elseif verbose > 1
    cplex_opt.display = 'on';
elseif verbose > 0
    cplex_opt.display = 'off';
end

if isempty(Ai) && isempty(Ae)
    unconstrained = 1;
    Ae = sparse(1, nx);
    be = 0;
else
    unconstrained = 0;
end

%% call the solver
if verbose
    alg_names = {
        'default',
        'primal simplex',
        'dual simplex',
        'network simplex',
        'barrier',
        'sifting',
        'concurrent'
    };
end
if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I' | ...
        vtype == 'S' | vtype == 'N'))
    mi = 0;
else
    mi = 1;
    %% expand vtype to nx elements if necessary
    if length(vtype) == 1 && nx > 1
        vtype = char(vtype * ones(1,nx));
    end
end

if mi
    if isempty(H) || ~any(any(H))
        if verbose
            fprintf('CPLEX Version %s -- %s MILP solver\n', ...
                vstr, alg_names{cplex_opt.lpmethod+1});
        end
        [x, f, eflag, output] = ...
            cplexmilp(c, Ai, bi, Ae, be, [], [], [], xmin, xmax, vtype, x0, cplex_opt);
        lam = [];
    else
        if verbose
            fprintf('CPLEX Version %s --  %s MIQP solver\n', ...
                vstr, alg_names{cplex_opt.qpmethod+1});
        end
        %% ensure H is numerically symmetric
        if ~isequal(H, H')
            H = (H + H')/2;
        end
        [x, f, eflag, output] = ...
            cplexmiqp(H, c, Ai, bi, Ae, be, [], [], [], xmin, xmax, vtype, x0, cplex_opt);
    end
    lam = [];
else
    if isempty(H) || ~any(any(H))
        if verbose
            fprintf('CPLEX Version %s -- %s LP solver\n', ...
                vstr, alg_names{cplex_opt.lpmethod+1});
        end
        [x, f, eflag, output, lam] = ...
            cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
    else
        if verbose
            fprintf('CPLEX Version %s --  %s QP solver\n', ...
                vstr, alg_names{cplex_opt.qpmethod+1});
        end
        %% ensure H is numerically symmetric
        if ~isequal(H, H')
            H = (H + H')/2;
        end
        [x, f, eflag, output, lam] = ...
            cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
    end
end

%% workaround for eflag == 5 (which we have seen return infeasible results)
%%          cplexstatus: 6
%%    cplexstatusstring: 'non-optimal'
%%              message: 'Solution with numerical issues'
if eflag > 1
    warning('qps_cplex: Undocumented ''exitflag'' value (%d)\n          cplexstatus: %d\n    cplexstatusstring: ''%s''\n              message: ''%s''', eflag, output.cplexstatus, output.cplexstatusstring, output.message);
    if eflag == 5 && mi
        eflag = 1;      %% give it a try for the MI phase
    else
        eflag = -100 - eflag;
    end
end

%% check for empty results (in case optimization failed)
if isempty(x)
    x = NaN(nx, 1);
end
if isempty(f)
    f = NaN;
end
if isempty(lam)
    lam.ineqlin = NaN(length(bi), 1);
    lam.eqlin   = NaN(length(be), 1);
    lam.lower   = NaN(nx, 1);
    lam.upper   = NaN(nx, 1);
    mu_l        = NaN(nA, 1);
    mu_u        = NaN(nA, 1);
else
    mu_l        = zeros(nA, 1);
    mu_u        = zeros(nA, 1);
end
if unconstrained
    lam.eqlin = [];
end

%% negate prices depending on version
if vnum < 12.003
    lam.eqlin   = -lam.eqlin;
    lam.ineqlin = -lam.ineqlin;
end

%% repackage lambdas
kl = find(lam.eqlin < 0);   %% lower bound binding
ku = find(lam.eqlin > 0);   %% upper bound binding

mu_l(ieq(kl)) = -lam.eqlin(kl);
mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));

mu_u(ieq(ku)) = lam.eqlin(ku);
mu_u(ilt) = lam.ineqlin(1:nlt);
mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));

lambda = struct( ...
    'mu_l', mu_l, ...
    'mu_u', mu_u, ...
    'lower', lam.lower, ...
    'upper', lam.upper ...
);

if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices)
    if verbose
        fprintf('--- Integer stage complete, starting price computation stage ---\n');
    end
    if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol)
        tol = opt.price_stage_warn_tol;
    else
        tol = 1e-7;
    end
    k = find(vtype == 'I' | vtype == 'B' | vtype == 'N' | ...
            (vtype == 'S' & x' == 0));
    x(k) = round(x(k));
    xmin(k) = x(k);
    xmax(k) = x(k);
    x0 = x;
    opt.cplex_opt.lpmethod = 1;     %% primal simplex
    opt.cplex_opt.qpmethod = 1;     %% primal simplex

    [x_, f_, eflag_, output_, lambda] = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt);
    if eflag ~= eflag_
        error('miqps_cplex: EXITFLAG from price computation stage = %d', eflag_);
    end
    if abs(f - f_)/max(abs(f), 1) > tol
        warning('miqps_cplex: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1));
    end
    xn = x;
    xn(abs(xn)<1) = 1;
    [mx, k] = max(abs(x - x_) ./ xn);
    if mx > tol
        warning('miqps_cplex: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k));
    end
    output.price_stage = output_;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%INTERSECTCONVEXCONVEX  Finds intersection of two convex polygons
%
%     [p12, area] = cv.intersectConvexConvex(p1, p2)
%     [...] = cv.intersectConvexConvex(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __p1__ first polygon, stored in numeric array (Nx2/Nx1x2/1xNx2) or cell
%   array of 2-element vectors (`{[x,y], ...}`).
% * __p2__ second polygon, stored in numeric array (Nx2/Nx1x2/1xNx2) or cell
%   array of 2-element vectors (`{[x,y], ...}`).
%
% ## Output
% * __p12__ intersection polygon points, a cell array of 2-element vectors
%   `{[x,y], ...}`
% * __a__ area
%
% ## Options
% * __HandleNested__ default true
%
% See also: cv.rotatedRectangleIntersection, cv.Rect.intersect, rectint
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef MimMemoryCacheItem < handle
    % MimMemoryCacheItem. Part of the internal framework of the TD MIM Toolkit.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the framework of the TD MIM Toolkit.
    %
    %     Used to hold a cached value and metadata in the memory cache.
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        Value
        Info
    end

    methods

        function obj = MimMemoryCacheItem(value, info)
            obj.Value = value;
            obj.Info = info;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%
%		 imgOut = bilateralFilterS(img, imgEdges, sigma_s, sigma_r)
%
%        This function implements a bilateral filter without
%        approximations. Note this function is very slow!
%
%		 Input:
%			-img: is an image to be filtered.
%           -imgEdges: is an edge image, where its edges will be transfered
%           to img. Note set this to [] if you do not want to transfer
%           edges.
%           -sigma_s: is the spatial sigma value.
%           -sigma_r: is the range sigma value.
%
%		 Output:
%			-imgOut: is the filtered image.
%
%     Copyright (C) 2014-18  Francesco Banterle
%
%     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 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

author: banterle


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

% You need to set these values correctly
X_norm = X;
mu = zeros(1, size(X, 2));
sigma = zeros(1, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: First, for each feature dimension, compute the mean
%               of the feature and subtract it from the dataset,
%               storing the mean value in mu. Next, compute the
%               standard deviation of each feature and divide
%               each feature by it's standard deviation, storing
%               the standard deviation in sigma.
%
%               Note that X is a matrix where each column is a
%               feature and each row is an example. You need
%               to perform the normalization separately for
%               each feature.
%
% Hint: You might find the 'mean' and 'std' functions useful.
%

mu = mean(X);
sigma = std(X);

for i = 1:size(X,2)
	X_norm(:,i) = (X(:,i) - mu(i)) / sigma(i);
end

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

end

author: gopaczewski

repo: coursera-ml


text:
% Do a reflection in the fourier domain for a 2-dimensional signal

function xf_reflected = reflect_spectrum2(xf)
% use fliplr and flipud to replace the flip function or the rot90 func
% xf_reflected = circshift(flip(flip(xf, 1), 2), [1 1 0]);
xf_reflected = circshift(rot90(xf, 2), [1 1 0]);

author: flyers

repo: drone-tracking