text:
function An = slnorm(A, p, d)
%SLNORM Compute the Lp-norms
%
% $ Syntax $
%   - An = slnorm(A)
%   - An = slnorm(A, p)
%   - An = slnorm(A, [], d)
%   - An = slnorm(A, p, d)
%
% $ Arguments $
%   - A:        the input array
%   - p:        the order of the norm (default = 2)
%   - d:        the dimension of subarrays (default = 1, column vectors)
%   - An:       the resultant array
%
% $ Description $
%   - An = slnorm(A) computes the L2 norm of column vectors of A. If A is
%     an array of size d1 x d2 x ... dn, then An would be of size
%     1 x d2 x ... dn. Each value of An is the L2 norm of corresponding
%     column vector.
%
%   - An = slnorm(A, p) computes the Lp norm of column vectors of A. If A
%     is an array of size d1 x d2 x ... dn, then An would be of size
%     1 x d2 x ... dn. Each value of An is the Lp norm of corresponding
%     column vector.
%
%   - An = slnorm(A, [], d) computes the L2 norm of sub-arrays along the
%     dimensions specified by d. For example, if A is a d1 x d2 x ... dn
%     array, and d = [1, 2], then An would be an array of size
%     1 x 1 x d3 x ... x dn. Each value of An is the square root of square-
%     sum of values in corresponding plane.
%
%   - An = slnorm(A, p, d) computes the Lp norm of sub-arrays along the
%     dimensions specified by d.
%
% $ Remarks $
%   # Lp norm is the pth order-root of the sum of p-th power of all values
%     in a subspace.
%   # p can be inf or -inf. If p = inf, then the Lp norm is simply the
%     maximum magnitude value; while if p = -inf, then the Lp norm is the
%     minimum magnitude value.
%
% $ History $
%   - Created by Dahua Lin on Nov 19th, 2005
%

%% parse and verify input arguments
if nargin < 2 || isempty(p)
    p = 2;
end
if nargin < 3 || isempty(d)
    d = 1;
end
if ~isscalar(p)
    error('sltoolbox:notscalar', 'p should be a scalar');
end
if p == 0
    error('sltoolbox:zeroerr', 'p should not be zero');
end

%% compute
An = abs(A);
if p == 1
    An = slsum(An, d);
elseif p == 2
    An = sqrt(slsum(An .* An ,d));
elseif p == inf
    An = slmax(An, d);
elseif p == (-inf)
    An = slmin(An, d);
else
    An = slsum(An.^p, d) .^ (1/p);
end

author: lmthang


text:
function cS = kyanite

%
cs_Ky = crystalSymmetry('-1',[0.907 1 0.71], [90, 101.2, 106]*degree);
N = Miller({1,0,0},{0,1,0},{0,0,1},{1,-1,0},{1,1,0},cs_Ky);
dist = [0.4, 1, 2.4, 1.05, 1.2];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

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


text:
% readneurolocs() - read neuroscan polar location file (.asc)
%
% Usage:
%   >> CHANLOCS = readneurolocs( filename );
%   >> CHANLOCS = readneurolocs( filename, method, 'key1', val1, ...);
%
% Inputs:
%   filename  - file name or matlab cell array { names x_coord y_coord }
%   method    - [1 2, 3, 4 or 5] different import methods
%
% Optional inputs:
%   same as caliblocs()
%   note that if no optional input are provided, re-centering will be
%   performed automatically and re-scaling of coordinates will be
%   performed for '.asc' files (not '.dat' files).
%
% Outputs:
%   CHANLOCS       - EEGLAB channel location data structure. See
%                    help readlocs()
%
% Author: Arnaud Delorme, CNL / Salk Institute, 4 March 2003
%
% See also: readlocs()

% Copyright (C) 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function chanlocs = readneurolocs( filename, method, varargin)

if nargin < 1
    help readneurolocs;
    return;
end
if nargin < 2
    plottag = 0;
end

if nargin < 2
    disp('There are 5 methods to import Neuroscan files, if one does not work, try others')
    try
        chanlocs = readneurolocs( filename, 1);
        return
    catch
        disp('Import method 1 failed, trying method 2');
        try
            chanlocs = readneurolocs( filename, 2);
            return
        catch
            disp('Import method 2 failed, trying method 3');
            try
                chanlocs = readneurolocs( filename, 3);
                return
            catch
                disp('Import method 3 failed, trying method 4');
                try
                    chanlocs = readneurolocs( filename, 4);
                    return
                catch
                    disp('Import method 3 failed, trying method 5');
                    try
                        chanlocs = readneurolocs( filename, 5);
                        return
                    catch
                        error('No import method worked');
                    end
                end
            end
        end
    end
end

% try new method
if method == 1
    locs = readtable(filename, 'filetype', 'text', 'Delimiter', { ';' ' ' '\t' }, 'ConsecutiveDelimitersRule', 'join');
    locs = table2cell(locs);

    if mod(size(locs,1),2) == 1
        error('Issue with file format')
    end
    halfHeight = size(locs,1)/2;
    loc1 = locs(1:halfHeight,:);
    loc2 = locs(halfHeight+1:end,:);
    x = [loc2{:,3}]; x = x-0.5;
    y = [loc2{:,4}]; y = y-0.5;
    radius = sqrt(x.^2 + y.^2);
    theta  = atan2d(x, y);

    for iChan = 1:halfHeight
        chanlocs(iChan).label = loc1{iChan,2};
        chanlocs(iChan).theta = theta(iChan);
        chanlocs(iChan).radius = radius(iChan);
    end
    chanlocs = convertlocs( chanlocs, 'topo2all');
    return
elseif method == 2
    % read location file
    % ------------------
    if ischar(filename)
        locs  = loadtxt( filename, 'delim', 9 );
    end

    if ~ischar(filename) || locs{1,1}(1) == ';' || size(locs,2) < 5
        if ~ischar(filename)
            names = filename{1};
            x     = filename{2};
            y     = filename{3};
        else
            if locs{1,1}(1) == ';'
                % remove trailing control channels
                % --------------------------------
                while isnumeric( locs{end,1} ) & locs{end,1} ~= 0
                    locs  = locs(1:end-1,:);
                end

                % find first numerical index
                % --------------------------
                index = 1;
                while ischar( locs{index,1} ) && index < size(locs,1)
                    index = index + 1;
                end

                % extract location array
                % ----------------------
                nchans = size( locs, 1 ) - index +1;
                chans  = [locs{end-nchans+1:end, 1:5}];
                chans  = reshape(chans,nchans,5);               %% Added this line in order to get x = chans(:,3)
                names  = locs(end-nchans*2+1: end-nchans, 2);
                for index = 1:length(names)
                    if ~ischar(names{index})
                        names{index} = int2str(names{index});
                    end
                end
                x = chans(:,3);
                y = -chans(:,4);
            else
                [tmp2 tmpind] = sort( [ locs{:,1} ]);
                locs = locs(tmpind,:);
                y      = [ locs{:,end} ];
                x      = [ locs{:,end-1} ];
                x      = x/513.1617*44;
                y      = y/513.1617*44;
                names = locs(:,end-2);
            end
        end

        % second solution using angle
        % ---------------------------
        [phi,theta] = cart2pol(x, y);
        phi = phi/pi*180;

        % convert to other types of coordinates
        % -------------------------------------
        labels = names';
        labels = cellfun(@num2str, labels, 'UniformOutput', false);
        chanlocs = struct('labels', labels', 'sph_theta_besa', mattocell(theta)', 'sph_phi_besa', mattocell(phi)');      %% labels instead of labels(:)
        chanlocs = convertlocs( chanlocs, 'sphbesa2all');

        for index = 1:length(chanlocs)
            chanlocs(index).labels = num2str(chanlocs(index).labels);
        end

        % re-calibration
        % --------------
        chanlocs = adjustlocs(chanlocs, 'autoscale', 'on', 'autorotate', 'off', varargin{:});
    end
elseif method == 3
    chanlocs = readneurodat(  filename);
elseif method == 4
    % read location file
    % ------------------
    if ischar(filename)
        locs  = loadtxt( filename, 'delim', [9 ' ']);
    end

    if size(locs,2) == 5
        % 5 rows, xyz positions
        for index = 1:size(locs,1)
            locs{index,3} = - locs{index,3};
        end
        chanlocs = struct('labels', locs(:,1), 'type', locs(:,2), 'X', locs(:,4), 'Y', locs(:,3), 'Z', locs(:,5));
        chanlocs = convertlocs( chanlocs, 'cart2all');
    elseif size(locs,2) == 4
        chanlocs = readlocs(filename, 'filetype', 'custom', 'format', { 'labels' '-Y' 'X' 'Z' });
    end
elseif method == 5
    chanlocs = readlocs(filename, 'filetype', 'custom', 'format', { 'labels' 'ignore' '-Y' 'X' 'Z' });
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function curveDisplay(x,y,varargin)
%curveDisplay Display of 2-D curve.
%   curveDisplay(X,Y,LINESPEC) displays the coordinates (x,y) of a CURVE.
%
%   LineSpec is a character string made from one element from any or all
%   the following 3 columns:
%
%            b     blue          .     point           -     solid
%            g     green         o     circle          :     dotted
%            r     red           x     x-mark          -.    dashdot
%            c     cyan          +     plus            --    dashed
%            m     magenta       *     star          (none)  no line
%            y     yellow        s     square
%            k     black         d     diamond
%            w     white         v     triangle (down)
%                                ^     triangle (up)
%                                <     triangle (left)
%                                >     triangle (right)
%                                p     pentagram
%                                h     hexagram
%
%
%   A symbol from the first column gives the color, the second is the
%   symbol used in the plot, and the third specifies the type of line
%   used to join the points in the plot. For example, to plot red circles
%   joined by straight lines we speficify the string 'ro-'. To plot just
%   red circles without any connecting lines we specify the string 'ro'.
%   The default for curveDisplay is black dots with no lines.
%
%   curveDisplay(X,Y,LINESPEC,NAME,VALUE) further modifies the lines and
%   markers using Name,Value pairs. Any Name,Value pairs supported by
%   the PLOT function are allowed, including:
%
%           Name                          Value
%
%        LineWidth         Width (in points) of the line and border
%                          of filled markers (circle, square, diamond,
%                          pentagram, hexagram, and the four triangles
%                          points). The default is 0.5 pt.
%
%          Color           The color of the line. The color can be
%                          specified as an RGB triplet (such as [0.5 1.0
%                          0.8]), a hexadecimal color code (such as
%                          '#FF8800'), a color name (such as 'red'), or
%                          a short color name (listed above in the
%                          description for LineSpec), or 'none'. The
%                          default is 'k' (black).
%
%     MarkerEdgeColor      Color of the marker or the color of the edge
%                          of the marker for filled markers, specified
%                          as described above for the Color parameter.
%                          The default is black, or, if a line is
%                          specified, the same color as the line joining
%                          the markers.
%
%    MarkerFaceColor       The fill color for the filled markers,
%                          specified as described above for the Color
%                          parameter. The default is 'none'.
%
%      MarkerSize          The size of the marker in points. The default
%                          is 7 pt.
%
%  Example: curveDisplay(x,y,'ro-','MarkerFaceColor','g','MarkerSize',4)
%  displays red circles of size 4 pt, connected by red solid line of 0.5
%  pt thick, with the circles filled in green.
%
%	To superimpose the curve on an image, f, use the following syntax:
%	figure, imshow(f)
%	hold on
%	curveDisplay(x,y,varargin)
%	hold off
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

narginchk(2,Inf)
if nargin == 2
   % Default.
   plot(y,x,'.')
else
   if ~isodd(length(varargin))
      error('Wrong number of inputs.');
   end
   plot(y,x,varargin{:});
end

author: dipum

repo: dipum-toolbox

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


text:
function res = fnamedat(this)
% Method for getting the name of the data file
% FORMAT res = fnamedat(this)
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: fnamedat.m 5025 2012-10-31 14:44:13Z vladimir $

if islinked(this)
    res = this.data.fname;
else
    res = [];
end

author: spm


text:
function C = competitiveHebb2(C,bmu,secbmu,y,Theta,alpha)
% This function applies type 2 Competitive Hebbian Learning,
% as described by relation (2.4) of reference [1].
% [1] Bruske J., Sommer G., "Dynamic Cell Structure Learns Perfectly Topology Preserving Map",
%                                                Neural Computation, vol. 7, Issue 4, July 1995, pp. 845-865.
c =  max(y(bmu)*y(secbmu),C(bmu,secbmu));

% Reduce all connections by multiplying by alpha:
C = alpha*C;

% Find those connections with: 0<strength<Theta:
[rows cols] = find(C<Theta & C>0);
... and erase them!
NumZeroConnections = length(rows);
for i=1:NumZeroConnections
     C(rows(i),cols(i)) = 0;
end

C(bmu,secbmu) = c;
C(secbmu,bmu) =  C(bmu,secbmu);

author: Sable

repo: mcbench-benchmarks


text:
function [dat, label, time, cfg] = preproc(dat, label, time, cfg, begpadding, endpadding)

% PREPROC applies various preprocessing steps on a single piece of EEG/MEG data
% that has been read from a data file.
%
% This low-level function serves as a subfunction for all FieldTrip modules that want
% to preprocess the data, such as FT_PREPROCESSING, FT_ARTIFACT_XXX,
% FT_TIMELOCKANALYSIS, etc. It ensures consistent handling of both MEG and EEG data
% and consistency in the use of all preprocessing configuration options.
%
% Use as
%   [dat, label, time, cfg] = preproc(dat, label, time, cfg, begpadding, endpadding)
%
% The required input arguments are
%   dat         Nchan x Ntime data matrix
%   label       Nchan x 1 cell-array with channel labels
%   time        Ntime x 1 vector with the latency in seconds
%   cfg         configuration structure, see below
% and the optional input arguments are
%   begpadding  number of samples that was used for padding (see below)
%   endpadding  number of samples that was used for padding (see below)
%
% The output is
%   dat         Nchan x Ntime data matrix
%   label       Nchan x 1 cell-array with channel labels
%   time        Ntime x 1 vector with the latency in seconds
%   cfg         configuration structure, optionally with extra defaults set
%
% Note that the number of input channels and the number of output channels can be
% different, for example when the user specifies that he/she wants to add the
% implicit EEG reference channel to the data matrix.
%
% The filtering of the data can introduce artifacts at the edges, hence it is better
% to pad the data with some extra signal at the begin and end. After filtering, this
% padding is removed and the other preprocessing steps are applied to the remainder
% of the data. The input fields begpadding and endpadding should be specified in
% samples. You can also leave them empty, which implies that the data is not padded.
%
% The configuration can contain
%   cfg.lpfilter      = 'no' or 'yes'  lowpass filter
%   cfg.hpfilter      = 'no' or 'yes'  highpass filter
%   cfg.bpfilter      = 'no' or 'yes'  bandpass filter
%   cfg.bsfilter      = 'no' or 'yes'  bandstop filter
%   cfg.dftfilter     = 'no' or 'yes'  line noise removal using discrete fourier transform
%   cfg.medianfilter  = 'no' or 'yes'  jump preserving median filter
%   cfg.lpfreq        = lowpass  frequency in Hz
%   cfg.hpfreq        = highpass frequency in Hz
%   cfg.bpfreq        = bandpass frequency range, specified as [low high] in Hz
%   cfg.bsfreq        = bandstop frequency range, specified as [low high] in Hz
%   cfg.dftfreq       = line noise frequencies for DFT filter, default [50 100 150] Hz
%   cfg.lpfiltord     = lowpass  filter order (default set in low-level function)
%   cfg.hpfiltord     = highpass filter order (default set in low-level function)
%   cfg.bpfiltord     = bandpass filter order (default set in low-level function)
%   cfg.bsfiltord     = bandstop filter order (default set in low-level function)
%   cfg.medianfiltord = length of median filter
%   cfg.lpfilttype    = digital filter type, 'but' (default) or 'firws' or 'fir' or 'firls'
%   cfg.hpfilttype    = digital filter type, 'but' (default) or 'firws' or 'fir' or 'firls'
%   cfg.bpfilttype    = digital filter type, 'but' (default) or 'firws' or 'fir' or 'firls'
%   cfg.bsfilttype    = digital filter type, 'but' (default) or 'firws' or 'fir' or 'firls'
%   cfg.lpfiltdir     = filter direction, 'twopass' (default), 'onepass' or 'onepass-reverse' or 'onepass-zerophase' (default for firws) or 'onepass-minphase' (firws, non-linear!)
%   cfg.hpfiltdir     = filter direction, 'twopass' (default), 'onepass' or 'onepass-reverse' or 'onepass-zerophase' (default for firws) or 'onepass-minphase' (firws, non-linear!)
%   cfg.bpfiltdir     = filter direction, 'twopass' (default), 'onepass' or 'onepass-reverse' or 'onepass-zerophase' (default for firws) or 'onepass-minphase' (firws, non-linear!)
%   cfg.bsfiltdir     = filter direction, 'twopass' (default), 'onepass' or 'onepass-reverse' or 'onepass-zerophase' (default for firws) or 'onepass-minphase' (firws, non-linear!)
%   cfg.lpinstabilityfix = deal with filter instability, 'no', 'reduce', 'split' (default  = 'no')
%   cfg.hpinstabilityfix = deal with filter instability, 'no', 'reduce', 'split' (default  = 'no')
%   cfg.bpinstabilityfix = deal with filter instability, 'no', 'reduce', 'split' (default  = 'no')
%   cfg.bsinstabilityfix = deal with filter instability, 'no', 'reduce', 'split' (default  = 'no')
%   cfg.lpfiltdf      = lowpass transition width (firws, overrides order, default set in low-level function)
%   cfg.hpfiltdf      = highpass transition width (firws, overrides order, default set in low-level function)
%   cfg.bpfiltdf      = bandpass transition width (firws, overrides order, default set in low-level function)
%   cfg.bsfiltdf      = bandstop transition width (firws, overrides order, default set in low-level function)
%   cfg.lpfiltwintype = lowpass window type, 'hann' or 'hamming' (default) or 'blackman' or 'kaiser' (firws)
%   cfg.hpfiltwintype = highpass window type, 'hann' or 'hamming' (default) or 'blackman' or 'kaiser' (firws)
%   cfg.bpfiltwintype = bandpass window type, 'hann' or 'hamming' (default) or 'blackman' or 'kaiser' (firws)
%   cfg.bsfiltwintype = bandstop window type, 'hann' or 'hamming' (default) or 'blackman' or 'kaiser' (firws)
%   cfg.lpfiltdev     = lowpass max passband deviation (firws with 'kaiser' window, default 0.001 set in low-level function)
%   cfg.hpfiltdev     = highpass max passband deviation (firws with 'kaiser' window, default 0.001 set in low-level function)
%   cfg.bpfiltdev     = bandpass max passband deviation (firws with 'kaiser' window, default 0.001 set in low-level function)
%   cfg.bsfiltdev     = bandstop max passband deviation (firws with 'kaiser' window, default 0.001 set in low-level function)
%   cfg.dftreplace    = 'zero' or 'neighbour', method used to reduce line noise, 'zero' implies DFT filter, 'neighbour' implies spectrum interpolation (default = 'zero')
%   cfg.dftbandwidth  = bandwidth of line noise frequencies, applies to spectrum interpolation, in Hz (default = [1 2 3])
%   cfg.dftneighbourwidth = bandwidth of frequencies neighbouring line noise frequencies, applies to spectrum interpolation, in Hz (default = [2 2 2])
%   cfg.plotfiltresp  = 'no' or 'yes', plot filter responses (firws, default = 'no')
%   cfg.usefftfilt    = 'no' or 'yes', use fftfilt instead of filter (firws, default = 'no')
%   cfg.demean        = 'no' or 'yes'
%   cfg.baselinewindow = [begin end] in seconds, the default is the complete trial
%   cfg.detrend       = 'no' or 'yes', this is done on the complete trial
%   cfg.polyremoval   = 'no' or 'yes', this is done on the complete trial
%   cfg.polyorder     = polynome order (default = 2)
%   cfg.derivative    = 'no' (default) or 'yes', computes the first order derivative of the data
%   cfg.hilbert       = 'no', 'abs', 'complex', 'real', 'imag', 'absreal', 'absimag' or 'angle' (default = 'no')
%   cfg.rectify       = 'no' or 'yes'
%   cfg.precision     = 'single' or 'double' (default = 'double')
%   cfg.absdiff       = 'no' or 'yes', computes absolute derivative (i.e.first derivative then rectify)
%
% Preprocessing options that you should only use for EEG data are
%   cfg.reref         = 'no' or 'yes' (default = 'no')
%   cfg.refchannel    = cell-array with new EEG reference channel(s)
%   cfg.refmethod     = 'avg', 'median', 'rest', 'bipolar' or 'laplace' (default = 'avg')
%   cfg.groupchans    = 'yes' or 'no', should channels be rereferenced in separate groups
%                       for bipolar and laplace methods, this requires channnels to be
%                       named using an alphanumeric code, where letters represent the
%                       group and numbers represent the order of the channel whithin
%                       its group (default = 'no')
%   cfg.leadfield     = matrix or cell-array, this is required when refmethod is 'rest'
%                       The leadfield can be a single matrix (channels X sources) which
%                       is calculated by using the forward theory, based on the
%                       electrode montage, head model and equivalent source model.
%                       It can also be the output of FT_PREPARE_LEADFIELD based on a
%                       realistic head model.
%   cfg.implicitref   = 'label' or empty, add the implicit EEG reference as zeros (default = [])
%   cfg.montage       = 'no' or a montage structure (default = 'no')
%
% See also FT_READ_DATA, FT_READ_HEADER

% Copyright (C) 2004-2012, 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$

% compute fsample
fsample = 1./nanmean(diff(time));

if nargin<5 || isempty(begpadding)
  begpadding = 0;
end
if nargin<6 || isempty(endpadding)
  endpadding = 0;
end

if iscell(cfg)
  % recurse over the subsequent preprocessing stages
  if begpadding>0 || endpadding>0
    ft_error('multiple preprocessing stages are not supported in combination with filter padding');
  end
  for i=1:length(cfg)
    tmpcfg = cfg{i};
    if nargout==1
      [dat                     ] = preproc(dat, label, time, tmpcfg, begpadding, endpadding);
    elseif nargout==2
      [dat, label              ] = preproc(dat, label, time, tmpcfg, begpadding, endpadding);
    elseif nargout==3
      [dat, label, time        ] = preproc(dat, label, time, tmpcfg, begpadding, endpadding);
    elseif nargout==4
      [dat, label, time, tmpcfg] = preproc(dat, label, time, tmpcfg, begpadding, endpadding);
      cfg{i} = tmpcfg;
    end
  end
  % ready with recursing over the subsequent preprocessing stages
  return
end

% set the defaults for the rereferencing options
cfg.reref                = ft_getopt(cfg, 'reref',       'no');
cfg.refchannel           = ft_getopt(cfg, 'refchannel',  {});
cfg.refmethod            = ft_getopt(cfg, 'refmethod',   'avg');
cfg.implicitref          = ft_getopt(cfg, 'implicitref', []);
cfg.leadfield            = ft_getopt(cfg, 'leadfield',   []);
cfg.groupchans           = ft_getopt(cfg, 'groupchans', 'no');
% set the defaults for the signal processing options
cfg.polyremoval          = ft_getopt(cfg, 'polyremoval', 'no');
cfg.polyorder            = ft_getopt(cfg, 'polyorder',   2);
cfg.detrend              = ft_getopt(cfg, 'detrend',     'no');
cfg.demean               = ft_getopt(cfg, 'demean',      'no');
cfg.baselinewindow       = ft_getopt(cfg, 'baselinewindow', 'all');
cfg.dftfilter            = ft_getopt(cfg, 'dftfilter', 'no');
cfg.dftfreq              = ft_getopt(cfg, 'dftfreq', [50 100 150]);
cfg.lpfilter             = ft_getopt(cfg, 'lpfilter',  'no');
cfg.hpfilter             = ft_getopt(cfg, 'hpfilter',  'no');
cfg.bpfilter             = ft_getopt(cfg, 'bpfilter',  'no');
cfg.bsfilter             = ft_getopt(cfg, 'bsfilter',  'no');
cfg.lpfiltord            = ft_getopt(cfg, 'lpfiltord', []);
cfg.hpfiltord            = ft_getopt(cfg, 'hpfiltord', []);
cfg.bpfiltord            = ft_getopt(cfg, 'bpfiltord', []);
cfg.bsfiltord            = ft_getopt(cfg, 'bsfiltord', []);
cfg.lpfilttype           = ft_getopt(cfg, 'lpfilttype', 'but');
cfg.hpfilttype           = ft_getopt(cfg, 'hpfilttype', 'but');
cfg.bpfilttype           = ft_getopt(cfg, 'bpfilttype', 'but');
cfg.bsfilttype           = ft_getopt(cfg, 'bsfilttype', 'but');
if strcmp(cfg.lpfilttype, 'firws'), cfg.lpfiltdir = ft_getopt(cfg, 'lpfiltdir', 'onepass-zerophase'); else, cfg.lpfiltdir = ft_getopt(cfg, 'lpfiltdir', 'twopass'); end
if strcmp(cfg.hpfilttype, 'firws'), cfg.hpfiltdir = ft_getopt(cfg, 'hpfiltdir', 'onepass-zerophase'); else, cfg.hpfiltdir = ft_getopt(cfg, 'hpfiltdir', 'twopass'); end
if strcmp(cfg.bpfilttype, 'firws'), cfg.bpfiltdir = ft_getopt(cfg, 'bpfiltdir', 'onepass-zerophase'); else, cfg.bpfiltdir = ft_getopt(cfg, 'bpfiltdir', 'twopass'); end
if strcmp(cfg.bsfilttype, 'firws'), cfg.bsfiltdir = ft_getopt(cfg, 'bsfiltdir', 'onepass-zerophase'); else, cfg.bsfiltdir = ft_getopt(cfg, 'bsfiltdir', 'twopass'); end
cfg.lpinstabilityfix     = ft_getopt(cfg, 'lpinstabilityfix', 'no');
cfg.hpinstabilityfix     = ft_getopt(cfg, 'hpinstabilityfix', 'no');
cfg.bpinstabilityfix     = ft_getopt(cfg, 'bpinstabilityfix', 'no');
cfg.bsinstabilityfix     = ft_getopt(cfg, 'bsinstabilityfix', 'no');
cfg.lpfiltdf             = ft_getopt(cfg, 'lpfiltdf',  []);
cfg.hpfiltdf             = ft_getopt(cfg, 'hpfiltdf', []);
cfg.bpfiltdf             = ft_getopt(cfg, 'bpfiltdf', []);
cfg.bsfiltdf             = ft_getopt(cfg, 'bsfiltdf', []);
cfg.lpfiltwintype        = ft_getopt(cfg, 'lpfiltwintype', 'hamming');
cfg.hpfiltwintype        = ft_getopt(cfg, 'hpfiltwintype', 'hamming');
cfg.bpfiltwintype        = ft_getopt(cfg, 'bpfiltwintype', 'hamming');
cfg.bsfiltwintype        = ft_getopt(cfg, 'bsfiltwintype', 'hamming');
cfg.lpfiltdev            = ft_getopt(cfg, 'lpfiltdev', []);
cfg.hpfiltdev            = ft_getopt(cfg, 'hpfiltdev', []);
cfg.bpfiltdev            = ft_getopt(cfg, 'bpfiltdev', []);
cfg.bsfiltdev            = ft_getopt(cfg, 'bsfiltdev', []);
cfg.plotfiltresp         = ft_getopt(cfg, 'plotfiltresp', 'no');
cfg.usefftfilt           = ft_getopt(cfg, 'usefftfilt', 'no');
cfg.medianfilter         = ft_getopt(cfg, 'medianfilter', 'no');
cfg.medianfiltord        = ft_getopt(cfg, 'medianfiltord', 9);
cfg.hilbert              = ft_getopt(cfg, 'hilbert', 'no');
cfg.derivative           = ft_getopt(cfg, 'derivative', 'no');
cfg.rectify              = ft_getopt(cfg, 'rectify', 'no');
cfg.boxcar               = ft_getopt(cfg, 'boxcar', 'no');
cfg.absdiff              = ft_getopt(cfg, 'absdiff', 'no');
cfg.precision            = ft_getopt(cfg, 'precision', []);
cfg.conv                 = ft_getopt(cfg, 'conv', 'no');
cfg.montage              = ft_getopt(cfg, 'montage', 'no');
cfg.dftinvert            = ft_getopt(cfg, 'dftinvert', 'no');
cfg.standardize          = ft_getopt(cfg, 'standardize', 'no');
cfg.denoise              = ft_getopt(cfg, 'denoise', '');
cfg.subspace             = ft_getopt(cfg, 'subspace', []);
cfg.custom               = ft_getopt(cfg, 'custom', '');
cfg.resample             = ft_getopt(cfg, 'resample', '');

% test whether the MATLAB signal processing toolbox is available
if strcmp(cfg.medianfilter, 'yes') && ~ft_hastoolbox('signal')
  ft_error('median filtering requires the MATLAB signal processing toolbox');
end

% do a sanity check on the filter configuration
if strcmp(cfg.bpfilter, 'yes') && ...
    (strcmp(cfg.hpfilter, 'yes') || strcmp(cfg.lpfilter,'yes'))
  ft_error('you should not apply both a bandpass AND a lowpass/highpass filter');
end

% do a sanity check on the hilbert transform configuration
if strcmp(cfg.hilbert, 'yes') && ~strcmp(cfg.bpfilter, 'yes')
  ft_warning('Hilbert transform should be applied in conjunction with bandpass filter')
end

% do a sanity check on hilbert and rectification
if strcmp(cfg.hilbert, 'yes') && strcmp(cfg.rectify, 'yes')
  ft_error('Hilbert transform and rectification should not be applied both')
end

% do a sanity check on the rereferencing/montage
if ~strcmp(cfg.reref, 'no') && ~strcmp(cfg.montage, 'no')
  ft_error('cfg.reref and cfg.montage are mutually exclusive')
end

% lnfilter is no longer used
if isfield(cfg, 'lnfilter') && strcmp(cfg.lnfilter, 'yes')
  ft_error('line noise filtering using the option cfg.lnfilter is not supported any more, use cfg.bsfilter instead')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do the rereferencing in case of EEG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(cfg.implicitref) && ~any(match_str(cfg.implicitref,label))
  label = {label{:} cfg.implicitref}';
  dat(end+1,:) = 0;
end

if strcmp(cfg.reref, 'yes')
  switch cfg.refmethod
    case {'avg', 'median'}
      % mean or median based derivation of specified or all channels
      cfg.refchannel = ft_channelselection(cfg.refchannel, label);
      refindx = match_str(label, cfg.refchannel);
      if isempty(refindx),ft_error('reference channel was not found');end
      dat = ft_preproc_rereference(dat, refindx, cfg.refmethod);

    case {'rest'}
      cfg.refchannel = ft_channelselection(cfg.refchannel, label);
      refindx = match_str(label, cfg.refchannel);
      if isempty(refindx)
        ft_error('reference channel was not found')
      end

      if isempty(cfg.leadfield)
        ft_error('A leadfield is required to re-refer to REST');
      end
      % check the leadfield, apparently the code contributor here wants to
      % support 3 case, either a matrix, a struct, or a cell-array. Yet the
      % original code is almost impossible to parse, with a lot of try and
      % catch statements, so the current version below is an attempt to clean
      % this up
      if isnumeric(cfg.leadfield)
        Nchann_lf = size(cfg.leadfield,1); % No. of channels in leadfield
        lf_label  = {};
        G         = cfg.leadfield;
      elseif isstruct(cfg.leadfield)
        Nchann_lf = numel(cfg.leadfield.label);
        lf_label  = cfg.leadfield.label;
        G         = cat(2, cfg.leadfield.leadfield{:});
      elseif iscell(cfg.leadfield)
        Nchann_lf = size(cfg.leadfield{1},1); % No. of channels in leadfield
        lf_label  = {};
        G         = cat(2, cfg.leadfield{:});
      end

      if Nchann_lf ~= length(label)
        ft_error('channels in the leadfield are not equal to those in the data');
      end

      if ~isempty(lf_label)
        [indx1, indx2] = match_str(lf_label(refindx), label(refindx));
        if ~isequal(indx1, indx2)
          ft_error('The order in leadfield may be NOT the same as in the data, please check the leadfield!');
        end
      else
        ft_warning('There is no label info in the leadfield, there is no guarantee that the order of the channels in leadfield is the same as in the data');
      end

      dat   = ft_preproc_rereference(dat, refindx, cfg.refmethod, [], G); % re-referencing
      label = label(refindx); % re-referenced channel labels

    case {'bipolar', 'laplace', 'doublebanana', 'longitudinal', 'circumferential', 'transverse'}
      % this is implemented as a montage that the user does not get to see
      tmpcfg = keepfields(cfg, {'refmethod', 'implicitref', 'refchannel', 'channel', 'groupchans'});
      tmpcfg.showcallinfo = 'no';
      montage = ft_prepare_montage(tmpcfg);

      % convert the data temporarily to a raw structure
      tmpdata.trial = {dat};
      tmpdata.time  = {time};
      tmpdata.label = label;

      % apply the montage to the data
      tmpdata = ft_apply_montage(tmpdata, montage, 'feedback', 'none');
      dat     = tmpdata.trial{1}; % the number of channels can have changed
      label   = tmpdata.label;    % the output channels can be different than the input channels
      clear tmpdata

    otherwise
      ft_error('unsupported value for cfg.refmethod');
  end % switch refmethod
end % if reref

if ~strcmp(cfg.montage, 'no') && ~isempty(cfg.montage)
  % convert the data temporarily to a raw structure
  tmpdata.trial = {dat};
  tmpdata.time  = {time};
  tmpdata.label = label;
  % apply the montage to the data
  tmpdata = ft_apply_montage(tmpdata, cfg.montage, 'feedback', 'none');
  dat   = tmpdata.trial{1}; % the number of channels can have changed
  label = tmpdata.label;    % the output channels can be different than the input channels
  clear tmpdata
end

if any(isnan(dat(:)))
  % filtering is not possible for at least a selection of the data
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaNs, not all processing methods are robust to NaNs, so the NaNs might spread');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do the filtering on the padded data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(cfg.denoise)
  hflag    = isfield(cfg.denoise, 'hilbert') && strcmp(cfg.denoise.hilbert, 'yes');
  datlabel = match_str(label, cfg.denoise.channel);
  reflabel = match_str(label, cfg.denoise.refchannel);
  tmpdat   = ft_preproc_denoise(dat(datlabel,:), dat(reflabel,:), hflag);
  dat(datlabel,:) = tmpdat;
end

% The filtering should in principle be done prior to the demeaning to
% ensure that the resulting mean over the baseline window will be
% guaranteed to be zero (even if there are filter artifacts).
% However, the filtering benefits from the data being pulled towards zero,
% causing fewer edge artifacts. That is why we start by removing the slow
% drift, then filter, and then repeat the demean/detrend/polyremove. Note
% that this might interfere with a situation where the requested baseline
% window is not present in the local time axis of the trial, e.g. when the
% trials are of different length
if any(strcmp({cfg.polyremoval cfg.detrend cfg.demean}, 'yes'))
  if strcmp(cfg.polyremoval, 'yes')
    polyorder = cfg.polyorder;
  elseif strcmp(cfg.detrend, 'yes')
    polyorder = 1;
  elseif strcmp(cfg.demean, 'yes')
    polyorder = 0;
  end
  datorig   = dat;
  nsamples  = size(dat,2);
  begsample = 1        + begpadding;
  endsample = nsamples - endpadding;
  dat       = ft_preproc_polyremoval(dat, polyorder, begsample, endsample); % this will also demean and detrend
  datdiff   = datorig - dat;
end

if strcmp(cfg.medianfilter, 'yes'), dat = ft_preproc_medianfilter(dat, cfg.medianfiltord); end
if strcmp(cfg.lpfilter, 'yes'),     dat = ft_preproc_lowpassfilter(dat, fsample, cfg.lpfreq, cfg.lpfiltord, cfg.lpfilttype, cfg.lpfiltdir, cfg.lpinstabilityfix, cfg.lpfiltdf, cfg.lpfiltwintype, cfg.lpfiltdev, cfg.plotfiltresp, cfg.usefftfilt); end
if strcmp(cfg.hpfilter, 'yes'),     dat = ft_preproc_highpassfilter(dat, fsample, cfg.hpfreq, cfg.hpfiltord, cfg.hpfilttype, cfg.hpfiltdir, cfg.hpinstabilityfix, cfg.hpfiltdf, cfg.hpfiltwintype, cfg.hpfiltdev, cfg.plotfiltresp, cfg.usefftfilt); end
if strcmp(cfg.bpfilter, 'yes'),     dat = ft_preproc_bandpassfilter(dat, fsample, cfg.bpfreq, cfg.bpfiltord, cfg.bpfilttype, cfg.bpfiltdir, cfg.bpinstabilityfix, cfg.bpfiltdf, cfg.bpfiltwintype, cfg.bpfiltdev, cfg.plotfiltresp, cfg.usefftfilt); end
if strcmp(cfg.bsfilter, 'yes')
  for i=1:size(cfg.bsfreq,1)
    % apply a bandstop filter for each of the specified bands, i.e. cfg.bsfreq should be Nx2
    dat = ft_preproc_bandstopfilter(dat, fsample, cfg.bsfreq(i,:), cfg.bsfiltord, cfg.bsfilttype, cfg.bsfiltdir, cfg.bsinstabilityfix, cfg.bsfiltdf, cfg.bsfiltwintype, cfg.bsfiltdev, cfg.plotfiltresp, cfg.usefftfilt);
  end
end
if strcmp(cfg.polyremoval, 'yes')
  % the begin and endsample of the polyremoval period correspond to the complete data minus padding
  nsamples  = size(dat,2);
  begsample = 1        + begpadding;
  endsample = nsamples - endpadding;
  dat = ft_preproc_polyremoval(dat, cfg.polyorder, begsample, endsample);
end
if strcmp(cfg.detrend, 'yes')
  % the begin and endsample of the detrend period correspond to the complete data minus padding
  nsamples  = size(dat,2);
  begsample = 1        + begpadding;
  endsample = nsamples - endpadding;
  dat = ft_preproc_detrend(dat, begsample, endsample);
end
if strcmp(cfg.demean, 'yes')
  if ischar(cfg.baselinewindow) && strcmp(cfg.baselinewindow, 'all')
    % the begin and endsample of the baseline period correspond to the complete data minus padding
    nsamples  = size(dat,2);
    begsample = 1        + begpadding;
    endsample = nsamples - endpadding;
    dat       = ft_preproc_baselinecorrect(dat, begsample, endsample);
  else
    % determine the begin and endsample of the baseline period and baseline correct for it
    begsample = nearest(time, cfg.baselinewindow(1));
    endsample = nearest(time, cfg.baselinewindow(2));
    if begsample==endsample && ...
        ((begsample==1 && time(begsample)>cfg.baselinewindow(1)) || (begsample==numel(time) && time(begsample)<cfg.baselinewindow(2)))
      ft_warning('requested baselinewindow does not have any samples in the time axis, no baseline correction applied in this trial');
      dat = dat + datdiff; % add back the previously subtracted stuff
    else
      if begsample==endsample
        ft_warning('requested baselinewindow just has a single sample in the time axis')
      end
      dat = ft_preproc_baselinecorrect(dat, begsample, endsample);
    end
  end
end
if strcmp(cfg.dftfilter, 'yes')
  datorig = dat;
  optarg = {};
  if isfield(cfg, 'dftreplace')
    optarg = cat(2, optarg, {'dftreplace', cfg.dftreplace});
    if strcmp(cfg.dftreplace, 'neighbour') && (begpadding>0 || endpadding>0)
      ft_error('Padding by data mirroring is not supported for spectrum interpolation.');
    end
  end
  if isfield(cfg, 'dftbandwidth')
    optarg = cat(2, optarg, {'dftbandwidth', cfg.dftbandwidth});
  end
  if isfield(cfg, 'dftneighbourwidth')
    optarg = cat(2, optarg, {'dftneighbourwidth', cfg.dftneighbourwidth});
  end
  dat     = ft_preproc_dftfilter(dat, fsample, cfg.dftfreq, optarg{:});
  if strcmp(cfg.dftinvert, 'yes')
    dat = datorig - dat;
  end
end
if ~strcmp(cfg.hilbert, 'no')
  dat = ft_preproc_hilbert(dat, cfg.hilbert);
end
if strcmp(cfg.rectify, 'yes')
  dat = ft_preproc_rectify(dat);
end
if isnumeric(cfg.boxcar)
  numsmp = round(cfg.boxcar*fsample);
  if ~rem(numsmp,2)
    % the kernel should have an odd number of samples
    numsmp = numsmp+1;
  end
  % kernel = ones(1,numsmp) ./ numsmp;
  % dat    = convn(dat, kernel, 'same');
  dat = ft_preproc_smooth(dat, numsmp); % better edge behavior
end
if isnumeric(cfg.conv)
  kernel = (cfg.conv(:)'./sum(cfg.conv));
  if ~rem(length(kernel),2)
    kernel = [kernel 0];
  end
  dat = convn(dat, kernel, 'same');
end
if strcmp(cfg.derivative, 'yes')
  dat = ft_preproc_derivative(dat, 1);
end
if strcmp(cfg.absdiff, 'yes')
  % this implements abs(diff(data), which is required for jump detection
  dat = abs([diff(dat, 1, 2) zeros(size(dat,1),1)]);
end
if strcmp(cfg.standardize, 'yes')
  dat = ft_preproc_standardize(dat, 1, size(dat,2));
end
if ~isempty(cfg.subspace)
  dat = ft_preproc_subspace(dat, cfg.subspace);
end
if ~isempty(cfg.custom)
  if ~isfield(cfg.custom, 'nargout')
    cfg.custom.nargout = 1;
  end
  if cfg.custom.nargout==1
    dat = feval(cfg.custom.funhandle, dat, cfg.custom.varargin);
  elseif cfg.custom.nargout==2
    [dat, time] = feval(cfg.custom.funhandle, dat, cfg.custom.varargin);
  end
end
if strcmp(cfg.resample, 'yes')
  if ~isfield(cfg, 'resamplefs')
    cfg.resamplefs = fsample./2;
  end
  if ~isfield(cfg, 'resamplemethod')
    cfg.resamplemethod = 'resample';
  end
  [dat               ] = ft_preproc_resample(dat,  fsample, cfg.resamplefs, cfg.resamplemethod);
  [time, dum, fsample] = ft_preproc_resample(time, fsample, cfg.resamplefs, cfg.resamplemethod);
end
if ~isempty(cfg.precision)
  % convert the data to another numeric precision, i.e. double, single or int32
  dat = cast(dat, cfg.precision);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% remove the filter padding and do the preprocessing on the remaining trial data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if begpadding~=0 || endpadding~=0
  dat = ft_preproc_padding(dat, 'remove', begpadding, endpadding);
  if strcmp(cfg.demean, 'yes') || nargout>2
    time = ft_preproc_padding(time, 'remove', begpadding, endpadding);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ i_lt, i_gt ] = r8vec_sorted_split ( n, a, split )

%*****************************************************************************80
%
%% R8VEC_SORTED_SPLIT "splits" a sorted R8VEC, given a splitting value.
%
%  Discussion:
%
%    Given SPLIT, the routine seeks indices I_LT and I_GT so that
%
%      A(I_LT) < SPLIT < A(I_GT),
%
%    and if there are intermediate index values between I_LT and
%    I_GT, then those entries of A are exactly equal to SPLIT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters
%
%    Input, integer N, the number of entries in A.
%
%    Input, real A(N), a sorted array.
%
%    Input, real SPLIT, a value to which the entries in A are
%    to be compared.
%
%    Output, integer I_LT:
%    0 if no entries are less than SPLIT;
%    N if all entries are less than SPLIT;
%    otherwise, the index of the last entry in A less than SPLIT.
%
%    Output, integer I_GT:
%    1 if all entries are greater than SPLIT;
%    N+1 if no entries are greater than SPLIT;
%    otherwise the index of the first entry in A greater than SPLIT.
%
  if ( n < 1 )
    i_lt = -1;
    i_gt = -1;
    return
  end

  if ( split < a(1) )
    i_lt = 0;
    i_gt = 1;
    return
  end

  if ( a(n) < split )
    i_lt = n;
    i_gt = n + 1;
    return
  end

  lo = 1;
  hi = n;

  while ( 1 )

    if ( lo + 1 == hi )
      i_lt = lo;
      break
    end

    mid = round ( ( lo + hi ) / 2 );

    if ( split <= a(mid) )
      hi = mid;
    else
      lo = mid;
    end

  end

  for i = i_lt + 1 : n
    if ( split < a(i) )
      i_gt = i;
      return
    end
  end

  i_gt = n + 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function GCO_SetLabelOrder(Handle,Order)
% GCO_SetLabelOrder   Set label order for Expansion/Swap moves.
%   GCO_SetLabelOrder(Handle,Order) tells Expansion/Swap to select labels
%   in a specific order when Order contains integers from 1..NumLabels.
%   By default, Expansion/Swap use a consistent, prescribed order of labels
%   in until convergence.
%   Example:
%      GCO_SetLabelOrder(Handle,5:10); % only operate on labels 5..10
%      GCO_SetLabelOrder(Handle,randperm(NumLabels)); % random label order
%

GCO_LoadLib();
gco_matlab('gco_setlabelorder',Handle,int32(Order));
end

author: andrewssobral

repo: lrslibrary


text:
%IISUM Sum of integral image
%
% S = IISUM(II, U1, V1, U2, V2) is the sum of pixels in the rectangular image
% region defined by its top-left (U1,V1) and bottom-right (U2,V2).  II is
% a precomputed integral image.
%
% See also INTGIMAGE.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function s = iisum(ii, c1, r1, c2, r2)

    r1 = r1 - 1;
    if r1 < 1
        sA = 0;
        sB = 0;
    else
        sB = ii(r1,c2);
    end
    c1 = c1 - 1;
    if c1 < 1
        sA = 0;
        sC = 0;
    else
        sC = ii(r2,c1);
    end
    if (r1 >= 1) && (c1 >= 1)
        sA = ii(r1,c1);
    end

    s = ii(r2,c2) + sA -sB - sC;

author: petercorke

repo: machinevision-toolbox-matlab


text:
function color=gen_color(curve_idx)

colors=[];
colors{end+1}='r';
colors{end+1}='b';
colors{end+1}='g';
colors{end+1}='m';
colors{end+1}='c';
colors{end+1}='black';
colors{end+1}=[0.7 0 0.7 ];
colors{end+1}=[0 0.7 0.7 ];
colors{end+1}=[ 0.83 0.33 0];

sel_idx=mod(curve_idx-1, length(colors))+1;
color=colors{sel_idx};

end

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function [favar]=favar_irfdisp(favar,IRFperiods,endo,IRFt,strctident,pref)

% we arranged everything before including re-transformation
% pick the relevant shocks (plotXshock) and we are ready to plot
favar_irf_estimates=favar.IRF.favar_irf_estimates(:,favar.IRF.plotXshock_index);

% plot the IRFs of the factors
if favar.IRFplot==1
    % one window for each shock, as we possibly plot a lot of variables in X
    for px=1:favar.IRF.npltXshck
        irf_favar=figure('Tag','BEARresults');
        % shock label
        if IRFt==1||IRFt==2||IRFt==3
            printlabels=endo{favar.IRF.plotXshock_index(1,px),1};
        elseif IRFt==4||IRFt==5||IRFt==6
            printlabels=strctident.signreslabels{favar.IRF.plotXshock_index(1,px)};
        end
        set(irf_favar,'name',['approximate impulse response functions (FAVAR) to shock ',printlabels]);

        numcol=ceil(sqrt(favar.npltX)); % rounded square root of npltX, make the plot quadratic
        numrow=ceil(favar.npltX/numcol); % and the number of rows we need
        startIRF=(px*favar.npltX)-favar.npltX+1;
        stopIRF=px*favar.npltX;
        count=0;
        for ii=startIRF:stopIRF
            count=count+1;
            subplot(numrow,numcol,count);
            hold on
            Xpatch=[(1:IRFperiods) (IRFperiods:-1:1)];
            Ypatch=[favar_irf_estimates{ii}(1,:) fliplr(favar_irf_estimates{ii}(3,:))];
            IRFpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
            set(IRFpatch,'facealpha',0.5);
            set(IRFpatch,'edgecolor','none');
            plot(favar_irf_estimates{ii}(2,:),'Color',[0.4 0.4 1],'LineWidth',2);
            plot([1,IRFperiods],[0 0],'k--');
            hold off
            %         minband=min(favar.IRF.irf_factors_final_plotX{ii}(1,:));
            %         maxband=max(favar.IRF.irf_factors_final_plotX{ii}(3,:));
            %         space=maxband-minband;
            %         Ymin=minband-0.2*space;
            %         Ymax=maxband+0.2*space;
            set(gca,'XLim',[1 IRFperiods],'FontName','Times New Roman'); %,'YLim',[Ymin Ymax]
            % title of subplot is variable name
            title(favar.informationvariablestrings{1,favar.plotX_index(count)},'FontWeight','normal','interpreter','none');
        end

        % top supertitle
        ax=axes('Units','Normal','Position',[.11 .075 .85 .88],'Visible','off');
        set(get(ax,'Title'),'Visible','on')
        title(['Shock: ',printlabels],'FontSize',11,'FontName','Times New Roman','FontWeight','normal','interpreter','none');
    end
end

%% save in excel
if pref.results==1
    bear.data.favar_excelrecord4fcn(favar, IRFperiods, endo, IRFt, strctident, favar_irf_estimates, pref)
end

author: european-central-bank

repo: BEAR-toolbox


text:
function y=num2bit(th_element,b,thetmin,thetmax)
%This function takes a real number and converts it to a b-bit
%representation. The value b is chosen in GAbit_roulette to satisfy
%the decimal accuracy requirements given in M (as in step 1 of encoding
%process of Subsection 9.3.2 of ISSO).
%
d=(thetmax-thetmin)/(2^b-1);
rnd=round((th_element-thetmin)/d);
y=zeros(1,b);
for i=b:-1:1
   ratio=rnd/(2^(i-1));
   if ratio >= 1
      y(b-i+1)=1;
      rnd=rnd-2^(i-1);
   else
   end
end


text:
function WJ_DistortionCorrection_3D_SItriggered(src,evt,varargin)
            global filename
            filename=[src.hSI.hScan2D.logFilePath,'\', src.hSI.hScan2D.logFileStem,'_',num2str(src.hSI.hScan2D.logFileCounter-1,'%05d'),'.tif'];
            disp(['Just saved raw data: ',filename]);
            disp('Start the distortion correction....');
            disp('Reloading the raw data....');

%             filename='97045_20210308_ML-400_AL-400_1Openfiled_00001.tif';

            [header, RawImage, imgInfo] = scanimage.util.opentif(filename);
            filename1=imgInfo.filename;
            disp(['Rawe 2P imaging data: "' filename1  '" is loaded']);
            save([filename(1:end-4),'_tifHeader.mat'],'header');% save the SI information, only necessasy for SI information saving. No need for DJ analysis
            disp([filename(1:end-4),'_tifHeader.mat',' is saved']);
            save([filename(1:end-4),'_imgInfo.mat'],'imgInfo');% save the imaging information,only necessasy for SI information saving. No need for DJ analysis
            disp([filename(1:end-4),'_imgInfo.mat',' is saved']);
            %
            Dimention=size(RawImage);
            Zs=header.SI.hStackManager.zs;
            % load in 3D transform matrx according to different frame szie
            if Dimention(1)==512
                % read 515 3D Tranform matrix
                TransformMatrix_tem=load('\\forskning.it.ntnu.no\ntnu\mh-kin\moser\open2pmini\FOV callibration\RubenWeijianNienke\MINI2P_L_WJ001\Report\TransMatrix_512_3D.mat');
                % read 512 3D FOV information
                FOV=csvread('\\forskning.it.ntnu.no\ntnu\mh-kin\moser\open2pmini\FOV callibration\RubenWeijianNienke\MINI2P_L_WJ001\Report\FOVreport_512_20211111.csv');
                Depth=csvread('\\forskning.it.ntnu.no\ntnu\mh-kin\moser\open2pmini\FOV callibration\RubenWeijianNienke\MINI2P_L_WJ001\Report\Depth20211111.csv');
                %
            elseif Dimention(1)==256
                % read 256 3D Tranform matrix
                TransformMatrix_tem=load('\\forskning.it.ntnu.no\ntnu\mh-kin\moser\open2pmini\FOV callibration\RubenWeijianNienke\MINI2P_L_WJ001\Report\TransMatrix_256_3D.mat');
                % read 256 3D FOV information
                FOV=csvread('\\forskning.it.ntnu.no\ntnu\mh-kin\moser\open2pmini\FOV callibration\RubenWeijianNienke\MINI2P_L_WJ001\Report\FOVreport_256_20211111.csv');
                Depth=csvread('\\forskning.it.ntnu.no\ntnu\mh-kin\moser\open2pmini\FOV callibration\RubenWeijianNienke\MINI2P_L_WJ001\Report\Depth20211111.csv');

                %
            else
            end
            %
            TransformMatrix=struct2cell(TransformMatrix_tem);

            PlaneMark=[];
            if length(Dimention)<=3
                ImageStack_corrected=RawImage;
                PlaneMark=ones(size(ImageStack_corrected(:)),1);
                disp('Image Stack Size:');
                disp(num2str(Dimention));
                disp('No HyperStack Merging is applied');
            elseif length(Dimention)==4
                disp('Image Stack Size:');
                disp(num2str(Dimention));
                disp('2D Timelapse imaging was infered.  Channel and Frame merging is applied');
                ImageStack_corrected=zeros(Dimention(1),Dimention(2),Dimention(3)*Dimention(4));
                k=1;
                for j=1:1:Dimention(4)
                    for i=1:1:Dimention(3)
                        ImageStack_corrected(:,:,k)=RawImage(:,:,i,j);
                        PlaneMark(k)=1;
                        k=k+1;
                    end
                end
            elseif length(Dimention)==5
                disp('Image Stack Size:');
                disp(num2str(Dimention));
                disp('Multi-layer sinlge TimePoint imaging was infered.  Channel, Frame and Plane merging is applied');
                ImageStack_corrected=zeros(Dimention(1),Dimention(2),Dimention(3)*Dimention(4)*Dimention(5));
                k=1;
                    for m=1:1:Dimention(5)
                        for j=1:1:Dimention(4)
                            for i=1:1:Dimention(3)
                                ImageStack_corrected(:,:,k)=RawImage(:,:,i,j,m);
                                PlaneMark(k)=m;
                                k=k+1;
                            end
                        end
                    end
            else
                disp('Image Stack Size:');
                disp(num2str(Dimention));
                disp('Multi-layer Timelapse imaging was infered.  Channel, Frame, Plane and Volume merging is applied');
                ImageStack_corrected=zeros(Dimention(1),Dimention(2),Dimention(3)*Dimention(4)*Dimention(5)*Dimention(6));
                k=1;

                for n=1:1:Dimention(6)
                    for m=1:1:Dimention(5)
                        for j=1:1:Dimention(4)
                            for i=1:1:Dimention(3)
                                ImageStack_corrected(:,:,k)=RawImage(:,:,i,j,m,n);
                                PlaneMark(k)=m;
                                k=k+1;
                            end
                        end
                    end
                end
            end
            clear RawImage
            % find which matrix to use
            PlaneClosest=PlaneMark;
            for i=1:1:size(PlaneMark,2)
                PlaneClosest(i)=Zs(PlaneMark(i));
                [minValue,closestIndex]=min(abs(Depth-PlaneClosest(i)));
                PlaneClosest(i)= closestIndex;
            end

            fTIF = DataIO.Fast_BigTiff_Write([filename(1:end-4),'_wrappiing-corrected.tif'],1*10000/FOV(3),0);
            for i=1:1:size(ImageStack_corrected,3)
                ImageStack_corrected(:,:,i)=imwarp(ImageStack_corrected(:,:,i),TransformMatrix{PlaneClosest(i),1},'OutputView',imref2d(size(ImageStack_corrected(:,:,1))));
                fTIF.WriteIMG(int16(ImageStack_corrected(:,:,i)'));
                disp(['Unwrapping and saving image number ',num2str(i)]);
            end
            fTIF.close;
            disp(['Distortion correction finished, data saved: ',filename(1:end-4),'_wrappiing corrected.tif']);
%             clear ImageStack_corrected
end

author: kavli-ntnu


text:
function [stiffness]=assemble(stiffness,k,index)
%----------------------------------------------------------
%  Purpose:
%     Assembly of element matrices into the system matrix
%
%  Synopsis:
%     [stiffness]=assemble(stiffness,k,index)
%
%  Variable Description:
%     stiffness - system matrix
%     k  - element matri
%     index - d.o.f. vector associated with an element
%-----------------------------------------------------------

 edof = length(index);
 for i=1:edof
   ii=index(i);
     for j=1:edof
       jj=index(j);
         stiffness(ii,jj)=stiffness(ii,jj)+k(i,j);
     end
 end


text:
function model = leReconstruct(leInfo, y)

% LERECONSTRUCT Reconstruct an LE form component parts.
% FORMAT
% DESC takes component parts of an LE model and reconstructs the
% LE model. The component parts are normally retrieved from a
% saved file.
% ARG leInfo : the active set and other information stored in a structure.
% ARG y : the output target training data for the LE.
% RETURN model : an LE model structure that combines the component
% parts.
%
% SEEALSO : leDeconstruct, leCreate
%
% COPYRIGHT : Neil D. Lawrence, 2009

% MLTOOLS

  model = leInfo;
  model.Y = y;

end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [bounds, best_cover] = pascal_overlap_bound(orig_sp, sp_parts, gt, gt_ind)
%function bounds = pascal_overlap_bound(orig_sp, sp_parts, gt, gt_ind)
% Computes theoretical upper bounds of Pascal overlap score for given
% initial superpixelation sp and ground truth segmentation gt.

assert(max(gt_ind) == length(gt_ind)); % This function assumes that gt_ind = 1:n, where n >= 1.

% remove empty superpixels
% for q = length(sp):-1:1
%     if size(sp{q}.pixels,1) == 0
%         sp(q) = [];
%     end
% end

scores = zeros(max(gt_ind),length(sp_parts));

for i = 1:length(sp_parts) % for each sp
    %pix = double(sp{i}.pixels); % using double() here is extremely important. Otherwise sub2ind below caps variable s at 65536, breaking everything.
    %s{i} = sub2ind([h,w], pix(:,1), pix(:,2)); % linear indices of pixels of current superpixel in the image

    indl = 0;
    for sus = sp_parts{i}
        indl = indl + length(orig_sp{sus}.spind);
    end

    s{i} = zeros(1, indl);
    indl = 0;
    for sus = sp_parts{i} % for each original sp part of the current sp
        s{i}(indl+1:indl+length(orig_sp{sus}.spind)) = orig_sp{sus}.spind;
        indl = indl + length(orig_sp{sus}.spind);
    end

    sp_classes = gt(s{i}); % classes of pixels in current sp
    scores(:,i) = histc(sp_classes, gt_ind)/length(sp_classes);

end % for each superpixel

a = zeros(1, max(gt_ind)); % intersection pixel count
b = zeros(1, max(gt_ind)); % union pixel count
best_cover = [];

for ind = gt_ind % for each object
    q = find(gt == ind);
    gt_seg{ind} = q;
    b(ind) = length(q); % size of gt object, initial size of union
end

for cl = gt_ind % for each object
    [cl_scores, scoreperm] = sort(scores(cl,:), 'descend'); % cl_scores = a*/(a* + b*)
    cl_scoresx = 1./(1./cl_scores - 1); % cl_scoresx = a*/b*

    for i = 1:length(cl_scores)
        % condition for a/b <= (a + a*)/(b + b*) is that a/b <= a*/b*. Since
        % always a <= b, having a* >= b*, i.e., scores >= 0.5 is enough too.
        if cl_scores(i) >= 0.5 || cl_scoresx(i) >= double(a(cl))/b(cl)
            a(cl) = a(cl) + length(intersect(s{scoreperm(i)}, gt_seg{cl})); % size of intersection increases by this value
            b(cl) = b(cl) + length(setdiff(s{scoreperm(i)}, gt_seg{cl})); % size of union increases by this value
        else
            best_cover{cl} = scoreperm(1:i);
            break; % scores from now on are too low to increase the overlap ratio, stop.
        end

    end % for each sp
end % for each object

bounds = double(a)./b;

%toc

author: Cloud-CV

repo: object-proposals


text:
function write_structure_to_db(conn,patient_id)
%function write_structure_to_db(conn,patient_id)
%
%Input: Writes all structures from global planC to database
%
%APA, 02/27/2011

global planC

indexS = planC{end};

columnNames = {'structure_id','patient_id','scan_id','structure_name','structure_uid','structure_color'};

%Loop over all scans, find structures belonging to scan and add to DB.
numStructs = length(planC{indexS.structures});
allStructsV = 1:numStructs;

for scanNum = 1:length(planC{indexS.scan})

    %Find matching scan in DB
    scanUID = planC{indexS.scan}(scanNum).scanUID;
    sqlq_find_scan = ['Select scan_id from scan where scan_uid = ''', scanUID,''''];
    scan_raw = exec(conn, sqlq_find_scan);
    scan = fetch(scan_raw);
    scan = scan.Data;
    if ~isstruct(scan)
        warning('Scan does not exist in DB. Cannnot write structures. Proceed to next scan')
    else
        scan_id = scan.scan_id{1};
    end

    %Loop over structures belonging to this scan
    assocScansV = getStructureAssociatedScan(allStructsV, planC);

    matchingScansV = find(assocScansV == scanNum);

    indMatchStructsV = allStructsV(matchingScansV);

    for structNum = indMatchStructsV

        structS = planC{indexS.structures}(structNum);

        %Find matching structure in DB
        whereclause = {['where structure_uid = ''', structS.strUID,'''']};
        sqlq_find_str = ['Select structure_id from structure where structure_uid = ''', structS.strUID,''''];
        str_raw = exec(conn, sqlq_find_str);
        str = fetch(str_raw);
        str = str.Data;
        if ~isstruct(str)
            structure_id = char(java.util.UUID.randomUUID);
            isNewRecord = 1;
        else
            structure_id = str.structure_id{1};
            isNewRecord = 0;
        end

        %structure_id
        recC{1} = structure_id;

        %patient_id
        recC{2} = patient_id;

        %scan_id
        recC{3} = scan_id;

        %structureName
        recC{4} = structS.structureName;

        %structure_uid
        recC{5} = structS.strUID;

        %structure_color
        recC{6} = mat2str(structS.structureColor);

        if isNewRecord
            insert(conn,'structure',columnNames,recC);
        else
            update(conn,'structure',columnNames,recC,whereclause);
        end

        pause(0.05)

    end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%MDL_FANUC10L  Create kinematic model of Fanuc AM120iB/10L robot
%
% MDL_FANUC10L is a script that creates the workspace variable R which
% describes the kinematic characteristics of a Fanuc AM120iB/10L robot
% using standard DH conventions.
%
% Also defines the workspace vector:
%   q0   mastering position.
%
% Notes::
% - SI units of metres are used.
%
% Author::
%  Wynand Swart,
%  Mega Robots CC, P/O Box 8412, Pretoria, 0001, South Africa
%  wynand.swart@gmail.com
%
% See also SerialLink, mdl_irb140, mdl_m16, mdl_motomanHP6, mdl_puma560.

% MODEL: Fanuc, AM120iB/10L, 6DOF, standard_DH

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.

%Cell: 073-1555-430
%30 Sep 2007
%Fanuc AM120iB/10L robot

%            theta    d      a    alpha
clear L
L(1) = Link([0        0     0.15  -pi/2   0  ]);
L(2) = Link([0        0     0.77   pi     0  ]);
L(3) = Link([0        0     0.1   -pi/2   0  ]);
L(4) = Link([0       -0.96  0      pi/2   0  ]);
L(5) = Link([0        0     0     -pi/2   0  ]);
L(6) = Link([0       -0.1   0      0      0  ]);
%##########################################################
%Pose 0; At MASTERING position;
%##########################################################
q0 =[0   -pi/2   0   0   0   0];
R=SerialLink(L, 'name', 'Fanuc AM120iB/10L');
%##########################################################

author: Allopart

repo: rbpf-gmapping


text:
function H = discriminant(f, x, y, varargin)
%DISCRIMINANT the determinant of Hessian of a CHEBFUN2 at (x,y)
%   H = DISCRIMINANT(F,x,y) returns the determinant of the Hessian of F at
%   (x,y). The gradient of F should be zero at (x,y).
%
%   H = DISCRIMINANT(F,G,x,y) returnes the determinant of the 'border' Hessian
%   of F at (x,y).
%
%   Note that we cannot represent the Hessian matrix because we do not allow
%   horizontal concatenation of CHEBFUN2 objects.
%
% See also JACOBIAN.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty( f ) )
    H = [];
    return
end

% Mixed second partial derivatives:
fxx = diff(f, 2, 2);
fyy = diff(f, 2, 1);
fxy = diff(f, [1 1]);

if ( nargin == 3 )      % Standard Hessian
    % Evaluate at (x,y):
    fxx = feval(fxx, x, y);
    fxy = feval(fxy, x, y);
    fyy = feval(fyy, x, y);
    H = fxx.*fyy - fxy.^2;

elseif ( nargin == 4 )  % Bordered Hessian
    % Parse user inputs:
    g = x;
    x = y;
    y = varargin{1};

    % Evaluate at (x,y):
    fxx = feval(fxx, x, y);
    fyy = feval(fxy, x, y);
    fxy = feval(fyy, x, y);

    % Partial diff and evaluate:
    gx = diff(g, 1, 2);
    gy = diff(g, 1, 1);
    gx = feval(gx, x, y);
    gy = feval(gy, x, y);

    % Determinant of bordered Hessian.
    H = -gx.*(gx.*fyy - gy.*fxy) + gy.*(gx.*fxy - gy.*fxx);

else
    error('CHEBFUN:CHEBFUN2:discriminant:badInput', 'Invalid input arguments.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function nlp = ZMPCost(nlp, varargin)

    w = varargin{1};

    plant = nlp.Plant;
    force = plant.Inputs.ConstraintWrench.ffoot;

    ws = SymVariable('w');

    zmp = -force(3)./force(2);
    constr = ws*(zmp+0.1).^2;
    m_fun = SymFunction(['zmp_cost_' plant.Name],tovector(constr),{force},{ws});
    addRunningCost(nlp,m_fun,{'ffoot'},{w});
end

author: ayonga

repo: frost-dev

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


text:
function  a = randint(m,n,a,b)
        %RANDINT  Randomly generated integral matrix.
        %         randint(m,n) returns an m-by-n such matrix with entries
        %         between 0 and 9.
        %         rand(m,n,a,b) return entries between integers  a  and  b .
        if nargin < 3, a = 0; b = 9; end
        a = floor((b-a+1)*rand(m,n)) + a;


text:
function c = rdivide(a,b)
%RDIVIDE      Interval elementwise right division a ./ b
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump  modified for infinity
% modified 06/06/98     S.M. Rump  modified for NaN+Nan*i
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    remove check for 'double'
%                                    take care of Matlab sparse Inf/NaN bug
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/03/05     S.M. Rump  sparse flag corrected
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 02/11/06     S.M. Rump  SparseInfNanFlag removed
%                                    improved performance
% modified 05/23/06     S.M. Rump  sparse Inf/NaN bug corrected in Version 7.2+
% modified 12/03/06     S.M. Rump  Sparse Bug global flag (thanks to Arnold)
% modified 10/23/07     S.M. Rump  complex numbers
% modified 02/18/09     S.M. Rump  NaN performance improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
  end

  % no need to take care about huge matrices: result would be almost full anyway
  % also no care necessary about previous Matlab sparse NaN bug (would be helpful, in fact)
  % make sure a full except b is scalar, b full anyway
  if isa(b,'intval')
    bcomplex = b.complex;
    if bcomplex
      b.mid = full(b.mid);
      b.rad = full(b.rad);
      makefull = ( prod(size(b.mid))~=1 );
      nanindex = isnan(b.mid) | isnan(b.rad);
    else
      b.inf = full(b.inf);
      b.sup = full(b.sup);
      makefull = ( prod(size(b.inf))~=1 );
      nanindex = isnan(b.inf) | isnan(b.sup);
    end
  else
    bcomplex = ~isreal(b);
    b = full(b);
    makefull = ( prod(size(b))~=1 );
    nanindex = isnan(b);
  end
  nanindex = sparse(nanindex);      % careful: full(True) | sparse = full
  if isa(a,'intval')
    acomplex = a.complex;
    if acomplex
      if makefull
        a.mid = full(a.mid);
        a.rad = full(a.rad);
      end
      nanindex = nanindex | isnan(a.mid) | isnan(a.rad);
    else
      if makefull
        a.inf = full(a.inf);
        a.sup = full(a.sup);
      end
      nanindex = nanindex | isnan(a.inf) | isnan(a.sup);
    end
  else
    acomplex = ~isreal(a);
    if makefull
      a = full(a);
    end
    nanindex = nanindex | isnan(a);
  end
  anynanindex = any(nanindex);
  anynanindex = any(anynanindex(:));

  ws = warning;
  warning off

  if acomplex | bcomplex                % numerator complex
    b = intval(b);                      % make sure b is interval
    if ~bcomplex                        % denominator is real
      c = a.*(1./b);
      return
    end
    x = real(b.mid);                    % denominator is complex
    y = imag(b.mid);
    setround(-1)
    Ninf = x.*x + y.*y + (-b.rad).*b.rad;
    index = ( Ninf<=0 );
    setround(1)
    Nsup = x.*x + y.*y + (-b.rad).*b.rad;
    x2 = max( x./Ninf , x./Nsup );
    y2 = max( y./Ninf , y./Nsup );
    setround(-1)
    x1 = max( x./Ninf , x./Nsup );
    y1 = max( y./Ninf , y./Nsup );
    c1 = x1 - j*y2;
    setround(1)
    c2 = x2 - j*y1;
    binv.complex = 1;
    binv.inf = [];
    binv.sup = [];
    binv.mid = c1 + 0.5*(c2-c1);
    binv.rad = abs( binv.mid - c1 ) + b.rad./Ninf;
    index = index | ( binv.rad<0 );
    if any(index(:))                          % division by zero
      binv.mid(index) = complex(NaN,NaN);
      binv.rad(index) = NaN;
    end
    binv = class(binv,'intval');
    c = a.*binv;
    if anynanindex
      c.mid(nanindex) = NaN;
      c.rad(nanindex) = NaN;                  % radius for sparse cannot be 0
    end
  else                                        % both a and b real
    c.complex = 0;
    if ~isa(a,'intval')                       % R ./ IR
      % be sure min/max works correct for zero upper bounds in b
      b.sup(b.sup==0) = -0;
      setround(-1)
      c.inf = min( a./b.inf , a./b.sup );
      setround(1)
      c.sup = max( a./b.inf , a./b.sup );
      index = ( b.inf<=0 ) & ( b.sup>=0 );    % 0 in b
      if ~isempty(find(index))
        c.inf(index) = -inf;
        c.sup(index) =  inf;
        index = index & ( a==0 );             % 0/0
        if any(index(:))
          c.inf(index) = NaN;
          c.sup(index) = NaN;
        end
        if anynanindex
          c.inf(nanindex) = NaN;
          c.sup(nanindex) = NaN;
        end
      end
    elseif ~isa(b,'intval')                     % IR ./ R
      setround(-1)
      c.inf = min( a.inf./b , a.sup./b );
      setround(1)
      c.sup = max( a.inf./b , a.sup./b );
      index = ( b==0 );                         % numerator/0
      if ~isempty(find(index))
        c.inf(index) = -inf;
        c.sup(index) =  inf;
        index = index & ( a.inf<=0 ) & ( 0<=a.sup );
        if any(index(:))
          c.inf(index) = NaN;
          c.sup(index) = NaN;
        end
        if anynanindex
          c.inf(nanindex) = NaN;
          c.sup(nanindex) = NaN;
        end
      end
    else                                        % IR ./ IR
      % be sure min/max works correct for zero upper bounds in b
      b.sup(b.sup==0) = -0;
      setround(-1)
      c.inf = min( a.inf./b.inf , a.inf./b.sup );
      c.inf = min( c.inf , a.sup./b.inf );
      c.inf = min( c.inf , a.sup./b.sup );
      setround(1)
      c.sup = max( a.inf./b.inf , a.inf./b.sup );
      c.sup = max( c.sup , a.sup./b.inf );
      c.sup = max( c.sup , a.sup./b.sup );
      index = ( b.inf<=0 ) & ( b.sup>=0 );      % 0 in b
      if ~isempty(find(index))
        if prod(size(b.inf))==1
          c.inf = -inf*ones(size(c.inf));
          c.sup = -c.inf;
        else
          c.inf(index) = -inf;
          c.sup(index) =  inf;
        end
        index = index & ( a.inf<=0 ) & ( a.sup>=0 );   % 0./0
        if any(index(:))
          c.inf(index) = NaN;
          c.sup(index) = NaN;
        end
      end
      if anynanindex
        c.inf(nanindex) = NaN;
        c.sup(nanindex) = NaN;
      end
    end
    c.mid = [];
    c.rad = [];
    c = class(c,'intval');
  end

  if issparse(b) & ( prod(size(c))~=1 )
    c = sparse(c);
  end

  warning(ws)
  setround(rndold)

repo: OpenMAS


text:
function [AF,BF]=framebounds(F,varargin);
%FRAMEBOUNDS  Frame bounds
%   Usage: fcond=framebounds(F);
%          [A,B]=framebounds(F);
%          [...]=framebounds(F,Ls);
%
%   `framebounds(F)` calculates the ratio $B/A$ of the frame bounds of the
%   frame given by *F*. The length of the system the frame bounds are
%   calculated for is given by `L=framelength(F,1)`.
%
%   `framebounds(F,Ls)` additionally specifies a signal length for which
%   the frame should work. The actual length used is `L=framelength(F,Ls)`.
%
%   `[A,B]=framebounds(F)` returns the frame bounds *A* and *B* instead of
%   just their ratio.
%
%
%   'framebounds` accepts the following optional parameters:
%
%     'fac'        Use a factorization algorithm. The function will throw
%                  an error if no algorithm is available.
%
%     'iter'       Call `eigs` to use an iterative algorithm.
%
%     'full'       Call `eig` to solve the full problem.
%
%     'auto'       Choose the `fac` method if possible, otherwise
%                  use the `full` method for small problems and the
%                  `iter` method for larger problems.
%                  This is the default.
%
%     'crossover',c
%                  Set the problem size for which the 'auto' method
%                  switches between `full` and `iter`. Default is 200.
%
%   The following parameters specifically related to the `iter` method:
%
%     'tol',t      Stop if relative residual error of eighs is less than the
%                  specified tolerance. Default is 1e-9
%
%     'maxit',n    Do at most n iterations in eigs. Default is 100.
%
%     'pcgtol',t   Stop if relative residual error of pcg is less than the
%                  specified tolerance. Default is 1e-6
%
%     'pcgmaxit',n Do at most n iterations in pcg. Default is 150.
%
%     'p',p        The number of Lanzcos basis vectors to use.  More vectors
%                  will result in faster convergence, but a larger amount of
%                  memory.  The optimal value of `p` is problem dependent and
%                  should be less than *L*.  The default value chosen
%                  automatically by eigs.
%
%     'print'      Display the progress.
%
%     'quiet'      Don't print anything, this is the default.
%
%   See also: frame, framered

complainif_notenoughargs(nargin,1,'FRAMEBOUNDS');
complainif_notvalidframeobj(F,'FRAMEBOUNDS');

% We handle the container frames first
  if strcmp(F.type,'fusion')
      AF=0;
      BF=0;
      for ii=1:F.Nframes
          [A,B]=framebounds(F.frames{ii},varargin{:});
          AF=AF+(A*F.w(ii)).^2;
          BF=BF+(B*F.w(ii)).^2;
      end;
      AF=sqrt(AF);
      BF=sqrt(BF);

      return;
  end;

  if strcmp(F.type,'tensor')
    AF=1;
    BF=1;
    for ii=1:F.Nframes
      [A,B]=framebounds(F.frames{ii},varargin{:});
      AF=AF*A;
      BF=BF*B;
    end;

    return;
  end;

  definput.keyvals.Ls=1;
  definput.keyvals.maxit=100;
  definput.keyvals.tol=1e-9;
  definput.keyvals.pcgmaxit=150;
  definput.keyvals.pcgtol=1e-6;
  definput.keyvals.crossover=200;
  definput.keyvals.p=[];
  definput.flags.print={'quiet','print'};
  definput.flags.method={'auto','fac','iter','full'};

  [flags,kv]=ltfatarghelper({'Ls'},definput,varargin);

  F=frameaccel(F,kv.Ls);
  L=F.L;

  % Default values, works for the pure frequency transforms.
  AF=1;
  BF=1;

  % Simple heuristic: If F.g is defined, the frame uses windows.
  if isfield(F,'g')
      if isempty(F.g)
          error('%s: No analysis frame is defined.', upper(mfilename));
      end;
      g=F.g;
      op    = @frana;
      opadj = @frsyn;
  end;

  F_isfac = isfield(F,'isfac') && F.isfac;

  if flags.do_fac && ~F_isfac
    error('%s: The type of frame has no factorization algorithm.',upper(mfilename));
  end;

  if (flags.do_auto && F_isfac) || flags.do_fac
    switch(F.type)
     case 'gen'
      V=svd(g);
      AF=min(V)^2;
      BF=max(V)^2;
     case {'dgt','dgtreal'}
      [AF,BF]=gabframebounds(g,F.a,F.M,L);
     case {'dwilt','wmdct'}
      [AF,BF]=wilbounds(g,F.M,L);
     case {'filterbank','ufilterbank'}
      [AF,BF]=filterbankbounds(g,F.a,L);
     case {'filterbankreal','ufilterbankreal'}
      [AF,BF]=filterbankrealbounds(g,F.a,L);
     case 'fwt'
      [AF,BF]=wfbtbounds({g,F.J,'dwt'},L);
     case 'wfbt'
      [AF,BF]=wfbtbounds(g,L);
     case 'ufwt'
      [AF,BF]=wfbtbounds({g,F.J,'dwt'},L,F.flags.scaling);
     case 'uwfbt'
      [AF,BF]=wfbtbounds(g,L,F.flags.scaling);
     case 'wpfbt'
      [AF,BF]=wpfbtbounds(g,L,F.flags.interscaling);
     case 'uwpfbt'
      [AF,BF]=wpfbtbounds(g,L,F.flags.interscaling,F.flags.scaling);
    end;
  end;

  if (flags.do_auto && ~F_isfac && F.L>kv.crossover) || flags.do_iter

    if flags.do_print
      opts.disp=1;
    else
      opts.disp=0;
    end;
    opts.isreal = F.realinput;
    opts.maxit  = kv.maxit;
    opts.tol    = kv.tol;
    opts.issym  = 0;
    if ~isempty(kv.p)
       opts.p      = kv.p;
    end

    pcgopts.maxit = kv.pcgmaxit;
    pcgopts.tol = kv.pcgtol;

    % Upper frame bound
    frameop = @(x) F.frsyn(F.frana(x));
    BF = real(eigs(frameop,L,1,'LM',opts));

    % Lower frame bound
    frameop2 = @(x) F.frsyn(F.frana(x));
    invfrop = @(x) pcgwrapper(frameop2,x,pcgopts.tol,pcgopts.maxit);

    % Test convergence of pcg
    test = randn(L,1);
    if ~F.realinput, test = test +1i*randn(L,1); end
    [~,flag] = invfrop(test);

    % If PCG converges, estimate the smallest eigenvalue, otherwise assume
    % AF = 0;
    if ~flag
        AF = real(eigs(invfrop,L,1,'SM',opts));
    else
        AF = 0;
    end

  end;

  if (flags.do_auto && ~F_isfac && F.L<=kv.crossover) || flags.do_full
    % Compute thee transform matrix.
    bigM=opadj(F,op(F,eye(L)));

    D=eig(bigM);

    % Clean the eigenvalues, we know they are real
    D=real(D);
    AF=min(D);
    BF=max(D);
  end;

  if nargout<2
    % Avoid the potential warning about division by zero.
    if AF==0
      AF=Inf;
    else
      AF=BF/AF;
    end;
  end;

end

% In order to mute pprinting out pcg progress
function [y,flag]=pcgwrapper(varargin)
[y,flag,dummyrelres,dummyiter] = pcg(varargin{:});
end

% The function has been written in this way, because Octave (at the time
% of writing) does not accept additional parameters at the end of the
% line of input arguments for eigs
function y=afun(x,F_in,op_in,opadj_in)
  persistent F;
  persistent op;
  persistent opadj;

  if nargin>1
    F     = F_in;
    op    = op_in;
    opadj = opadj_in;
  else
    y=opadj(F,op(F,x));
  end;

end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [en,U,lam] = kentropy(X, A1, A2, A3, A4)
% Quadratic Renyi Entropy for a kernel based estimator
%
% Given the eigenvectors and the eigenvalues of the kernel matrix, the entropy is computed by
%
% >> H = kentropy(X, U, lam)
%
% The eigenvalue decomposition can also be computed (or
% approximated) implicitly:
%
% >> H = kentropy(X, kernel, sig2)
%
%
% Full syntax
%
% >> H = kentropy(X, kernel, kernel_par)
% >> H = kentropy(X, kernel, kernel_par, type)
% >> H = kentropy(X, kernel, kernel_par, type, nb)
%
%       Outputs
%         H          : Quadratic Renyi entropy of the kernel matrix
%       Inputs
%         X          : N x d matrix with the training data
%         kernel     : Kernel type (e.g. 'RBF_kernel')
%         kernel_par : Kernel parameter (bandwidth in the case of the 'RBF_kernel')
%         type(*)    : 'eig'(*), 'eigs', 'eign'
%         nb(*)      : Number of eigenvalues/eigenvectors used in the eigenvalue decomposition approximation
%
%
% >> H = kentropy(X, U, lam)
%
%       Outputs
%         H   : Quadratic Renyi entropy of the kernel matrix
%       Inputs
%         X   : N x d matrix with the training data
%         U   : N x nb matrix with principal eigenvectors
%         lam : nb x 1 vector with eigenvalues of principal components
%
% See also:
%   kernel_matrix, RBF_kernel, demo_fixedsize

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

n= size(X,1);

if isstr(A1), % kernel_matrix

  kernel = A1;
  kernel_par = A2;
  eval('etype = A3;','etype =''eig'';');
  if ~(strcmp(etype, 'eig') |strcmp(etype, 'eigs') |strcmp(etype,'eign')),
    error('type has to be ''eig'', ''eigs'' or ''eign''...');
  end
  eval('nb = A4;',' ');

  if strcmp(etype,'eign'),
    eval('[U,lam] = eign(X,kernel,kernel_par,nb);','[U,lam] = eign(X,kernel,kernel_par);');
  else
    omega = kernel_matrix(X, kernel, kernel_par);
    eval('[U,lam] = feval(etype,omega,nb);','[U,lam] = feval(etype,omega);');
    if size(lam,1)==size(lam,2), lam = diag(lam); end
    %onen = ones(n,1)./n; en = -log(onen'*omega*onen);
  end

else
  U = A1;
  lam = A2;
end
en = -log((sum(U,1)/n).^2 * lam);

author: peterhcharlton

repo: RRest


text:
function bc = specific_bc(xbd,ybd)
%solution3_bc   Reference problem 1.3  boundary condition
%   bc = specific_bc(xbd,ybd);
%   input
%          xbd          x boundary coordinate vector
%          ybd          y boundary coordinate vector
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bc=2*(1+ybd)./((ybd+1).*(ybd+1) + (xbd+3).*(xbd+3));
return


text:
%Program for one sided image reflection along a Line
function R = reflection(I,L)
%Author : Jeny Rajan
%I - Image to be Reflected
%L - Line position
%R - Reflected image
%eg : R = reflection(I,128);
% L should be between 1 and number of columns
[x y z]=size(I);
R=zeros(x,y,z);
R(:,L+1:y,:)=I(:,1:y-L,:);
R(:,1:L,:)=I(:,L:-1:1,:);
R=uint8(R);
end


text:
% Authors: A. Iscen, G. Tolias, Y. Avrithis, T. Furon, O. Chum. 2017.
% create image and feature ids for features of an image collection
% [image_ids, feature_ids] = imgfeatids (nof)
% nof: Nx1 vector with number of features per image
% image_ids: Mx1 image id per feature. M = sum(nof)
% feautre_ids: Mx1 unique feature ids
function [image_ids, feature_ids] = imgfeatids (nof)

% image ids for each database descriptor
cs = cumsum(double (nof));
[~, image_ids] = histc (1: cs (end), [1 cs+1]); %image ids here

% to create feature ids
feature_ids = 1:sum (nof);
rng = cs (1) + 1: length (feature_ids); %range such that values of first image are left unchanged
feature_ids (rng) = feature_ids (rng) - cs ( image_ids (rng) - 1); % subtract number of features of previous images such that fids becomes number of feature id per image

author: ahmetius

repo: diffusion-retrieval


text:
function p=getvarfrom(ws,name)
%
% p=getvarfrom(ws,name)
%
% get variable value by name from specified work-space
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%    ws: name of the work-space, for example, 'base'
%    name: name string of the variable
%
% output:
%    p: the value of the specified variable, if the variable does not
%       exist, return empty array
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

wsname=ws;
if(~iscell(ws))
   wsname=cell(1);
   wsname{1}=ws;
end

p=[];
for i=1:length(wsname)
    isdefined=evalin(wsname{i},['exist(''' name ''')']);
    if(isdefined==1)
        p=evalin(wsname{i},name);
        break;
    end
end


text:
function sp = spdim(n, d, options)
% SPDIM   Computes the number of sparse grid points
%    SP = SPDIM(N,D) Computes the number of points of the sparse
%    grid of dimension D and level N.
%
%    SP = SPDIM(N,D,OPTIONS) Computes the number of points as
%    above, but with default grid type replaced by the grid type
%    specified in OPTIONS, an argument created with the SPSET
%    function. See SPSET for details.
%
%    See also SPINTERP, SPGRID, SPVALS.

% Author : Andreas Klimke
% Version: 1.2
% Date   : November 18, 2007

% Change log:
% V1.0   : September 24, 2003
%          Initial version
% V1.1   : June 15, 2004
%          Added new grid type : Chebyshev distributed nodes
%          (at the extrema of the Chebyshev polynomials)
% V1.2   : November 18, 2007
%          Added new grid type : Gauss-Patterson

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

if nargin < 3, options = []; end

if d == 0
	% By definition; useful when computing recurrence formulae such
  % as Bungartz, "Finite Elements of Higher Order on Sparse
  % Grids", p.35, 1998.
	sp = 1;
	return;
end

gridtype = spget(options, 'GridType', 'Clenshaw-Curtis');

if strcmpi(gridtype, 'clenshaw-curtis')
	sp = spdimcc(n,d);
elseif strcmpi(gridtype, 'maximum')
	sp = spdimm(n,d);
elseif strcmpi(gridtype, 'noboundary')
	sp = spdimm(n,d,0);
elseif strcmpi(gridtype, 'chebyshev')
	% number of nodes same as CC-grid
	sp = spdimcc(n,d);
elseif strcmpi(gridtype, 'gauss-patterson')
	% number of nodes same as NB-grid
	sp = spdimm(n,d,0);
else
	error('MATLAB:spinterp:badopt',['Unknown grid type ''' gridtype '''.']);
end


text:
function [ ] = display_graph(x_category, y_category, algorithm_list, w_list, info_list, scale, line, width, xlim_range_in, ylim_range_in)
% SHow graphs of optimizations
%
% Inputs:
%       x_category          "numofgrad" or "iter" or "epoch" or "grad_calc_count"
%       y_category          "cost" or "optimality_gap" or "gnorm" or "subgnorm"
%       algorithms_list     algorithms to be evaluated
%       w_list              solution produced by each algorithm
%       info_list           statistics produced by each algorithm
%
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created by H.Kasai on Oct. 23, 2016
% Modified by H.Kasai on Apr. 24, 2017

    if nargin < 6
        scale_type = 'semilogy';
    else
        scale_type = scale;
    end

    if nargin < 7
        line_type = 'line';
    else
        line_type = line;
    end

    if nargin < 8
        linewidth = 2;
    else
        linewidth = width;
    end

    if nargin < 9
        xlim_range = [];
    else
        xlim_range = xlim_range_in;
    end

    if nargin < 10
        ylim_range = [];
    else
        ylim_range = ylim_range_in;
    end

    % for plotting
    if strcmp(line_type, 'line')
        linetype = {'r','b','m','g','c','y','r--','b--','c--','g--','m--','y--','r:','b:','c:','g:','m:','y:','r.','b.','c.','g.','m.','y.'};
    elseif strcmp(line_type, 'line-with-mark')
        linetype = {'ro-','bo-','mo-','go-','co-','yo-','r*-','b*-','m*-','g*-','c*-','y*-','r+--','b+--','m+--','g+--','c+--','y+--','rs:','bs:','ms:','gs:','cs:','ys:','r.','b.','c.','g.','m.','y.'};
    else
    end
    fontsize = 16;
    markersize = 5;

    % initialize
    legend_str = cell(1);
    alg_num = 0;

    % plot
    figure;
    for alg_idx=1:length(algorithm_list)
        if ~isempty(info_list{alg_idx})
            alg_num = alg_num + 1;

            if strcmp(x_category, 'numofgrad')
                x_plot_data = info_list{alg_idx}.grad_calc_count;
            elseif strcmp(x_category, 'iter')
                if isfield(info_list{alg_idx}, 'iter')
                    x_plot_data = info_list{alg_idx}.iter;
                elseif isfield(info_list{alg_idx}, 'total_iter')
                    x_plot_data = info_list{alg_idx}.total_iter;
                end
            elseif strcmp(x_category, 'epoch')
                x_plot_data = info_list{alg_idx}.epoch;
            elseif strcmp(x_category, 'inner_iter')
                x_plot_data = info_list{alg_idx}.subinfos.inner_iter;
            elseif strcmp(x_category, 'grad_calc_count')
                x_plot_data = info_list{alg_idx}.grad_calc_count;
            elseif strcmp(x_category, 'time')
                x_plot_data = info_list{alg_idx}.time;
            elseif strcmp(x_category, 'lambda') || strcmp(x_category, 'l1-norm')
                x_plot_data = w_list;
            elseif strcmp(x_category, 'coeff_pos')
                x_plot_data = [1:length(w_list{alg_idx})];
            else
            end

            if strcmp(y_category, 'cost')
                y_plot_data = info_list{alg_idx}.cost;
            elseif strcmp(y_category, 'best_cost')
                y_plot_data = info_list{alg_idx}.best_cost;
            elseif strcmp(y_category, 'cost_lag')
                y_plot_data = info_list{alg_idx}.cost_lag;
            elseif strcmp(y_category, 'optimality_gap')
                y_plot_data = info_list{alg_idx}.optgap;
            elseif strcmp(y_category, 'abs_optimality_gap')
                y_plot_data = info_list{alg_idx}.absoptgap;
            elseif strcmp(y_category, 'best_optimality_gap')
                y_plot_data = info_list{alg_idx}.best_optgap;
            elseif strcmp(y_category, 'sol_optimality_gap')
                y_plot_data = info_list{alg_idx}.sol_optgap;
            elseif strcmp(y_category, 'dual_gap')
                y_plot_data = info_list{alg_idx}.dual_gap;
            elseif strcmp(y_category, 'const_norm')
                y_plot_data = info_list{alg_idx}.const_norm;
            elseif strcmp(y_category, 'econst_norm')
                y_plot_data = info_list{alg_idx}.econst_norm;
            elseif strcmp(y_category, 'ineconst_norm')
                y_plot_data = info_list{alg_idx}.ineconst_norm;
            elseif strcmp(y_category, 'inv_rho')
                y_plot_data = info_list{alg_idx}.inv_rho;
            elseif strcmp(y_category, 'eta')
                y_plot_data = info_list{alg_idx}.eta;
            elseif strcmp(y_category, 'gradL_norm')
                y_plot_data = info_list{alg_idx}.gradL_norm;
            elseif strcmp(y_category, 'gnorm')
                y_plot_data = info_list{alg_idx}.gnorm;
            elseif strcmp(y_category, 'subgnorm')
                y_plot_data = info_list{alg_idx}.subgnorm;
            elseif strcmp(y_category, 'K')
                y_plot_data = info_list{alg_idx}.K;
            elseif strcmp(y_category, 'reg') || strcmp(y_category, 'l1-norm') || strcmp(y_category, 'trace_norm')
                y_plot_data = info_list{alg_idx}.reg;
            elseif strcmp(y_category, 'coeffs') || strcmp(y_category, 'aprox_err')
                y_plot_data = info_list{1};
            elseif strcmp(y_category, 'coeff_amp') || strcmp(y_category, 'aprox_err')
                y_plot_data = info_list{alg_idx};
            elseif strcmp(y_category, 'dnorm')
                y_plot_data = info_list{alg_idx}.subinfos.dnorm;
            end

            if strcmp(scale_type, 'semilogy')
                semilogy(x_plot_data, y_plot_data, linetype{alg_num}, 'MarkerSize', markersize, 'Linewidth', linewidth); hold on;
            elseif strcmp(scale_type, 'loglog')
                loglog(x_plot_data, y_plot_data, linetype{alg_num}, 'MarkerSize', markersize, 'Linewidth', linewidth); hold on;
            elseif strcmp(scale_type, 'linear')
                if strcmp(y_category, 'coeffs')
                    plot(x_plot_data, y_plot_data, 'Linewidth', linewidth); hold on;
                else
                    plot(x_plot_data, y_plot_data, linetype{alg_num}, 'MarkerSize', markersize, 'Linewidth', linewidth); hold on;
                end
            else
                error('Invalid scale type');
            end

            if ~strcmp(y_category, 'coeff')
                legend_str{alg_num} = algorithm_list{alg_idx};
            end
        else
            %
        end
    end
    hold off;

    % X label
    if strcmp(x_category, 'numofgrad')
        xlabel('Number of gradient evaluations', 'FontSize', fontsize);
    elseif strcmp(x_category, 'iter')
        xlabel('Iteration', 'FontSize', fontsize);
    elseif strcmp(x_category, 'epoch')
        xlabel('Epoch', 'FontSize', fontsize);
    elseif strcmp(x_category, 'grad_calc_count')
        xlabel('Number of gradient calculations', 'FontSize', fontsize);
    elseif strcmp(x_category, 'time')
        xlabel('Time', 'FontSize', fontsize);
    elseif strcmp(x_category, 'lambda')
        xlabel('$$\lambda$$', 'FontSize', fontsize,'Interpreter', 'Latex');
    elseif strcmp(x_category, 'l1-norm')
        xlabel('$$\ell$$-1 norm', 'FontSize', fontsize,'Interpreter', 'Latex');
    elseif strcmp(x_category, 'coeff_pos')
        xlabel('Coefficient position', 'FontSize', fontsize);
    elseif strcmp(x_category, 'inner_iter')
        xlabel('Inner Iteration', 'FontSize', fontsize);
    else
    end

    % Y label
    if strcmp(y_category, 'cost')
        ylabel('Cost', 'FontSize', fontsize);
    elseif strcmp(y_category, 'best_cost')
        ylabel('Best cost', 'FontSize', fontsize);
    elseif strcmp(y_category, 'cost_lag')
        ylabel('Lagrangian cost', 'FontSize', fontsize);
    elseif strcmp(y_category, 'optimality_gap')
        ylabel('Optimality gap', 'FontSize', fontsize);
    elseif strcmp(y_category, 'abs_optimality_gap')
        ylabel('Absolute optimality gap', 'FontSize', fontsize);
    elseif strcmp(y_category, 'best_optimality_gap')
        ylabel('Best optimality gap', 'FontSize', fontsize);
    elseif strcmp(y_category, 'sol_optimality_gap')
        ylabel('Solution optimality gap', 'FontSize', fontsize);
    elseif strcmp(y_category, 'dual_gap')
        ylabel('Dual gap', 'FontSize', fontsize);
    elseif strcmp(y_category, 'const_norm')
        ylabel('Norm of constraints', 'FontSize', fontsize);
    elseif strcmp(y_category, 'econst_norm')
        ylabel('Norm of equality constraints', 'FontSize', fontsize);
    elseif strcmp(y_category, 'ineconst_norm')
        ylabel('Norm of inequality constraints', 'FontSize', fontsize);
    elseif strcmp(y_category, 'gradL_norm')
        ylabel('Norm of gradient of Lagrangian', 'FontSize', fontsize);
    elseif strcmp(y_category, 'inv_rho')
        ylabel('Inverse of \rho', 'FontSize', fontsize);
    elseif strcmp(y_category, 'eta')
        ylabel('\eta', 'FontSize', fontsize);
    elseif strcmp(y_category, 'gnorm')
        ylabel('Norm of gradient', 'FontSize', fontsize);
    elseif strcmp(y_category, 'subgnorm')
        ylabel('Norm of subgradient', 'FontSize', fontsize);
    elseif strcmp(y_category, 'K')
        ylabel('Batch size', 'FontSize', fontsize);
    elseif strcmp(y_category, 'reg')
        ylabel('Regularizer', 'FontSize', fontsize);
    elseif strcmp(y_category, 'trace_norm')
        ylabel('Trace (nuclear) norm', 'FontSize', fontsize);
    elseif strcmp(y_category, 'l1-norm')
        ylabel('$$\ell$$-1 norm', 'FontSize', fontsize, 'Interpreter', 'Latex', 'FontName','Arial');
    elseif strcmp(y_category, 'coeffs')
        ylabel('Coefficient', 'FontSize', fontsize);
    elseif strcmp(y_category, 'aprox_err')
        ylabel('Approximation error', 'FontSize', fontsize);
    elseif strcmp(y_category, 'coeff_amp')
        ylabel('Coefficient amplitude', 'FontSize', fontsize);
    elseif strcmp(y_category, 'dnorm')
        ylabel('Norm of direction', 'FontSize', fontsize);
    end

    % range
    if ~isempty(xlim_range)
        xlim([xlim_range])
    end

    if ~isempty(ylim_range)
        ylim([ylim_range])
    end

    % legend
    if ~strcmp(y_category, 'coeffs')
        legend(legend_str);
    end

    set(gca, 'FontSize', fontsize);
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function [y] = spm_gx_hdm_sck(x,u,P,M)
% Simulated BOLD response to input.
% FORMAT [y] = spm_gx_hdm_sck(x,u,P,M)
% y    - BOLD response (%)
% x    - state vector     (see spm_fx_fmri)
% P    - Parameter vector (see spm_fx_fmri)
%__________________________________________________________________________
%
% This function implements the BOLD signal model described in:
%
% Stephan KE, Weiskopf N, Drysdale PM, Robinson PA, Friston KJ (2007)
% Comparing hemodynamic models with DCM. NeuroImage 38: 387-401.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston & Klaas Enno Stephan
% $Id: spm_gx_hdm_sck.m 6263 2014-11-17 13:48:36Z karl $

% biophysical constants for 1.5 T:
%==========================================================================

% echo time (seconds)
%--------------------------------------------------------------------------
try
    TE = M(1).TE;
catch
    TE = 0.04;
end

% resting venous volume
%--------------------------------------------------------------------------
V0    = 100*0.08;

% slope r0 of intravascular relaxation rate R_iv as a function of oxygen
% saturation Y:  R_iv = r0*[(1-Y)-(1-Y0)]
%--------------------------------------------------------------------------
r0    = 25;

% frequency offset at the outer surface of magnetized vessels
%--------------------------------------------------------------------------
nu0   = 40.3;

% region-specific resting oxygen extraction fractions
%--------------------------------------------------------------------------
E0    = P(5,:);

% region-specific ratios of intra- to extravascular components of
% the gradient echo signal (prior mean = 1, log-normally distributed
% scaling factor)
%--------------------------------------------------------------------------
epsi  = exp(P(6,:));

% coefficients in BOLD signal model
%--------------------------------------------------------------------------
k1    = 4.3.*nu0.*E0.*TE;
k2    = epsi.*r0.*E0.*TE;
k3    = 1 - epsi;

% exponentiation of hemodynamic state variables
%--------------------------------------------------------------------------
x     = exp(x);

% BOLD signal
%--------------------------------------------------------------------------
v     = x(3,:);
q     = x(4,:);
y     = V0.*(k1.*(1 - q) + k2.*(1 - (q./v)) + k3.*(1 - v));

author: spm


text:
function rgb = convertYuvToRgb(yuv)
% convert row vector YUV [0, 255] in row vector RGB [0, 255]

load conversion.mat; % load conversion matrices

yuv = double(yuv);

yuv(:, 2 : 3) = yuv(:, 2 : 3) - 127;
rgb = (yuvToRgb *yuv.').';

rgb = uint8(clipValue(rgb, 0, 255));

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL