text:
function varargout = sqrt(varargin)
%SQRT (overloaded)
%
% t = sqrt(x)
%
% The variable t can only be used in concavity preserving
% operations such as t>=1, max t etc.
%
% When SQRT is used in a problem, the domain constraint
% (x>=0) is automatically added to the problem.
%
% In nonconvex cases, use sqrtm instead.
%
% See also CPOWER

switch class(varargin{1})

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        X = varargin{1};
        [n,m] = size(X);
        if is(varargin{1},'real') %& (n*m==1)
            varargout{1} = InstantiateElementWise(mfilename,varargin{:});
        else
            error('SQRT can only be applied to real scalars');
        end

    case 'char' % YALMIP send 'model' when it wants the epigraph or hypograph
        switch varargin{1}
            case 'graph'

                t = varargin{2}; % Second arg is the extended operator variable
                X = varargin{3}; % Third arg and above are the args user used when defining t.
                if is(X,'linear')
                    varargout{1} = (cone([(X-1)/2;t],(X+1)/2));
                    varargout{2} = struct('convexity','concave','monotonicity','increasing','definiteness','positive');
                    varargout{3} = X;
                elseif is(X,'quadratic')
                    [F,x] = check_for_special_cases(X,t);
                    if isempty(F)
                        varargout{1} = [];
                        varargout{2} = [];
                        varargout{3} = [];
                    else
                        varargout{1} = F;
                        varargout{2} = struct('convexity','convex','monotonicity','none','definiteness','positive');
                        varargout{3} = x;
                    end

                else
                    varargout{1} = [];
                    varargout{2} = [];
                    varargout{3} = [];
                end

            otherwise
                varargout{1} = [];
                varargout{2} = [];
                varargout{3} = [];
        end
    otherwise
end

function [F,x] = check_for_special_cases(q,t)
% Check if user is constructing sqrt(quadratic). If that is the case,
% return norm(linear)
F = [];
x = [];
if length(q)>1
    return
end
[Q,c,f,x,info] = quaddecomp(q);
if info==0 & nnz(Q)>0
    index = find(any(Q,2));
    if length(index) < length(Q)
        Qsub = Q(index,index);
        [Rsub,p]=chol(Qsub);
        if p==0
            [i,j,k] = find(Rsub);
            R = sparse((i),index(j),k,length(Qsub),length(Q));
        else
            R = [];
        end
    else
        [R,p]=chol(Q);
        if p & min(eig(full(Q)))>=-1e-12
            [U,S,V] = svd(full(Q));
            r = max(find(diag(S)));
            R = sqrtm(S(1:r,1:r))*V(:,1:r)';
            p = 0;
        end
    end
    d = 0.5*(R'\c);
    if p==0 &  f-d'*d>-1e-12
        F = cone([R*x+d;sqrt(f-d'*d)],t);
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function events = in_events_kdf(sFile, EventFile)
% IN_EVENTS_KDF: Read the events descriptions from a KRISS .trg file.
%
% USAGE:  events = in_events_kdf(sFile, EventFile)

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

% Initialize returned variable
events = repmat(db_template('event'), 0);
% Load trigger file
EventDat = load(EventFile, '-ASCII');
if isempty(EventDat)
    return;
end

% Consider both the stim and responses as events of the same level
AllEvt = [EventDat(:,1:2); EventDat(EventDat(:,4)~=0, 3:4)];
% Get all the unique events
uniqueEvt = unique(AllEvt(:,2));
% Create events structures: one per category of event
for iEvt = 1:length(uniqueEvt)
    % Find all the occurrences of event #iEvt
    iMrk = find(AllEvt(:,2) == uniqueEvt(iEvt));
    % Add event structure
    events(iEvt).label      = num2str(uniqueEvt(iEvt));
    events(iEvt).epochs     = ones(1, length(iMrk));
    events(iEvt).times      = AllEvt(iMrk,1)' ./ sFile.prop.sfreq;
    events(iEvt).reactTimes = [];
    events(iEvt).select     = 1;
    events(iEvt).channels   = [];
    events(iEvt).notes      = [];
end

author: brainstorm-tools


text:
function fiff_write_float(fid,kind,data)
%
% fiff_write_float(fid,kind,data)
%
% Writes a single-precision floating point tag to a fif file
%
%     fid           An open fif file descriptor
%     kind          Tag kind
%     data          The data
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.3  2006/04/27 22:38:37  msh
%   Splitting an empty list now results in an empty output.
%   Added fiff_write_double and fiff_write_short
%   Write an array of floats, ints, and doubles instead of just one value.
%   Fixed help text of fiff_pick_channels.
%
%   Revision 1.2  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.1  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%
%

me='MNE:fiff_write_float';

if nargin ~= 3
        error(me,'Incorrect number of arguments');
end

FIFFT_FLOAT=4;
FIFFV_NEXT_SEQ=0;
nel=numel(data);
datasize=nel*4;
count = fwrite(fid,int32(kind),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFT_FLOAT),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(datasize),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFV_NEXT_SEQ),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,single(data),'single');
if count ~= nel
    error(me,'write failed');
end

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function x = relaxdouble(x)
%RELAXDOUBLE (overloaded sdpvar/relaxdouble on double)

% Do nothing, just return input.


text:
load('pend_data', 'S', 'n');

lmi = lmistruct(S, n);
lmi = lmi_asym_decay(lmi, 0.5); % |x(t)| < C * exp(-0.5 t)
lmi = lmi_input(lmi, 30, 0.3);  % max 30N force when |x| < 0.3
K = lmi_solve(lmi);

save('pend_data', '-append', 'K');

author: Sable

repo: mcbench-benchmarks


text:
function SaveRippleEvents(filename,ripples,channelID)

%SaveRippleEvents - Save hippocampal ripple (~200Hz oscillations) events.
%
%  USAGE
%
%    SaveRippleEvents(filename,ripples,channelID)
%
%    filename       file to save to
%    ripples        ripple info as provided by FindRipples
%    channelID      channel ID (appended to the event description)
%
%  SEE
%
%    See also FindRipples, RippleStats, PlotRippleStats, SaveEvents.

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

if nargin < 3,
  error('Incorrect number of parameters (type ''help SaveRippleEvents'' for details).');
end

n = size(ripples,1);
r = ripples(:,1:3)';
events.time = r(:);
for i = 1:3:3*n,
	events.description{i,1} = ['Ripple start ' int2str(channelID)];
	events.description{i+1,1} = ['Ripple peak ' int2str(channelID)];
	events.description{i+2,1} = ['Ripple stop ' int2str(channelID)];
end

SaveEvents(filename,events);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% PLOTTOPO - plot concatenated multichannel data epochs in a topographic
% or
%              rectangular array. Uses a channel location file with the same
%              format as TOPOPLOT, or else plots data on a rectangular grid.
%              If data are all positive, they are assumed to be spectra.
% Usage:
%    >> plottopo(data, 'key1', 'val1', 'key2', 'val2')
% Or
%    >> plottopo(data,'chan_locs',frames,limits,title,channels,...
%                      axsize,colors,ydir,vert) % old function call
% Inputs:
%   data       = data consisting of consecutive epochs of (chans,frames)
%                or (chans,frames,n)
%
% Optional inputs:
%  'chanlocs'  = [struct] channel structure or file plot ERPs at channel
%                locations. See help READLOCS for data channel format.
%  'geom'      = [rows cols] plot ERP in grid (overwrite previous option).
%                Grid size for rectangular matrix. Example: [6 4].
%  'frames'    = time frames (points) per epoch {def|0 -> data length}
%  'limits'    = [xmin xmax ymin ymax]  (x's in ms or Hz) {def|0
%                 (or both y's 0) -> use data limits)
%  'ylim'      = [ymin ymax] y axis limits. Overwrite option above.
%  'title'     = [string] plot title {def|'' -> none}
%  'chans'     = vector of channel numbers to plot {def|0 -> all}
%  'axsize'    = [x y] axis size {default [.05 .08]}
%  'legend'    = [cell array] cell array of string for the legend. Note
%                the last element can be an integer to set legend
%                position.
%  'showleg'   = ['on'|'off'] show or hide legend.
%  'colors'    = [cell array] cell array of plot aspect. E.g. { 'k' 'k--' }
%                for plotting the first curve in black and the second one
%                in black dashed. Can also contain additional formatting.
%                { { 'k' 'linewidth' 2 } 'k--' } same as above but
%                the first line is bolded.
%  'ydir'      = [1|-1] y-axis polarity (pos-up = 1; neg-up = -1) {def -> -1}
%  'vert'      = [vector] of times (in ms or Hz) to plot vertical lines
%                {def none}
%  'hori'      = [vector] plot horizontal line at given ordinate values.
%  'regions'   = [cell array] cell array of size nchan. Each cell contains a
%                float array of size (2,n) each column defining a time region
%                [low high] to be highlighted.
%  'plotfunc'  = [cell] use different function for plotting data. The format
%                is { funcname arg2 arg3 arg2 ... }. arg1 is taken from the
%                data.
%
% Author: Scott Makeig and Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 3-2-98
%
% See also: PLOTDATA, TOPOPLOT

% Copyright (C) 3-2-98 from PLOTDATA Scott Makeig, SCCN/INC/UCSD,
% scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 5-11-98 added channels arg -sm
% 7-15-98 added ydir arg, made pos-up the default -sm
% 7-23-98 debugged ydir arg and pos-up default -sm
% 12-22-99 added grid size option, changed to SBPLOT order -sm
% 03-16-00 added AXCOPY feature -sm & tpj
% 08-21-00 debugged axheight/axwidth setting -sm
% 01-25-02 reformated help & license, added links -ad
% 03-11-02 change the channel names plotting position and cutomize pop-up -ad
% 03-15-02 added readlocs and the use of eloc input structure -ad
% 03-15-02 debugging chanlocs structure -ad & sm

%  'chan_locs' = file of channel locations as in >> topoplot example   {grid}
%                ELSE: [rows cols] grid size for rectangular matrix. Example: [6 4]
%   frames     = time frames (points) per epoch {def|0 -> data length}
%  [limits]    = [xmin xmax ymin ymax]  (x's in ms or Hz) {def|0
%                 (or both y's 0) -> use data limits)
%  'title'     = plot title {def|0 -> none}
%   channels   = vector of channel numbers to plot & label {def|0 -> all}
%                   else, filename of ascii channel-name file
%   axsize     = [x y] axis size {default [.07 .07]}
%  'colors'    = file of color codes, 3 chars per line
%                ( '.' = space) {0 -> default color order}
%   ydir       = y-axis polarity (pos-up = 1; neg-up = -1) {def -> pos-up}
%   vert       = [vector] of times (in ms or Hz) to plot vertical lines {def none}
%   hori        = [vector] of amplitudes (in uV or dB) to plot horizontal lines {def none}
%

function plottopo(data, varargin);

%
%%%%%%%%%%%%%%%%%%%%% Graphics Settings - can be customized %%%%%%%%%%%%%%%%%%
%
LINEWIDTH     = 0.7;     % data line widths (can be non-integer)
FONTSIZE      = 10;      % font size to use for labels
CHANFONTSIZE  = 7;       % font size to use for channel names
TICKFONTSIZE  = 8;       % font size to use for axis labels
TITLEFONTSIZE = 12;      % font size to use for the plot title
PLOT_WIDTH    = 0.95;     % 0.75, width and height of plot array on figure
PLOT_HEIGHT   = 0.88;    % 0.88
gcapos = get(gca,'Position'); axis off;
PLOT_WIDTH    = gcapos(3)*PLOT_WIDTH; % width and height of gca plot array on gca
PLOT_HEIGHT   = gcapos(4)*PLOT_HEIGHT;
curfig = gcf;            % learn the current graphic figure number
%
%%%%%%%%%%%%%%%%%%%% Default settings - use commandline to override %%%%%%%%%%%
%
DEFAULT_AXWIDTH  = 0.05; %
DEFAULT_AXHEIGHT = 0.08; %
DEFAULT_SIGN = -1;                        % Default - plot positive-up
ISRECT = 0;                               % default

if nargin < 1
    help plottopo
    return
end

if length(varargin) > 0
    if length(varargin) == 1 || ~ischar(varargin{1}) || isempty(varargin{1}) || ...
        (length(varargin)>2 &  ~ischar(varargin{3}))
        options = { 'chanlocs' varargin{1} };
        if nargin > 2, options = { options{:} 'frames' varargin{2} }; end
        if nargin > 3, options = { options{:} 'limits' varargin{3} }; end
        if nargin > 5, options = { options{:} 'chans'  varargin{5} }; end
        if nargin > 6, options = { options{:} 'axsize' varargin{6} }; end
        if nargin > 7, options = { options{:} 'colors' varargin{7} }; end
        if nargin > 8, options = { options{:} 'ydir'   varargin{8} }; end
        if nargin > 9, options = { options{:} 'vert'   varargin{9} }; end
        if nargin > 10,options = { options{:} 'hori'  varargin{10} }; end
        if nargin > 4 && ~isequal(varargin{4}, 0), options = {options{:} 'title'  varargin{4} }; end
        %    , chan_locs,frames,limits,plottitle,channels,axsize,colors,ydr,vert)
    else
        options = varargin;
    end
else
    options = varargin;
end
g = finputcheck(options, { 'chanlocs'  ''    []          '';
                    'frames'    'integer'               [1 Inf]     size(data,2);
                    'chans'     { 'integer','string' }  { [1 Inf] [] }    0;
                    'geom'      'integer'               [1 Inf]     [];
                    'channames' 'string'                []          '';
                    'limits'    'float'                 []          0;
                    'ylim'      'float'                 []          [];
                    'title'     'string'                []          '';
                    'plotfunc'  'cell'                  []          {};
                    'axsize'    'float'                 [0 1]       [nan nan];
                    'regions'   'cell'                  []          {};
                    'colors'    { 'cell','string' }     []          {};
                    'legend'    'cell'                  []          {};
                    'showleg'   'string'                {'on','off'} 'on';
                    'ydir'      'integer'               [-1 1]      DEFAULT_SIGN;
                    'vert'      'float'                 []          [];
                    'hori'      'float'                 []          []});
if ischar(g), error(g); end
data = reshape(data, size(data,1), size(data,2), size(data,3));
%if length(g.chans) == 1 & g.chans(1) ~= 0, error('can not plot a single ERP'); end

[chans,framestotal]=size(data);           % data size

%
%%%%%%%%%%%%%%% Substitute defaults for missing parameters %%%%%%%%%%%%%%%%
%

axwidth  = g.axsize(1);
if length(g.axsize) < 2
    axheight = NaN;
else
    axheight = g.axsize(2);
end
if isempty(g.chans) || g.chans(1) == 0
   g.chans = 1:size(data,1);
elseif ~ischar(g.chans)
   g.chans = g.chans;
end

nolegend = 0;
if isempty(g.legend), nolegend = 1; end

if ~isempty(g.ylim)
    g.limits(3:4) = g.ylim;
end
plotgrid = 0;
if isempty(g.chanlocs) % plot in a rectangular grid
    plotgrid = 1;
elseif ~isfield(g.chanlocs, 'theta')
    plotgrid = 1;
end
if length(g.chans) < 4 && ~plotgrid
    disp('Not enough channels, does not use channel coordinate to plot axis');
    plotgrid = 1;
end
if plotgrid && isempty(g.geom)
  n = ceil(sqrt(length(g.chans)));
  g.geom = [n ceil(length(g.chans)/n)];
end
if ~isempty(g.geom)
    plotgrid = 1;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Test parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  icadefs; % read BACKCOLOR, MAXPLOTDATACHANS constant from icadefs.m
  if g.frames <=0,
    g.frames = framestotal;    % default
    datasets=1;
  elseif g.frames==1,
    fprintf('plottopo: cannot plot less than 2 frames per trace.\n');
    return
    datasets=1;
  else
    datasets = fix(framestotal/g.frames);        % number of traces to overplot
  end

  if max(g.chans) > chans
    fprintf('plottopo(): max channel index > %d channels in data.\n',...
                       chans);
    return
  end
  if min(g.chans) < 1
    fprintf('plottopo(): min channel index (%g) < 1.\n',...
                       min(g.chans));
    return
  end
  if length(g.chans)>MAXPLOTDATACHANS,
    fprintf('plottopo(): not set up to plot more than %d traces.\n',...
                       MAXPLOTDATACHANS);
    return
  end

  if datasets>MAXPLOTDATAEPOCHS
      fprintf('plottopo: not set up to plot more than %d epochs.\n',...
                       MAXPLOTDATAEPOCHS);
    return
  end
  if datasets<1
      fprintf('plottopo: cannot plot less than 1 epoch!\n');
      return
  end

  if ~isempty(g.geom)
      if isnan(axheight) % if not specified
          axheight = gcapos(4)/(g.geom(1)+1);
          axwidth  = gcapos(3)/(g.geom(2)+1);
      end
      % if chan_locs(2) > 5
      %     axwidth = 0.66/(chan_locs(2)+1);
      % end
  elseif isnan(axwidth)
      axheight = DEFAULT_AXHEIGHT;
      axwidth =  DEFAULT_AXWIDTH;
  end
    fprintf('Plotting data using axis size [%g,%g]\n',axwidth,axheight);

    %
    %%%%%%%%%%%%% Extend the size of the plotting area in the window %%%%%%%%%%%%
    %
    curfig = gcf;
    h=figure(curfig);
    set(h,'PaperUnits','normalized'); % use percentages to avoid US/A4 difference
    set(h,'PaperPosition',[0.0235308 0.0272775 0.894169 0.909249]); % equivalent
    orient portrait
    axis('normal');

    set(gca,'Color',BACKCOLOR);               % set the background color

    axcolor= get(0,'DefaultAxesXcolor'); % find what the default x-axis color is
    vertcolor = 'k';
    horicolor = vertcolor;

%     %
%     %%%%%%%%%%%%%%%%%%%% Read the channel names %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %
%     if isempty(g.channames)
%         if ~ischar(g.chans)
%             % g.channames = zeros(MAXPLOTDATACHANS,4);
%             % for c=1:length(g.chans),
%             %     g.channames(c,:)= sprintf('%4d',g.chans(c));
%             % end
%             if length(g.chans) > 1 | g.chans(1) ~= 0
%                 g.channames = num2str(g.chans(:));                   %%CJH
%             end
%         else % ischar(g.chans)
%             chid = fopen(g.chans,'r');
%             if chid <3,
%                 fprintf('PLOTTOPO: cannot open file %s.\n',g.chans);
%                 return
%             else
%                 fprintf('PLOTTOPO: opened file %s.\n',g.chans);
%             end
%
%             %%%%%%%
%             % fid=fopen('fgetl.m');
%             % while 1
%             %   line = fgetl(fid);
%             %   if ~ischar(line), break, end
%             %     disp(line)
%             %   end
%             % end
%             % fclose(fid);
%             %%%%%%%
%
%             g.channames = fscanf(chid,'%s',[4 MAXPLOTDATACHANS]);
%             g.channames = g.channames';
%             [r c] = size(g.channames);
%             for i=1:r
%                 for j=1:c
%                     if g.channames(i,j)=='.',
%                         g.channames(i,j)=' ';
%                     end
%                 end
%             end
%         end; % setting g.channames
%     end
%
    %
    %%%%%%%%%%%%%%%%%%%%%%%%% Plot and label specified channels %%%%%%%%%%%%%%%%%%
    %
    data = data(g.chans,:);
    chans = length(g.chans);

    %
    %%%%%%%%%%%%%%%%%%%%%%%%% Read the color names %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if ischar(g.colors) % filename for backward compatibility but not documented
        cid = fopen(g.colors,'r');
        % fprintf('cid = %d\n',cid);
        if cid <3,
            fprintf('plottopo: cannot open file %s.\n',g.colors);
            return
        end
        g.colors = fscanf(cid,'%s',[3 MAXPLOTDATAEPOCHS]);
        g.colors = g.colors';
        [r c] = size(g.colors);
        for i=1:r
            for j=1:c
                if g.colors(i,j)=='.',
                    g.colors(i,j)=' ';
                end
            end
        end
        g.colors = cellstr(g.colors);
        for c=1:length(g.colors)   % make white traces black unless axis color is white
            if g.colors{c}(1)=='w' && axcolor~=[1 1 1]
                g.colors{c}(1)='k';
            end
        end
    else % use default color order (no yellow!)
        tmpcolors = { 'b' 'r' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' ...
                      'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm' 'r' 'b' 'g' 'c' 'm'};
        g.colors = {g.colors{:} tmpcolors{:}};  % make > 64 available
    end
    %
    %%%%%%%%%%%%%%%%%%%%%%% Read and adjust limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if g.limits==0,      % == 0 or [0 0 0 0]
        xmin=0;
        xmax=g.frames-1;
        % for abs max scaling:
        ymax=max(max(abs(data)));
        ymin=ymax*-1;
        % for data limits:
        %ymin=min(min(data));
        %ymax=max(max(data));
    else
        if length(g.limits)~=4,
            error('plottopo: limits should be 0 or an array [xmin xmax ymin ymax].\n');
        end
        xmin = g.limits(1);
        xmax = g.limits(2);
        ymin = g.limits(3);
        ymax = g.limits(4);
    end

    if xmax == 0 && xmin == 0,
        x = (0:1:g.frames-1);
        xmin = 0;
        xmax = g.frames-1;
    else
        dx = (xmax-xmin)/(g.frames-1);
        x=xmin*ones(1,g.frames)+dx*(0:g.frames-1); % compute x-values
    end
    if xmax<=xmin,
        fprintf('plottopo() - xmax must be > xmin.\n')
        return
    end

    if ymax == 0 && ymin == 0,
        % for abs max scaling:
        ymax=max(max(abs(data)));
        ymin=ymax*-1;
        % for data limits:
        %ymin=min(min(data));
        %ymax=max(max(data));
    end
    if ymax<=ymin,
        fprintf('plottopo() - ymax must be > ymin.\n')
        return
    end

    xlabel = 'Time (ms)';
    %
    %%%%%%%%%%%%%%%%%%%%%% Set up plotting environment %%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % h = gcf;
    % set(h,'YLim',[ymin ymax]);       % set default plotting parameters
    % set(h,'XLim',[xmin xmax]);
    % set(h,'FontSize',18);
    % set(h,'DefaultLineLineWidth',1); % for thinner postscript lines
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%% Print plot info %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % clf;   % clear the current figure

    % print plottitle over (left) subplot 1
    figure(curfig); h=gca;title(g.title,'FontSize',TITLEFONTSIZE); % title plot
    hold on
    msg = ['Plotting %d traces of %d frames with colors: '];

    for c=1:datasets
        cind = mod(c-1, length(g.colors))+1;
        if iscell(g.colors{cind})
            msg = [msg  '''' g.colors{cind}{1} ''' ' ];
        else
            msg = [msg  '''' g.colors{cind} ''' ' ];
        end
    end
    msg = [msg '\n'];    % print starting info on screen . . .
    fprintf('limits: [xmin,xmax,ymin,ymax] = [%4.1f %4.1f %4.2f %4.2f]\n',...
            xmin,xmax,ymin,ymax);
    fprintf(msg,datasets,g.frames);

    set(h,'FontSize',FONTSIZE);           % choose font size
    set(h,'YLim',[ymin ymax]);            % set default plotting parameters
    set(h,'XLim',[xmin xmax]);

    axis('off')
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%% Read chan_locs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if plotgrid
        ISRECT = 1;
        ht = g.geom(1);
        wd = g.geom(2);
        if chans > ht*wd
            fprintf('plottopo(): (%d) channels to be plotted > grid size [%d %d]\n',...
                    chans,ht,wd);
            return
        end
        xvals = 0; yvals = 0;
        if isempty(g.channames)
            if isfield(g.chanlocs,'labels') && ~iscellstr({g.chanlocs.labels})
                g.channames = strvcat(g.chanlocs.labels);
            else
                g.channames = repmat(' ',ht*wd,4);
                for i=1:ht*wd
                    channum = num2str(i);
                    g.channames(i,1:length(channum)) = channum;
                end
            end
        end

    else % read chan_locs file
         % read the channel location file
         % ------------------------------
        if isstruct(g.chanlocs)
            nonemptychans = cellfun('isempty', { g.chanlocs.theta });
            nonemptychans = find(~nonemptychans);
            [tmp g.channames Th Rd] = readlocs(g.chanlocs(nonemptychans));
            g.channames = strvcat({ g.chanlocs.labels });
        else
            [tmp g.channames Th Rd] = readlocs(g.chanlocs);
            g.channames = strvcat(g.channames);
            nonemptychans = [1:length(g.channames)];
        end
        Th = pi/180*Th;                 % convert degrees to radians
        Rd = Rd;

        if length(g.chans) > length(g.chanlocs),
            error('plottopo(): data channels must be <= ''chanlocs'' channels')
        end

        [yvalstmp,xvalstmp] = pol2cart(Th,Rd); % translate from polar to cart. coordinates
        xvals(nonemptychans) = xvalstmp;
        yvals(nonemptychans) = yvalstmp;

        % find position for other channels
        % --------------------------------
        totalchans = length(g.chanlocs);
        emptychans = setdiff_bc(1:totalchans, nonemptychans);
        totalchans = floor(sqrt(totalchans))+1;
        for index = 1:length(emptychans)
            xvals(emptychans(index)) = 0.7+0.2*floor((index-1)/totalchans);
            yvals(emptychans(index)) = -0.4+mod(index-1,totalchans)/totalchans;
        end
        g.channames = g.channames(g.chans,:);
        xvals     = xvals(g.chans);
        yvals     = yvals(g.chans);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % xvals = 0.5+PLOT_WIDTH*xvals;   % controls width of  plot array on page!
    % yvals = 0.5+PLOT_HEIGHT*yvals;  % controls height of plot array on page!
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if length(xvals) > 1
        if length(unique(xvals)) > 1
            xvals = (xvals-mean([max(xvals) min(xvals)]))/(max(xvals)-min(xvals)); % recenter
            xvals = gcapos(1)+gcapos(3)/2+PLOT_WIDTH*xvals;   % controls width of plot
                                                              % array on current axes
        end
    end
    yvals = gcapos(2)+gcapos(4)/2+PLOT_HEIGHT*yvals;  % controls height of plot
                                                      % array on current axes
                                                      %
                                                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot traces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                                      %

    xdiff=xmax-xmin;
    ydiff=ymax-ymin;

    Axes = [];
    for P=0:datasets-1, %  for each data epoch
        fprintf('trace %d: ',P+1);

        for c=1:chans, %%%%%%%% for each data channel %%%%%%%%%%%%%%%%%%%%%%%%%%

            if P>0 % subsequent pages (Axes specified)
                axes(Axes(c))
                hold on;                      % plot down left side of page first
                axis('off')
            else   % first page, specify axes
                if plotgrid
                    Axes = [ Axes sbplot(g.geom(1), g.geom(2), c)];
                else
                    xcenter = xvals(c);
                    ycenter = yvals(c);
                    Axes = [Axes axes('Units','Normal','Position', ...
                                      [xcenter-axwidth/2 ycenter-axheight/2 axwidth axheight])];
                end
                %axes(Axes(c))
                axis('off')

                hold on;                      % plot down left side of page first
                                              % set(h,'YLim',[ymin ymax]);    % set default plotting parameters
                                              % set(h,'XLim',[xmin xmax]);

                axislcolor = get(gca,'Xcolor');   %%CJH

                axis('off');

                % secondx = 200;                             % draw second vert axis
                % axis('off');plot([secondx secondx],[ymin ymax],'color',axislcolor);
                %
                %%%%%%%%%%%%%%%%%%%%%%% Print channel names %%%%%%%%%%%%%%%%%%%%%%%%%%
                %
                NAME_OFFSET = -.25;
                NAME_OFFSETY = .2;
                if ymin <= 0 && ymax >= 0,
                    yht = 0;
                else
                    yht = mean(data(c,1+P*g.frames:1+P*g.frames+g.frames-1));
                end
                if ~ISRECT    % print before traces
                    xt = double(xmin-NAME_OFFSET*xdiff);
                    yt = double(yht-NAME_OFFSETY*ydiff);
                    str = [g.channames(c,:)];
                    h=text(xt,yt,str);
                    set(h,'HorizontalAlignment','right');
                    %set(h,'FontSize',CHANFONTSIZE);           % choose font size
                else % ISRECT
                    xmn = xdiff/2+xmin;
                    h=text(double(xmn),double(ymax+0.05*ymax),[g.channames(c,:)]);
                    set(h,'HorizontalAlignment','right');
                    %set(h,'FontSize',CHANFONTSIZE);            % choose font size
                end % ISRECT

                %
                %%%%%%%%%%%%%%%%%%%%%%% Highlight regions %%%%%%%%%%%%%%%%%%%%%%%%%%
                %
                if ~isempty(g.regions)
                    for index=1:size(g.regions{c},2)
                        tmpreg = g.regions{c}(:,index);
                        if tmpreg(1) ~= tmpreg(2)
                            tmph = patch([tmpreg(1) tmpreg(2) tmpreg(2) tmpreg(1)], ...
                                         [-100 -100 100 100], [0.9 0.9 0.9]); hold on;
                            set(tmph, 'edgecolor', [0.9 0.9 0.9]); %,'facealpha',0.5,'edgealpha',0.5);
                        end
                    end
                end

            end; % P=0

            %
            %%%%%%%%%%%%%%%%%%%%%%% Plot data traces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            Pind = mod(P+1-1, length(g.colors))+1;
            if ~iscell( g.colors{Pind} ), tmpcolor = { g.colors{Pind} 'linewidth' LINEWIDTH };
            else                          tmpcolor = g.colors{Pind};
            end
            ymn = min([ymax ymin]);
            ymx = max([ymax ymin]);
            if isempty(g.plotfunc)
                if ischar(tmpcolor{1}) && length(tmpcolor) > 1
                    plot(x,data(c,1+P*g.frames:1+P*g.frames+g.frames-1), tmpcolor{1}, tmpcolor{2:end});
                else
                    plot(x,data(c,1+P*g.frames:1+P*g.frames+g.frames-1), 'color', tmpcolor{:});
                end;
                if g.ydir == -1
                    set(gca, 'ydir', 'reverse');
                end
                axis([xmin xmax ymn ymx]);          % set axis bounds
            elseif P == 1
                func = eval( [ '@' g.plotfunc{1} ] );
                feval(func, data(c,:), g.plotfunc{2:end});
            end

            if P == datasets-1 % last pass
                               %
                               %%%%%%%%%%%%%%%%%%%%%%% Plot lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                               %
                plot([0 0],[ymin ymax],'color',axislcolor); % draw vert axis at time 0
                axis('off');
                plot([xmin xmax],[0 0],'color',axislcolor);  % draw horizontal axis

                %
                %%%%%%%%%%%%%%%%%%%% plot vertical lines (optional) %%%%%%%%%%%%%%%%%
                %

                if isempty(g.vert)
                    g.vert = [xmin xmax];
                    ymean = (ymin+ymax)/2;
                    vmin = ymean-0.1*(ymean-ymin);
                    vmax = vmin*-1;  %ymean+0.2*(ymax-ymean);
                elseif ~isnan(g.vert)
                    ymean = (ymin+ymax)/2;
                    vmin = ymean-0.1*(ymean-ymin);
                    vmax = vmin*-1;  %ymean+0.2*(ymax-ymean);
                    for v = g.vert
                        plot([v v],[vmin vmax],'color',vertcolor); % draw vertical lines
                    end
                end

                %
                %%%%%%%%%%%%%%%%%%%% plot horizontal lines (optional) %%%%%%%%%%%%%%%
                %
                if isempty(g.hori)
                    g.hori = [ymin ymax];
                end
                if ~isnan(g.hori)
                    xmean = 0;
                    hmin = xmean-0.2*(xmean-xmin);
                    hmax = hmin*-1; %xmean+0.3*(xmax-xmean);
                    for v = g.hori
                        plot([hmin hmax],[v v], 'color',horicolor); % draw horizontal lines
                    end
                end

            end

            fprintf(' %d',c); % finished with channel plot
        end; % c, chans / subplot
             % handle legend
        if nolegend, g.legend{P+1} = ['Data ' int2str(P) ]; end

        fprintf('\n');
    end; % P / epoch

    %
    %%%%%%%%%%%%%%%%%%%%% Make time and amp cal bar %%%%%%%%%%%%%%%%%%%%%%%%%
    %
    ax = axes('Units','Normal','Position', ...
              [0.85 0.1 axwidth axheight]); % FIX!!!!
    axes(ax)
    axis('off');
    if xmin <=0
        figure(curfig);p=plot([0 0],[ymn ymx],'color','k'); % draw vert axis at zero
    else
        figure(curfig);p=plot([xmin xmin],[ymn ymx],'color','k'); % draw vert axis at zero
    end
    if g.ydir == -1
        set(gca, 'ydir', 'reverse');
    end
    axis([xmin xmax ymn ymx]);        % set axis values
    hold on
    %set(p, 'Clipping','off');        % center text
    figure(curfig);p=plot([xmin xmax],[0 0],'color',axislcolor); % draw horizontal axis
    axis([xmin xmax ymin ymax]);        % set axis values
                                        %
                                        %%%%%%%%%%%%%%%%%%%% plot vertical lines (optional) %%%%%%%%%%%%%%%%%
                                        %
    if ~isnan(g.vert)
        for v = g.vert
            figure(curfig);plot([v v],[vmin vmax],'color',vertcolor); % draw vertical lines
        end
    end
    %
    %%%%%%%%%%%%%%%%%%%% plot horizontal lines (optional) %%%%%%%%%%%%%%%%%
    %
    if ~isnan(g.hori)
        xmean = 0;
        hmin = xmean-0.2*(xmean-xmin);
        hmax = hmin*-1; %xmean+0.3*(xmax-xmean);
        for v = g.hori
            figure(curfig);plot([hmin hmax],[v v], 'color',horicolor); % draw horizontal lines
        end
    end

    % secondx = 200;                    % draw second vert axis
    % axis('off');plot([secondx secondx],[ylo ymax],'color',axislcolor);
    %
    %%%%%%%%%%%%%%%%%%%%% Plot negative-up %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    signx = xmin-0.15*xdiff;
    figure(curfig);axis('off');h=text(double(signx),double(ymin),num2str(ymin,3)); % text ymin
    set(h,'FontSize',TICKFONTSIZE);               % choose font size
    set(h,'HorizontalAlignment','right','Clipping','off');

    figure(curfig);axis('off');h=text(double(signx), double(ymax),['+' num2str(ymax,3)]);  % text +ymax
    set(h,'FontSize',TICKFONTSIZE);         % choose font size
    set(h,'HorizontalAlignment','right','Clipping','off');

    ytick = g.ydir*(-ymax-0.3*ydiff);
    figure(curfig);tick = [int2str(xmin)]; h=text(double(xmin),double(ytick),tick);
    set(h,'FontSize',TICKFONTSIZE);         % choose font size
    set(h,'HorizontalAlignment','center',...
          'Clipping','off');  % center text

    tick = [xlabel]; figure(curfig);h=text(double(xmin+xdiff/2),double(ytick-0.5*g.ydir*ydiff),tick);
    set(h,'FontSize',TICKFONTSIZE);         % choose font size
    set(h,'HorizontalAlignment','center',...
          'Clipping','off');  % center text

    tick = [int2str(xmax)]; figure(curfig);h=text(double(xmax),double(ytick),tick);
    set(h,'FontSize',TICKFONTSIZE);         % choose font size
    set(h,'HorizontalAlignment','center',...
          'Clipping','off');  % center text

    if length(g.legend) > 1 && strcmpi(g.showleg, 'on')
        tmpleg = vararg2str(g.legend);
        quotes = find(tmpleg == '''');
        for index = length(quotes):-1:1
            tmpleg(quotes(index)+1:end+1) = tmpleg(quotes(index):end);
            tmpleg(quotes(index)) = '''';
        end
        tmpleg = [ 'legend(' tmpleg ');' ];
    else tmpleg = '';
    end
    com = [ 'axis on;' ...
            'clear xlabel ylabel;' tmpleg ...
            'xlabel(''''Time (ms)'''');' ...
            'ylabel(''''Potential (\muV)'''');' ];
    axcopy(gcf, com); % turn on popup feature
                      %
                      %%%%%%%%%%%%%%%%%% Make printed figure fill page %%%%%%%%%%%%%%%%%%%%%%%%%%%
                      %
    orient tall
    % curfig = gcf;
    % h=figure(curfig);
    % set(h,'PaperPosition',[0.2 0.3 7.6 10]); % stretch out the plot on the page

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%FeatureMatch Feature correspondence object
%
% This class represents the correspondence between two PointFeature
% objects.  A vector of FeatureMatch objects can represent the
% correspondence between sets of points.
%
% Methods::
% plot       Plot corresponding points
% show       Show summary statistics of corresponding points
%
% ransac     Determine inliers and outliers
% inlier     Return inlier matches
% outlier    Return outlier matches
% subset     Return a subset of matches
% remove
%
% display    Display value of match
% char       Convert value of match to string
%
% Properties::
% p1         Point coordinates in view 1 (2x1)
% p2         Point coordinates in view 2 (2x1)
% p          Point coordinates in view 1 and 2 (4x1)
% distance   Match strength between the points
%
% Properties of a vector of FeatureMatch objects are returned as a vector.
% If F is a vector (Nx1) of FeatureMatch objects then F.p1 is a 2xN matrix
% with each column the corresponding view 1 point coordinate.
%
% Note::
%  - FeatureMatch is a reference object.
%  - FeatureMatch objects can be used in vectors and arrays
%  - Operates with all objects derived from PointFeature, such as
%    ScalePointFeature, SurfPointFeature and SiftPointFeature.
%
% See also PointFeature, SurfPointFeature, SiftPointFeature.

% TODO:
% distance, strength should be converted to similarity
% p should be p12

% 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/>.

classdef FeatureMatch < handle

    properties
        % the trailing underscore is to distinguish these properties from the methods
        % of almost similar name.  Finding the property of a vector of objects
        % F.x results in a list rather than a vector.  You need to write [F.x] so we
        % create methods to provide this.  Using dependent properties does not work.
        xy_          % x1 y1 x2 y2 of corresponding points
        distance_    % strength of match
        inlier_      % NaN - indeterminate
                    % true - inlier
                    % false - outlier
    end

    methods

        function m = FeatureMatch(f1, f2, s)
        %FeatureMatch.FeatureMatch Create a new FeatureMatch object
        %
        % M = FeatureMatch(F1, F2, S) is a new FeatureMatch object describing a
        % correspondence between point features F1 and F2 with a strength of S.
        %
        % M = FeatureMatch(F1, F2) as above but the strength is set to NaN.
        %
        % Notes::
        % - Only the coordinates of the PointFeature are kept.
        %
        % See also PointFeature, SurfPointFeature, SiftPointFeature.

            if nargin == 0
                m.xy_ = [];
                m.distance_ = [];
                m.inlier_ = [];
                return;
            end

            m.xy_ = [f1.u_ f1.v_ f2.u_ f2.v_]';
            if nargin < 3
                m.distance = NaN;
            else
                m.distance_ = s;
            end
            m.inlier_ = NaN;
        end

        function v = inlier(m)
        %FeatureMatch.inlier Inlier features
        %
        % M2 = M.inlier() is a subset of the FeatureMatch vector M that are
        % considered to be inliers.
        %
        % Notes::
        % - Inliers are not determined until after RANSAC is run.
        %
        % See also FeatureMatch.outlier, FeatureMatch.ransac.
            v = m([m.inlier_] == true);
        end

        function v = outlier(m)
        %FeatureMatch.outlier Outlier features
        %
        % M2 = M.outlier() is a subset of the FeatureMatch vector M that are
        % considered to be outliers.
        %
        % Notes::
        % - Outliers are not determined until after RANSAC is run.
        %
        % See also FeatureMatch.inlier, FeatureMatch.ransac.
            v = m([m.inlier_] == false);
        end

        function v = distance(m)
            v = [m.distance_];
        end

        function v = inlierx(m)
            v = find([m.inlier_]);
        end

        function display(m)
        %FeatureMatch.display Display value
        %
        % M.display() displays a compact human-readable representation of the
        % feature pair.  If M is a vector then the elements are printed one per line.
        %
        % Notes::
        % - This method is invoked implicitly at the command line when the result
        %   of an expression is a FeatureMatch object and the command has no trailing
        %   semicolon.
        %
        % See also FeatureMatch.char.

            disp(' ');
            disp([inputname(1), ' = '])
            disp(' ');
            if length(m) > 20
                fprintf('%d corresponding points (listing suppressed)\n', length(m));
            else
                disp( char(m) );
            end
        end % display()

        function s = char(matches)
        %FeatureMatch.char Convert to string
        %
        % S = M.char() is a compact string representation of the match object.
        % If M is a vector then the string has multiple lines, one per element.

            s = '';
            for m=matches
                ss = sprintf('(%g, %g) <-> (%g, %g), dist=%f', ...
                    m.xy_, m.distance_);
                switch m.inlier_
                case true
                    ss = [ss ' +'];
                case false
                    ss = [ss ' -'];
                end
                s = strvcat(s, ss);
            end
        end

        function s = show(m)
        %FeatureMatch.show Display summary statistics of the FeatureMatch vector
        %
        % M.show() is a compact summary of the FeatureMatch vector M that gives
        % the number of matches, inliers and outliers (and their percentages).
            s = sprintf('%d corresponding points\n', length(m));
            in = [m.inlier_];
            s = [s sprintf('%d inliers (%.1f%%)\n', ...
                sum(in==true), sum(in==true)/length(m)*100)];
            s = [s sprintf('%d outliers (%.1f%%)\n', ...
                sum(in==false), sum(in==false)/length(m)*100) ];
        end

        function v = subset(m, n, varargin)
        %FeatureMatch.subset Subset of matches
        %
        % M2 = M.subset(N) is a FeatureMatch vector with no more than N elements
        % sampled uniformly from M.

            opt.random = false;
            opt = tb_optparse(opt, varargin);

            if opt.random
                i = randi(length(m), n);
            else
                i = round(linspace(1, length(m), n));
            end
            v = m(i);
        end

        function s = p1(m, k)
        %FeatureMatch.p1 Feature point coordinates from view 1
        %
        % P = M.p1() is a 2xN matrix containing the feature points coordinates
        % from view 1.  These are the (u,v) properties of the feature F1 passed
        % to the constructor.
        %
        % See also FeatureMatch.FeatureMatch, FeatureMatch.p2, FeatureMatch.p.
            xy = [m.xy_];
            s = xy(1:2,:);
        end

        function s = p2(m, k)
        %FeatureMatch.p2 Feature point coordinates from view 2
        %
        % P = M.p2() is a 2xN matrix containing the feature points coordinates
        % from view 1.  These are the (u,v) properties of the feature F2 passed
        % to the constructor.
        %
        % See also FeatureMatch.FeatureMatch, FeatureMatch.p1, FeatureMatch.p.

            xy = [m.xy_];
            xy = [m.xy_];
            s = xy(3:4,:);
        end

        function s = p(m, k)
        %FeatureMatch.p Feature point coordinate pairs
        %
        % P = M.p() is a 4xN matrix containing the feature point coordinates.
        % Each column contains the coordinates of a pair of corresponding
        % points [u1,v1,u2,v2].
        %
        % See also FeatureMatch.p1, FeatureMatch.p2.
            s = [m.xy_];
        end

        function plot(m, varargin)
        %FeatureMatch.plot Show corresponding points
        %
        % M.plot() overlays the correspondences in the FeatureMatch vector M
        % on the current figure.  The figure must comprise views 1 and 2 side
        % by side, for example by:
        %
        %      idisp({im1,im2})
        %      m.plot()
        %
        % M.plot(LS) as above but the optional line style arguments LS are
        % passed to plot.
        %
        % Notes::
        % - Using IDISP as above adds UserData to the figure, and an error is
        %   created if this UserData is not found.
        % See also IDISP.

            opt.offset = [];
            [opt,args] = tb_optparse(opt, varargin);

            if isempty(opt.offset)
            try
                ud = get(gca, 'UserData');
                u0 = ud.u0;
            catch
                error('Current image is not a pair displayed by idisp');
            end
            opt.offset = [0 u0(2) 0 0];
            end

            xy = [m.xy_];
            hold on
            for k=1:numcols(xy),
                plot([xy(1,k)+opt.offset(1) xy(3,k)+opt.offset(2)], ...
                    [xy(2,k)+opt.offset(3), xy(4,k)+opt.offset(4)],  args{:});
            end
            hold off
            figure(gcf);
        end % plot

        function [MM,rr] = ransac(m, func, varargin)
        %FeatureMatch.ransac Apply RANSAC
        %
        % M.ransac(FUNC, OPTIONS) applies the RANSAC algorithm to fit the point
        % correspondences to the model described by the function FUNC.  The
        % OPTIONS are passed to the RANSAC() function.  Elements of the
        % FeatureMatch vector have their status updated in place to indicate
        % whether they are inliers or outliers.
        %
        % Example::
        %      f1 = isurf(im1);
        %      f2 = isurf(im2);
        %      m = f1.match(f2);
        %      m.ransac( @fmatrix, 1e-4);
        %
        % See also FMATRIX, HOMOGRAPHY, RANSAC.

            opt.retry = 1;
            [opt,args] = tb_optparse(opt, varargin);
            if opt.verbose
                args = [args 'verbose'];
            end

            while true
                try
                    [M,in,resid] = ransac(func, [m.xy_], args{:});
                    break;
                catch err
                    opt.retry = opt.retry - 1;
                    if opt.retry > 0
                        continue;
                    else
                        rethrow(err);
                    end
                end
            end

            % mark all as outliers
            for i=1:length(m)
                m(i).inlier_ = false;
            end
            for i=in
                m(i).inlier_ = true;
            end

            if nargout >= 1
                MM = M;
            end
            if nargout >= 2
                rr = resid;
            end
        end
    end

end

author: petercorke

repo: machinevision-toolbox-matlab


text:
classdef FCP1 < PROBLEM
% <multi> <real> <constrained>
% Benchmark constrained MOP proposed by Jiawei Yuan

%------------------------------- Reference --------------------------------
% J. Yuan, H. Liu, Y. Ong, and Z. He, Indicator-based evolutionary
% algorithm for solving constrained multi-objective optimization problems,
% IEEE Transactions on Evolutionary Computation, 2022, 26(2): 379-391.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Jiawei Yuan

    methods
        %% Initialization
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D)
                obj.D = 30;
            end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            g = 1 + 9*mean(PopDec(:,2:end),2);
            PopObj(:,1) = PopDec(:,1).*g;
            PopObj(:,2) = (1-PopDec(:,1)).*g;
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            g = 1 + 9*mean(PopDec(:,2:end),2);
            Dis    = abs(9-g);
            %%%%% Type-I constraints
            y1     = Dis.^2-0.25;
            y2     = 1./(Dis+1e-6).*(1.2+sin(Dis*pi));
            PopCon = min([y1,y2],[],2);
        end
        %% Sample reference points on Pareto front
        function P = GetOptimum(obj,N)
            P=8.5*UniformPoint(N,obj.M);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
clear
close all
PAR_DEFAULT_NPROC =  0;
PAR_DEFAULT_TAU    = 0.25;
PAR_DEFAULT_LAMBDA = 0.15;
PAR_DEFAULT_THETA  = 0.3;
PAR_DEFAULT_NSCALES =5;
PAR_DEFAULT_ZFACTOR =0.5;
PAR_DEFAULT_NWARPS = 5;
PAR_DEFAULT_EPSILON = 0.01;
useOracle = 'a';
lrdist = 0;
hBin = 0;
h02 = 0;
[width, height, nframes, nchannels, sequence] = ReadInputData('traffic.avi');
fparams.tau     = PAR_DEFAULT_TAU;
fparams.theta   = PAR_DEFAULT_THETA;
fparams.nscales = PAR_DEFAULT_NSCALES;
fparams.zfactor = PAR_DEFAULT_ZFACTOR;
fparams.warps  = PAR_DEFAULT_NWARPS;
fparams.epsilon = PAR_DEFAULT_EPSILON;
fparams.verbose = 0;
fparams.iflagMedian = 0;
N = 1 + log(hypot(width, height)/16.0) / log(1/fparams.zfactor);
%
   options.osigma = 	{"s:", 0, "5.0", [], "noise standard deviation"};

    options.oiFrame = {"i:", 0, "-1", [], "frame to denoise (-1: denoise all frames)"};

    %%Parameters
    options.oibloc = 	{"b:", 0, "12", [], "radius of search region"};

    options.oiwin = 	{"w:", 0, "2", [], "radius of patch"};

   options.oitemp = 	{"t:", 0, "7", [], "radius of temporal neighborhood"};

   options.oiknn = 	{"k:", 0, [], [], "minimum number of patches (recommended: 55 gray images, 95 color images)"};

   options.oflat = 	{"f:", 0, "0.85", [], "flat parameter"};

   options.olrdist = {"c:", 0, "1.0", [], "threshold for left-right coherence in occlusions mask"};

    options.ohbin = 	{"h:", 0, "0.5", [], "occlusion binarization threshold"};

    options.ofocc = 	{"o:", 0, "5.5", [], "occlusion factor"};

    options.ofpca1 = 	{"p:", 0, "1.8",[], "PCA factor 1st step"};

   options.ofpca2 = 	{"q:", 0, "1.45", [], "PCA factor 2nd step"};

   options.odist1 = 	{"d:", 0, "0.0", [], "3D blocks distances 1st step"};

    options.odist2 = 	{"e:", 0, "2.0",[], "3D blocks distances 2nd step"};

   options.olambda1 = 	{"l:", 0, "0.075", [], "optical flow lambda 1st step"};

   options.olambda2 = 	{"m:", 0, "0.15", [], "optical flow lambda 2nd step"};
  parameter.pinput = {"input", [], "input file"};

    parameter.pout = {"out", [], "output file"};

    %video Denoising Parameters not setted yet
    dparams.fSigma = str2double(options.osigma{3});
    iFrame = str2double(options.oiFrame{3});
    iscolor= nchannels;
    iBloc = 2*str2double(options.oibloc{3})+1;
    iWin  = 2*str2double(options.oiwin{3})+1;
    iTemp = 2*str2double(options.oitemp{3})+1;
    iKnn = 0;

   if(iscolor == 1)
       iKnn = 95;
   else
       iKnn = 55;
   end
    hBinTh =  str2double(options.ohbin{3});
    factorocc =  str2double(options.ofocc{3});
    factorflat = str2double(options.oflat{3});
    lrdistTh =  str2double(options.olrdist{3});
    factorPCA_1st = str2double(options.ofpca1{3});
    thdist_1st =  str2double(options.odist1{3});
    lambda_1st = str2double(options.olambda1{3});
    factorPCA_2nd = str2double(options.ofpca2{3});
    thdist_2nd = str2double(options.odist2{3});
    lambda_2nd = str2double(options.olambda2{3});
    % Global parameters denoising part %
    dparams.iFrames=nframes;
    dparams.iBloc = iBloc;
    dparams.iWin  = iWin;
    dparams.iTemp = iTemp;
    dparams.useFlatPar=factorflat;
    dparams.ifK = iKnn;
    % first iteration parameters %
    dparams.useOracle = 0;

    dparams.fRMult =  factorPCA_1st;
    dparams.fFixedThrDist = thdist_1st;
    fparams.lambda = lambda_1st;

    useOracle=dparams.useOracle;
    % denoising function %
    denoise_function(fparams, dparams, -1,nframes,useOracle,sequence);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% %Question No: 5 (a)
% Determine the weights of a network with 4 input and 2 output units using
% Perceptron Learning Law for the following input-output pairs:

% Input: [1100]' [1001]' [0011]' [0110]'
% output: [11]' [10]' [01]' [00]'
% Discuss your results for different choices of the learning rate
% parameters.
% Use suitable values for the initial weights.

in=[1 1 0 0 -1;1 0 0 1 -1; 0 0 1 1 -1; 0 1 1 0 -1];
out=[1 1; 1 0; 0 1; 0 0];
eta=input('Enter the learning rate value = ');
it=input('Enter the number of iterations required = ');
wgt=input('Enter the weights,2 by 5 matrix(including weight for bias):\n');
for x=1:it
    for i=1:4
        s1=0;
        s2=0;
        for j=1:5
          s1=s1+in(i,j)*wgt(1,j);
          s2=s2+in(i,j)*wgt(2,j);
        end
        wi=eta*(out(i,1)-sign(s1))*in(i,:);
        wgt(1,:)=wgt(1,:)+wi;
        wi=eta*(out(i,2)-sign(s2))*in(i,:);
        wgt(2,:)=wgt(2,:)+wi;
    end
end
wgt


text:
function [xDwCentr,yDwCentr,zDwCentr] = calc_doseWeightedCentroid(structNum,doseNum,normalizeFlag,planC)
%function [xDwCentr,yDwCentr,zDwCentr] = calc_doseWeightedCentroid(structNum,doseNum,normalizeFlag,planC)
%
% This function returns the dose-weighted centroid for the input structure and dose.
%
% INPUTS:
% structNum: structure index within the planC object
% doseNum: dose index within the planC object
% normalizeFlag: Flag to normalize the centroid between 0 and 1.
% If normalizeFlag == 1, the structure coordinates are scaled between and 0
% and 1 based on min/max coordinates along that dimension. x,y,z are based
% on RTOG coordinate system: https://github.com/cerr/CERR/wiki/Coordinate-system
% planC: can be accessed from Viewer by using 'global planC' or by loading
% the CERR .mat file in memory.
%
% OUTPUTS:
% x,y,z: dose weighted centroid coordinates
%
% Example:
%
% structNum = 10;
% doseNum = 1;
% normFlg = 1;
% global planC
% [x,y,z] = calc_doseWeightedCentroid(structNum,doseNum,normFlg,planC)
%
%APA, 05/15/2018
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Check if plan passed, if not use global.
if ~exist('planC')
    global planC;
end
indexS = planC{end};

[iV,jV,kV] = getUniformStr(structNum, planC);
assocScanNum = getStructureAssociatedScan(structNum,planC);
[xV,yV,zV] = getUniformScanXYZVals(planC{indexS.scan}(assocScanNum));
xCalcV = xV(jV);
yCalcV = yV(iV);
zCalcV = zV(kV);

%Get scan transM
transM_scan = getTransM('scan',assocScanNum,planC);
if isempty(transM_scan)
    transM_scan = eye(4);
end

%Get dose transM
transM_dose = getTransM('dose',doseNum,planC);
if isempty(transM_dose)
    transM_dose = eye(4);
end

%Apply transM to dose calc pts
if ~isequal(transM_scan,transM_dose)
    [xCalcV, yCalcV, zCalcV] = applyTransM(transM_scan*inv(transM_dose), xCalcV, yCalcV, zCalcV);
    [jnk,jnk,zV] = applyTransM(transM_scan*inv(transM_dose), zV*0, zV*0 ,zV);
end

dosesV = getDoseAt(doseNum, xCalcV, yCalcV, zCalcV, planC);

%Normalize x,y,z coordinates
if normalizeFlag
    xCalcV = (xCalcV - min(xCalcV)) ./ (max(xCalcV) - min(xCalcV));
    yCalcV = (yCalcV - min(yCalcV)) ./ (max(yCalcV) - min(yCalcV));
    zCalcV = (zCalcV - min(zCalcV)) ./ (max(zCalcV) - min(zCalcV));
end

xDwCentr = sum(dosesV.*xCalcV)/sum(dosesV);
yDwCentr = sum(dosesV.*yCalcV)/sum(dosesV);
zDwCentr = sum(dosesV.*zCalcV)/sum(dosesV);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [beta_gibbs,F_gibbs,L_gibbs,phi_gibbs,sigma_gibbs,lambda_t_gibbs,sigma_t_gibbs,sbar,favar,It,Bu]=...
    favar_stvol3gibbs(Xbart,Xt,yt,B0,phi0,alpha0,delta0,f0,upsilon0,betahat,sigmahat,gamma,G,I_o,omega,T,n,k,It,Bu,pick,pickf,favar,data_endo,lags)

%% preliminaries
% initialise variables
nfactorvar=favar.nfactorvar;
numpc=favar.numpc;
favarX=favar.X(:,favar.plotX_index);
onestep=favar.onestep;
Sigma=bear.nspd(favar.Sigma);
Ll=favar.L;
favar_X=favar.X;
% load priors
L0=favar.L0*eye(n);
a0=favar.a0;
b0=favar.b0;
% sigmahat=(1/T)*(EPS'*EPS);

% preallocation
Ll_gibbs=zeros(size(Ll(:),1),It-Bu);
R2_gibbs=zeros(size(favarX,2),It-Bu);

if onestep==0 %static factors in this case
    FY=data_endo;
    pbstring='two-step'; %string for the progress bar
    % elseif onestep==1
    %     pbstring='one-step'; %string for the progress bar
end

% preliminary elements for the algorithm
% compute the product G'*I_gamma*G (to speed up computations of deltabar)
GIG=G'*I_o*G;
% compute alphabar
alphabar=T+alpha0;

% initiate the Gibbs sampler
% initiate the counting of iterations
count=1;
pickcount=1;
% initiate the record matrices and cells
beta_gibbs=[];
F_gibbs=[];
L_gibbs=[];
phi_gibbs=[];
sigma_gibbs=[];
lambda_t_gibbs={};
sigma_t_gibbs={};

% step 1: determine initial values for the algorithm

% initial value for beta
beta=betahat;
B=reshape(beta,k,n);
% initial value for f_2,...,f_n
% obtain the triangular factorisation of sigmahat
[Fhat Lambdahat]=bear.triangf(sigmahat);
% obtain the initial value for F
F=Fhat;
% obtain the inverse of Fhat
[invFhat]=bear.invltod(Fhat,n);
% create the cell storing the different vectors of invF
Finv=cell(n,1);
% store the vectors
for ii=2:n
    Finv{ii,1}=invFhat(ii,1:ii-1);
end
% initial values for L
L=zeros(T,1);
% initial values for phi
phi=1;

% step 2: determine the sbar values and Lambda
sbar=diag(Lambdahat);
Lambda=sparse(diag(sbar));
% then determine sigma^(0)
sigma=F*Lambda*F';

% step 3: recover the series of initial values for lambda_1,...,lambda_T and sigma_1,...,sigma_T
lambda_t=repmat(diag(sbar),1,1,T);
sigma_t=repmat(sigmahat,1,1,T);

% create a progress bar
hbar = bear.parfor_progressbar(It,['Progress of the Gibbs sampler (',pbstring,').']);

%% run the Gibbs sampler
while count<=It

    % step 4: draw beta from its conditional posterior
    % first compute the summations required for omegabar and betabar
    summ1=zeros(k,k);
    summ2=zeros(k,n);
    % run the summation
    for jj=1:T
        prodt=Xt{jj,1}'*exp(-L(jj,1));
        summ1=summ1+prodt*Xt{jj,1};
        summ2=summ2+prodt*yt(:,:,jj)';
    end
    % then obtain the inverse of phi0
    invphi0=diag(1./diag(phi0));
    % obtain the inverse of phibar
    invphibar=summ1+invphi0;
    % recover phibar
    C=chol(bear.nspd(invphibar),'Lower')';
    invC=C\speye(k);
    phibar=invC*invC';
    % recover Bbar
    Bbar=phibar*(summ2+invphi0*B0);
    % draw B from its posterior
    B=bear.matrixndraw(Bbar,sigma,phibar,k,n);
    % finally recover beta by vectorising
    beta=B(:);

    % step 5: draw the series f_2,...,f_n from their conditional posteriors
    % recover first the residuals
    for jj=1:T
        epst(:,:,jj)=yt(:,:,jj)-Xbart{jj,1}*beta;
    end
    % then draw the vectors in turn
    for jj=2:n
        % first compute the summations required for upsilonbar and fbar
        summ1=zeros(jj-1,jj-1);
        summ2=zeros(jj-1,1);
        % run the summation
        for kk=1:T
            prodt=epst(1:jj-1,1,kk)*exp(-L(kk,1));
            summ1=summ1+prodt*epst(1:jj-1,1,kk)';
            summ2=summ2+prodt*epst(jj,1,kk)';
        end
        summ1=(1/sbar(jj,1))*summ1;
        summ2=(-1/sbar(jj,1))*summ2;
        % then obtain the inverse of upsilon0
        invupsilon0=diag(1./diag(upsilon0{jj,1}));
        % obtain upsilonbar
        invupsilonbar=summ1+invupsilon0;
        C=chol(bear.nspd(invupsilonbar));
        invC=C\speye(jj-1);
        upsilonbar=full(invC*invC');
        % recover fbar
        fbar=upsilonbar*(summ2+invupsilon0*f0{jj,1});
        % finally draw f_i^(-1)
        Finv{jj,1}=fbar+chol(bear.nspd(upsilonbar),'lower')*randn(jj-1,1);
    end
    % recover the inverse of F
    invF=eye(n);
    for jj=2:n
        invF(jj,1:jj-1)=Finv{jj,1};
    end
    % eventually recover F
    F=bear.invltod(invF,n);
    % update sigma
    sigma=F*Lambda*F';

    % step 6: draw phi from its conditional posterior
    % estimate deltabar
    deltabar=L'*GIG*L+delta0;
    % draw the value phi_i
    phi=bear.igrandn(alphabar/2,deltabar/2);

    % step 7: draw the series lambda_t from their conditional posteriors, t=1,...,T
    % consider periods in turn
    for kk=1:T
        % a candidate value will be drawn from N(lambdabar,phibar)
        % the definitions of lambdabar and phibar varies with the period, thus define them first
        % if the period is the first period
        if kk==1
            lambdabar=(gamma*L(2,1))/(1/omega+gamma^2);
            phibar=phi/(1/omega+gamma^2);
            % if the period is the final period
        elseif kk==T
            lambdabar=gamma*L(T-1,1);
            phibar=phi;
            % if the period is any period in-between
        else
            lambdabar=(gamma/(1+gamma^2))*(L(kk-1,1)+L(kk+1,1));
            phibar=phi/(1+gamma^2);
        end
        % now draw the candidate
        cand=lambdabar+phibar^0.5*randn;
        % compute the acceptance probability
        prob=bear.mhprob3(cand,L(kk,1),sbar,epst(:,1,kk),Finv,n);
        % draw a uniform random number
        draw=rand;
        % keep the candidate if the draw value is lower than the prob
        if draw<=prob
            L(kk,1)=cand;
            % if not, just keep the former value
        end
    end
    % then recover the series of matrices lambda_t and sigma_t
    for kk=1:T
        lambda_t(:,:,kk)=exp(L(kk,1))*diag(sbar);
        sigma_t(:,:,kk)=F*lambda_t(:,:,kk)*F';
    end

    %% Sample Sigma and Ll (static)
    [Sigma,Ll]=bear.favar_SigmaL(Sigma,Ll,nfactorvar,numpc,onestep,n,favar_X,FY,a0,b0,T,lags,L0);

    %% record phase
    % if the burn-in sample phase is not yet over
    if count<=Bu
        % simply add 1 to the iteration count
        count=count+1;
        % on the other hand, if the burn-in sample phase is over
    elseif count>Bu
        % adding one iteration to the count will depend on wether post-burn selection applies
        % if there is no post burn selection
        if pick==0
            % record the results
            beta_gibbs(:,count-Bu)=beta;
            F_gibbs(:,:,count-Bu)=F;
            L_gibbs(:,:,count-Bu)=L;
            phi_gibbs(count-Bu,1)=phi;
            sigma_gibbs(:,count-Bu)=sigma(:);
            % save the factors and loadings (keep the notation in the code consistent, although - except L - they don't change)
            Ll_gibbs(:,count-Bu)=Ll(:);

            % compute R2 (Coefficient of Determination) for plotX variables (keep the notation in the code consistent, although they don't change)
            R2=bear.favar_R2(favarX,FY);
            R2_gibbs(:,count-Bu)=R2(:);
            for jj=1:T
                lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
                sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
            end
            % then add one to the count
            count=count+1;
            % if there is post burn selection, only one draw over 'fpick' draws will be retained
        elseif pick==1
            % if the iteration does not correspond to fpick, don't record the results, don't increase the regular count, but do increase pickcount by 1, and do record the acceptance rate of the Metropolis-Hastings step
            if pickcount~=pickf
                pickcount=pickcount+1;
                % on the other hand, if the iteration does correspond to fpick
            elseif pickcount==pickf
                % do record the results
                beta_gibbs(:,count-Bu)=beta;
                F_gibbs(:,:,count-Bu)=F;
                L_gibbs(:,:,count-Bu)=L;
                phi_gibbs(count-Bu,1)=phi;
                sigma_gibbs(:,count-Bu)=sigma(:);
                % save the factors and loadings (keep the notation in the code consistent, although - except L - they don't change)
                Ll_gibbs(:,count-Bu)=Ll(:);

                % compute R2 (Coefficient of Determination) for plotX variables (keep the notation in the code consistent, although they don't change)
                R2=bear.favar_R2(favarX,FY);
                R2_gibbs(:,count-Bu)=R2(:);
                for jj=1:T
                    lambda_t_gibbs{jj,1}(:,:,count-Bu)=lambda_t(:,:,jj);
                    sigma_t_gibbs{jj,1}(:,:,count-Bu)=sigma_t(:,:,jj);
                end
                % then increase the regular count by 1 and re-initialise pickcount
                count=count+1;
                pickcount=1;
            end
        end
    end

    % update progress by one iteration
    hbar.iterate(1);

end

% in case we have thinning of the draws,
thin=abs(round(favar.thin)); % should be a positive integer
if thin~=1
    beta_gibbs=beta_gibbs(:,thin:thin:end);
    F_gibbs=F_gibbs(:,thin:thin:end);
    L_gibbs=L_gibbs(:,thin:thin:end);
    phi_gibbs=phi_gibbs(:,thin:thin:end);
    sigma_gibbs=sigma_gibbs(:,thin:thin:end);
    Ll_gibbs=Ll_gibbs(:,thin:thin:end);
    R2_gibbs=R2_gibbs(:,thin:thin:end);
    for jj=1:T
        lambda_t_gibbs=lambda_t_gibbs(:,:,thin:thin:end);
        sigma_t_gibbs=sigma_t_gibbs(:,:,thin:thin:end);
    end
    It=(1/thin)*It;
    Bu=(1/thin)*Bu;
end

% save in favar structure
favar.L_gibbs=Ll_gibbs;
favar.R2_gibbs=R2_gibbs;

close(hbar);   %close progress bar

author: european-central-bank

repo: BEAR-toolbox


text:
function [i,j] = xy2ind(ebsd,x,y)
% convert x,y coordinates into indeces of ebsd
%
% Syntax
%
%   ind = xy2ind(ebsd,x,y)
%   ebsd(ind)
%
%   [i,j] = ind = xy2ind(ebsd,x,y)
%   ebsd(i,j)
%
% Input
%  ebsd - @EBSDsquare
%  x,y  - spatial coordinates
%
% Output
%  ind  - index to @EBSDsquare
%  i,j  - indeces to @EBSDsquare
%

if nargin == 2
  y = x(:,2);
  x = x(:,1);
end

i = 1+round((y - ebsd.prop.y(1))./ebsd.dy);
j = 1+round((x - ebsd.prop.x(1))./ebsd.dx);

if nargout == 1, i = sub2ind(size(ebsd),i,j); end

author: mtex-toolbox

repo: mtex

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


text:
% pop_comments() - edit comments
%
% Usage:
%   >> newcomments = pop_comments( oldcomments);
%   >> newcomments = pop_comments( oldcomments, title, newcomments, concat);
%
% Inputs:
%   oldcomments - old comments (string or cell array of strings)
%   title       - optional window title (string)
%   newcomments - new comments (string or cell array of strings)
%                 to assign (during commandline calls only)
%   concat      - [0|1] 1 concatenate the newcomments to the old one.
%                 Default is 0.
%
% Outputs:
%   newcomments - new comments, string
%
% Note: if new comments are given as input, there are simply
%       converted and returned by the function; otherwise a
%       window pops up.
%
% Example
%  EEG.comments = pop_comments( { 'This is the first line.' ' ' ...
%               'This is the third line.' }, 'Editing');
% EEG.comments = pop_comments(EEG.comments,'','This is the fourth line.",1);
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: eeglab()

% Copyright (C) 2001 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

% 01-25-02 reformated help & license -ad
% 03-16-02 text interface editing -sm & ad

function [newcomments, com] = pop_comments( comments, plottitle, newcomments, concat );

com = '';
if exist('comments') ~=1, comments = '';
elseif iscell(comments), comments = strvcat(comments{:});
end;

% remove trailing blanks and make multiline
comments = strmultiline( comments, 53);

if nargin < 3
    newcomments = comments;
	try, icadefs;
	catch,
		BACKCOLOR  =  [.8 .8 .8];
		GUIBUTTONCOLOR   = [.8 .8 .8];
	end;
	figure('menubar', 'none', 'tag', 'comment', 'color', BACKCOLOR, 'userdata', 0, ...
		   'numbertitle', 'off', 'name', 'Read/Enter comments -- pop_comments()');
	pos = get(gca,'position'); % plot relative to current axes
	q = [pos(1) pos(2) 0 0];
	s = [pos(3) pos(4) pos(3) pos(4)]./100;
	if exist('plottitle') ~=1, plottitle = ''; end;

	h = title(plottitle);
	set(h, 'fontname','Helvetica','fontweight', 'bold', 'interpreter', 'none');

	axis off;

	% create the buttons
	% ------------------
  	uicontrol('Parent',gcf, ...
  	'Units','Normalized', ...
	'Position', [0 -5 20 10].*s+q, ...
	'backgroundcolor', GUIBUTTONCOLOR, ...
	'string','CANCEL', 'callback', 'close(findobj(''tag'', ''comment''));' );

  	uicontrol('Parent',gcf, ...
  	'Units','Normalized', ...
	'Position', [80 -5 20 10].*s+q, ...
	'backgroundcolor', GUIBUTTONCOLOR, ...
	'string','SAVE', 'callback', ...
		[ 'set(gcbf, ''userdata'', ' ...
		'get(findobj(''parent'', gcbf, ''tag'', ''edit''), ''string''));' ]);

	%hh = text( q(1), 100*s(2)+q(2), comments, 'tag', 'edit');
	%set( hh, 'editing', 'on', 'verticalalignment', 'top');

    %hh = uicontrol('Parent',gcf, ...
  	%'Units','Normalized', ...
  	%'style', 'text', ...
	%'Position', [0 100 105 5].*s+q, ...
	%'string', 'Warning: each blank line must contain at least a ''space'' character', ...
	%'horizontalalignment', 'left', ...
    %'backgroundcolor', BACKCOLOR );

    hh = uicontrol('Parent',gcf, ...
  	'Units','Normalized', ...
  	'style', 'edit', ...
  	'tag', 'edit', ...
	'Position', [0 10 105 85].*s+q, ...
	'string', comments, ...
	'backgroundcolor', [ 1 1 1], ...
	'horizontalalignment', 'left', ...
	'max', 3, ...
	'fontsize', 12);

    % Try to use 'courier' since it has constant character size
    lf = listfonts;
    tmppos = strmatch('Courier', lf);
    if ~isempty(tmppos)
        set(hh, 'fontname', lf{tmppos(1)}, 'fontsize', 10);
    end;

    waitfor(gcf, 'userdata');

    % find return mode
    if isempty(get(0, 'currentfigure')), return; end;
    tmp = get(gcf, 'userdata');
    if ~isempty(tmp) & isstr(tmp)
        newcomments = tmp; % ok button
    else return;
    end;

	close(findobj('tag', 'comment'));
else
    if iscell(newcomments)
        newcomments = strvcat(newcomments{:});
    end;
    if nargin > 3 & concat == 1
        newcomments = strvcat(comments, newcomments);
    end;
    return;
end;

I = find( comments(:) == '''');
comments(I) = ' ';
if nargout > 1
        if ~strcmp( comments, newcomments)
          allsame = 1;
            for index = 1:size(comments, 1)
                if ~strcmp(comments(index,:), newcomments(index,:)), allsame = 0; end;
            end;
        else
            allsame = 0;
        end;
        if allsame & ~isempty(comments)
             com =sprintf('EEG.comments = pop_comments(EEG.comments, '''', %s, 1);', vararg2str(newcomments(index+1:end,:)));
        else
            com =sprintf('EEG.comments = pop_comments('''', '''', %s);', vararg2str(newcomments));
        end;
end;
return;

author: PatternRecognition

repo: OpenBMI


text:
% DEMSWISSROLLFULLLLE5 Demonstrate LLE on the oil data.

% MLTOOLS

[Y, lbls] = lvmLoadData('swissRollFull');
%Y = Y(1:20, :);
options = lleOptions(4);
options.acyclic=true;
options.isNormalised = false;
model = lleCreate(2, size(Y, 2), Y, options);
model = lleOptimise(model, 2);

%lvmScatterPlotColor(model, model.Y(:, 2));

%if exist('printDiagram') & printDiagram
%  lvmPrintPlot(model, model.Y(:, 2), 'SwissRollFull', 1, true);
%end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = solver_apply( ndxs, func, varargin )
global tfocs_count___

%SOLVER_APPLY	Internal TFOCS routine.
%	A wrapper function used to facilitate the counting of function calls, and
%   to standardize the number of output arguments for the algorithms.

[ varargout{1:max(nargout,1)} ] = func( varargin{:} );
tfocs_count___(ndxs) = tfocs_count___(ndxs) + 1;

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

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
% IndexToAssignment Convert index to variable assignment.
%
%   A = IndexToAssignment(I, D) converts an index, I, into the .val vector
%   into an assignment over variables with cardinality D. If I is a vector,
%   then the function produces a matrix of assignments, one assignment
%   per row.
%
%   See also AssignmentToIndex.m and SampleFactors.m

function A = IndexToAssignment(I, D)

D = D(:)'; % ensure that D is a row vector
A = mod(floor(repmat(I(:) - 1, 1, length(D)) ./ repmat(cumprod([1, D(1:end - 1)]), length(I), 1)), ...
        repmat(D, length(I), 1)) + 1;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function y = im2jpeg(x, quality)
%IM2JPEG Compresses an image using a JPEG approximation.
%   Y = IM2JPEG(X, QUALITY) compresses image X based on 8 x 8 DCT
%   transforms, coefficient quantization, and Huffman symbol
%   coding. Input QUALITY determines the amount of information that
%   is lost and compression achieved.  Y is an encoding structure
%   containing fields:
%
%      Y.size      Size of X
%      Y.numblocks Number of 8-by-8 encoded blocks
%      Y.quality   Quality factor (as percent)
%      Y.huffman   Huffman encoding structure, as returned by
%                  MAT2HUFF
%
%   See also JPEG2IM.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/10/26 18:37:41 $

error(nargchk(1, 2, nargin));             % Check input arguments
if ndims(x) ~= 2 | ~isreal(x) | ~isnumeric(x) | ~isa(x, 'uint8')
   error('The input must be a UINT8 image.');
end
if nargin < 2
   quality = 1;   % Default value for quality.
end

m = [16 11  10  16  24  40  51  61        % JPEG normalizing array
     12  12  14  19  26  58  60  55       % and zig-zag redordering
     14  13  16  24  40  57  69  56       % pattern.
     14  17  22  29  51  87  80  62
     18  22  37  56  68  109 103 77
     24  35  55  64  81  104 113 92
     49  64  78  87  103 121 120 101
     72  92  95  98  112 100 103 99] * quality;

order = [1 9  2  3  10 17 25 18 11 4  5  12 19 26 33  ...
        41 34 27 20 13 6  7  14 21 28 35 42 49 57 50  ...
        43 36 29 22 15 8  16 23 30 37 44 51 58 59 52  ...
        45 38 31 24 32 39 46 53 60 61 54 47 40 48 55  ...
        62 63 56 64];

[xm, xn] = size(x);                % Get input size.
x = double(x) - 128;               % Level shift input
t = dctmtx(8);                     % Compute 8 x 8 DCT matrix

% Compute DCTs of 8x8 blocks and quantize the coefficients.
y = blkproc(x, [8 8], 'P1 * x * P2', t, t');
y = blkproc(y, [8 8], 'round(x ./ P1)', m);

y = im2col(y, [8 8], 'distinct');  % Break 8x8 blocks into columns
xb = size(y, 2);                   % Get number of blocks
y = y(order, :);                   % Reorder column elements

eob = max(x(:)) + 1;               % Create end-of-block symbol
r = zeros(numel(y) + size(y, 2), 1);
count = 0;
for j = 1:xb                       % Process 1 block (col) at a time
   i = max(find(y(:, j)));         % Find last non-zero element
   if isempty(i)                   % No nonzero block values
      i = 0;
   end
   p = count + 1;
   q = p + i;
   r(p:q) = [y(1:i, j); eob];      % Truncate trailing 0's, add EOB,
   count = count + i + 1;          % and add to output vector
end

r((count + 1):end) = [];           % Delete unusued portion of r

y.size      = uint16([xm xn]);
y.numblocks = uint16(xb);
y.quality   = uint16(quality * 100);
y.huffman   = mat2huff(r);

repo: weiyanmin


text:
function [sqrtx,sqrtinvx,di] = randiwishart(sigma,df,di)
%RANDIWISHART Generate inverse Wishart random matrix
%   W=RANDIWISHART(SIGMA,DF) generates a random matrix W from the inverse
%   Wishart distribution with parameters SIGMA and DF.  The inverse of W
%   has the Wishart distribution with covariance matrix inv(SIGMA) and DF
%   degrees of freedom.
%
%   W=RANDIWISHART(SIGMA,DF,DI) expects DI to be the Cholesky factor of
%   the inverse of SIGMA.
%
%   [W,DI]=RANDIWISHART(SIGMA,DF) returns DI so it can be used again in
%   future calls to RANDIWISHART.

n = size(sigma,1);
if (df<n) % require this to ensure invertibility
   error('randiwish:BadDf',...
         'Degrees of freedom must be no smaller than the dimension of SIGMA.');
end

% Get Cholesky factor for inv(sigma) unless that's already done
if nargin<3
    %     [d,p] = chol(sigma,0);
    %     if p~=0
    %         error('stats:iwishrnd:BadCovariance',...
    %             'Covariance matrix must be symmetric and positive definite.');
    %     end
    d = chol(sigma);
    di = d'\eye(size(d));  % either take inverse here and scale chol of
    %randwishart sample and then take inverse of sample, or take inverse of
    %sample and then scale after w/o the inverse.
end

a = randwishart(df/2,n);
sqrtinvx = sqrt(2)*a*di;
sqrtx = (sqrtinvx\eye(size(sqrtinvx)))';

% x = 2*(a'*a);
% x = x\eye(size(x));
% x = d'*(x*d);

% sqrtx = sqrt(2)*a;
% sqrtx = (sqrtx\eye(size(sqrtx)))';
% sqrtx = sqrtx*d;

author: michaelchughes

repo: NPBayesHMM


text:
%% Check rate of convergence for 3D H1 interface problem
%     -div(A grad u)  = f,    x\in \Omega
%      where A is a piecewise constant on Omega^+ and Omega^-.
%
% Domain: Rectangular domain: [xmin,xmax] X [ymin,ymax] X [zmin,zmax]
% Mesh: Cartesian triangular mesh.
% Method: VIFE

%% Geometry and Boundary Conditions
%clear
%close all
%clc

%path(pathdef)
addpath(genpath(pwd),'-begin');
rmpath(genpath('./.git'));
rmpath(genpath('./docs'));
savepath;

domain = [-1,1,-1,1,-1,1];
bc = [1,1,1,1,1,1]; % Dirichelet BC

%% Finite Element Type
femtype = 'P1';
disp(['FEM Type =  Conforming ', femtype]);

%% Initial Partition
nx0 = 10;
ny0 = nx0;
nz0 = nx0;

%% Task
showErr = 0;
showMesh = 0;
computErr = 1;

%% PDE
test = 5;
switch test
    case 0
        pde = poissonNonPoly3D;
    case 1 % circular interface
        %r = pi/5; bm = 1; bp = pi/(2*r^2);
        r = sqrt(pi/2); bm = 1; bp = pi/(2*r^2);domain = [-2,2,-2,2,-2,2];
        x0 = 0; y0 = 0; z0 = 0; rx = r; ry = r; rz = r;
        pde = elli3DcircIntf(bm,bp,r,x0,y0,z0,rx,ry,rz);
    case 2 % orthotorus interface
        domain = [-1.2,1.2,-1.2,1.2,-1.2,1.2];
        bm = 1; bp = 1;
        rx = 1; ry = 0.075; rz = 3;
        pde = elli3DorthocircIntf(bm,bp,rx,ry,rz);
    case 3 % line interface
        bm = 1; bp = 10;
        rx = 1; ry = 0; rz = 0; cx = pi/10; cy = 0; cz = 0;
        cm = 1; cp = 1; a = 1;
        pde = elli3DlinIntf(bm,bp,cx,cy,cz,rx,ry,rz,a,cm,cp);
    case 4 % line interface but with zero boundary conditions
        bm = 1; bp = 10^(2); delta = pi/10; % delta can not be -1,1
        pde = elli3DlinIntf2(bm,bp,delta);
    case 5 % circular interface
        r = pi/4; bm = 1; bp = 10; domain = [-1,1,-1,1,-1,1];
        x0 = 0; y0 = 0; z0 = 0;
        pde = elli3DcircIntf5(bm,bp,r,x0,y0,z0);
    case 6 % circular interface
        bm = 1; bp = 10; domain = [-1.3,1.3,-1.3,1.3,-1.3,1.3];
        x1 = -0.3; y1 = 0; z1 = 0; r11 = pi/5; r12 = 0.2;
        x2 =  0.3; y2 = 0; z2 = 0; r21 = pi/5; r22 = 0.2;
        pde = elli3DtorusTwin(bm,bp,x1,y1,z1,r11,r12,x2,y2,z2,r21,r22);
end

%% Max Iteration
maxIt = 4;
time = zeros(9,maxIt);
error = zeros(maxIt,4);
ratio = zeros(maxIt,4);

for i = 1:maxIt

    %% 1. Generate Mesh
    tic
    nx = nx0 + 10*(i-1); h = (domain(2) - domain(1))/nx;
    ny = ny0 + 10*(i-1);
    nz = nz0 + 10*(i-1);
    time(1,i) = nx;
    disp(' ')
    disp('**************************************************************************************')
    disp(['Partition =  ',int2str(nx),' X ',int2str(ny),' X ',int2str(nz)]);
    disp(' ')

    mesh = genMesh3D(domain, nx, ny, nz);
    mesh = enrichMesh3D(mesh,1); % Mesh detail level = 2 (for PPIFE).
    disp(['number of element =  ', int2str(length(mesh.t))]);
    mesh = genIntfMesh3D(mesh,pde.intf);
    disp(['number of interface element =  ', int2str(-min(mesh.tLoc)),...
        ', is ', num2str(100*-min(mesh.tLoc)/length(mesh.t)), '% of all elements']);
    time(2,i) = toc;

    tic
    fem = genFEM3D(mesh,femtype,bc);
    disp(['number of DoF =  ', int2str(length(fem.p))]);
    time(3,i) = toc;

    %% 2. Generate FEMI data including IFE functions and their matrices
    tic
    meshI = genIVmesh(mesh);
    time(4,i) = toc;
    length(unique(union(union(meshI.tface(:,1),meshI.tface(:,2)),meshI.tface(:,3))))

    tic
    option.mass = 0; % generate mass matrix
    option.rhs = 1;
    femI = genP1VIFEM3D(meshI,mesh,pde,option);
    time(5,i) = toc;

    %% 3. Assemble Matrix
    tic
    disp(' '); disp('Start Assembling Matrix');
    option.mass = 0; % generate mass matrix
    matrix = genMatVIFE3D(pde,mesh,fem,meshI,femI);
    time(6,i) = toc;

    %% 4. Solve the linear system Au = f
    tic
    disp(' ')

    disp('Start Solving Linear System: using PCG with ichol precond');
    %L = ichol(matrix.A);
    alpha = 10;
    L = ichol(matrix.A,struct('type','ict','droptol',1e-3,'diagcomp',alpha));
    [u,flag,relres,iter,resvec] = pcg(matrix.A,matrix.f,1e-8,1000,L,L');

    time(7,i) = toc;
    tu = matrix.tu;
    uh = tu; %uh(matrix.mapper) = u;

    %% 5. Postprocess: Calculating Errors
    tic
    errND = max(abs(uh - tu)); % Error on nodes
    err.nd = errND; err.inf = 0; err.l2 = 0; err.h1 = 0;
    if computErr == 1
        disp(' ');  disp('Start computing error in Inf norm');
        eNorm = 'inf';
        errInfN = getErrVIFE3D(uh, pde, mesh, meshI, femI, fem, eNorm);
        eNorm = 'L2'; disp(['Start computing error in ',eNorm,'  norm']);
        errL2 = getErrVIFE3D(uh, pde, mesh, meshI, femI, fem, eNorm);
        eNorm = 'H1x'; disp(['Start computing error in ',eNorm,' norm']);
        errH1x = getErrVIFE3D(uh, pde, mesh, meshI, femI, fem, eNorm);
        eNorm = 'H1y'; disp(['Start computing error in ',eNorm,' norm']);
        errH1y = getErrVIFE3D(uh, pde, mesh, meshI, femI, fem, eNorm);
        eNorm = 'H1z'; disp(['Start computing error in ',eNorm,' norm']);
        errH1z = getErrVIFE3D(uh, pde, mesh, meshI, femI, fem, eNorm);

        disp(' ')
        err.inf = max([errND,errInfN]);
        err.l2 = errL2;
        err.h1 = sqrt(errH1x^2+errH1y^2+errH1z^2);
    end
    time(8,i) = toc;
    time(9,i) = 1e6*sum(time(2:8,i))/length(mesh.t);
    error(i,1) = errND; error(i,2) = err.inf;
    error(i,3) = err.l2; error(i,4) = err.h1;

    %% 6: Display Output
    disp(' ')
    disp('Errors')
    disp('Node        Inf norm    L2 norm     H1 norm')
    formatSpec = '%6.4e  %6.4e  %6.4e  %6.4e\n';
    fprintf(formatSpec, err.nd, err.inf, err.l2, err.h1)

    if i > 1
        format short
        rNd = log(err0.nd/err.nd)./log(h0/h);
        rInf = log(err0.inf/err.inf)./log(h0/h);
        rL2 = log(err0.l2/err.l2)./log(h0/h);
        rH1 = log(err0.h1/err.h1)./log(h0/h);
        disp(' ')
        disp('Convergence Rate')
        disp('Node        Inf norm    L2 norm     H1 norm')
        formatSpec = '%6.4f      %6.4f      %6.4f      %6.4f\n';
        fprintf(formatSpec, rNd, rInf, rL2, rH1)
        ratio(i,:) = [rNd,rInf,rL2,rH1];
    end
    err0 = err; h0 = h;

    disp(' '); disp('CPU Time')
    disp('   N     Mesh     FEM      MeshI    FemI    Matrix   Solve    Error    Time/1M cell')
    formatSpec = '%4i  %7.2f  %7.2f  %7.2f  %7.2f  %7.2f  %7.2f  %7.2f   %7.2f\n';
    fprintf(formatSpec, time(:,1:i))

    %% 6. Plot Solution and Error
    if showErr == 1
    end

end

repo: ifem


text:
function [s, d, S_rf, S_sf, S_k, S_l] = projHmgLinIntoPinHoleOnRob(Rf, Sf, Sk, l)

% PROJHMGLININTOPINHOLEONROB Project Hmg line into pinhole on robot.
%   [s, d] = PROJHMGLININTOPINHOLEONROB(Rf, Sf, Sk, l) projects the
%   inverse-depth line l into a pin/hole camera with intrinsic parameters
%   Sk mounted on a robot. The robot frame is Rf and the sensor frame in
%   the robot is Sf.
%
%   The results are a 2D segment S and a depths vector D.
%
%   [s, d, S_rf, S_sf, S_k, S_l] = (...) returns the Jacobians wrt all
%   input parameters.

%   Copyright 2009 Teresa Vidal.

if nargout <= 2

    sw = hmgLin2seg(l);
    s   = projSegLinIntoPinHoleOnRob(Rf, Sf, Sk, sw);
    d   = 1./l([6 9]);

else

    [sw, SW_l] = hmgLin2seg(l);
    [s, d, S_rf, S_sf, S_k, S_sw] = projSegLinIntoPinHoleOnRob(Rf, Sf, Sk, sw);
    S_l  = S_sw*SW_l;

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function hyperball_integrals_test ( )

%*****************************************************************************80
%
%% HYPERBALL_INTEGRALS_TEST tests the HYPERBALL_INTEGRALS library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HYPERBALL_INTEGRALS_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the HYPERBALL_INTEGRALS library.\n' );

  hyperball_integrals_test01 ( );
  hyperball_integrals_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HYPERBALL_INTEGRALS_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return;
end

author: johannesgerer

repo: jburkardt-m


text:
function hm=plotsurf(node,face,varargin)
%
% hm=plotsurf(node,face,opt)
%
% plot 3D surface meshes (2d manifold) or polylines (1d manifold)
%
% author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
%
% input:
%      node: node coordinates, dimension (nn,3); if node has a
%            4th column, it will be used to set the color at each node.
%      face: triangular surface face list; if face has a 4th column,
%            it will be used to separate the surface into
%            sub-surfaces and display them in different colors;
%            face can be a cell array, each element of the array represents
%            a polyhedral facet of the mesh, if an element is an array with
%            two array subelements, the first one is the node index, the
%            second one is a scalar as the group id of the facet.
%      opt:  additional options for the plotting, see plotmesh
%
% output:
%   hm: handle or handles (vector) to the plotted surfaces
%
% example:
%
%   h=plotsurf(node,face);
%   h=plotsurf(node,face,'facecolor','r');
%   h=plotsurf(node,edges,'linestyle','-','linewidth',2,'color','r');
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

rngstate = rand ('state');

if(nargin>=2)
  randseed=hex2dec('623F9A9E'); % "U+623F U+9A9E"
  if(isoctavemesh) randseed=randseed+3; end
  if(~isempty(getvarfrom({'caller','base'},'ISO2MESH_RANDSEED')))
        randseed=getvarfrom({'caller','base'},'ISO2MESH_RANDSEED');
  end
  rand('state',randseed);

  if(iscell(face))
    sc=sparse(10,3); % face colormap
    sc(1:10,:)=rand(3,10)';
    len=length(face);
    newsurf=cell(1);
    % reorganizing each labeled surface into a new cell
    for i=1:len
        fc=face{i};
        if(iscell(fc) && length(fc)>=2)
            if(fc{2}+1>10)
                sc(fc{2}+1,:)=rand(1,3);
            end
            if(fc{2}+1>length(newsurf))
                newsurf{fc{2}+1}={};
            end
            newsurf{fc{2}+1}{end+1}=fc{1};
        else % unlabeled facet is tagged by 0
            if(iscell(fc))
                newsurf{1}{end+1}=cell2mat(fc);
            else
                newsurf{1}{end+1}=fc;
            end
        end
    end
    hold on;
	h=[];
    newlen=length(newsurf);

    for i=1:newlen
        if(isempty(newsurf{i})); continue; end
        try
            subface=cell2mat(newsurf{i}')';
            if(size(subface,1)>1 && ismatrix(subface))
               subface=subface';
            end
            h=[h patch('Vertices',node,'Faces',subface,'facecolor',sc(i,:),varargin{:})];
        catch
            for j=1:length(newsurf{i})
                h=[h patch('Vertices',node,'Faces',newsurf{i}{j},'facecolor',sc(i,:),varargin{:})];
            end
        end
    end
  else
    if(size(face,2)==4)
        tag=face(:,4);
		types=unique(tag);
        hold on;
		h=[];
	    for i=1:length(types)
            if(size(node,2)==3)
                h=[h plotasurf(node,face(tag==types(i),1:3),'facecolor',rand(3,1),varargin{:})];
            else
                h=[h plotasurf(node,face(tag==types(i),1:3),varargin{:})];
            end
        end
    else
        h=plotasurf(node,face,varargin{:});
    end
  end
end
if(~isempty(h))
  axis equal;
  if(all(get(gca,'view')==[0 90]))
      view(3);
  end
end
if(~isempty(h) && nargout>=1)
  hm=h;
end

rand ('state',rngstate);

%-------------------------------------------------------------------------
function hh=plotasurf(node,face,varargin)
isoct=isoctavemesh;
if(size(face,2)<=2)
    h=plotedges(node,face,varargin{:});
else
  if(size(node,2)==4)
	if(isoct && ~exist('trisurf','file'))
	    h=trimesh(face(:,1:3),node(:,1),node(:,2),node(:,3),node(:,4),'edgecolor','k',varargin{:});
	else
	    h=trisurf(face(:,1:3),node(:,1),node(:,2),node(:,3),node(:,4),varargin{:});
	end
  else
	if(isoct && ~exist('trisurf','file'))
	    h=trimesh(face(:,1:3),node(:,1),node(:,2),node(:,3),'edgecolor','k',varargin{:});
	else
	    h=trisurf(face(:,1:3),node(:,1),node(:,2),node(:,3),varargin{:});
    end
  end
end
if(exist('h','var')) hh=h; end

author: fieldtrip

repo: fieldtrip


text:
function [paramsOut, okay] = MasterGUI(hObject, callbackdata, userData, EEG) %#ok<INUSL>
geometry = {1, [1, 1], [1, 1], [1, 1]};
geomvert = [1,1,1,1];
title = 'PREP pipeline control panel';
inputData = struct('signal', EEG, 'name', title, 'userData', userData);
closeOpenWindows(inputData.name);
uilist= {{'style', 'text', 'string', 'Override default parameters for processing step:'}...
    {'style', 'pushbutton', 'string', 'Boundary', ...
    'Callback', {@boundaryGUI, inputData}} ...
     {'style', 'pushbutton', 'string', 'Reference', ...
    'Callback', {@referenceGUI, inputData}} ...
    {'style', 'pushbutton', 'string', 'Detrend', ...
    'Callback', {@detrendGUI, inputData}} ...
     {'style', 'pushbutton', 'string', 'Report', ...
    'Callback', {@reportGUI, inputData}} ...
    {'style', 'pushbutton', 'string', 'Line noise', ...
    'Callback', {@lineNoiseGUI, inputData}}...
    {'style', 'pushbutton', 'string', 'Post process', ...
    'Callback', {@postProcessGUI, inputData}}};
[~, userData, buttonPressed, ~] = inputgui('geometry', geometry, 'geomvert', geomvert, ...
    'uilist', uilist, 'title', title, ...
    'helpcom', 'pophelp(''pop_prepPipeline'')');

okay = okayPressed(buttonPressed);
paramsOut = struct();
if ~isempty(userData)
    fNames = fieldnames(userData);
    for k = 1:length(fNames)
        nextStruct = userData.(fNames{k});
        nextNames = fieldnames(nextStruct);
        for j = 1:length(nextNames)
            paramsOut.(nextNames{j}) = nextStruct.(nextNames{j});
        end
    end
end

    function okay = okayPressed(buttonPressed)
        % Returns true if the okay button is pressed, false if otherwise
        okay = false;
        if ~isempty(buttonPressed)
            okay = true;
        end % okayPressed
    end

end % MasterGUI

author: VisLab

repo: EEG-Clean-Tools


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

%%Calculate the force/torque needed to produce the force torque needed
%%at the end effector:
%% we now that [T]=[J]'*[F]
%% where F=force\torque of the end effector (in this case only torques)
%% T=forces\torques of the joints (in this case only torques)

function T=staticForce(J6,F)

T=J6'*F;

end

author: Sable

repo: mcbench-benchmarks