text:
function [norm] = electroderealign(cfg);

% ELECTRODEREALIGN rotates and translates electrode positions to
% template electrode positions or towards the head surface. It can
% either perform a rigid body transformation, in which only the
% coordinate system is changed, or it can apply additional deformations
% to the input electrodes.
%
% Use as
%   [elec] = electroderealign(cfg)
%
% Three different methods for aligning the input electrodes are implemented:
% based on a warping method, based on the fiducials or interactive with a
% graphical user interface. Each of these approaches is described below.
%
% 1) You can apply a spatial deformation method (i.e. 'warp') that
% automatically minimizes the distance between the electrodes and the
% averaged standard. The warping methods use a non-linear search to
% optimize the error between input and template electrodes or the
% head surface.
%
% 2) You can apply a rigid body realignment based on three fiducial locations.
% Realigning using the fiducials only ensures that the fiducials (typically
% nose, left and right ear) are along the same axes in the input electrode
% set as in the template electrode set.
%
% 3) You can display the electrode positions together with the skin surface,
% and manually (using the graphical user interface) adjust the rotation,
% translation and scaling parameters, so that the two match.
%
% The configuration can contain the following options
%   cfg.method         = different methods for aligning the electrodes
%                        'rigidbody'       apply a rigid-body warp
%                        'globalrescale'   apply a rigid-body warp with global rescaling
%                        'traditional'     apply a rigid-body warp with individual axes rescaling
%                        'nonlin1'         apply a 1st order non-linear warp
%                        'nonlin2'         apply a 2nd order non-linear warp
%                        'nonlin3'         apply a 3rd order non-linear warp
%                        'nonlin4'         apply a 4th order non-linear warp
%                        'nonlin5'         apply a 5th order non-linear warp
%                        'realignfiducial' realign the fiducials
%                        'interactive'     manually using graphical user interface
%   cfg.channel        = Nx1 cell-array with selection of channels (default = 'all'),
%                        see CHANNELSELECTION for details
%   cfg.fiducial       = cell-array with the name of three fiducials used for
%                        realigning (default = {'nasion', 'lpa', 'rpa'})
%   cfg.casesensitive  = 'yes' or 'no', determines whether string comparisons
%                        between electrode labels are case sensitive (default = 'yes')
%   cfg.feedback       = 'yes' or 'no' (default = 'no')
%
% The electrode set that will be realigned is specified as
%   cfg.elecfile       = string with filename, or alternatively
%   cfg.elec           = structure with electrode definition
%
% If you want to align the electrodes to a single template electrode set
% or to multiple electrode sets (which will be averaged), you should
% specify the template electrode sets as
%   cfg.template       = single electrode set that serves as standard
% or
%   cfg.template{1..N} = list of electrode sets that are averaged into the standard
% The template electrode sets can be specified either as electrode
% structures (i.e. when they are already read in memory) or as electrode
% files.
%
% If you want to align the electrodes to the head surface as obtained from
% an anatomical MRI (using one of the warping methods), you should specify
% the head surface
%   cfg.headshape      = a filename containing headshape, a structure containing a
%                        single triangulated boundary, or a Nx3 matrix with surface
%                        points
%
% In case you only want to realign the fiducials, the template electrode
% set only has to contain the three fiducials, e.g.
%   cfg.template.pnt(1,:) = [110 0 0]  % location of the nose
%   cfg.template.pnt(2,:) = [0  90 0]  % left ear
%   cfg.template.pnt(3,:) = [0 -90 0]  % right ear
%   cfg.template.label    = {''nasion', 'lpa', 'rpa'}
%
% See also READ_FCDC_ELEC, VOLUMEREALIGN

% Copyright (C) 2005-2006, Robert Oostenveld
%
% $Log: electroderealign.m,v $
% Revision 1.1  2009/01/30 04:02:02  arno
% *** empty log message ***
%
% Revision 1.6  2007/08/06 09:20:14  roboos
% added support for bti_hs
%
% Revision 1.5  2007/07/26 08:00:09  roboos
% also deal with cfg.headshape if specified as surface, set of points or ctf_hs file.
% the construction of the tri is now done consistently for all headshapes if tri is missing
%
% Revision 1.4  2007/02/13 15:12:51  roboos
% removed cfg.plot3d option
%
% Revision 1.3  2006/12/12 11:28:33  roboos
% moved projecttri subfunction into seperate function
%
% Revision 1.2  2006/10/04 07:10:07  roboos
% updated documentation
%
% Revision 1.1  2006/09/13 07:20:06  roboos
% renamed electrodenormalize to electroderealign, added "deprecated"-warning to the old function
%
% Revision 1.10  2006/09/13 07:09:24  roboos
% Implemented support for cfg.method=interactive, using GUI for specifying and showing transformations. Sofar only for electrodes+headsurface.
%
% Revision 1.9  2006/09/12 15:26:06  roboos
% implemented support for aligning electrodes to the skin surface, extended and improved documentation
%
% Revision 1.8  2006/04/20 09:58:34  roboos
% updated documentation
%
% Revision 1.7  2006/04/19 15:42:53  roboos
% replaced call to warp_pnt with new function name warp_optim
%
% Revision 1.6  2006/03/14 08:16:00  roboos
% changed function call to warp3d into warp_apply (thanks to Arno)
%
% Revision 1.5  2005/05/17 17:50:37  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
%
% Revision 1.4  2005/03/21 15:49:43  roboos
% added cfg.casesensitive for string comparison of electrode labels
% added cfg.feedback and cfg.plot3d option for debugging
% changed output: now ALL electrodes of the input are rerurned, after applying the specified transformation
% fixed small bug in feedback regarding distarnce prior/after realignfiducials)
% added support for various warping strategies, a.o. traditional, rigidbody, nonlin1-5, etc.
%
% Revision 1.3  2005/03/16 09:18:56  roboos
% fixed bug in fprintf feedback, instead of giving mean squared distance it should give mean distance before and after normalization
%
% Revision 1.2  2005/01/18 12:04:39  roboos
% improved error handling of missing fiducials
% added other default fiducials
% changed debugging output
%
% Revision 1.1  2005/01/17 14:56:06  roboos
% new implementation
%

% set the defaults
if ~isfield(cfg, 'channel'),       cfg.channel = 'all';       end
if ~isfield(cfg, 'feedback'),      cfg.feedback = 'no';       end
if ~isfield(cfg, 'casesensitive'), cfg.casesensitive = 'yes'; end
if ~isfield(cfg, 'headshape'),     cfg.headshape = [];        end
if ~isfield(cfg, 'template'),      cfg.template = [];         end

% this is a common mistake which can be accepted
if strcmp(cfg.method, 'realignfiducials')
  cfg.method = 'realignfiducial';
end

if strcmp(cfg.method, 'warp')
  % rename the default warp to one of the method recognized by the warping toolbox
  cfg.method = 'traditional';
end

if strcmp(cfg.feedback, 'yes')
  % use the global fb field to tell the warping toolbox to print feedback
  global fb
  fb = 1;
else
  global fb
  fb = 0;
end

usetemplate  = isfield(cfg, 'template')  && ~isempty(cfg.template);
useheadshape = isfield(cfg, 'headshape') && ~isempty(cfg.headshape);

if usetemplate
  % get the template electrode definitions
  if ~iscell(cfg.template)
    cfg.template = {cfg.template};
  end
  Ntemplate = length(cfg.template);
  for i=1:Ntemplate
    if isstruct(cfg.template{i})
      template(i) = cfg.template{i};
    else
      template(i) = read_fcdc_elec(cfg.template{i});
    end
  end
elseif useheadshape
  % get the surface describing the head shape
  if isstruct(cfg.headshape) && isfield(cfg.headshape, 'pnt')
    % use the headshape surface specified in the configuration
    headshape = cfg.headshape;
  elseif isnumeric(cfg.headshape) && size(cfg.headshape,2)==3
    % use the headshape points specified in the configuration
    headshape.pnt = cfg.headshape;
  elseif ischar(cfg.headshape) && filetype(cfg.headshape, 'ctf_shape')
    % read the headshape from file
    headshape = read_ctf_shape(cfg.headshape);
  elseif ischar(cfg.headshape) && filetype(cfg.headshape, '4d_hs')
    % read the headshape from file
    headshape     = [];
    headshape.pnt = read_bti_hs(cfg.headshape);
  else
    error('cfg.headshape is not specified correctly')
  end
  if ~isfield(headshape, 'tri')
    % generate a closed triangulation from the surface points
    headshape.tri = projecttri(headshape.pnt);
  end
else
  error('you should either specify template electrode positions, template fiducials or a head shape');
end

% get the electrode definition that should be warped
if isfield(cfg, 'elec')
  elec = cfg.elec;
else
  elec = read_fcdc_elec(cfg.elecfile);
end

% remember the original electrode locations and labels
orig = elec;

% convert all labels to lower case for string comparisons
% this has to be done AFTER keeping the original labels and positions
if strcmp(cfg.casesensitive, 'no')
  for i=1:length(elec.label)
    elec.label{i} = lower(elec.label{i});
  end
  for j=1:length(template)
    for i=1:length(template(j).label)
      template(j).label{i} = lower(template(j).label{i});
    end
  end
end

if strcmp(cfg.feedback, 'yes')
  % create an empty figure, continued below...
  figure
  axis equal
  axis vis3d
  hold on
  xlabel('x')
  ylabel('y')
  zlabel('z')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if usetemplate && any(strcmp(cfg.method, {'rigidbody', 'globalrescale', 'traditional', 'nonlin1', 'nonlin2', 'nonlin3', 'nonlin4', 'nonlin5'}))
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % determine electrode selection and overlapping subset for warping
  cfg.channel = channelselection(cfg.channel, elec.label);
  for i=1:Ntemplate
    cfg.channel = channelselection(cfg.channel, template(i).label);
  end

  % make subselection of electrodes
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label = elec.label(datsel);
  elec.pnt   = elec.pnt(datsel,:);
  for i=1:Ntemplate
    [cfgsel, datsel] = match_str(cfg.channel, template(i).label);
    template(i).label = template(i).label(datsel);
    template(i).pnt   = template(i).pnt(datsel,:);
  end

  % compute the average of the template electrode positions
  all = [];
  for i=1:Ntemplate
    all = cat(3, all, template(i).pnt);
  end
  avg    = mean(all,3);
  stderr = std(all, [], 3);

  fprintf('warping electrodes to template... '); % the newline comes later
  [norm.pnt, norm.m] = warp_optim(elec.pnt, avg, cfg.method);
  norm.label = elec.label;

  dpre  = mean(sqrt(sum((avg - elec.pnt).^2, 2)));
  dpost = mean(sqrt(sum((avg - norm.pnt).^2, 2)));
  fprintf('mean distance prior to warping %f, after warping %f\n', dpre, dpost);

  if strcmp(cfg.feedback, 'yes')
    % plot all electrodes before warping
    my_plot3(elec.pnt, 'r.');
    my_plot3(elec.pnt(1,:), 'r*');
    my_plot3(elec.pnt(2,:), 'r*');
    my_plot3(elec.pnt(3,:), 'r*');
    my_text3(elec.pnt(1,:), elec.label{1}, 'color', 'r');
    my_text3(elec.pnt(2,:), elec.label{2}, 'color', 'r');
    my_text3(elec.pnt(3,:), elec.label{3}, 'color', 'r');

    % plot all electrodes after warping
    my_plot3(norm.pnt, 'm.');
    my_plot3(norm.pnt(1,:), 'm*');
    my_plot3(norm.pnt(2,:), 'm*');
    my_plot3(norm.pnt(3,:), 'm*');
    my_text3(norm.pnt(1,:), norm.label{1}, 'color', 'm');
    my_text3(norm.pnt(2,:), norm.label{2}, 'color', 'm');
    my_text3(norm.pnt(3,:), norm.label{3}, 'color', 'm');

    % plot the template electrode locations
    my_plot3(avg,      'b.');
    my_plot3(avg(1,:), 'b*');
    my_plot3(avg(2,:), 'b*');
    my_plot3(avg(3,:), 'b*');
    my_text3(avg(1,:), norm.label{1}, 'color', 'b');
    my_text3(avg(2,:), norm.label{2}, 'color', 'b');
    my_text3(avg(3,:), norm.label{3}, 'color', 'b');

    % plot lines connecting the input/warped electrode locations with the template locations
    my_line3(elec.pnt, avg, 'color', 'r');
    my_line3(norm.pnt, avg, 'color', 'm');
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif useheadshape && any(strcmp(cfg.method, {'rigidbody', 'globalrescale', 'traditional', 'nonlin1', 'nonlin2', 'nonlin3', 'nonlin4', 'nonlin5'}))
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % determine electrode selection and overlapping subset for warping
  cfg.channel = channelselection(cfg.channel, elec.label);

  % make subselection of electrodes
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label = elec.label(datsel);
  elec.pnt   = elec.pnt(datsel,:);

  fprintf('warping electrodes to head shape... '); % the newline comes later
  [norm.pnt, norm.m] = warp_optim(elec.pnt, headshape, cfg.method);
  norm.label = elec.label;

  dpre  = warp_error([],     elec.pnt, headshape, cfg.method);
  dpost = warp_error(norm.m, elec.pnt, headshape, cfg.method);
  fprintf('mean distance prior to warping %f, after warping %f\n', dpre, dpost);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(cfg.method, 'realignfiducial')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % try to determine the fiducials automatically if not specified
  option1 = {'nasion' 'left' 'right'};
  option2 = {'nasion' 'lpa' 'rpa'};
  option3 = {'nz' 'lpa' 'rpa'};
  if ~isfield(cfg, 'fiducial')
    if length(match_str(elec.label, option1))==3
      cfg.fiducial = option1;
    elseif length(match_str(elec.label, option2))==3
      cfg.fiducial = option2;
    elseif length(match_str(elec.label, option3))==3
      cfg.fiducial = option3;
    else
      error('could not determine three fiducials, please specify cfg.fiducial')
    end
  end
  fprintf('using fiducials {''%s'', ''%s'', ''%s''}\n', cfg.fiducial{1}, cfg.fiducial{2}, cfg.fiducial{3});

  % determine electrode selection
  cfg.channel = channelselection(cfg.channel, elec.label);
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label = elec.label(datsel);
  elec.pnt   = elec.pnt(datsel,:);

  if length(cfg.fiducial)~=3
    error('you must specify three fiducials');
  end

  % do case-insensitive search for fiducial locations
  nas_indx = match_str(lower(elec.label), lower(cfg.fiducial{1}));
  lpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{2}));
  rpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{3}));
  if length(nas_indx)~=1 || length(lpa_indx)~=1 || length(rpa_indx)~=1
    error('not all fiducials were found in the electrode set');
  end
  elec_nas = elec.pnt(nas_indx,:);
  elec_lpa = elec.pnt(lpa_indx,:);
  elec_rpa = elec.pnt(rpa_indx,:);

  % find the matching fiducials in the template and average them
  templ_nas = [];
  templ_lpa = [];
  templ_rpa = [];
  for i=1:Ntemplate
    nas_indx = match_str(lower(template(i).label), lower(cfg.fiducial{1}));
    lpa_indx = match_str(lower(template(i).label), lower(cfg.fiducial{2}));
    rpa_indx = match_str(lower(template(i).label), lower(cfg.fiducial{3}));
    if length(nas_indx)~=1 || length(lpa_indx)~=1 || length(rpa_indx)~=1
      error(sprintf('not all fiducials were found in template %d', i));
    end
    templ_nas(end+1,:) = template(i).pnt(nas_indx,:);
    templ_lpa(end+1,:) = template(i).pnt(lpa_indx,:);
    templ_rpa(end+1,:) = template(i).pnt(rpa_indx,:);
  end
  templ_nas = mean(templ_nas,1);
  templ_lpa = mean(templ_lpa,1);
  templ_rpa = mean(templ_rpa,1);

  % realign both to a common coordinate system
  elec2common  = headcoordinates(elec_nas, elec_lpa, elec_rpa);
  templ2common = headcoordinates(templ_nas, templ_lpa, templ_rpa);

  % compute the combined transform and realign the electrodes to the template
  norm       = [];
  norm.m     = elec2common * inv(templ2common);
  norm.pnt   = warp_apply(norm.m, elec.pnt, 'homogeneous');
  norm.label = elec.label;

  nas_indx = match_str(lower(elec.label), lower(cfg.fiducial{1}));
  lpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{2}));
  rpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{3}));
  dpre  = mean(sqrt(sum((elec.pnt([nas_indx lpa_indx rpa_indx],:) - [templ_nas; templ_lpa; templ_rpa]).^2, 2)));
  nas_indx = match_str(lower(norm.label), lower(cfg.fiducial{1}));
  lpa_indx = match_str(lower(norm.label), lower(cfg.fiducial{2}));
  rpa_indx = match_str(lower(norm.label), lower(cfg.fiducial{3}));
  dpost = mean(sqrt(sum((norm.pnt([nas_indx lpa_indx rpa_indx],:) - [templ_nas; templ_lpa; templ_rpa]).^2, 2)));
  fprintf('mean distance between fiducials prior to realignment %f, after realignment %f\n', dpre, dpost);

  if strcmp(cfg.feedback, 'yes')
    % plot the first three electrodes before transformation
    my_plot3(elec.pnt(1,:), 'r*');
    my_plot3(elec.pnt(2,:), 'r*');
    my_plot3(elec.pnt(3,:), 'r*');
    my_text3(elec.pnt(1,:), elec.label{1}, 'color', 'r');
    my_text3(elec.pnt(2,:), elec.label{2}, 'color', 'r');
    my_text3(elec.pnt(3,:), elec.label{3}, 'color', 'r');

    % plot the template fiducials
    my_plot3(templ_nas, 'b*');
    my_plot3(templ_lpa, 'b*');
    my_plot3(templ_rpa, 'b*');
    my_text3(templ_nas, ' nas', 'color', 'b');
    my_text3(templ_lpa, ' lpa', 'color', 'b');
    my_text3(templ_rpa, ' rpa', 'color', 'b');

    % plot all electrodes after transformation
    my_plot3(norm.pnt, 'm.');
    my_plot3(norm.pnt(1,:), 'm*');
    my_plot3(norm.pnt(2,:), 'm*');
    my_plot3(norm.pnt(3,:), 'm*');
    my_text3(norm.pnt(1,:), norm.label{1}, 'color', 'm');
    my_text3(norm.pnt(2,:), norm.label{2}, 'color', 'm');
    my_text3(norm.pnt(3,:), norm.label{3}, 'color', 'm');
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(cfg.method, 'interactive')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % open a figure
  fig = figure;
  % add the data to the figure
  set(fig, 'CloseRequestFcn', @cb_close);
  setappdata(fig, 'elec', elec);
  setappdata(fig, 'transform', eye(4));
  if useheadshape
    setappdata(fig, 'surf', headshape);
  end
  if usetemplate
    % FIXME interactive realigning to template electrodes is not yet supported
    % this requires a consistent handling of channel selection etc.
    setappdata(fig, 'template', template);
  end
  % add the GUI elements
  cb_creategui(gca);
  cb_redraw(gca);
  rotate3d on
  waitfor(fig);
  % get the data from the figure that was left behind as global variable
  global norm
  tmp = norm;
  clear global norm
  norm = tmp;
  clear tmp

else
  error('unknown method');
end

% apply the spatial transformation to all electrodes, and replace the
% electrode labels by their case-sensitive original values
if any(strcmp(cfg.method, {'rigidbody', 'globalrescale', 'traditional', 'nonlin1', 'nonlin2', 'nonlin3', 'nonlin4', 'nonlin5'}))
  norm.pnt   = warp_apply(norm.m, orig.pnt, cfg.method);
else
  norm.pnt   = warp_apply(norm.m, orig.pnt, 'homogenous');
end
norm.label = orig.label;

% add version information to the configuration
try
  % get the full name of the function
  cfg.version.name = mfilename('fullpath');
catch
  % required for compatibility with Matlab versions prior to release 13 (6.5)
  [st, i] = dbstack;
  cfg.version.name = st(i);
end
cfg.version.id = '$Id: electroderealign.m,v 1.1 2009/01/30 04:02:02 arno Exp $';

% remember the configuration
norm.cfg = cfg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some simple SUBFUNCTIONs that facilitate 3D plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = my_plot3(xyz, varargin)
h = plot3(xyz(:,1), xyz(:,2), xyz(:,3), varargin{:});
function h = my_text3(xyz, varargin)
h = text(xyz(:,1), xyz(:,2), xyz(:,3), varargin{:});
function my_line3(xyzB, xyzE, varargin)
for i=1:size(xyzB,1)
  line([xyzB(i,1) xyzE(i,1)], [xyzB(i,2) xyzE(i,2)], [xyzB(i,3) xyzE(i,3)], varargin{:})
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION to layout a moderately complex graphical user interface
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = layoutgui(fig, geometry, position, style, string, value, tag, callback);
horipos  = geometry(1); % lower left corner of the GUI part in the figure
vertpos  = geometry(2); % lower left corner of the GUI part in the figure
width    = geometry(3); % width  of the GUI part in the figure
height   = geometry(4); % height of the GUI part in the figure
horidist = 0.05;
vertdist = 0.05;
options  = {'units', 'normalized', 'HorizontalAlignment', 'center'}; %  'VerticalAlignment', 'middle'
Nrow     = size(position,1);
h        = cell(Nrow,1);
for i=1:Nrow
  if isempty(position{i})
    continue;
  end
  position{i} = position{i} ./ sum(position{i});
  Ncol = size(position{i},2);
  ybeg = (Nrow-i  )/Nrow + vertdist/2;
  yend = (Nrow-i+1)/Nrow - vertdist/2;
  for j=1:Ncol
    xbeg    = sum(position{i}(1:(j-1))) + horidist/2;
    xend    = sum(position{i}(1:(j  ))) - horidist/2;
    pos(1) = xbeg*width  + horipos;
    pos(2) = ybeg*height + vertpos;
    pos(3) = (xend-xbeg)*width;
    pos(4) = (yend-ybeg)*height;
    h{i}{j} = uicontrol(fig, ...
      options{:}, ...
      'position', pos, ...
      'style',    style{i}{j}, ...
      'string',   string{i}{j}, ...
      'tag',      tag{i}{j}, ...
      'value',    value{i}{j}, ...
      'callback', callback{i}{j} ...
      );
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_creategui(hObject, eventdata, handles);
% define the position of each GUI element
position = {
  [2 1 1 1]
  [2 1 1 1]
  [2 1 1 1]
  [1]
  [1]
  [1]
  [1]
  [1 1]
  };

% define the style of each GUI element
style = {
  {'text' 'edit' 'edit' 'edit'}
  {'text' 'edit' 'edit' 'edit'}
  {'text' 'edit' 'edit' 'edit'}
  {'pushbutton'}
  {'pushbutton'}
  {'toggle'}
  {'toggle'}
  {'text' 'edit'}
  };

% define the descriptive string of each GUI element
string = {
  {'rotate'    0 0 0}
  {'translate' 0 0 0}
  {'scale'     1 1 1}
  {'redisplay'}
  {'apply'}
  {'toggle grid'}
  {'toggle axes'}
  {'alpha' 0.7}
  };

% define the value of each GUI element
value = {
  {[] [] [] []}
  {[] [] [] []}
  {[] [] [] []}
  {[]}
  {[]}
  {0}
  {0}
  {[] []}
  };

% define a tag for each GUI element
tag = {
  {'' 'rx' 'ry' 'rz'}
  {'' 'tx' 'ty' 'tz'}
  {'' 'sx' 'sy' 'sz'}
  {''}
  {''}
  {'toggle grid'}
  {'toggle axes'}
  {'' 'alpha'}
  };

% define the callback function of each GUI element
callback = {
  {[] @cb_redraw @cb_redraw @cb_redraw}
  {[] @cb_redraw @cb_redraw @cb_redraw}
  {[] @cb_redraw @cb_redraw @cb_redraw}
  {@cb_redraw}
  {@cb_apply}
  {@cb_redraw}
  {@cb_redraw}
  {[] @cb_redraw}
  };

fig = get(hObject, 'parent');
layoutgui(fig, [0.7 0.05 0.25 0.50], position, style, string, value, tag, callback);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_redraw(hObject, eventdata, handles);
fig = get(hObject, 'parent');
surf = getappdata(fig, 'surf');
elec = getappdata(fig, 'elec');
template = getappdata(fig, 'template');
% get the transformation details
rx = str2num(get(findobj(fig, 'tag', 'rx'), 'string'));
ry = str2num(get(findobj(fig, 'tag', 'ry'), 'string'));
rz = str2num(get(findobj(fig, 'tag', 'rz'), 'string'));
tx = str2num(get(findobj(fig, 'tag', 'tx'), 'string'));
ty = str2num(get(findobj(fig, 'tag', 'ty'), 'string'));
tz = str2num(get(findobj(fig, 'tag', 'tz'), 'string'));
sx = str2num(get(findobj(fig, 'tag', 'sx'), 'string'));
sy = str2num(get(findobj(fig, 'tag', 'sy'), 'string'));
sz = str2num(get(findobj(fig, 'tag', 'sz'), 'string'));
R = rotate   ([rx ry rz]);
T = translate([tx ty tz]);
S = scale    ([sx sy sz]);
H = S * T * R;
elec.pnt = warp_apply(H, elec.pnt);
axis vis3d; cla
xlabel('x')
ylabel('y')
zlabel('z')
if ~isempty(surf)
  triplot(surf.pnt, surf.tri,  [], 'faces_skin');
  alpha(str2num(get(findobj(fig, 'tag', 'alpha'), 'string')));
end
if ~isempty(template)
  triplot(template.pnt, [], [], 'nodes_blue')
end
triplot(elec.pnt, [], [], 'nodes');
if isfield(elec, 'line')
  triplot(elec.pnt, elec.line, [], 'edges');
end
if get(findobj(fig, 'tag', 'toggle axes'), 'value')
  axis on
else
  axis off
end
if get(findobj(fig, 'tag', 'toggle grid'), 'value')
  grid on
else
  grid off
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_apply(hObject, eventdata, handles);
fig = get(hObject, 'parent');
elec      = getappdata(fig, 'elec');
transform = getappdata(fig, 'transform');
% get the transformation details
rx = str2num(get(findobj(fig, 'tag', 'rx'), 'string'));
ry = str2num(get(findobj(fig, 'tag', 'ry'), 'string'));
rz = str2num(get(findobj(fig, 'tag', 'rz'), 'string'));
tx = str2num(get(findobj(fig, 'tag', 'tx'), 'string'));
ty = str2num(get(findobj(fig, 'tag', 'ty'), 'string'));
tz = str2num(get(findobj(fig, 'tag', 'tz'), 'string'));
sx = str2num(get(findobj(fig, 'tag', 'sx'), 'string'));
sy = str2num(get(findobj(fig, 'tag', 'sy'), 'string'));
sz = str2num(get(findobj(fig, 'tag', 'sz'), 'string'));
R = rotate   ([rx ry rz]);
T = translate([tx ty tz]);
S = scale    ([sx sy sz]);
H = S * T * R;
elec.pnt = warp_apply(H, elec.pnt);
transform = H * transform;
set(findobj(fig, 'tag', 'rx'), 'string', 0);
set(findobj(fig, 'tag', 'ry'), 'string', 0);
set(findobj(fig, 'tag', 'rz'), 'string', 0);
set(findobj(fig, 'tag', 'tx'), 'string', 0);
set(findobj(fig, 'tag', 'ty'), 'string', 0);
set(findobj(fig, 'tag', 'tz'), 'string', 0);
set(findobj(fig, 'tag', 'sx'), 'string', 1);
set(findobj(fig, 'tag', 'sy'), 'string', 1);
set(findobj(fig, 'tag', 'sz'), 'string', 1);
setappdata(fig, 'elec', elec);
setappdata(fig, 'transform', transform);
cb_redraw(hObject);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_close(hObject, eventdata, handles);
% make the current transformation permanent and subsequently allow deleting the figure
cb_apply(gca);
% get the updated electrode from the figure
fig    = hObject;
% hmmm, this is ugly
global norm
norm   = getappdata(fig, 'elec');
norm.m = getappdata(fig, 'transform');
set(fig, 'CloseRequestFcn', @delete);
delete(fig);

author: PatternRecognition

repo: OpenBMI


text:
function state = psoplotswarm(options,state,flag,ijk)
% Plots the positions of particle swarm.
%
% ijk is an additional parameter, in the form of a 1x3 vector. It lists the
% dimension of the problem to be plotted, for a multidimensional problem.
%
% Updated this to conform with MATLAB GA toolbox specifications.

% We'll only plot the first two dimensions by default. User can
% set different dimensions to be plotted using ijk.

if nargin < 4 && size(state.Population,2) > 1 % Go to defaults
    ijk = [1,2] ; % k doesn't exist
elseif size(state.Population,2) == 1 % For 1-D problems
    ijk = 1 ;
else % Robust input check
    ijk = reshape(ijk,1,[]) ;
end

% Input checking
if size(ijk,2) > 3
    warning('PSO:Plotting:tooManyDimensions',...
        'Unable to plot %g dimensions. Plotting first 3 only',...
        length(ijk))
end % if length

if strcmpi(flag(1:4),'init') % Initialize
    delete(findobj(gca,'-regexp','Tag','*Locations'))
    initLoc = line(state.Population(:,ijk(1)),state.Population(:,ijk(2)),...
        'Color',0.75*ones(1,3),...
        'Marker','.',...
        'LineStyle','none',...
        'Tag','Initial Locations') ;

    % Set reasonable axes limits
    % ---------------------------------------------------------------------
    xlim([options.PopInitRange(1,ijk(1)) options.PopInitRange(2,ijk(1))])
    if size(ijk,2) > 1
        ylim([options.PopInitRange(1,ijk(2)) ...
            options.PopInitRange(2,ijk(2))])
        if size(ijk,2) > 2
            zlim([options.PopInitRange(1,ijk(3)) ...
                options.PopInitRange(2,ijk(3))])
        end % if size
    end % if size
    % ---------------------------------------------------------------------

    title('Swarm positions')
    set(gca,'Tag','Swarm Plot','NextPlot','add')

    % Initialize plots
    % ---------------------------------------------------------------------
    if size(ijk,2) == 1 %  One dimensional
        currentLoc = line(state.Population(:,ijk(1)),...
            zeros(size(state.Population,1),1)) ;
    elseif size(ijk,2) == 2 % Two dimensional
        currentLoc = line(state.Population(:,ijk(1)),...
            state.Population(:,ijk(2))) ;
    elseif size(ijk,2) == 3 % Three dimensional
        currentLoc = line(state.Population(:,ijk(1)),...
            state.Population(:,ijk(2)),...
            state.Population(:,ijk(3))) ;
    end % if size
    set(currentLoc,...
            'LineStyle','none',...
            'Marker','.',...
            'Color','blue',...
            'Tag','Swarm Locations') ;
    % ---------------------------------------------------------------------
elseif strcmpi(flag(1:4),'iter') % Iterate
    currentLoc = findobj(gca,'Tag','Swarm Locations','Type','line') ;
    if size(ijk,2) == 1 %  One dimensional
        set(currentLoc,'XData',state.Population(:,ijk(1)))
    elseif size(ijk,2) == 2 % Two dimensional
        set(currentLoc,...
            'XData',state.Population(:,ijk(1)),...
            'YData',state.Population(:,ijk(2)))
    elseif size(ijk,2) == 3 % Three dimensional
        set(currentLoc,...
            'XData',state.Population(:,ijk(1)),...
            'YData',state.Population(:,ijk(2)),...
            'ZData',state.Population(:,ijk(3)))
    end
end

if strcmpi(flag(1:4),'init')
    legend([initLoc,currentLoc],{'Initial Positions','Current Positions'})
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function output = callsdpagmp(interfacedata)

% CALLSDPAGMP.m Call SDPA-GMP from YALMIP
% ----------------------------------------------------------------------- %
%        Author:    Giovanni Fantuzzi
%                   Department of Aeronautics
%                   Imperial College London
%       Created:    23/08/2016
%
%     Copyright (C) 2016  Giovanni Fantuzzi
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ----------------------------------------------------------------------- %

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

% Bounded variables converted to constraints
if ~isempty(ub)
    [F_struc,K] = addStructureBounds(F_struc,K,ub,lb);
end

% Convert from internal (sedumi) format
[mDIM,nBLOCK,bLOCKsTRUCT,c,F] = sedumi2sdpa(F_struc,c,K);

if options.verbose==0
    options.sdpa_gmp.print = 'no';
else
    options.sdpa_gmp.print = 'display';
end

if options.savedebug
    ops = options.sdpa_gmp;
    save sdpa_gmpdebug mDIM nBLOCK bLOCKsTRUCT c F ops
end

if options.showprogress
    showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALL SDPA-GMP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
solvertime = tic;
[objVal,x,X,Y,INFO] = sdpagmp(mDIM,nBLOCK,bLOCKsTRUCT,c,F,[],[],[],options.sdpa_gmp);
solvertime = toc(solvertime);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% From here onwards, like in YALMIP native callsdpa
% Create variables in YALMIP internal format
Primal = x;

Dual = [];
for i = 1:length(Y)
    Dual = [Dual;Y{i}(:)];
end

Slack = [];
if options.saveduals
    for i = 1:length(X)
        Slack = [Slack;X{i}(:)];
    end
end

switch (INFO.phasevalue)
    case 'pdOPT'
        problem = 0;
    case {'noINFO','pFEAS','dFEAS'}
        problem = 3;
    case {'pdFEAS'}
        problem = 4;
    case 'pFEAS_dINF'
        problem = 2;
    case 'pINF_dFEAS'
        problem = 1;
    case 'pUNBD'
        problem = 2;
    case 'dUNBD'
        problem = 1;
    case 'pdINF'
        problem = 12;
    otherwise
        problem = -1;
end
infostr = yalmiperror(problem,interfacedata.solver.tag);

if options.savesolveroutput
    solveroutput.objVal = objVal;
    solveroutput.x = x;
    solveroutput.X = X;
    solveroutput.Y = Y;
    solveroutput.INFO = INFO;
else
    solveroutput = [];
end

if options.savesolverinput
    solverinput.mDIM = mDIM;
    solverinput.nBLOCK=nBLOCK;
    solverinput.bLOCKsTRUCT=bLOCKsTRUCT;
    solverinput.c=c;
    solverinput.F=F;
else
    solverinput = [];
end

% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [c, ceq] = ConstraintFCN_models(u,uold,x,N,LBo,UBo,LBdu,UBdu,p,select_model)
%% Constraint function of nonlinear MPC for F8 system
%
% Inputs:
%   u:      optimization variable, from time k to time k+N-1
%   x:      current state at time k
%   Ts:     controller sample time
%   N:      prediction horizon length
%   uold:   latest applied control input
%   LBo:    Lower bound of output x
%   UBo:    Upper bound of output x
%   LBdu:   Lower bound for input difference uk - uk-1
%   UBdu:   Upper bound for input difference uk - uk-1
%   p:      Parameters for model
%   select_model: Selects model future-state prediction
%
% Output:
%   c:      inequality constraints applied across prediction horizon
%   ceq:    equality constraints  (empty)
%
Nvar = length(x);
%% Nonlinear MPC design parameters
% range of angle of attack
zMin = LBo(1);
zMax = UBo(1);

%% Integrate system
if strcmp(select_model,'DelayDMDc')
     [xk,~] = lsim(p.sys,[p.udelay(1:N);u'],[0:N-1].*p.dt,[p.xdelay(:,1); x]-[p.xmean; p.xmean]);
    xk = xk(:,Nvar+1:2*Nvar);
    xk = xk + repmat(p.xmean',[N 1]); xk = xk';
elseif strcmp(select_model,'DMDc')
    [xk,~] = lsim(p.sys,[u',0],[0:N].*p.dt,x-p.xmean);
    xk = xk(2:end,:) + repmat(p.xmean',[N 1]); xk = xk';
elseif strcmp(select_model,'SINDYc')
    Ns = size(x,1);
    xk = zeros(Ns,N+1); xk(:,1) = x;
    for ct=1:N
        % Obtain plant state at next prediction step.
        xk(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xk(:,ct),u(ct),p.dt,1,[],p);
    end
    xk = xk(:,2:N+1);
elseif strcmp(select_model,'NARX')
    Hu = [u',0];
    Hx = zeros(Nvar,length(Hu)); Hx(:,1) = x;
    [Us,Ui,Si] = preparets(p.net,con2seq(Hu),{},con2seq(Hx));
    xk = p.net(Us,Ui,Si);
    xk = cell2mat(xk);
end

%% Inequality constraints calculation
c = zeros(2*N,1);

% Apply N population size constraints across prediction horizon, from time
% k+1 to k+N
duk = u(1)-uold;
for ct=1:N
    c(2*ct-1) = -duk+LBdu;
    c(2*ct) = duk-UBdu;
    if ct<N
        duk = u(ct+1)-u(ct);
    end
end

c1 = zeros(2*N,1);
c2 = zeros(2*N,1);
c3 = zeros(2*N,1);
for ct=1:N
    c1(2*ct-1) = -xk(1,ct)+zMin;
    c2(2*ct-1) = -xk(2,ct)+LBo(2);
    c3(2*ct-1) = -xk(3,ct)+LBo(3);
    c1(2*ct) = xk(1,ct)-zMax;
    c2(2*ct) = xk(2,ct)-UBo(2);
    c3(2*ct) = xk(3,ct)-UBo(3);
end

c = [c;c1;c2;c3];

%% No equality constraints
ceq = [];

author: eurika-kaiser

repo: SINDY-MPC


text:
function output = precision(tp,fp)
%
%
% Precision = tp/(tp+fp)
% (see page 268 of Manning and Schutze)
%
% Inputs
%    tp: True Positive
%    fp: False Positive
% Outputs
%    precision measure

output = tp/(tp+fp);

end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function img = discreteCylinder(varargin)
%DISCRETECYLINDER Discretize a 3D cylinder
%
%   IMG = discreteCylinder(LX, LY, LZ, P1, P2, RADIUS)
%   LX, LY and LZ are row vectors specifying position of vertex centers
%   along each coordinate.
%   P1 is the starting point of the cylinder, given as a 1-by-3 row vector
%   P2 is the ending point of the cylinder, given as a 1-by-3 row vector
%   RADIUS is the cylinder radius.
%
%   IMG = discreteCylinder(LX, LY, LZ, CYLINDER)
%   send parameters in a row vector, where CYLINDER is a 1-by-7 row vector
%   containing coordinates of start point, of end point, and radius:
%   CYLINDER = [X1 Y1 Z1 X2 Y2 Z2 R];
%
%   Example
%     % Compute the union of three mutually orthogonal cylinders
%     p0 = [30 30 30];
%     p1 = [90 30 30];
%     p2 = [30 90 30];
%     p3 = [30 30 90];
%     cyl1 = discreteCylinder(1:100, 1:100, 1:100, [p0 p1 25]);
%     cyl2 = discreteCylinder(1:100, 1:100, 1:100, [p0 p2 25]);
%     cyl3 = discreteCylinder(1:100, 1:100, 1:100, [p0 p3 25]);
%     cylUnion = cyl1 | cyl2 | cyl3;
%     [f v] = isosurface(cylUnion, .5);
%     figure;
%     drawMesh(v, f, 'linestyle', 'none', 'facecolor', 'r');
%     l = light; view([120 20]);
%
%   See also
%   imShapes, discreteBall, discreteCapsule3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-10-21
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY
%   2010-10-21 create from discreteCube

%% Process input arguments

% compute coordinate of image voxels
[lx, ly, lz, varargin] = parseGridArgs3d(varargin{:});
[x, y, z] = meshgrid(lx, ly, lz);

% process input parameters
if length(varargin) == 1
    % input is a 1-by-7 row vector
    var = varargin{1};
    if length(var) ~= 7
        error('Should specify a row vector with 7 inputs');
    end

    % extract first and last point coordinates
    p1 = var(1:3);
    p2 = var(4:6);
    radius = var(7);

elseif length(varargin) == 3
    % inputs are P1, P2 and R
    p1 = varargin{1};
    p2 = varargin{2};
    radius = varargin{3};

else
    error('wrong number of arguments: should be 1 or 3');
end

%% Transform voxel coordinates

% compute cylinder direction angle
dirVect = p2-p1;

[theta, phi, height] = cart2sph2(dirVect);

% compute coordinate of image voxels in cylinder reference system
% (cylinder pointing upwards)
trans = composeTransforms3d(...
    createTranslation3d(-p1),...
    createRotationOz(-phi),...
    createRotationOy(-theta),...
    createScaling3d(1./[radius radius height]));
[x, y, z] = transformPoint3d(x, y, z, trans);

%% Compute final image

% create image: simple threshold over 3 dimensions, and test z axis
img = ((x.*x + y.*y) < 1) & (z>0) & (z<1);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function data = LoadBinary_Down_ss(filename,varargin)

% LoadBinary - Load data from a multiplexed binary file.
%
%  Reading a subset of the data can be done in two different manners: either
%  by specifying start time and duration (more intuitive), or by indicating
%  the position and size of the subset in terms of number of samples per
%  channel (more accurate).
%
%  USAGE
%
%    data = LoadBinary(filename,<options>)
%
%    filename       file to read
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'frequency'   sampling rate (in Hz, default = 20kHz)
%     'start'       position to start reading (in s, default = 0)
%     'duration'    duration to read (in s, default = Inf)
%     'offset'      position to start reading (in samples per channel,
%                   default = 0)
%     'samples'     number of samples (per channel) to read (default = Inf)
%     'nChannels'   number of data channels in the file (default = 1)
%     'channels'    channels to read (default = all)
%     'precision'   sample precision (default = 'int16')
%     'skip'        number of bytes to skip after each value is read
%                   (default = 0)
%     'downsample'  factor by which to downample by (default = 1)
%    =========================================================================

% Copyright (C) 2004-2011 by Michaël Zugaro
%Modified by DLevenstein 2016 to include downsampling
%
% 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.

% Default values
nChannels = 1;
precision = 'int16';
skip = 0;
frequency = 20000;
channels = [];
start = 0;
duration = Inf;
offset = 0;
nSamplesPerChannel = Inf;
time = false;
samples = false;
downsamplefactor = 1;

if nargin < 1 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help LoadBinary'' for details).');
end

% Parse options
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+3) ' is not a property (type ''help LoadBinary'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'frequency',
			frequency = varargin{i+1};
			if ~isdscalar_ss(frequency,'>0'),
				error('Incorrect value for property ''frequency'' (type ''help LoadBinary'' for details).');
			end
		case 'start',
			start = varargin{i+1};
			if ~isdscalar_ss(start),
				error('Incorrect value for property ''start'' (type ''help LoadBinary'' for details).');
			end
			if start < 0, start = 0; end
			time = true;
		case 'duration',
			duration = varargin{i+1};
			if ~isdscalar_ss(duration,'>=0'),
				error('Incorrect value for property ''duration'' (type ''help LoadBinary'' for details).');
			end
			time = true;
		case 'offset',
			offset = varargin{i+1};
			if ~isiscalar_ss(offset),
				error('Incorrect value for property ''offset'' (type ''help LoadBinary'' for details).');
			end
			if offset < 0, offset = 0; end
			samples = true;
		case 'samples',
			nSamplesPerChannel = varargin{i+1};
			if ~isdscalar_ss(nSamplesPerChannel,'>=0'),
				error('Incorrect value for property ''samples'' (type ''help LoadBinary'' for details).');
			end
			samples = true;
		case 'nchannels',
			nChannels = varargin{i+1};
			if ~isiscalar_ss(nChannels,'>0'),
				error('Incorrect value for property ''nChannels'' (type ''help LoadBinary'' for details).');
			end
		case 'channels',
			channels = varargin{i+1};
			if ~isivector_ss(channels,'>0'),
				error('Incorrect value for property ''channels'' (type ''help LoadBinary'' for details).');
			end
		case 'precision',
			precision = varargin{i+1};
			if ~isstring_ss(precision),
				error('Incorrect value for property ''precision'' (type ''help LoadBinary'' for details).');
			end
		case 'skip',
			skip = varargin{i+1};
			if ~isiscalar_ss(skip,'>=0'),
				error('Incorrect value for property ''skip'' (type ''help LoadBinary'' for details).');
            end
		case 'downsample',
			downsamplefactor = varargin{i+1};
			if ~isiscalar_ss(downsamplefactor,'>=0'),
				error('Incorrect value for property ''downsample'' (type ''help LoadBinary'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help LoadBinary'' for details).']);
	end
end

% Either start+duration, or offset+size
if time && samples,
	error(['Data subset can be specified either in time or in samples, but not both (type ''help LoadBinary'' for details).']);
end

% By default, load all channels
if isempty(channels),
	channels = 1:nChannels;
end

% Check consistency between channel IDs and number of channels
if any(channels>nChannels),
	error('Cannot load specified channels (listed channel IDs inconsistent with total number of channels).');
end

% Open file
if ~exist(filename),
	error(['File ''' filename ''' not found.']);
end
f = fopen(filename,'r');
if f == -1,
	error(['Cannot read ' filename ' (insufficient access rights?).']);
end

% Size of one data point (in bytes)
sampleSize = 0;
switch precision,
	case {'uchar','unsigned char','schar','signed char','int8','integer*1','uint8','integer*1'},
		sampleSize = 1;
	case {'int16','integer*2','uint16','integer*2'},
		sampleSize = 2;
	case {'int32','integer*4','uint32','integer*4','single','real*4','float32','real*4'},
		sampleSize = 4;
	case {'int64','integer*8','uint64','integer*8','double','real*8','float64','real*8'},
		sampleSize = 8;
end

% Position and number of samples (per channel) of the data subset
if time,
	dataOffset = floor(start*frequency)*nChannels*sampleSize;
	nSamplesPerChannel = floor((duration*frequency));
else
	dataOffset = offset*nChannels*sampleSize;
end

% Position file index for reading
status = fseek(f,dataOffset,'bof');
if status ~= 0,
	fclose(f);
	error('Could not start reading (possible reasons include trying to read past the end of the file).');
end

% Determine total number of samples in file
fileStart = ftell(f);
status = fseek(f,0,'eof');
if status ~= 0,
	fclose(f);
	error('Error reading the data file (possible reasons include trying to read past the end of the file).');
end
fileStop = ftell(f);
% (floor in case all channels do not have the same number of samples)
maxNSamplesPerChannel = floor(((fileStop-fileStart)/nChannels/sampleSize));
frewind(f);
status = fseek(f,dataOffset,'bof');
if status ~= 0,
	fclose(f);
	error('Could not start reading (possible reasons include trying to read past the end of the file).');
end

if isinf(nSamplesPerChannel) || nSamplesPerChannel > maxNSamplesPerChannel,
	nSamplesPerChannel = maxNSamplesPerChannel;
end

if downsamplefactor>1
    precision = [num2str(nChannels),'*',precision];
    skip = nChannels*(downsamplefactor-1)*sampleSize;
    nSamplesPerChannel = floor(nSamplesPerChannel./downsamplefactor);
end

% For large amounts of data, read chunk by chunk
maxSamplesPerChunk = 10000;
nSamples = nSamplesPerChannel*nChannels;
if nSamples <= maxSamplesPerChunk,
	data = LoadChunk(f,nChannels,channels,nSamples/nChannels,precision,skip);
else
	% Determine chunk duration and number of chunks
	nSamplesPerChunk = floor(maxSamplesPerChunk/nChannels)*nChannels;
	nChunks = floor(nSamples/nSamplesPerChunk);
	% Preallocate memory
	data = zeros(nSamplesPerChannel,length(channels));
	% Read all chunks
	i = 1;
	for j = 1:nChunks,
		d = LoadChunk(f,nChannels,channels,nSamplesPerChunk/nChannels,precision,skip);
		[m,n] = size(d);
		if m == 0, break; end
		data(i:i+m-1,:) = d;
		i = i+m;
	end
	% If the data size is not a multiple of the chunk size, read the remainder
	remainder = nSamples - nChunks*nSamplesPerChunk;
	if remainder ~= 0,
		d = LoadChunk(f,nChannels,channels,remainder/nChannels,precision,skip);
		[m,n] = size(d);
		if m ~= 0,
			data(i:i+m-1,:) = d;
		end
	end
end

fclose(f);

% ---------------------------------------------------------------------------------------------------------

function data = LoadChunk(fid,nChannels,channels,nSamples,precision,skip)

if skip ~= 0,
	data = fread(fid,[nChannels nSamples],precision,skip);
    %data = fread(fid,[nChannels nSamples],[num2str(nChannels),'*',precision],nChannels*skip);
else
	data = fread(fid,[nChannels nSamples],precision);
end
data=data';

if isempty(data),
	warning('No data read (trying to read past file end?)');
elseif ~isempty(channels),
	data = data(:,channels);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [Vx,Vy] = computeFlowBkgSup(im1curr,im2curr,params,gparams)
% function [Vx,Vy] = computeFlowBkgSup(im1curr,im2curr,params)
% Computes flow and then suppresses flow in the background.

%There is some flow towards the center that can't be estimated
% properly. The flow increases as we go away from the center.
% d_err is to account for that.

sz = round(size(im1curr));
bwimg = zeros(sz(1),sz(2));
ctr = [ceil( (sz(1)+1)/2),ceil( (sz(2)+1)/2)];
bwimg(ctr(1),ctr(2))=1;
dimg = bwdist(bwimg,'euclidean');
[xx,yy]= meshgrid(1:sz(2),1:sz(1));
aimg = atan2(-(yy-ctr(1)),-(xx-ctr(2)));

dd_err = dimg/size(im1curr,1);
flow_thres = gparams.flow_thres;

uv = estimate_flow_interface(im1curr,im2curr,...
  'hs-brightness',{'max_warping_iters',gparams.warping_iters});

if ~params.stationary,
  Vx = uv(:,:,1);
  Vy = uv(:,:,2);
  return;
end

cdx = params.dx;
cdy = params.dy;
curt = params.theta;
rotd = [cos(curt) sin(curt); -sin(curt) cos(curt)]*[cdx;cdy];
cdx = -rotd(1); cdy = -rotd(2);

ctheta = params.dtheta;
rotflowu = dimg.*(cos(aimg+ctheta)-cos(aimg));
rotflowv = dimg.*(sin(aimg+ctheta)-sin(aimg));

dd1 = sqrt( (uv(:,:,1)-cdx-rotflowu).^2 + (uv(:,:,2)-cdy-rotflowv).^2);

crdx = params.rdx;  crdy = params.rdy;
rotd = [cos(curt) sin(curt); -sin(curt) cos(curt)]*[crdx;crdy];
fly_flow = rotd;

dd2 = sqrt( (uv(:,:,1)-fly_flow(1)).^2 + (uv(:,:,2)-fly_flow(2)).^2);
dd = min(dd1,dd2);
for ndx = 1:2
  tt = uv(:,:,ndx);
  tt( dd<(dd_err+flow_thres)) = 0;
  uv(:,:,ndx) = tt;
end
Vx = uv(:,:,1);
Vy = uv(:,:,2);

author: kristinbranson

repo: JAABA


text:
function fmincon_test ( )

%*****************************************************************************80
%
%% FMINCON_TEST tests the FMINCON library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FMINCON_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the FMINCON library.\n' );

  fmincon_test01 ( );
  fmincon_test02 ( );
  fmincon_test03 ( );
  fmincon_test04 ( );
  fmincon_test05 ( );
  fmincon_test06 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FMINCON_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout=sde_decorrelate(d,r2)
%SDE_DECORRELATE  Decorrelate correlated values.
%   R1 = SDE_DECORRELATE(D, R2) returns the matrix R1 of M decorrelated values
%   given the N-by-N diffusion matrix D and the M-by-N matrix of (correlated)
%   values R2. The first column of R1 is always SQRT(D(1,1)) divided by the
%   first column of R2.
%
%   [R1, C] = SDE_DECORRELATE(D, R2) also returns the N-by-N correlation matrix,
%   C, decorrelated from the diffusion matrix, D.
%
%   C = SDE_DECORRELATE(D) without a second input argument returns the N-by-N
%   correlation matrix, C, decorrelated from the diffusion matrix, D.
%
%   Note:
%       R1 = SDE_DECORRELATE(D, SDE_CORRELATE(C, R1)) is only accurate to within
%       the floating point machine precision, EPS.
%
%   Example:
%       % Correlate and then decorrelate normally-distributed samples
%       r1 = randn(5,3)
%       c = [1 0.8 0.2;0.8 1 0.5;0.2 0.5 1];
%       [r2, d] = sde_correlate(c,r1)
%       r3 = sde_decorrelate(d,r2)
%
%   See also: SDE_CORRELATE, RAND, RANDN, RANDSTREAM, RANDSTREAM/RANDN, EPS

%   Andrew D. Horchler, horchler @ gmail . com, Created 5-20-13
%   Revision: 1.2, 7-12-13

if nargout > min(nargin,2)
    error('SDETools:sde_decorrelate:TooManyOutputs',...
          'Too many output arguments for number of supplied inputs.');
end

if isempty(d) || isempty(r2)
    if nargin == 1
        varargout{1} = [];
    else
        varargout{2} = [];
    end
else
    [m,n] = size(d);
    if ndims(d) ~= 2 || m ~= n              %#ok<ISMAT>
        error('SDETools:sde_decorrelate:NonSquareMatrix',...
              'The diffusion matrix must be square.');
    end

    if isscalar(d)
        c = d^2;
        isDiag = false;
	elseif sde_isdiag(d)
        d = diag(d);
        c = diag(d.^2);
        isDiag = true;
    else
        c = d^2;
        isDiag = false;
    end

    if nargin == 1
        varargout{1} = c;
    else
        if ndims(r2) ~= 2 || size(r2,2) ~= m	%#ok<ISMAT>
            error('SDETools:sde_decorrelate:DimensionMismatch',...
                 ['The number of columns in the matrix of normally '...
                  'distributed values must equal the dimension of the '...
                  'correlation matrix.']);
        end

        if isDiag
            varargout{1} = bsxfun(@mrdivide,r2,d);
        else
            varargout{1} = r2/d;
        end
        if nargout == 2
            varargout{2} = c;
        end
    end
end

author: horchler

repo: SDETools

save_path: github-repos/MATLAB/horchler-SDETools


text:
classdef BackgroundSubtractorLSBP < handle
    %BACKGROUNDSUBTRACTORLSBP  Background Subtraction using Local SVD Binary Pattern
    %
    % More details about the algorithm can be found at [LGuo2016].
    %
    % It is based on LSBP feature descriptors and achieves state-of-the-art
    % performance on the CDnet 2012 dataset. LSBP descriptors are particularly
    % good in regions with illumination variation, noise and shadows. So, this
    % algorithm has better performance in this kind of regions.
    %
    % After extraction of LSBP descriptors, the algorithm processes frames
    % pixel-wise (i.e independently). Thus the implementation is parallelized
    % and fast enough for real-time processing.
    %
    % ## References
    % [LGuo2016]:
    % > L. Guo, D. Xu, and Z. Qiang. "Background Subtraction using Local SVD
    % > Binary Pattern". In 2016 IEEE Conference on Computer Vision and
    % > Pattern Recognition Workshops (CVPRW), pages 1159-1167, June 2016.
    % > [PDF](http://www.cv-foundation.org/openaccess/content_cvpr_2016_workshops/w24/papers/Guo_Background_Subtraction_Using_CVPR_2016_paper.pdf)
    %
    % See also: cv.BackgroundSubtractorLSBP.BackgroundSubtractorLSBP,
    %  cv.BackgroundSubtractorLSBP.apply,
    %  cv.BackgroundSubtractorLSBP.getBackgroundImage
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    %% BackgroundSubtractor
    methods
        function this = BackgroundSubtractorLSBP(varargin)
            %BACKGROUNDSUBTRACTORLSBP  Creates a LSBP Background Subtractor
            %
            %     bs = cv.BackgroundSubtractorLSBP()
            %     bs = cv.BackgroundSubtractorLSBP('OptionName', optionValue, ...)
            %
            % ## Options
            % * __MotionCompensation__ Whether to use camera motion
            %   compensation. One of:
            %   * __None__ (default)
            %   * __LK__
            % * __NSamples__ Number of samples to maintain at each point of
            %   the frame. default 20
            % * __LSBPRadius__ LSBP descriptor radius. default 16
            % * __TLower__ Lower bound for T-values. See [LGuo2016] for
            %   details. default 2.0
            % * __TUpper__ Upper bound for T-values. See [LGuo2016] for
            %   details. default 32.0
            % * __TInc__ Increase step for T-values. See [LGuo2016] for
            %   details. default 1.0
            % * __TDec__ Decrease step for T-values. See [LGuo2016] for
            %   details. default 0.05
            % * __RScale__ Scale coefficient for threshold values. default 10.0
            % * __RIncDec__ Increase/Decrease step for threshold values.
            %   default 0.005
            % * __NoiseRemovalThresholdFacBG__ Strength of the noise removal
            %   for background points. default 0.0004
            % * __NoiseRemovalThresholdFacFG__ Strength of the noise removal
            %   for foreground points. default 0.0008
            % * __LSBPThreshold__ Threshold for LSBP binary string. default 8
            % * __MinCount__ Minimal number of matches for sample to be
            %   considered as foreground. default 2
            %
            % See also: cv.BackgroundSubtractorLSBP
            %
            this.id = BackgroundSubtractorLSBP_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     bs.delete()
            %
            % See also: cv.BackgroundSubtractorLSBP
            %
            if isempty(this.id), return; end
            BackgroundSubtractorLSBP_(this.id, 'delete');
        end

        function fgmask = apply(this, im, varargin)
            %APPLY  Updates the background model and computes the foreground mask
            %
            %     fgmask = bs.apply(im)
            %     fgmask = bs.apply(im, 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __im__ Next video frame.
            %
            % ## Output
            % * __fgmask__ The output foreground mask as an 8-bit binary image
            %   (0 for background, 255 for foregound).
            %
            % ## Options
            % * __LearningRate__ The value between 0 and 1 that indicates how
            %   fast the background model is learnt. Negative parameter value
            %   makes the algorithm to use some automatically chosen learning
            %   rate. 0 means that the background model is not updated at all,
            %   1 means that the background model is completely reinitialized
            %   from the last frame. default -1
            %
            % See also: cv.BackgroundSubtractorLSBP.getBackgroundImage
            %
            fgmask = BackgroundSubtractorLSBP_(this.id, 'apply', im, varargin{:});
        end

        function bgImg = getBackgroundImage(this)
            %GETBACKGROUNDIMAGE  Computes a background image
            %
            %     bgImg = bs.getBackgroundImage()
            %
            % ## Output
            % * __bgImg__ The output background image.
            %
            % See also: cv.BackgroundSubtractorLSBP.apply
            %
            bgImg = BackgroundSubtractorLSBP_(this.id, 'getBackgroundImage');
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.BackgroundSubtractorLSBP.empty
            %
            BackgroundSubtractorLSBP_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the algorithm is empty (e.g. in the very
            %   beginning or after unsuccessful read).
            %
            % See also: cv.BackgroundSubtractorLSBP.clear
            %
            b = BackgroundSubtractorLSBP_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.BackgroundSubtractorLSBP.save,
            %  cv.BackgroundSubtractorLSBP.load
            %
            name = BackgroundSubtractorLSBP_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in a file storage.
            %
            % See also: cv.BackgroundSubtractorLSBP.load
            %
            BackgroundSubtractorLSBP_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from a file storage.
            % The previous model state is discarded.
            %
            % See also: cv.BackgroundSubtractorLSBP.save
            %
            BackgroundSubtractorLSBP_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

    %% BackgroundSubtractorLSBPDesc
    methods (Static)
        function desc = computeLSBPDesc(frame, LSBPSamplePoints)
            %COMPUTELSBPDESC  This is for calculation of the LSBP descriptors
            %
            %     desc = cv.BackgroundSubtractorLSBP.computeLSBPDesc(frame, LSBPSamplePoints)
            %
            % ## Input
            % * __frame__ input frame
            % * __LSBPSamplePoints__ 32 sample points
            %
            % ## Output
            % * __desc__ LSBP descriptors
            %
            desc = BackgroundSubtractorLSBP_(0, 'computeLSBPDesc', frame, LSBPSamplePoints);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Fig. 9.17   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%script for Figure 5.59
sysG=tf([1],[1 0.2 1 0]);
sysH=tf(123*[1 .18 .81],[1 20 100]);
sysL=sysG*sysH;
rlocus(sysL)
title('Figure 9.17  Root locus for the system of Figure 9.18')
axis([-3 1 -1.5 1.5])
hold on
sysCL=feedback(sysG,sysH);
r=eig(sysCL);
plot(r,'*')
z=0:.1:.9;
wn= .5:.5:3;
sgrid(z, wn)
hold off

author: Sable

repo: mcbench-benchmarks


text:
%% ODF Analysis
%
%%

author: mtex-toolbox

repo: mtex

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


text:
clear all; close all; clc;

% Draws poles and zeros on a complex plane
% Coefficients of the rational filter
b = [1 0.7 0.6];
a = [1 -1.5 0.9];
spx.graphics.plot.rational_poles_zeros(b,a);

author: indigits

repo: sparse-plex


text:
function [p] = integrateSubIntervals(x, cdf)

[xLB, xUB] = subIntervals(x);

cdfUB = cdf(xUB);
cdfLB = cdf(xLB);

p = (cdfUB - cdfLB) ./ (xUB - xLB);

end


text:
function [satpos,satclk,vare,varc,stat]=pephpos(time,sat,nav)

global glc
stat=1; norder=10; %polynomial fitting order
tdiff=zeros(norder+1,1); p=zeros(3,norder+1);
satpos=zeros(3,1); s=zeros(3,1); satclk=0;
varc=0;vare=0;

if nav.np<norder||timediff(time,nav.peph(1).time)<-900||...
        timediff(time,nav.peph(end).time)>900
    stat=0; return;
end

%binary search
i=0;
j=nav.np-1;
while i<j
    k=fix((i+j)/2);
    if timediff(nav.peph(k+1).time,time)<0
        i=k+1;
    else
        j=k;
    end
end

if i<=0,index=0;  else,index=i-1;  end
index=index+1;

idx=index-norder/2; %start of fitting interval
if idx<1,idx=1;  elseif idx+norder>nav.np,idx=nav.np-norder;  end

nt=1;
for j=0:norder
    tdiff(nt)=timediff(nav.peph(idx+j).time,time); %time difference
    nt=nt+1;
    if dot(nav.peph(idx+j).pos(sat,1:3),nav.peph(idx+j).pos(sat,1:3))<=0
        stat=0; return;
    end
end

for j=0:norder
    pos=nav.peph(idx+j).pos(sat,:);

    %correciton for earh rotation
    sinl=sin(glc.OMGE*tdiff(j+1));
    cosl=cos(glc.OMGE*tdiff(j+1));

    p(1,j+1)=cosl*pos(1)-sinl*pos(2);
    p(2,j+1)=sinl*pos(1)+cosl*pos(2);
    p(3,j+1)=pos(3);

end

%polynomial interpolation using Neville's algorithm
for i=1:3
    satpos(i,1)=polyinter(tdiff,p(i,:));
end

%calculate satellite position variance
for i=1:3
    s(i)=nav.peph(index).std(sat,i);
end
std=sqrt(dot(s,s)); %error caused by ephemeris
if tdiff(1)>0
    std=std+(5E-7*tdiff(1)^2)/2; %error caused by fitting
elseif tdiff(end)<0
    std=std+(5E-7*tdiff(end)^2)/2;
end
vare=std^2;

%linear interpolation for clock
t0=timediff(time,nav.peph(index).time);
t1=timediff(time,nav.peph(index+1).time);
c0=nav.peph(index).pos(sat,4);
c1=nav.peph(index+1).pos(sat,4);

if t0<=0
    satclk=c0;
    if satclk~=0
        varc=(nav.peph(index).std(sat,4)*glc.CLIGHT-1E-3*t0)^2;
    end
elseif t1>=0
    satclk=c1;
    if satclk~=0
        varc=(nav.peph(index+1).std(sat,4)*glc.CLIGHT+1E-3*t1)^2;
    end
elseif c0~=0 && c1~=0
    satclk=(c1*t0-c0*t1)/(t0-t1);
    if abs(t0)-abs(t1)>0
        varc=(nav.peph(index+1).std(sat,4)*glc.CLIGHT+1E-3*abs(t1))^2;
    else
        varc=(nav.peph(index).std(sat,4)*glc.CLIGHT+1E-3*abs(t0))^2;
    end
else
    satclk=0;
    varc=0;
end

return

repo: GINav


text:
classdef prtUiManagerImage < prtUiManagerAxes

    properties (Dependent)

    end

    properties
        imageHandle = [];
        climSliderHandle = [];
    end
    properties (Hidden)
        axesPosition
    end

    methods

        function tickoff(self)
            set(self.managedHandle,'xtick',[]);
            set(self.managedHandle,'ytick',[]);
        end

        function self = prtUiManagerImage(varargin)
           if nargin
               self = prtUtilAssignStringValuePairs(self, varargin{:});
           end
        end

        function reImage(self, varargin)
            if isempty(self.imageHandle)
                self.imageHandle = image(varargin{:},'Parent',self.managedHandle);
            else
                self.updateCData(varargin{1});
            end
        end

        function reImagesc(self, varargin)
            if isempty(self.imageHandle)
                self.imageHandle = imagesc(varargin{:},'Parent',self.managedHandle);
            else
                self.updateCData(varargin{1});
            end
        end

        function reImshow(self, varargin)
            if isempty(self.imageHandle)
                self.imageHandle = imshow(varargin{:},'Parent',self.managedHandle);
            else
                self.updateCData(varargin{1});
            end
        end

        function image(self,varargin)
            self.imageHandle = image(varargin{:},'Parent',self.managedHandle);
        end

        function imagesc(self,varargin)
            self.imageHandle = imagesc(varargin{:},'Parent',self.managedHandle);
        end

        function imshow(self,varargin)
            self.imageHandle = imshow(varargin{:},'Parent',self.managedHandle);
        end

        function updateCData(self, cData)
            % updateData(self, cData)
            self.checkValid;

            prevCdata = get(self.imageHandle,'cdata');
            if ~isequal(size(prevCdata),size(cData))
                warning('prt:prtUiManagerImage:changeCdataSize','Attempt to change the size of the data in an image by setting the cData, try re-creating the image');
            end
            set(self.imageHandle,'CData',cData);
        end

        function ih = get.imageHandle(self)
            ih = self.imageHandle;
        end

        function shapePoly = getRoiPoly(self)
            % Output either a struct or an empty array
            self.checkValid;

            axes(self.managedHandle); %#ok<MAXES> %Because roipoly doesn't let you specify axes
            [~,xi,yi] = roipoly;
            shapePoly = [];
            if ~isempty(xi)
                shapePoly = struct('type','poly','xi',xi,'yi',yi);
            end
        end

        function shapeLine = getRoiLine(self)
            % Output either a struct or an empty array
            self.checkValid;

            axes(self.managedHandle); %#ok<MAXES> %Because roipoly doesn't let you specify axes
            [~,xi,yi] = roipoly;
            shapeLine = [];
            if ~isempty(xi)
                xi = xi(1:end-1);
                yi = yi(1:end-1);
                shapeLine = struct('type','line','xi',xi,'yi',yi);
            end
        end

        function shapeRect = getRoiRect(self)
            % Output either a struct or an empty array
            self.checkValid;

            axes(self.managedHandle); %#ok<MAXES>  %Because imcrop doesn't let you specify axes
            [~,rect] = imcrop;
            shapeRect = [];
            if ~isempty(rect)
                xi = [rect(1);rect(1);rect(1)+rect(3);rect(1)+rect(3);rect(1)];
                yi = [rect(2);rect(2)+rect(4);rect(2)+rect(4);rect(2);rect(2)];
                shapeRect = struct('type','rect','xi',xi,'yi',yi);
            end

        end

        function shapePoint = getRoiPoint(self)
            % Output either a struct or an empty array
            self.checkValid;

            axes(self.managedHandle); %#ok<MAXES>  %Because imcrop doesn't let you specify axes
            [xi, yi] = getpts;
            shapePoint = [];
            if ~isempty(xi)
                shapePoint = struct('type','point','xi',xi(end),'yi',yi(end));
            end

        end

        function is = checkValid(self)
            assert(ishandle(self.imageHandle),'prt:prtUiManagerImage:badHandle','Some or all of the requested handles are no longer valid. The axes may have been deleted.');
            is = true;
        end

        function h = plotShape(self,shapeArray)
            self.checkValid;

            axes(self.managedHandle); %#ok<MAXES>  %Because ishold doesn't let you specify axes
            holdState = ishold;
            h = nan(size(shapeArray));
            hold(self.managedHandle,'on');
            %Plot each shape:
            try
                for i = 1:length(shapeArray)
                    if length(shapeArray(i).xi) == 1
                        h(i) = plot(shapeArray(i).xi,shapeArray(i).yi,'bo'); %points
                    else
                        h(i) = plot(shapeArray(i).xi,shapeArray(i).yi,'b');
                    end
                end
            catch ME
                disp(ME)
                debugdlg('I had trouble plotting one of your objects.  Do you want to debug?');
            end
            if ~holdState
                hold(self.managedHandle,'off');
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [a,t]=rotqr2ax(q)
%ROTQR2AX converts a real quaternion to the corresponding rotation axis and angle
% Inputs:
%
%     Q(4,1)   real-valued quaternion (with magnitude = 1)
%
% Outputs:
%
%     A(3,1)   Unit vector in the direction of the rotation axis.
%     T        Rotation angle in radians (in range 0 to 2pi)
%
% In the quaternion representation of a rotation, and q(1) = cos(t/2)
% where t is the angle of rotation in the range 0 to 2pi
% and q(2:4)/sin(t/2) is a unit vector lying along the axis of rotation
% a positive rotation about [0 0 1] takes the X axis towards the Y axis.
%
%      Copyright (C) Mike Brookes 2007
%      Version: $Id: rotqr2ax.m,v 1.2 2007/11/21 12:42:36 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=q(2:4);
m=sqrt(a'*a);
t=2*atan2(m,q(1));      % avoids problems if unnormalized
if m>0
    a=a/m;
else
    a=[0 0 1]';
end

author: decouples


text:
function [p params] = glm_getHrfParams(params, select);
%
% subParams = glm_getHrfParams(params, [select=0]);
%
% If you're using a pre-defined HRF for GLM analyses,
% such as 'boynton', 'spm', or 'dale&buckner', return those
% arguments which should be called for those canned function.
% These arguments are kept in a params.glmHRF_params field.
% If this field isn't found, or isn't the right size, ask
% the user to specify these, with appropriate defaults.
% Will also return the params struct with the glmHRF_params field
% set properly.
%
% If the select flag is passed in as 1, will put up a dialog for the
% user to specify the parameters. Otherwise, will set to reasonable
% defaults.
%
% ras, 01/2007.

% EDIT HISTORY:
% ras, 02/2007 -- added select flag, GUI dialogs.
% DY, 02/20/2007 -- removed checks for number of parameters such that
% whenever select ~=1, the parameters are automatically set to the
% defaults, even if there do appear to be the right number of parameters
% (since odd bugs can result if there are the right number of parameters,
% but the values are off).
% ras, 02/21/2007 -- Davie, I get your point that there is this possibility;
% however, since this function is used to retrieve stored values, your
% change would cause only the defaults to ever be used (see glm_hrf). It
% may be possible to change the way this function is called, but this would
% require additional logic (since you also need a check that the right
% number of HRF params are specified at HRF creation); and the possibility
% you're concerned about could only occur if you manually muck with the
% parameters -- none of the accessor functions can set the wrong # of
% params. So, I'm reverting for now, and we'll see if we need to do more.

p = [];

if ~exist('select', 'var') | isempty(select), select = 0;   end

% this is only needed for predefined HRF options
if ~ismember(params.glmHRF, [2 3 4]), return; end

% init params field
if ~isfield(params, 'glmHRF_params')
    params.glmHRF_params = [];
end

% get params
switch params.glmHRF
    case 2,     % boynton
        if select==1
            params.glmHRF_params = boyntonHIRF_dialog;
        elseif length(params.glmHRF_params) ~= 3 % wrong size
            % [n tau delta]
            params.glmHRF_params = [3 1.08 2.05];
        end

    case 3,     % spm
        if select==1
            params.glmHRF_params = spmHRF_dialog;
            %	p(1) - delay of response (relative to onset)	   6
            %	p(2) - delay of undershoot (relative to onset)    16
            %	p(3) - dispersion of response			   1
            %	p(4) - dispersion of undershoot			   1
            %	p(5) - ratio of response to undershoot		   6
            %	p(6) - onset (seconds)				   0
            %	p(7) - length of kernel (seconds)		  32
		elseif length(params.glmHRF_params) ~= 6 % wrong size
            maxT = max(params.timeWindow);
            params.glmHRF_params = [6 16 1 1 6 0 maxT];
        end

    case 4,     % dale & buckner
        if select==1
            params.glmHRF_params = daleBucknerHIRF_dialog;
            % [delta tau]
%             params.glmHRF_params = [2.25 1.25]; % more event-friendly
        elseif length(params.glmHRF_params) ~= 2 % wrong size
            params.glmHRF_params = [1.25 2.5]; % vals used in er_runSelxavgBlock (old code)
        end
end

p = params.glmHRF_params;

return
% /----------------------------------------------------------------/ %

% /----------------------------------------------------------------/ %
function p = boyntonHIRF_dialog;
% dialog to set [n tau delta] for boyntonHIRF function.
dlg(1).fieldName = 'eqn';
dlg(1).style = 'text';
dlg(1).string = 'HRF equation:';
dlg(1).value = 'h(t) = [(t/tau) ^ (n-1) * exp(-t/tau)] / [tau(n-1)!]';

dlg(2).fieldName = 'n';
dlg(2).style = 'edit';
dlg(2).string = 'n (integer):';
dlg(2).value = '3';

dlg(3).fieldName = 'tau';
dlg(3).style = 'edit';
dlg(3).string = 'tau (decay), secs:';
dlg(3).value = '1.08';

dlg(4).fieldName = 'delta';
dlg(4).style = 'edit';
dlg(4).string = 'delay (delta), secs:';
dlg(4).value = '2.05';

resp = generalDialog(dlg, 'Boynton HRF');
if isempty(resp)
    error('User canceled.')
end

p = [str2num(resp.n) str2num(resp.tau) str2num(resp.delta)];
return
% /----------------------------------------------------------------/ %

% /----------------------------------------------------------------/ %
function p = spmHRF_dialog;
% dialog to set the params for spm_hrf:
%	p(1) - delay of response (relative to onset)	   6
%	p(2) - delay of undershoot (relative to onset)    16
%	p(3) - dispersion of response			   1
%	p(4) - dispersion of undershoot			   1
%	p(5) - ratio of response to undershoot		   6
%	p(6) - onset (seconds)				   0
%	p(7) - length of kernel (seconds)		  32
dlg(1).fieldName = 'responseDelay';
dlg(end).style = 'edit';
dlg(end).string = 'delay of response (relative to onset)';
dlg(end).value = '6';

dlg(end+1).fieldName = 'undershootDelay';
dlg(end).style = 'edit';
dlg(end).string = 'delay of undershoot (relative to onset)';
dlg(end).value = '16';

dlg(end+1).fieldName = 'responseDispersion';
dlg(end).style = 'edit';
dlg(end).string = 'dispersion of response';
dlg(end).value = '1';

dlg(end+1).fieldName = 'undershootDispersion';
dlg(end).style = 'edit';
dlg(end).string = 'dispersion of response';
dlg(end).value = '1';

dlg(end+1).fieldName = 'ratio';
dlg(end).style = 'edit';
dlg(end).string = 'ratio of response to undershoot';
dlg(end).value = '6';

dlg(end+1).fieldName = 'onset';
dlg(end).style = 'edit';
dlg(end).string = 'onset (seconds)';
dlg(end).value = '0';

dlg(end+1).fieldName = 'kernelLength';
dlg(end).style = 'edit';
dlg(end).string = 'length of kernel (seconds)';
dlg(end).value = '32';

resp = generalDialog(dlg, 'SPM HRF');
if isempty(resp)
    error('User canceled.')
end

p = [str2num(resp.responseDelay) str2num(resp.undershootDelay) ...
     str2num(resp.responseDispersion) str2num(resp.undershootDispersion) ...
     str2num(resp.ratio) str2num(resp.onset) str2num(resp.kernelLength)];

return
% /----------------------------------------------------------------/ %

% /----------------------------------------------------------------/ %
function p = daleBucknerHIRF_dialog;
% dialog to set [delta tau] for fmri_hemodyn function.
dlg(1).fieldName = 'eqn';
dlg(1).style = 'text';
dlg(1).string = 'HRF equation:';
dlg(1).value = {'h(t>delta)  = ((t-delta)/tau)^2 * exp(-(t-delta)/tau)'; ...
                 'h(t<=delta) = 0;'};

dlg(2).fieldName = 'delta';
dlg(2).style = 'edit';
dlg(2).string = 'delay (delta), secs:';
dlg(2).value = '1.25';

dlg(3).fieldName = 'tau';
dlg(3).style = 'edit';
dlg(3).string = 'tau (decay), secs:';
dlg(3).value = '2.5';

resp = generalDialog(dlg, 'Dale & Buckner HRF');
if isempty(resp)
    error('User canceled.')
end

p = [str2num(resp.delta) str2num(resp.tau)];

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [PR, PN, PC, PL] = subspaceProjector(A, printLevel, sub_space)
% Returns the matrix for projection onto the subspace of matrix A,
% specified by `sub_space`
%   * 'all'
%   * 'R' row space         i.e. range(A')
%   * 'N' nullspace         i.e. null(A)
%   * 'C' column space      i.e. range(A)
%   * 'L' left nullspace    i.e. null(A')
%
% Let M denote the Moore-Penrose pseudoinverse of the A and the subscripts are the following
% `_R` row space        i.e. range(A')
% `_N` nullspace        i.e. null(A)
% `_C` column space     i.e. range(A)
% `_L` left nullspace   i.e. null(A')
%
% Let
%
% .. math::
%    v   &= v_R + v_N \\
%    v_R &= M A v = PR v \\
%    v_N &= (I - M A) v = PN v
%
% Let
%
% .. math::
%    u   &= u_C + u_L \\
%    u_C &= A M u = PC u \\
%    u_L &= (I - A M) u = PL u
%
% Examples:
%
% Given :math:`A v = b`, then :math:`v_R = M b`
%
% Given :math:`A^Tu = q`, then :math:`u_C = M^T q`
%
% USAGE:
%
%    [PR, PN, PC, PL] = subspaceProjector(model, printLevel, sub_space)
%
% INPUT:
%    A          `m x n` matrix
%
% OPTIONAL INPUTS:
%    printLevel:          {(1), 0}, 1 = print diagnostics, 0 = silent
%    sub_space:           returns projection matrices onto all or one select
%
%                           * sub_space
%                           * 'all'
%                           * 'R' row space
%                           * 'N' nullspace
%                           * 'C' column space
%                           * 'L' left nullspace
%
% OUTPUTS:
%    [PR, PN, PC, PL]:    matrices for projection onto the row, null,
%                         column and left nullspace of A, respectively
%
% .. Author:
%       - 10 July 2009 : Ronan Fleming. First Version.
%       - 10 Aug  2009 : Changed to use Micheal Saunders faster approach
%       -    Jan  2018 : Changed to take a matrix A

if ~exist('printLevel','var')
    printLevel=1;
end

if ~exist('sub_space','var')
    sub_space='all';
end

[nMet,nRxn]=size(A);

archstr = computer('arch');
switch archstr
    case {'glnx86','glnxa64'}
        %A = U1*D1*V1'
        if printLevel
            fprintf('%s','Calculating SVD ...');
            tic
        end
        %Michael Saunders code
        [U1,D1,V1,r] = subspaceSVD(A);
        if printLevel
            fprintf('%s\n',[' finished. toc = ' num2str(toc)]);
        end
        PR=[];PN=[];PC=[];PL=[];
        if strcmp(sub_space,'R')
            PR=V1*V1';
        elseif strcmp(sub_space,'N')
            PN=eye(nRxn) - V1*V1';
        elseif strcmp(sub_space,'C')
            PC=U1*U1';
        elseif strcmp(sub_space,'L')
            PL=eye(nMet) - U1*U1';
        elseif strcmp(sub_space,'all')
            PR=V1*V1';
            PN=eye(nRxn) - V1*V1';
            PC=U1*U1';
            PL=eye(nMet) - U1*U1';
        end
    otherwise
        %for other architectures calculate the Moore-Penrose Pseudoinverse
        if printLevel
            fprintf('%s','Calculating the Moore-Penrose Pseudoinverse...');
            tic
        end
        M=pinv(full(A));
        if printLevel
            fprintf('%s\n',[' finished. toc = ' num2str(toc)]);
        end

        PR=[];PN=[];PC=[];PL=[];
        if strcmp(sub_space,'R')
            PR=M*A;
        elseif strcmp(sub_space,'N')
            PN=eye(nRxn)-M*A;
        elseif strcmp(sub_space,'C')
            PC=A*M;
        elseif strcmp(sub_space,'L')
            PL=eye(nMet)-A*M;
        elseif strcmp(sub_space,'all')
            PR=M*A;
            PN=eye(nRxn)-M*A;
            PC=A*M;
            PL=eye(nMet)-A*M;
        end
end

author: opencobra

repo: cobratoolbox


text:
function fx = p14_fun ( n, x )

%*****************************************************************************80
%
%% P14_FUN evaluates the integrand for problem 14.
%
%  Discussion:
%
%    S&S gives "exact" value as     1.0634618101...
%    Mathematica returns            1.0634618101722400407...
%    S&S gives Laguerre(16) as      1.0634713425...
%    S&S gives EXP_TRANSFORM(16) as 1.0634618101...
%
%    The FORTRAN version of this routine, compiled with G95, was getting
%    a floating point exception when evaluating the integrand
%    and using a Laguerre rule of order 64.  So I have had to truncate
%    the evaluation of the exponential.
%
%  Integral:
%
%    Integral ( 0 <= x < +oo ) sin ( exp ( - x ) + exp ( - 4 x ) ) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the function values.
%
  fx(1:n) = sin ( exp ( -x(1:n) ) + exp ( -4.0 * x(1:n) ) );

  return
end


text:
function perm0_uniform_test ( )

%*****************************************************************************80
%
%% PERM0_UNIFORM_TEST tests PERM0_UNIFORM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2015
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PERM0_UNIFORM_TEST\n' );
  fprintf ( 1, '  PERM0_UNIFORM randomly selects a permutation of 0,...,N-1.\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for test = 1 : 5

    [ p, seed ] = perm0_uniform ( n, seed );

    fprintf ( 1, '  ' );
    for i = 1 : n
      fprintf ( 1, '%4d', p(i) );
    end
    fprintf ( 1, '\n' );

  end

  return
end


text:
%

function [f, df] = sampleEijOpt(x)
% optimization code

% param
vert = x(4:7);
w_res = x(8);

% box case (4 lines), three free parameters
x0 = x(1); y0 = x(2); % camera center
yc = x(3); % corner

% assume vertical lines are fixed
a_aob = (vert(2) - vert (1))/w_res*2*pi;
a_boc = (vert(3) - vert (2))/w_res*2*pi;
a_cod = (vert(4) - vert (3))/w_res*2*pi;
a_doa = (vert(1) + w_res - vert (4))/w_res*2*pi;

% energy
v_ao  = [x0 y0]; v_bo = [x0-1 y0]; v_co = [x0-1 y0 - yc]; v_do = [x0 y0-yc];
n_v_ao = norm(v_ao); n_v_bo = norm(v_bo); n_v_co = norm(v_co); n_v_do = norm(v_do);
b_aob = acos(dot(v_ao, v_bo)/n_v_ao/n_v_bo);
if det([v_ao;v_bo]) < 0
    b_aob = 2*pi - b_aob;
end
b_boc = acos(dot(v_bo, v_co)/n_v_bo/n_v_co);
if det([v_bo;v_co]) < 0
    b_boc = 2*pi - b_boc;
end
b_cod = acos(dot(v_co, v_do)/n_v_co/n_v_do);
if det([v_co;v_do]) < 0
    b_cod = 2*pi - b_cod;
end
b_doa = acos(dot(v_do, v_ao)/n_v_do/n_v_ao);
if det([v_do;v_ao]) < 0
    b_doa = 2*pi - b_doa;
end

f = (b_aob - a_aob)^2 + (b_boc - a_boc)^2 + (b_cod - a_cod)^2 + (b_doa - a_doa)^2;

% gradient
% x0
d_aob_x0 = (2*x0-1)*n_v_ao*n_v_bo + dot(v_ao, v_bo) * (x0*n_v_bo/n_v_ao + (x0-1)*n_v_ao/n_v_bo);
d_aob_x0 = d_aob_x0 * (-1/(sqrt(1-cos(b_aob)*cos(b_aob))+eps))/n_v_ao/n_v_ao/n_v_bo/n_v_bo;
if det([v_ao;v_bo]) < 0
    d_aob_x0 = -d_aob_x0;
end
d_aob_x0 = 2*(b_aob - a_aob) * d_aob_x0;

d_boc_x0 = 2*(x0-1)*n_v_bo*n_v_co + dot(v_bo, v_co) * ((x0-1)*n_v_co/n_v_bo + (x0-1)*n_v_bo/n_v_co);
d_boc_x0 = d_boc_x0 * (-1/(sqrt(1-cos(b_boc)*cos(b_boc))+eps))/n_v_bo/n_v_bo/n_v_co/n_v_co;
if det([v_bo;v_co]) < 0
    d_boc_x0 = -d_boc_x0;
end
d_boc_x0 = 2*(b_boc - a_boc) * d_boc_x0;

d_cod_x0 = (2*x0-1)*n_v_co*n_v_do + dot(v_co, v_do) * ((x0-1)*n_v_do/n_v_co + x0*n_v_co/n_v_do);
d_cod_x0 = d_cod_x0 * (-1/(sqrt(1-cos(b_cod)*cos(b_cod))+eps))/n_v_co/n_v_co/n_v_do/n_v_do;
if det([v_co;v_do]) < 0
    d_cod_x0 = -d_cod_x0;
end
d_cod_x0 = 2*(b_cod - a_cod) * d_cod_x0;

d_doa_x0 = 2*x0*n_v_do*n_v_ao +dot(v_do, v_ao) * (x0*n_v_ao/n_v_do + x0*n_v_do/n_v_ao);
d_doa_x0 = d_doa_x0 * (-1/(sqrt(1-cos(b_doa)*cos(b_doa))+eps))/n_v_do/n_v_do/n_v_ao/n_v_ao;
if det([v_do;v_ao]) < 0
    d_doa_x0 = -d_doa_x0;
end
d_doa_x0 = 2*(b_doa - a_doa) * d_doa_x0;

% y0
d_aob_y0 = 2*y0*n_v_ao*n_v_bo + dot(v_ao, v_bo) * (y0*n_v_bo/n_v_ao + y0*n_v_ao/n_v_bo);
d_aob_y0 = d_aob_y0 * (-1/(sqrt(1-cos(b_aob)*cos(b_aob))+eps))/n_v_ao/n_v_ao/n_v_bo/n_v_bo;
if det([v_ao;v_bo]) < 0
    d_aob_y0 = -d_aob_y0;
end
d_aob_y0 = 2*(b_aob - a_aob) * d_aob_y0;

d_boc_y0 = (2*y0-yc)*n_v_bo*n_v_co + dot(v_bo, v_co) * (y0*n_v_co/n_v_bo + (y0-yc)*n_v_bo/n_v_co );
d_boc_y0 = d_boc_y0/n_v_bo/n_v_bo/n_v_co/n_v_co;
if det([v_bo;v_co]) < 0
    d_boc_y0 = -d_boc_y0;
end
d_boc_y0 = 2*(b_boc - a_boc) * d_boc_y0;

d_cod_y0 = 2*(y0-yc)*n_v_co*n_v_do + dot(v_co, v_do) * ((y0-yc)*n_v_do/n_v_co + (y0-yc)*n_v_co/n_v_do);
d_cod_y0 = d_cod_y0 * (-1/(sqrt(1-cos(b_cod)*cos(b_cod))+eps))/n_v_co/n_v_co/n_v_do/n_v_do;
if det([v_co;v_do]) < 0
    d_cod_y0 = -d_cod_y0;
end
d_cod_y0 = 2*(b_cod - a_cod) * d_cod_y0;

d_doa_y0 = (2*y0-yc)*n_v_do*n_v_ao + dot(v_do, v_ao) * ((y0-yc)*n_v_ao/n_v_do + y0*n_v_do/n_v_ao);
d_doa_y0 = d_doa_y0 * (-1/(sqrt(1-cos(b_doa)*cos(b_doa))+eps))/n_v_do/n_v_do/n_v_ao/n_v_ao;
if det([v_do;v_ao]) < 0
    d_doa_y0 = -d_doa_y0;
end
d_doa_y0 = 2*(b_doa - a_doa) * d_doa_y0;

% yc
d_boc_yc = (-y0)*n_v_bo*n_v_co + dot(v_bo, v_co) * n_v_bo* (yc-y0)/n_v_co;
d_boc_yc = d_boc_yc/n_v_bo/n_v_bo/n_v_co/n_v_co;
if det([v_bo;v_co]) < 0
    d_boc_yc = -d_boc_yc;
end
d_boc_yc = 2*(b_boc - a_boc) * d_boc_yc;

d_cod_yc = 2*(yc-y0)*n_v_co*n_v_do + dot(v_co, v_do) * ((yc-y0)*n_v_do/n_v_co+(yc-y0)*n_v_co/n_v_do);
d_cod_yc = d_cod_yc * (-1/(sqrt(1-cos(b_cod)*cos(b_cod))+eps))/n_v_co/n_v_co/n_v_do/n_v_do;
if det([v_co;v_do]) < 0
    d_cod_yc = -d_cod_yc;
end
d_cod_yc = 2*(b_cod - a_cod) * d_cod_yc;

d_doa_yc = (-y0)*n_v_do*n_v_ao + dot(v_do, v_ao) * n_v_ao * (yc-y0)/n_v_do;
d_doa_yc = d_doa_yc * (-1/(sqrt(1-cos(b_doa)*cos(b_doa))+eps))/n_v_do/n_v_do/n_v_ao/n_v_ao;
if det([v_do;v_ao]) < 0
    d_doa_yc = -d_doa_yc;
end
d_doa_yc = 2*(b_doa - a_doa) * d_doa_yc;

d_x0 = d_aob_x0 + d_boc_x0 + d_cod_x0 + d_doa_x0;
d_y0 = d_aob_y0 + d_boc_y0 + d_cod_y0 + d_doa_y0;
d_yc = d_boc_yc + d_cod_yc + d_doa_yc;

df = [d_x0; d_y0; d_yc; 0; 0; 0; 0; 0];

author: zouchuhang

repo: LayoutNet


text:
function enum = deranged_enum ( n )

%*****************************************************************************80
%
%% DERANGED_ENUM returns the number of derangements of N objects.
%
%  Discussion:
%
%    A derangement of N objects is a permutation with no fixed
%    points.  If we symbolize the permutation operation by "P",
%    then for a derangment, P(I) is never equal to I.
%
%  Recursion:
%
%      D(0) = 1
%      D(1) = 0
%      D(2) = 1
%      D(N) = (N-1) * ( D(N-1) + D(N-2) )
%
%    or
%
%      D(0) = 1
%      D(1) = 0
%      D(N) = N * D(N-1) + (-1)**N
%
%  Formula:
%
%    D(N) = N! * ( 1 - 1/1! + 1/2! - 1/3! ... 1/N! )
%
%    Based on the inclusion/exclusion law.
%
%    D(N) is the number of ways of placing N non-attacking rooks on
%    an N by N chessboard with one diagonal deleted.
%
%    Limit ( N -> Infinity ) D(N)/N! = 1 / e.
%
%    The number of permutations with exactly K items in the right
%    place is COMB(N,K) * D(N-K).
%
%  First values:
%
%     N         D(N)
%     0           1
%     1           0
%     2           1
%     3           2
%     4           9
%     5          44
%     6         265
%     7        1854
%     8       14833
%     9      133496
%    10     1334961
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of objects to be permuted.
%
%    Output, integer ENUM, the number of derangements of N objects.
%
  if ( n < 0 )

    dn = 0;

  elseif ( n == 0 )

    dn = 1;

  elseif ( n == 1 )

    dn = 0;

  elseif ( n == 2 )

    dn = 1;

  else

    dnm1 = 0;
    dn = 1;

    for i = 3 : n
      dnm2 = dnm1;
      dnm1 = dn;
      dn = ( i - 1 ) * ( dnm1 + dnm2 );
    end

  end

  enum = dn;

  return
end


text:
 function r1d = ir_mri_dce_r1d(r10, Ct, varargin)
%function r1d = ir_mri_dce_r1d(r10, Ct, varargin)
%|
%| Make time series of R1=1/T1 values given baseline value(s) R10
%| and contrast agent time-concentration curves Ct.
%|
%| in
%|	r10	[1] | [Nc]	[1/s] baseline R1 values
%|				for Nc tissue classes
%|	Ct	[Nc Nt]		[mMmol] tissue constrast concentration curves
%|				for Nt time points
%|
%| option
%|	'r1'	[1]		relaxivity [1/mMol 1/s]; default: 4.5 (for 3T)
%|
%| out
%|	r1d	[Nc Nt]		[1/s] dynamic R1=1/T1 curves
%|
%| 2014-08-21 Jeff Fessler and Mai Le, University of Michigan

if nargin == 1 && streq(r10, 'test'), ir_mri_dce_r1d_test, return, end
if nargin < 2, ir_usage, end

arg.r1 = 4.5; % [mMol^-1 sec^-1] for Gd at 3T, Table 2 of sasaki:05:eea
arg = vararg_pair(arg, varargin);

if numel(r10) ~= 1 && numel(r10) ~= nrow(Ct)
	fail 'Nc mismatch'
end
r1d = ir_mri_dce_r1d_do(r10(:), Ct, arg.r1);

function r1d = ir_mri_dce_r1d_do(r10, Ct, r1)

if numel(r10) ~= 1
	r10 = repmat(r10, [1 ncol(Ct)]);
end
r1d = r10 + r1 * Ct;

% ir_mri_dce_r1d_test()
function ir_mri_dce_r1d_test

TR = 5e-3; % 5 msec
duration = 4; % [min]
ti = linspace(0, duration, 401); % [min]
Ktrans = [3.0 0.6 0.2 0];
kep = [6.0 2.0 1.3 1];
leg = {'rapid', 'moderate', 'slow', 'none'};
Ct = ir_mri_dce_aif1(ti, Ktrans, kep);

f = mri_brainweb_params('grey-matter');
t10 = f.t1 / 1000; % [s]
r10 = 1 / t10; % [1/s]

r1d = ir_mri_dce_r1d(r10, Ct);
if im
	clf
	subplot(121)
	plot(ti, r1d, '.-'), ylabelf 'R1 [1/s]'
	legend(leg{:})
	t1d = 1 ./ r1d;
	subplot(122)
	plot(ti, 1000*t1d, '.-'), ylabelf 'T1 [ms]'
	xlabel 't [min]'
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function R = p2p_rcu(a, n_vals, eps)
% Compute the achievable sum rate (at the symmetrical rate point) from the point-to-point RCU bound
% for the joint source distribution: [a 1/3*(1-a); 1/3*(1-a) 1/3*(1-a)].

% Input:
% a: parameter to define the joint source distribution, should be in range [0.25, 1)
% n_vals: a vector containing the blocklengths
% eps: target error probability
% Output:
% R: the achievable rates for the blocklengths specified in n_vals.

    p = [a 1/3*(1-a); 1/3*(1-a) 1/3*(1-a)];
    H = sum(sum(p.*log(1./p)));  % Source joint entropy

    % Compute log binomial coefficients
    log_b = binomial_coeff(max(n_vals));

    R = zeros(length(n_vals),1);

    for i = 1:length(n_vals)
        n = n_vals(i);
        m = 0:n;  % Number of occurrence of the largest joint probability mass
        % log_Pr stores all the log probabilities in descending order
        log_Pr = m*log(a)+(n-m)*log(1/3*(1-a));

        % [x, y]: initial range for the bisection algorithm
        x = H;
        y = H + 2;
        %%%%%%%%%% Use this part to check if the range is valid for bisection algorithm %%%%%%%%%%
        c = x;
        log_indicator = zeros(1,n+1);
        for j = 0:n
            log_indicator(j+1) = min(0, my_logsumexp(log_b{n+1}(j+1:n+1)+(n-(j:n))*log(3))-n*c);
        end
        err = exp(my_logsumexp(log_Pr + log_b{n+1} + (n-m)*log(3) + log_indicator));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % Start of bisection iteration
        iter = 0;
        while abs(x - y) > 0.00001 && iter <= 25
            c = (x+y)/2;
            %%%%%%%%%%%%%%%% body of iteration %%%%%%%%%%%%%%%%
            log_indicator = zeros(1,n+1);
            for j = 0:n
                log_indicator(j+1) = min(0, my_logsumexp(log_b{n+1}(j+1:n+1)+(n-(j:n))*log(3))-n*c);
            end
            err = exp(my_logsumexp(log_Pr + log_b{n+1} + (n-m)*log(3) + log_indicator));
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if err - eps < 0
                y = c;
            else
                x = c;
            end
            iter = iter + 1;
        end
        R(i) = c;
    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre