text:
function varargout = contourf(varargin)
%CONTOURF   Filled contour plot of a CHEBFUN2.
%   CONTOURF(...) is the same as CONTOUR(...) except that the areas between
%   contours are filled with colors according to the Z-value for each level.
%   Contour regions with data values at or above a given level are filled with
%   the color that maps to the interval.
%
%   NaN's in the Z-data leave white holes with black borders in the contour
%   plot.
%
%   When you use the CONTOURF(Z, V) syntax to specify a vector of contourf
%   levels (V must increase monotonically), contourf regions with Z-values less
%   than V(1) are not filled (are rendered in white). To fill such regions with
%   a color, make V(1) less than or equal to the minimum Z-data value.
%
%   CONTOURF(F, 'NUMPTS', N) computes the contourf lines on a N by N grid. If N
%   is larger than 200 then the contourf lines are drawn with more detail.
%
%   [C, H] = CONTOURF(...) also returns a handle H to a CONTOURGROUP object.
%
% See also CONTOUR.

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

[varargout{1:nargout}] = contourf@separableApprox(varargin{:});

end


text:
function [ a, b] = p04_interval ( )

%*****************************************************************************80
%
%% P04_INTERVAL returns a starting interval for optimization for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, two points defining an interval in which
%    the local minimizer should be sought.
%
  a = 0.0001;
  b = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [h,T,perm] = polardendrogram(Z,varargin)
%POLARDENDROGRAM plots a polar dendrogram plot, taking same options as
%dendrogram and giving same outputs.
% Example:
% X= rand(100,2);
% Y= pdist(X,'cityblock');
% Z= linkage(Y,'average');
% [H,T] = polardendrogram(Z,'colorthreshold','default');

%Plot a normal dendrogram
[h,T,perm] = dendrogram(Z,varargin{:});

%Get x and y ranges
xlim = get(gca,'XLim');
ylim = get(gca,'YLim');
minx = xlim(1);
maxx = xlim(2);
miny = ylim(1);
maxy = ylim(2);
xrange = maxx-minx;
yrange = maxy-miny;

%Reshape into a polar plot
for i=1:size(h)
    xdata = get(h(i),'XData');
    ydata = get(h(i),'YData');
    %Rescale xdata to go from pi/12 to 2pi - pi/12
    xdata = (((xdata-minx)/xrange)*(pi*11/6))+(pi/12);
    %Rescale ydata to go from 1 to 0, cutting off lines
    %which drop below the axis limit
    ydata = max(ydata,miny);
    ydata = 1-((ydata-miny)/yrange);
    %To make horizontal lines look more circular,
    %insert ten points into the middle of the line before
    %polar transform
    newxdata = [xdata(1), linspace(xdata(2),xdata(3),10), xdata(4)];
    newydata = [ydata(1), repmat(ydata(2),1,10), ydata(4)];
    %Transform to polar coordinates
    [xdata,ydata]=pol2cart(newxdata,newydata);
    %Reset line positions to new polar positions
    set(h(i),'XData',xdata);
    set(h(i),'YData',ydata);
end

%Relabel leaves
for i=minx+1:maxx-1
    [x,y]=pol2cart((((i-minx)/xrange)*(pi*11/6))+(pi*1/12),1.1);
    text(x,y,num2str(perm(i)));
end

%Add and label gridlines
hold on
lineh(1) = polar([0,0],[0,1],'-');
lineh(2) = polar(linspace(0,2*pi,50),ones(1,50),':');
lineh(3) = polar(linspace(0,2*pi,50),ones(1,50)*0.75,':');
lineh(4) = polar(linspace(0,2*pi,50),ones(1,50)*0.5,':');
lineh(5) = polar(linspace(0,2*pi,50),ones(1,50)*0.25,':');
set(lineh,'Color',[0.5,0.5,0.5]);
for i=1:4
    [x,y]=pol2cart(0,i/4);
    str = sprintf('%2.1f',((maxy-miny)*(4-i)/4)+miny);
    text(x,y,str,'VerticalAlignment','bottom');
end

%Prettier
set(gca,'XLim',[-1.5,1.5],'YLim',[-1.5,1.5],'Visible','off');
view(3)
axis fill
daspect([1,1,100]);
zoom(2.8);

author: Sable

repo: mcbench-benchmarks


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

% FT_REJECTARTIFACT removes data segments containing artifacts. It returns a
% configuration structure with a modified trial definition which can be used for
% preprocessing of only the clean data.
%
% You should start by detecting the artifacts in the data using the function
% FT_ARTIFACT_xxx where xxx is the type of artifact. Subsequently FT_REJECTARTIFACT
% looks at the detected artifacts and removes them from the trial definition or from
% the data. In case you wish to replace bad parts by nans, you have to specify data
% as an input parameter.
%
% Use as
%   cfg = ft_rejectartifact(cfg)
% with the cfg as obtained from FT_DEFINETRIAL, or as
%   data = ft_rejectartifact(cfg, data)
% with the data as obtained from FT_PREPROCESSING
%
% The following configuration options are supported:
%   cfg.artfctdef.reject          = 'none', 'partial', 'complete', 'nan', or 'value' (default = 'complete')
%   cfg.artfctdef.minaccepttim    = when using partial rejection, minimum length
%                                   in seconds of remaining trial (default = 0.1)
%   cfg.artfctdef.crittoilim      = when using complete rejection, reject trial only when artifacts occur within
%                                   this time window (default = whole trial). This only works with in-memory data,
%                                   since trial time axes are unknown for data on disk.
%   cfg.artfctdef.feedback        = 'yes' or 'no' (default = 'no')
%   cfg.artfctdef.invert          = 'yes' or 'no' (default = 'no')
%   cfg.artfctdef.value           = scalar value to replace the data in the artifact segments (default = nan)
%   cfg.artfctdef.eog.artifact    = Nx2 matrix with artifact segments, this is added to the cfg by using FT_ARTIFACT_EOG
%   cfg.artfctdef.jump.artifact   = Nx2 matrix with artifact segments, this is added to the cfg by using FT_ARTIFACT_JUMP
%   cfg.artfctdef.muscle.artifact = Nx2 matrix with artifact segments, this is added to the cfg by using FT_ARTIFACT_MUSCLE
%   cfg.artfctdef.zvalue.artifact = Nx2 matrix with artifact segments, this is added to the cfg by using FT_ARTIFACT_ZVALUE
%   cfg.artfctdef.visual.artifact = Nx2 matrix with artifact segments, this is added to the cfg by using FT_DATABROWSER
%   cfg.artfctdef.xxx.artifact    = Nx2 matrix with artifact segments, this could be added by your own artifact detection function
%
% A trial that contains an artifact can be rejected completely or partially. In case
% of partial rejection, a minimum length of the resulting sub-trials can be
% specified using minaccepttim.
%
% Output:
%   If cfg is used as the only input parameter, the output is a cfg structure with an updated trl.
%   If cfg and data are both input parameters, the output is an updated raw data structure with only the clean data segments.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_ARTIFACT_EOG, FT_ARTIFACT_MUSCLE, FT_ARTIFACT_JUMP, FT_ARTIFACT_THRESHOLD,
% FT_ARTIFACT_CLIP, FT_ARTIFACT_ECG, FT_DATABROWSER, FT_REJECTVISUAL

% Undocumented local options:
% cfg.headerfile
% cfg.reject
% cfg.trl
% cfg.trlold
% cfg.version
%
% These old configuration options are still supported
% cfg.rejectmuscle      = 'no' or 'yes'
% cfg.rejecteog         = 'no' or 'yes'
% cfg.rejectjump        = 'no' or 'yes'
% cfg.rejectfile        = string with filename
% cfg.artfctdef.writerej = filename of rejection file
% cfg.artfctdef.type    = cell-array with strings, e.g. {'eog', 'muscle' 'jump'}

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

% FIXME this function contains a lot of lines of code that pertain to backward
% compatibility support that dates back to 2004/2005. It would be good to strip
% that code and only keep the relevant parts

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

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

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

% ft_checkdata is done further down

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');

% set the defaults
cfg.artfctdef              = ft_getopt(cfg, 'artfctdef');
cfg.artfctdef.type         = ft_getopt(cfg.artfctdef, 'type', {});
cfg.artfctdef.reject       = ft_getopt(cfg.artfctdef, 'reject', 'complete');
cfg.artfctdef.minaccepttim = ft_getopt(cfg.artfctdef, 'minaccepttim', 0.1);
cfg.artfctdef.crittoilim   = ft_getopt(cfg.artfctdef, 'crittoilim', []);
cfg.artfctdef.feedback     = ft_getopt(cfg.artfctdef, 'feedback', 'no');
cfg.artfctdef.invert       = ft_getopt(cfg.artfctdef, 'invert', 'no');
cfg.artfctdef.value        = ft_getopt(cfg.artfctdef, 'value', nan);

% convert from old-style to new-style configuration
if isfield(cfg,'reject')
  ft_warning('converting from old-style artifact configuration to new-style');
  cfg.artfctdef.reject = cfg.reject;
  cfg = rmfield(cfg, 'reject');
end

% convert from old-style to new-style configuration
if isfield(cfg.artfctdef,'common')
  ft_warning('converting from old-style artifact configuration to new-style');
  if isfield(cfg.artfctdef.common,'minaccepttim')
    cfg.artfctdef.minaccepttim = cfg.artfctdef.common.minaccepttim;
    cfg.artfctdef = rmfield(cfg.artfctdef, 'common');
  end
end

% ensure that it is a cell-array
if ischar(cfg.artfctdef.type)
  cfg.artfctdef.type = {cfg.artfctdef.type};
end

% support the rejectXXX cfg settings for backward compatibility
if isfield(cfg, 'rejectmuscle')
  dum = find(strcmp('muscle', cfg.artfctdef.type));
  if strcmp(cfg.rejectmuscle,'yes') && isempty(dum)
    % this overrules the other setting, add it to the type-list
    cfg.artfctdef.type = cat(1, {'muscle'}, cfg.artfctdef.type(:));
  elseif strcmp(cfg.rejectmuscle,'no') && ~isempty(dum)
    % this overrules the other setting, remove it from the type-list
    cfg.artfctdef.type(dum) = [];
  end
  cfg = rmfield(cfg, 'rejectmuscle');
end

% support the rejectXXX cfg settings for backward compatibility
if isfield(cfg, 'rejecteog')
  dum = find(strcmp('eog', cfg.artfctdef.type));
  if strcmp(cfg.rejecteog,'yes') && isempty(dum)
    % this overrules the other setting, add it to the type-list
    cfg.artfctdef.type = cat(1, {'eog'}, cfg.artfctdef.type(:));
  elseif strcmp(cfg.rejecteog,'no') && ~isempty(dum)
    % this overrules the other setting, remove it from the type-list
    cfg.artfctdef.type(dum) = [];
  end
  cfg = rmfield(cfg, 'rejecteog');
end

% support the rejectXXX cfg settings for backward compatibility
if isfield(cfg, 'rejectjump')
  dum = find(strcmp('jump', cfg.artfctdef.type));
  if strcmp(cfg.rejectjump,'yes') && isempty(dum)
    % this overrules the other setting, add it to the type-list
    cfg.artfctdef.type = cat(1, {'jump'}, cfg.artfctdef.type(:));
  elseif strcmp(cfg.rejectjump,'no') && ~isempty(dum)
    % this overrules the other setting, remove it from the type-list
    cfg.artfctdef.type(dum) = [];
  end
  cfg = rmfield(cfg, 'rejectjump');
end

% support the rejectXXX cfg settings for backward compatibility
if isfield(cfg, 'rejectfile')
  % this is slightly different to the ones above, since rejectfile is either 'no' or contains the filename
  dum = find(strcmp('file', cfg.artfctdef.type));
  if ~strcmp(cfg.rejectfile,'no') && isempty(dum)
    % this overrules the other setting, add it to the type-list
    cfg.artfctdef.type = cat(1, {'file'}, cfg.artfctdef.type(:));
  elseif strcmp(cfg.rejectfile,'no') && ~isempty(dum)
    % this overrules the other setting, remove it from the type-list
    cfg.artfctdef.type(dum) = [];
  end
end

% the data can be specified as input variable or through cfg.inputfile
hasdata = exist('data', 'var');

if hasdata
  % check if the input data is valid for this function
  data = ft_checkdata(data, 'hassampleinfo', 'yes');

  if isfield(data, 'sampleinfo')
    trl = zeros(numel(data.trial), 3);
    trl(:,[1 2]) = data.sampleinfo;

    % recreate offset vector (artifact functions depend on this)
    % TODO: the artifact rejection stuff should be rewritten to avoid
    % needing this workaround
    for ntrl = 1:numel(data.trial)
      trl(ntrl,3) = time2offset(data.time{ntrl}, data.fsample);
    end

    if isfield(data, 'trialinfo')
      if istable(data.trialinfo)
        % convert table into normal array, keep the column labels
        VariableNames = data.trialinfo.Properties.VariableNames;
        data.trialinfo = table2array(data.trialinfo);
      end
      trl = [trl data.trialinfo];
    end
  else
    trl = [];
  end

elseif isfield(cfg, 'trl')
  trl = cfg.trl;
end

% ensure the crittoilim input argument is valid
if ~isempty(cfg.artfctdef.crittoilim)

  if (size(cfg.artfctdef.crittoilim,2) ~= 2 ...
      || (size(cfg.artfctdef.crittoilim,1) ~= size(trl,1) ...
      && size(cfg.artfctdef.crittoilim,1) ~= 1))
    ft_error('if specified, cfg.artfctdef.crittoilim should be a 1x2 or Nx2 vector');
  end

  % if specified as 1x2 vector, expand into Nx2
  if (size(cfg.artfctdef.crittoilim,1) == 1)
    cfg.artfctdef.crittoilim = repmat(cfg.artfctdef.crittoilim,size(trl,1),1);
  end

  checkCritToi = 1; % flag for convenience
else
  checkCritToi = 0;
end

% ensure that there are trials that can be scanned for artifacts and/or rejected
if isempty(trl)
  ft_error('no trials were selected, cannot perform artifact detection/rejection');
end

% prevent double occurences of artifact types, ensure that the order remains the same
[dum, i] = unique(cfg.artfctdef.type);
cfg.artfctdef.type = cfg.artfctdef.type(sort(i));
% ensure that it is a row vector
cfg.artfctdef.type = cfg.artfctdef.type(:)';

% If bad parts are to be filled with nans, make sure data is available
if ~hasdata && (strcmp(cfg.artfctdef.reject, 'nan') || strcmp(cfg.artfctdef.reject, 'value'))
  ft_error('If bad parts are to be filled with nans or another value, the input data has to be specified');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% call the appropriate function for each of the artifact types
% this will produce a Nx2 matrix with the begin and end sample of artifacts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for type=1:length(cfg.artfctdef.type)
  funhandle = ft_getuserfun(cfg.artfctdef.type{type}, 'artifact');
  fprintf('evaluating %s\n', func2str(funhandle));
  % each call to artifact_xxx adds cfg.artfctdef.xxx.artifact
  if hasdata
    cfg = feval(funhandle, cfg, data);
  else
    cfg = feval(funhandle, cfg);
  end
end

% collect the artifacts that have been detected from cfg.artfctdef.xxx.artifact
dum = fieldnames(cfg.artfctdef);
sel = false(size(dum));
artifact = {};
for i=1:length(dum)
  sel(i) = issubfield(cfg.artfctdef, dum{i}) && issubfield(cfg.artfctdef, [dum{i} '.artifact']);
  if sel(i)
    artifact{end+1} = getsubfield(cfg.artfctdef, [dum{i} '.artifact']);
    num = size(artifact{end}, 1);
    if isempty(num)
      num = 0;
    end
    fprintf('detected %3d %s artifacts\n', num, dum{i});
  end
end

% update the configuration to reflect the artifacts types that were scanned
cfg.artfctdef.type = dum(sel);

% combine all trials into a single boolean vector
trialall = convert_event(trl(:,1:2), 'boolvec');

% combine all artifacts into a single vector
rejectall = zeros(1, max(trl(:,2)));
for i=1:length(cfg.artfctdef.type)
  dum = artifact{i};
  for j=1:size(dum,1)
    rejectall(dum(j,1):dum(j,2)) = i;  % the artifact type is coded here
  end
end

% ensure that both vectors are the same length
if length(trialall)>length(rejectall)
  rejectall(length(trialall)) = 0;
elseif length(trialall)<length(rejectall)
  trialall(length(rejectall)) = 0;
end

% make header, needed only for sampling frequency
if hasdata
  hdr = ft_fetch_header(data);
else
  if isfield(cfg, 'headerformat')
    hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat);
  else
    hdr = ft_read_header(cfg.headerfile);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% give visual feedback on the trial definition and the artifacts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(cfg.artfctdef.feedback, 'yes')
  COLOR = 'grcmykgrcmykgrcmykgrcmyk';
  % use the trial definition present in the local configuration
  trl  = trl;
  % compute the time axis that corresponds with each trial
  time = {};
  for i=1:size(trl,1)
    time{i} = offset2time(trl(i,3), hdr.Fs, trl(i,2)-trl(i,1)+1);
  end

  figure
  title('linear display of the continuous data')
  xlabel('sample number');
  hold on
  x = 1:length(trialall);
  y = ones(size(trialall)) * 0.53;
  y(~trialall) = nan;
  plot(x, y, 'b');
  for i=1:length(cfg.artfctdef.type)
    x = 1:length(rejectall);
    y = ones(size(rejectall)) * (0.5 - i*0.03);
    y(~(rejectall==i)) = nan;
    plot(x, y, COLOR(i));
  end
  set(gca, 'ytick', []);
  dum = (trialall~=0) | (rejectall~=0);
  [dum2, dum] = find(dum);
  smpbeg = dum(1)   - hdr.Fs;
  smpend = dum(end) + hdr.Fs;
  axis([smpbeg smpend 0 1]);
  legend({'defined trials', cfg.artfctdef.type{:}});

  figure
  title('individual trials after alignment')
  xlabel('time (s)');
  ylabel('trial number');
  hold on
  for i=1:size(trl,1)
    x = [time{i}(1) time{i}(end)];
    y = [i i];
    plot(x, y, 'b')
    for j=1:length(cfg.artfctdef.type)
      x = time{i};
      y = rejectall(trl(i,1):trl(i,2));
      y(y~=j) = nan;
      y(y==j) = i + j*0.03;
      plot(x, y, COLOR(j));
    end
    timebeg(i) = time{i}(1);
    timeend(i) = time{i}(end);
  end
  axis([min(timebeg)-0.1 max(timeend)+0.1 0.5 size(trl,1)+0.5]);
  axis ij
  legend({'defined trials', cfg.artfctdef.type{:}});
end % feedback

% convert to logical, this is required for the subsequent code
rejectall = (rejectall~=0);

% invert the artifact selection
if istrue(cfg.artfctdef.invert)
  fprintf('inverting selection of clean/artifactual data\n');
  rejectall = ~rejectall;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write the rejection to an EEP format file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfield(cfg.artfctdef, 'writerej') && ~isempty(cfg.artfctdef.writerej)
  fid = fopen(cfg.artfctdef.writerej, 'w');
  if fid<0
    ft_error('could not open rejection file for writing');
  else
    % determine the begin and end of each rejected period (in samples)
    rejectonset = find(diff([0 rejectall])== 1);
    rejectofset = find(diff([rejectall 0])==-1);
    % determine the begin and end of each rejected period (in seconds)
    rejectonset = (rejectonset-1)/hdr.Fs;
    rejectofset = (rejectofset-1)/hdr.Fs;
    for rejlop=1:length(rejectonset)
      fprintf(fid, '%f-%f\n', rejectonset(rejlop), rejectofset(rejlop));
    end
    fclose(fid);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% remove the trials that (partially) coincide with a rejection mark
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if any(strcmp(cfg.artfctdef.reject, {'partial', 'complete', 'nan', 'value'}))
  trialok = [];

  count_complete_reject = 0;
  count_partial_reject  = 0;
  count_nan             = 0;
  count_value           = 0;
  count_outsidecrit     = 0;

  trlCompletelyRemovedInd = [];
  trlPartiallyRemovedInd  = [];

  for trial=1:size(trl,1)
    % cut out the part of the rejection-axis corresponding with this trial
    rejecttrial = rejectall(trl(trial,1):trl(trial,2));

    if all(not(rejecttrial))
      % the whole trial is good
      trialok = [trialok; trl(trial,:)];

    elseif all(rejecttrial) && strcmp(cfg.artfctdef.reject, 'nan')
      % the whole trial is bad, but it is requested to be replaced with nans
      data.trial{trial}(:,rejecttrial) = nan;
      count_nan = count_nan + 1;
      trialok = [trialok; trl(trial,:)]; % Mark the trial as good as nothing will be removed

    elseif all(rejecttrial) && strcmp(cfg.artfctdef.reject, 'value')
      % the whole trial is bad, but it is requested to be replaced with a specific value
      data.trial{trial}(:,rejecttrial) = cfg.artfctdef.value;
      count_value = count_value + 1;
      trialok = [trialok; trl(trial,:)]; % Mark the trial as good as nothing will be removed

    elseif all(rejecttrial)
      % the whole trial is bad
      count_complete_reject = count_complete_reject + 1;
      trlCompletelyRemovedInd = [trlCompletelyRemovedInd trial];

    elseif any(rejecttrial) && strcmp(cfg.artfctdef.reject, 'complete')
      % some part of the trial is bad, check if within crittoilim?
      if (checkCritToi)
        critInd = (data.time{trial} >= cfg.artfctdef.crittoilim(trial,1) ...
          & data.time{trial} <= cfg.artfctdef.crittoilim(trial,2));
        if (any(critInd & rejecttrial))
          count_complete_reject = count_complete_reject + 1;
          trlCompletelyRemovedInd = [trlCompletelyRemovedInd trial];
          continue;
        else
          trialok = [trialok; trl(trial,:)];
          count_outsidecrit = count_outsidecrit + 1;
        end
      else % no crittoilim checking required
        count_complete_reject = count_complete_reject + 1;
        trlCompletelyRemovedInd = [trlCompletelyRemovedInd trial];
        continue;
      end

    elseif any(rejecttrial) && strcmp(cfg.artfctdef.reject, 'partial')
      % some part of the trial is bad, reject only the bad part
      trialnew = [];
      rejecttrial = [0 not(rejecttrial) 0];
      % the first column is the begin sample, the second the end sample and the third is the offset
      trialnew(:,1) = find(diff(rejecttrial(:))== 1) + trl(trial,1) - 1;
      trialnew(:,2) = find(diff(rejecttrial(:))==-1) + trl(trial,1) - 2;
      trialnew(:,3) = find(diff(rejecttrial(:))== 1) - 1 + trl(trial,3);
      % some people use additional columns in the trl matrix to store trigger values and/or reaction times
      % these should remain linked to the original trial, i.e. they should be copied for each new fragment
      for i=4:size(trl,2)
        trialnew(:,i) = trl(trial,i);
      end
      minacceptnumsmp = round(cfg.artfctdef.minaccepttim .* hdr.Fs);
      trialnew((trialnew(:,2)-trialnew(:,1))<minacceptnumsmp,:) = [];
      count_partial_reject = count_partial_reject + 1;
      trialok = [trialok; trialnew];
      trlPartiallyRemovedInd = [trlPartiallyRemovedInd trial];

    elseif any(rejecttrial) && strcmp(cfg.artfctdef.reject, 'nan')
      % Some part of the trial is bad, replace bad part with nans
      data.trial{trial}(:,rejecttrial) = nan;
      count_nan = count_nan + 1;
      trialok = [trialok; trl(trial,:)]; % Mark the trial as good as nothing will be removed

    elseif any(rejecttrial) && strcmp(cfg.artfctdef.reject, 'value')
      % Some part of the trial is bad, replace bad part with specified value
      data.trial{trial}(:,rejecttrial) = cfg.artfctdef.value;
      count_value = count_value + 1;
      trialok = [trialok; trl(trial,:)]; % Mark the trial as good as nothing will be removed

    end
  end % for each trial

  fprintf('rejected  %3d trials completely\n', count_complete_reject);
  fprintf('rejected  %3d trials partially\n', count_partial_reject);
  fprintf('filled parts of  %3d trials with nans\n', count_nan);
  fprintf('filled parts of  %3d trials with the specified value\n', count_value);
  if (checkCritToi)
    fprintf('retained  %3d trials with artifacts outside critical window\n', count_outsidecrit);
  end
  fprintf('resulting %3d trials\n', size(trialok,1));
  cfg.trlold = trl;      % return the original trial definition in the configuration
  cfg.trl    = trialok;  % return the cleaned trial definition in the configuration

  if strcmp(cfg.artfctdef.feedback, 'yes')
    fprintf('the following trials were completely removed: ');
    for k = trlCompletelyRemovedInd
      fprintf('%d ', k);
    end
    fprintf('\nthe following trials were partially removed: ');
    for k = trlPartiallyRemovedInd
      fprintf('%d ', k);
    end
    fprintf('\n');
  end

else
  fprintf('not rejecting any data, only marking the artifacts\n');
end

if isempty(cfg.trl)
  ft_error('No trials left after artifact rejection.')
else
  if hasdata && ~strcmp(cfg.artfctdef.reject, 'nan') % Skip this step to avoid removing parts that should be filled with nans
    % apply the updated trial definition on the data
    tmpcfg      = keepfields(cfg, {'trl', 'showcallinfo'});
    data        = removefields(data, {'trialinfo'});
    data        = ft_redefinetrial(tmpcfg, data);
    % restore the provenance information
    [cfg, data] = rollback_provenance(cfg, data);
  end
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble provenance
if hasdata
  ft_postamble previous data
  ft_postamble history data
  ft_postamble savevar data
  % return the data, the output variable is called cfg instead of data
  cfg = data;
end

author: spm


text:
function gX = multiKernGradX(kern, x, x2)

% MULTIKERNGRADX Gradient of MULTI kernel with respect to a point x.
% FORMAT
% DESC computes the gradient of the multiple output block
% kernel with respect to the input positions.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x : locations against which gradients are being computed.
% RETURN gX : the returned gradients. The gradients are returned in
% a matrix which is numData x numInputs x numData. Where numData is
% the number of data points and numInputs is the number of input
% dimensions in X.
%
% FORMAT
% DESC computes the gradident of the multiple output block
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO : multiKernParamInit, kernGradX, multiKernDiagGradX
%
% COPYRIGHT : Mauricio A. Alvarez and Neil D. Lawrence, 2008

% KERN

if iscell(x)
    if nargin > 3 && ~iscell(x2)
        error('Time course information is not matched in Cell format!');
    end
    % Collate arguments.
    for i=1:kern.numBlocks
        arg{i}{1} = x{i};
        if nargin > 2
            arg{i}{2} = x2{i};
        else
            arg{i}{2} = arg{i}{1};
        end
    end
    for i = 1:kern.numBlocks
        if nargin > 2
            gX = multiKernGradientBlockX(kern, arg{i}{:}, i, i);
        else
            gX = multiKernGradientBlockX(kern, arg{i}{1}, i, i);
        end
        for j = 1:i-1
            if ~isempty(kern.block{i}.cross{j})
                gX = multiKernGradientBlockX(kern, arg{i}{1}, arg{j}{2}, i, j);
            end
        end
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function output = write_image(pfObj, imageFormat, filename, height, width, scale)

% Function to write plotly figures to a supported image format, which are the following: "png", "jpg", "jpeg", "webp", "svg", "pdf", "eps", "json"

debug=0;
if nargin < 2
    imageFormat="png";
    filename='figure.png';
    height=pfObj.layout.height;
    width=pfObj.layout.width;
    scale=1;
elseif nargin < 3
    filename=['figure.',char(imageFormat)];
    height=pfObj.layout.height;
    width=pfObj.layout.width;
    scale=1;
elseif nargin < 4
    height=pfObj.layout.height;
    width=pfObj.layout.width;
    scale=1;
elseif nargin < 5
    width=pfObj.layout.width;
    scale=1;
elseif nargin < 6
    scale=1;
end

if strcmpi(imageFormat,'jpg')
    imageFormat = "jpeg";
end

wd=fileparts(fileparts(mfilename('fullpath')));
output=[];

if ~isa(pfObj,'plotlyfig')
    fprintf('\nError: Input is not a plotlyfig object.\n\n');
    return
end

if isunix()
    kExec = string(fullfile(wd,'kaleido','kaleido'));
    cc="cat";
else
    kExec = string(fullfile(wd,'kaleido','kaleido.cmd'));
    cc="type";
end
plyJsLoc = string(fullfile(wd,'kaleido','plotly-latest.min.js'));

if ~isfile(kExec) || ~isfile(plyJsLoc)
    status=getKaleido();
else
    status=1;
end

if status == 0
    return
end

mjLoc = replace(string(fullfile(wd,'kaleido','etc','mathjax','MathJax.js')),'\','/');
scope="plotly";

% Prepare input plotly object for Kaleido
q=struct();
q.data.data = pfObj.data;
q.data.layout = pfObj.layout;
q.data.layout = rmfield(q.data.layout,'height');
q.data.layout = rmfield(q.data.layout,'width');
q.format = string(imageFormat);
q.height = height;
q.scale = scale;
q.width = width;

pfJson = native2unicode(jsonencode(q),'UTF-8');
tFile = string(fullfile(wd,'kaleido','temp.txt'));
f=fopen(tFile,'w');
fprintf(f,"%s",pfJson);
fclose(f);

cmd=[cc," ",tFile," | ",kExec," ",scope," --plotlyjs='",plyJsLoc,"' ","--mathjax='file:///",mjLoc,"' --no-sandbox --disable-gpu --allow-file-access-from-files --disable-breakpad --disable-dev-shm-usage"];

if debug
    inputCmd=char(join(cmd,''));
    fprintf('\nDebug info:\n%s\n\n',inputCmd);
end

[code,out]=system(char(join(cmd,'')));
if debug
    disp(out);
end

if code ~= 0
    fprintf('\nFatal: Failed to run Kaleido.\n\n');
    return;
else
    a=string(split(out,newline));
    if a(end)==""
        a(end)=[];
    end
    output = jsondecode(a(end));
end

if output.code ~= 0
    fprintf('\nError: %s\n',output.message);
else
    out=unicode2native(output.result,'UTF-8');
    out=base64decode(out);
    f=fopen(char(filename),'wb');
    fwrite(f,out);
    fclose(f);
end

author: plotly


text:
function newTargets = getImgOverPlanes(ops1)
nFOV = size(ops1,2);
Ly = ops1{1}.Ly;
Lx = ops1{1}.Lx;
planesToInterpolate = ops1{1}.planesToInterpolate;
ds = zeros(length(planesToInterpolate), 2, nFOV);
% planesToInterpolate should be indices of non-flyback planes
% (only those should be aligned)
for i = 2:length(planesToInterpolate) % align consecutive planes
    pl = planesToInterpolate(i);
    for j = 1:nFOV
        ds(i,:,j) = registration_offsets(ops1{pl,j}.mimg, ops1{pl-1,j}, 0);
    end
end
ds = cumsum(ds,1);
ds = bsxfun(@minus, ds, mean(ds,1));
images = zeros(Ly, Lx, length(planesToInterpolate), nFOV);
for i = 1:length(planesToInterpolate)
    for j = 1:nFOV
        images(:,:,i,j) = ops1{planesToInterpolate(i),j}.mimg;
    end
end
ds = reshape(permute(ds, [1 3 2]), [], 2);
images = reshape(images, Ly, Lx, []);
newTargets = register_movie(images, ops1{1}, ds);
newTargets = reshape(newTargets, Ly, Lx, length(planesToInterpolate), nFOV);

author: cortex-lab


text:
function c8_uniform_01_test ( )

%*****************************************************************************80
%
%% C8_UNIFORM_01_TEST tests C8_UNIFORM_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2015
%
%  Author:
%
%    John Burkardt
%
  seed = 123456678;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_UNIFORM_01_TEST\n' );
  fprintf ( 1, '  C8_UNIFORM_01 returns a uniformly random "unit" C8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       C1=C8_UNIFORM_01\n' );
  fprintf ( 1, '     ---------------------\n' );
  fprintf ( 1, '\n' );

  for test = 1 : 10

    [ c1, seed ] = c8_uniform_01 ( seed );

    fprintf ( 1, '  (%12f  %12f)\n', real ( c1 ), imag ( c1 ));

  end

  return
end


text:
function [Grad,t]=StdBit(tStart,tEnd,gAmp)

%create a minimum signal pulse starting from tStart
%tStart signal pulse start time
%tEnd   signal pulse end time

[Grad,t]=StdRect(tStart,tEnd,gAmp,3); % signal bit
Grad(1)=0;
Grad(end)=0;

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function x = simplex_grid_index_to_point ( m, n, ng, g, v )

%*****************************************************************************80
%
%% SIMPLEX_GRID_INDEX_TO_POINT returns  points corresponding to simplex indices.
%
%  Discussion:
%
%    The M-dimensional simplex is defined by M+1 vertices.
%
%    Given a regular grid that uses N subintervals along the edge between
%    each pair of vertices, a simplex grid index G is a set of M+1 values
%    each between 0 and N, and summing to N.
%
%    This function determines the coordinates X of the point corresponding
%    to the index G.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of subintervals.
%
%    Input, integer NG, the number of grid indices to be converted.
%
%    Input, integer G(M+1,NG), the grid indices of 1 or more points.
%
%    Input, real V(M,M+1), the coordinates of the vertices of the simplex.
%
%    Output, real X(M,NG), the coordinates of one or more points.
%
  x = v(1:m,1:m+1) * g(1:m+1,1:ng) / n;

  return
end


text:
function [E,H] = Maxwell1D(E,H,eps,mu,FinalTime);

% function [E,H] = Maxwell1D(E,H,eps,mu,FinalTime)
% Purpose  : Integrate 1D Maxwell's until FinalTime starting with conditions (E(t=0),H(t=0))
%            and materials (eps,mu).

Globals1D;
time = 0;

% Runge-Kutta residual storage
resE = zeros(Np,K); resH = zeros(Np,K);

% compute time step size
xmin = min(abs(x(1,:)-x(2,:)));
CFL=1.0;  dt = CFL*xmin;
Nsteps = ceil(FinalTime/dt); dt = FinalTime/Nsteps

% outer time step loop
for tstep=1:Nsteps
   for INTRK = 1:5
      [rhsE, rhsH] = MaxwellRHS1D(E,H,eps,mu);

      resE = rk4a(INTRK)*resE + dt*rhsE;
      resH = rk4a(INTRK)*resH + dt*rhsH;

      E = E+rk4b(INTRK)*resE;
      H = H+rk4b(INTRK)*resH;
   end
   % Increment time
   time = time+dt;
end
return

repo: Aero-matlab


text:
errData = data2plot;

clear ph
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
symbols = {'o', 'd', 's'};
xx = 1:0.1:11;
yy = zeros(3,length(xx));
for iM = 1:3
    if iM == 2
        xm = zeros(1,N_ETA);
        for k = 1:N_ETA
            TF = isnan(errData(1:Nr,k,iM));
            xm(k) = median(errData(TF==0,k,iM),1);
        end
    else
        xm = median(errData(1:Nr,:,iM),1);
    end
    yy(iM,:) = spline(1:11,xm,xx);
end
figure; hold on
if shaded_region == 1
    fillh1 = fill([0.1 11.9 11.9 0.1], [10 10 ylimregion.*ones(1,2)],0.9*ones(1,3));
    fillh1.EdgeColor = 0.9*ones(1,3); fillh1.FaceAlpha = 0.5;
end
ph(1)=plot(xx,yy(1,:),'-r','LineWidth',1); hold on,
ph(2)=plot(xx(1:end),yy(2,1:end),'-g','LineWidth',1);
ph(3)=plot(xx(1:end),yy(3,1:end),'-','Color',[0.7,0.7,1],'LineWidth',1);
for iM = 1:3
    boxplot(errData(:,:,iM),'PlotStyle','compact','Colors',ccolors(iM,:), 'Symbol', symbols{iM}, ...
        'Widths',0.1); hold on
    delete(findobj(gca,'Type','text'))
end
set(gca,'xtick',[1:2:11],'xticklabel',num2str(eta_vec(1:2:end)'),'Position',[0.2 0.22 0.75 0.73])
xlim([0 12]), ylim(yaxlim)

if exist('LOG_SCALE')
    if LOG_SCALE == 1
        set(gca,'YScale','log')
        ylim([10^-3 yaxlim(2)])
        set(gca,'ytick',[10.^[-3:1:1]])
    end
end
xt = xlabel('Eta'); xt.Position = [115 -20 -0.1];
ylabel(ytext)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_noise_',NARXtraining,'_',PostName,'_noleg.eps']);

l1 = legend(ph,'DMDc','SINDYc','NARX');
if exist('LegLoc')==1
    set(l1,'Location',LegLoc)
else
    set(l1,'Location','NorthWest')%'NorthWest')SouthEast
    l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
end
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_PREDPERF_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_noise_',NARXtraining,'_',PostName,'.eps']);

author: eurika-kaiser

repo: SINDY-MPC


text:
function colorsPlot(red, green, blue, justGreen)
%CREATEFIGURE(red,green,blue,justGreen)
%  red:  image cdata
%  green:  image cdata
%  blue:  image cdata
%  justGreen:  image cdata

%  Auto-generated by MATLAB on 12-Nov-2007 09:55:06

% Copyright 2011 The MathWorks, Inc.

% Create figure
figure1 = figure('Colormap',...
    [0 0 0;0.01587 0.01587 0.01587;0.03175 0.03175 0.03175;0.04762 0.04762 0.04762;0.06349 0.06349 0.06349;0.07937 0.07937 0.07937;0.09524 0.09524 0.09524;0.1111 0.1111 0.1111;0.127 0.127 0.127;0.1429 0.1429 0.1429;0.1587 0.1587 0.1587;0.1746 0.1746 0.1746;0.1905 0.1905 0.1905;0.2063 0.2063 0.2063;0.2222 0.2222 0.2222;0.2381 0.2381 0.2381;0.254 0.254 0.254;0.2698 0.2698 0.2698;0.2857 0.2857 0.2857;0.3016 0.3016 0.3016;0.3175 0.3175 0.3175;0.3333 0.3333 0.3333;0.3492 0.3492 0.3492;0.3651 0.3651 0.3651;0.381 0.381 0.381;0.3968 0.3968 0.3968;0.4127 0.4127 0.4127;0.4286 0.4286 0.4286;0.4444 0.4444 0.4444;0.4603 0.4603 0.4603;0.4762 0.4762 0.4762;0.4921 0.4921 0.4921;0.5079 0.5079 0.5079;0.5238 0.5238 0.5238;0.5397 0.5397 0.5397;0.5556 0.5556 0.5556;0.5714 0.5714 0.5714;0.5873 0.5873 0.5873;0.6032 0.6032 0.6032;0.619 0.619 0.619;0.6349 0.6349 0.6349;0.6508 0.6508 0.6508;0.6667 0.6667 0.6667;0.6825 0.6825 0.6825;0.6984 0.6984 0.6984;0.7143 0.7143 0.7143;0.7302 0.7302 0.7302;0.746 0.746 0.746;0.7619 0.7619 0.7619;0.7778 0.7778 0.7778;0.7937 0.7937 0.7937;0.8095 0.8095 0.8095;0.8254 0.8254 0.8254;0.8413 0.8413 0.8413;0.8571 0.8571 0.8571;0.873 0.873 0.873;0.8889 0.8889 0.8889;0.9048 0.9048 0.9048;0.9206 0.9206 0.9206;0.9365 0.9365 0.9365;0.9524 0.9524 0.9524;0.9683 0.9683 0.9683;0.9841 0.9841 0.9841;1 1 1]);

% Create subplot
subplot1 = subplot(2,2,1,'Parent',figure1,'YDir','reverse','Layer','top');
% Uncomment the following line to preserve the X-limits of the axes
 xlim([0.5 320.5]);
% Uncomment the following line to preserve the Y-limits of the axes
 ylim([0.5 240.5]);
box('on');
hold('all');

% Create image
image(red,'Parent',subplot1,'CDataMapping','scaled');

% Create title
title('Red');

% Create subplot
subplot2 = subplot(2,2,2,'Parent',figure1,'YDir','reverse','Layer','top');
% Uncomment the following line to preserve the X-limits of the axes
 xlim([0.5 320.5]);
% Uncomment the following line to preserve the Y-limits of the axes
 ylim([0.5 240.5]);
box('on');
hold('all');

% Create title
title('Green');

% Create image
image(green,'Parent',subplot2,'CDataMapping','scaled');

% Create subplot
subplot3 = subplot(2,2,3,'Parent',figure1,'YDir','reverse','Layer','top');
% Uncomment the following line to preserve the X-limits of the axes
 xlim([0.5 320.5]);
% Uncomment the following line to preserve the Y-limits of the axes
 ylim([0.5 240.5]);
box('on');
hold('all');

% Create title
title('Blue');

% Create image
image(blue,'Parent',subplot3,'CDataMapping','scaled');

% Create subplot
subplot4 = subplot(2,2,4,'Parent',figure1,'YDir','reverse','Layer','top');
% Uncomment the following line to preserve the X-limits of the axes
 xlim([0.5 320.5]);
% Uncomment the following line to preserve the Y-limits of the axes
 ylim([0.5 240.5]);
box('on');
hold('all');

% Create title
title('Just Green');

% Create image
image(justGreen,'Parent',subplot4,'CDataMapping','scaled');


text:
function test3
%TEST3 test sparse on int8, int16, and logical
% Example:
%   test3
% See also cholmod_test

% Copyright 2006-2007, Timothy A. Davis, University of Florida

fprintf ('=================================================================\n');
fprintf ('test3: test sparse on int8, int16, and logical\n') ;

clear all
c =  ['a' 'b' 0 'd']							    %#ok
sparse(c)
sparse2(c)
sparse(c')
sparse2(c')
whos
nzmax(ans)  %#ok

try % this will fail
    sparse(int8(c))
catch
    fprintf ('sparse(int8(c)) fails in MATLAB\n') ;
end
sparse2(int8(c))

sparse2 (int16(c))
whos
s = logical(rand(4) > .5)						    %#ok
sparse (s)
whos
sparse2(s)
whos

x = rand(4)								    %#ok
sparse (x > .5)								    %#ok
whos
sparse2 (x > .5)							    %#ok
whos

fprintf ('test3 passed\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ a, b ] = p26_lim ( dim_num )

%*****************************************************************************80
%
%% P26_LIM returns the integration limits for problem 26.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end


text:
function y = perform_windowed_fourier_transform(x,w,q,n, options)

% perform_windowed_fourier_transform - compute a local Fourier transform
%
% Forward transform:
%   MF = perform_windowed_fourier_transform(M,w,q,n, options);
% Backward transform:
%   M  = perform_windowed_fourier_transform(MF,w,q,n, options);
%
%   w is the width of the window used to perform local computation.
%   q is the spacing betwen each window.
%
%   MF(:,:,,i,j) contains the spectrum around point ((i-1)*q,(j-1)*q)+1.
%
%   A typical use, for an redundancy of 2x2 could be w=2*q+1
%
%   options.bound can be either 'per' or 'sym'
%
%   options.normalization can be set to
%       'tightframe': tight frame transform, with energy conservation.
%       'unit': unit norm basis vectors, usefull to do thresholding
%
%   Copyright (c) 2006 Gabriel Peyre

options.null = 0;
if size(x,3)>1
    dir = -1;
    if nargin<4
        % assume power of 2 size
        n = q*size(x,3);
        n = 2^floor(log2(n));
    end
else
    dir = 1;
    n = size(x,1);
end

if isfield(options, 'bound')
    bound = options.bound;
else
    bound = 'sym';
    bound = 'per';
end

if isfield(options, 'transform_type')
    transform_type = options.transform_type;
else
    transform_type = 'fourier';
end

if isfield(options, 'normalization')
    normalization = options.normalization;
else
    normalization = 'tightframe';
end

% perform sampling
t = 1:q:n+1;
[Y,X] = meshgrid(t,t);
p = size(X,1);

if mod(w,2)==1
% w = ceil((w-1)/2)*2+1;
    w1 = (w-1)/2;
    t = -w1:w1;
else
    t = -w/2+1:w/2;
end
[dY,dX] = meshgrid(t,t);

X = reshape(X,[1 1 p p]);
Y = reshape(Y,[1 1 p p]);
X = repmat( X, [w w 1 1] );
Y = repmat( Y, [w w 1 1] );
dX = repmat( dX, [1 1 p p] );
dY = repmat( dY, [1 1 p p] );

X1 = X+dX;
Y1 = Y+dY;

switch lower(bound)
    case 'sym'
        X1(X1<1) = 1-X1(X1<1);
        X1(X1>n) = 2*n+1-X1(X1>n);
        Y1(Y1<1) = 1-Y1(Y1<1);
        Y1(Y1>n) = 2*n+1-Y1(Y1>n);
    case 'per'
        X1 = mod(X1-1,n)+1;
        Y1 = mod(Y1-1,n)+1;
end

% build a weight function
if isfield(options, 'window_type')
    window_type = options.window_type;
else
    window_type = 'sin';
end
if isfield(options, 'eta')
    eta = options.eta;
else
    eta = 1;
end

if strcmp(window_type, 'sin')
    t = linspace(0,1,w);
    W = sin(t(:)*pi).^2;
    W = W * W';
elseif strcmp(window_type, 'constant')
    W = ones(w);
else
    error('Unkwnown winow.');
end

I = X1 + (Y1-1)*n;

%% renormalize the windows
weight = zeros(n);
for i=1:p
    for j=1:p
        weight(I(:,:,i,j)) = weight(I(:,:,i,j)) + W.^2;
    end
end
weight = sqrt(weight);
Weight = repmat(W, [1 1 p p]);
for i=1:p
    for j=1:p
        Weight(:,:,i,j) = Weight(:,:,i,j) ./ weight(I(:,:,i,j));
    end
end

if strcmp(normalization, 'unit')
    if strcmp(transform_type, 'fourier')
        % for Fourier it is easy
        Renorm = sqrt( sum( sum( Weight.^2, 1 ), 2 ) )/w;
    else
        % for DCT it is less easy ...
        % take a typical window in the middle of the image
        weight = Weight(:,:,round(end/2),round(end/2));
        % compute diracs
        [X,Y,fX,fY] = ndgrid(0:w-1,0:w-1,0:w-1,0:w-1);
        A = 2 * cos( pi/w * ( X+1/2 ).*fX ) .* cos( pi/w * ( Y+1/2 ).*fY ) / w;
        A(:,:,1,:) = A(:,:,1,:) / sqrt(2); % scale zero frequency
        A(:,:,:,1) = A(:,:,:,1) / sqrt(2);
        A = A .* repmat( weight, [1 1 w w] );
        Renorm = sqrt( sum( sum( abs(A).^2, 1 ),2  ) );
    end
    Renorm = reshape( Renorm, w,w );
end

%% compute the transform
if dir==1
    y = zeros(eta*w,eta*w,p,p);
    if mod(w,2)==1
        m = (eta*w+1)/2; w1 = (w-1)/2;
        sel = m-w1:m+w1;
    else
        m = (eta*w)/2+1; w1 = w/2;
        sel = m-w1:m+w1-1;
    end
    y(sel,sel,:,:) = x(I) .* Weight;
    % perform the transform
    y = my_transform( y, +1, transform_type );
    % renormalize if necessary
    if strcmp(normalization, 'unit')
        y = y ./ repmat( Renorm, [1 1 p p] );
    end
else
    if strcmp(normalization, 'unit')
        x = x .* repmat( Renorm, [1 1 p p] );
    end
    x = my_transform( x, -1, transform_type );
    x = real( x.*Weight );
    y = zeros(n);
    for i=1:p
        for j=1:p
            y(I(:,:,i,j)) = y(I(:,:,i,j)) + x(:,:,i,j);
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = my_transform(x,dir,transform_type)

% my_transform - perform either FFT or DCT with energy conservation.
%   Works on array of size (w,w,a,b) on the 2 first dimensions.

w = size(x,1);
if strcmp(transform_type, 'fourier')
    % normalize for energy conservation
    if dir==1
        y = fftshift( fft2(x) ) / w;
    else
        y = ifft2( ifftshift(x*w) );
    end
elseif strcmp(transform_type, 'dct')
    for i=1:size(x,3)
        for j=1:size(x,4)
            y(:,:,i,j) = perform_dct_transform(x(:,:,i,j),dir);
        end
    end
else
    error('Unknown transform');
end

author: gpeyre

repo: matlab-toolboxes


text:
function f = fix(f)
%FIX   Round a CHEBFUN pointwise toward zero.
%   G = FIX(F) returns the CHEBFUN G such that G(x) = FIX(F(x)) for each x in
%   F.domain.
%
%   If F is complex, then the G = FIX(REAL(F)) + 1i*FIX(IMAG(F)).
%
% See also ROUND, CEIL, FLOOR.

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

% Deal with the trivial empty case:
if ( isempty(f) )
    return
end

% Deal with unbounded functions:
if ( ~isfinite(f) )
    error('CHEBFUN:CHEBFUN:fix:inf', ...
        'fix() is not defined for functions which diverge to infinity.');
end

% Deal with complex-valued functions:
if ( ~isreal(f) )
    if ( isreal(1i*f) )
        f = 1i*fix(imag(f));
    else
        f = fix(real(f)) + 1i*fix(imag(f));
    end
    return
end

% Find all the integer crossings for f:
mmvals = minandmax(f); % [TODO]: Only need a good bound?
minf = min(mmvals(1, :));
maxf = max(mmvals(2, :));
range = floor([minf, maxf]);
for k = (range(1)+1):range(2)
    f = addBreaksAtRoots(f - k) + k;
end

% Loop over the FUNs:
for k = 1:numel(f.funs)
    f.funs{k} = fix(f.funs{k});
end

% Fix the pointValues:
f.pointValues = fix(f.pointValues);

% Simplify the result:
f = merge(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [Y, X] = max2( f )
%MAX2   Global maximum of a SEPARABLEAPPROX.
%   Y = MAX2(F) returns the global maximum of F over its domain.
%
%   [Y, X] = MAX2(F) returns the global maximum in Y and its location X.
%
%  This command may be faster if the OPTIMIZATION TOOLBOX is installed.
%
% See also MIN2, MINANDMAX2.

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

% Check for empty:
if ( isempty( f ) )
    return
end

% Call MINANDMAX2:
[Y, X] = minandmax2(f);

% Extract out maximum:
Y = Y(2);
X = X(2,:);

end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   q = solve_spherical_wrist2(robot, q, T, wrist)
%   Solves the inverse kinematic problem for a spherical wrist. This is for
%   the particular RELATIVE orientation of the last three reference systems in
%   ABB robots
%
%   robot: robot structure.
%   q: vector containing the values of the joints 1, 2 and 3.
%   T: orientation of the last reference system.
%   wrist: select -1 or 1 for two possible solutions (wrist up, wrist down)
%
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

function q = solve_spherical_wrist2(robot, q, T, wrist, method)

switch method

    %algebraic solution
    case 'algebraic'
        A01=dh(robot, q, 1);
        A12=dh(robot, q, 2);
        A23=dh(robot, q, 3);

        Q=inv(A23)*inv(A12)*inv(A01)*T;

        %detect the degenerate case when q(5)=0, this leads to zeros
        % in Q13, Q23, Q31 and Q32 and Q33=1
        thresh=1e-12;
        %detect if q(5)==0
        % this happens when cos(q5) in the matrix Q is close to 1
        if abs(Q(3,3)-1)>thresh
            %normal solution
            if wrist==1 %wrist up
                q(4)=atan2(Q(2,3),Q(1,3));
                q(6)=atan2(-Q(3,2),Q(3,1));
            else %wrist down
                q(4)=atan2(Q(2,3),Q(1,3))+pi;
                q(6)=atan2(-Q(3,2),Q(3,1))+pi;
            end
            if abs(cos(q(6)+q(4)))>thresh
                cq5=(-Q(1,1)-Q(2,2))/cos(q(4)+q(6))-1;
            end
            if abs(sin(q(6)+q(4)))>thresh
                cq5=(Q(1,2)-Q(2,1))/sin(q(4)+q(6))-1;
            end
            if abs(sin(q(6)))>thresh
                sq5=Q(3,2)/sin(q(6));
            end
            if abs(cos(q(6)))>thresh
                sq5=-Q(3,1)/cos(q(6));
            end
            q(5)=atan2(sq5,cq5);
        else %degenerate solution, in this case, q4 cannot be determined,
            % so q(4)=0 is assigned
            if wrist==1 %wrist up
                q(4)=0;
                q(5)=0;
                q(6)=atan2(Q(1,2)-Q(2,1),-Q(1,1)-Q(2,2));
            else %wrist down
                q(4)=-pi;
                q(5)=0;
                q(6)=atan2(Q(1,2)-Q(2,1),-Q(1,1)-Q(2,2))+pi;
            end

        end

        %algebraic solution
    case 'geometric'

        % Obtain the position and orientation of the system 3
        % using the already computed joints q1, q2 and q3
        T01=dh(robot, q, 1);
        T12=dh(robot, q, 2);
        T23=dh(robot, q, 3);
        T03=T01*T12*T23;

        x3=T03(1:3,1);
        y3=T03(1:3,2);
        z3=T03(1:3,3);

        % T= [ nx ox ax Px;
        %     ny oy ay Py;
        %     nz oz az Pz];
        a=T(1:3,3);

        % find z4 normal to the plane formed by z3 and a
        z4=cross(z3, a);	% end effector's vector a: T(1:3,3)

        % in case of degenerate solution,
        % when z3 and z6 are parallel, choose q(4)=0 as solution
        if norm(z4) <= 0.000001
            if wrist == 1 %wrist up
                q(4)=0;
            else
                q(4)=-pi;
            end
        else
            cq4=wrist*dot(z4, -y3);
            sq4=wrist*dot(z4, x3);
            q(4)=atan2(sq4, cq4);
        end

        % solve for q5
        T34=dh(robot, q, 4);
        T04=T03*T34;
        x4=T04(1:3, 1);
        y4=T04(1:3, 2);

        z5=T(1:3, 3); % The vector a T(1:3,3) is coincident with z5

        cq5=dot(z5, y4);
        sq5=dot(z5, -x4);
        q(5)=atan2(sq5, cq5);

        % solve for q6
        x6=T(1:3, 1);

        T45=dh(robot, q, 5);
        T05=T04*T45;
        x5=T05(1:3, 1);
        y5=T05(1:3, 2);

        cq6=dot(x6, -x5);
        sq6=dot(x6, -y5);
        q(6)=atan2(sq6, cq6);

    otherwise
        disp('no method specified in solve_spherical_wrist');
end

repo: ARTE


text:
function [fusion,w0] = qfuser_v6(w,scores,wfuse)

if nargin==0
    test_this();
    return;
end

% block 1
f1 = linear_fuser([],scores.scores);
w1 = wfuse;
[whead,wtail] = splitvec_fh(length(w1));
f1 = f1(whead);

% block 2
modelQ = scores.modelQ;
[q,n1] = size(modelQ);
modelQ = [modelQ;ones(1,n1)];
segQ = scores.segQ;
[q2,n2] = size(segQ);
segQ = [segQ;ones(1,n2)];
assert(q==q2);
q = q + 1;

wq = q*(q+1)/2;
r = AWB_fh(modelQ',segQ,tril_to_symm_fh(q));
[whead,wtail] = splitvec_fh(wq,wtail);
r = r(whead);
w2 = zeros(wq,1);w2(end) = -5;

% block 3
s = AWB_fh(modelQ',segQ,tril_to_symm_fh(q,wtail));
w3 = w2;

% assemble
rs = stack([],r,s);
fusion = scalibration_fh(stack(w,f1,rs));
w0 = [w1;w2;w3];

end

function test_this()

m = 3;
k = 2;
n1 = 4;
n2 = 5;

scores.scores = randn(m,n1*n2);
scores.modelQ = randn(k,n1);
scores.segQ = randn(k,n2);

wfuse = [1,2,3,4]';

[fusion,w0] = qfuser_v6([],scores,wfuse);

test_MV2DF(fusion,w0);

[fusion(w0),linear_fuser(wfuse,scores.scores)]

%fusion(w0)

end

author: nesl


text:
s = 256;
mkdir flickr30k-images-256
p = dir('./flickr30k-images/*.jpg');
wp = sprintf('./flickr30k-images-%d/',s);
h1 = [];
w1 = [];
for i = 1:numel(p)
    %disp(i);
    str = [ './flickr30k-images/',p(i).name];
    im = imread(str);
    [h,w,~] = size(im);
    h1 = cat(1,h1,h);
    w1 = cat(1,w1,w);

    if(h<w)
        rate = s/h;
        imr = imresize(im,[s,w*rate]);
    else
        rate = s/w;
        imr = imresize(im,[h*rate,s]);
    end
    imwrite(imr,[ wp,p(i).name]);
    %}
end

disp( mean(h1) );  % 395
disp( mean(w1) );  % 460

author: layumi

repo: Image-Text-Embedding


text:
function Test_trustregions

    % We attempt to compute an intrinsic mean of subspaces, that is, of
    % points on the Grassmann manifold. Let there be m subspaces of
    % dimension p embedded in R^n. We generate random data for our tests:
    % X is a random nxpxm matrix such that each slice of size nxp is an
    % orthonormal basis of a subspace.
    n = 50;
    p = 3;
    m = 100;
    Gr_multi = grassmannfactory(n, p, m);
    X = Gr_multi.rand();

    % Our search space is the Grassmann manifold: we want to locate one
    % point on the Grassmannian that "averages" the m given points.
    Gr = grassmannfactory(n, p);

    % The cost is the sum of squared Riemannian distances to the data
    % points.
    function f = cost(x)
        f = 0;
        for i = 1 : m
            xi = X(:, :, i);
            f = f + Gr.dist(x, xi)^2;
        end
        f = f/(2*m);
    end

    % The gradient is based on the Riemannian logarithmic map at the
    % current point, which gives tangent vectors at x pointing towards the
    % data points.
    function g = grad(x)
        g = Gr.zerovec(x);
        for i = 1 : m
            xi = X(:, :, i);
            g = g - Gr.log(x, xi);
        end
        g = g / m;
    end

    % Setup the problem structure, with the manifold M and the cost and its
    % gradient. Notice that we do not provide a Hessian, because it's
    % rather tricky to compute.
    problem.M = Gr;
    problem.cost = @cost;
    problem.grad = @grad;
%     problem.hess = @(x, xdot) Gr.zerovec(x) ;% * Gr.norm(x, xdot);

    % For peace of mind, check that the gradient is correct.
    % checkgradient(problem);
    % pause;

    % As a simple minded initial guess, choose one of the data points.
    x0 = X(:, :, 1);

    % Setup some options for the trustregions algorihm
    options.tolgradnorm = 1e-16;
    options.maxtime = 30;
    options.maxiter = 200;
    options.verbosity = 2;
    options.debug = 0;
    options.rho_regularization = 1e3;

    % We did not specify a Hessian, but use trustregions anyway. Hence, the
    % Hessian will be approximated, and we should be warned about it. To
    % disable the warning, you may execute this command:
    warning('off', 'manopt:getHessian:approx');

    [x, cost_x, info] = trustregions(problem, x0, options); %#ok<ASGLU>

    xdata = [info.time];
    ydata = [info.gradnorm];
    semilogy(xdata, ydata, '.-');
    xlabel('Time [s]');
    ylabel('Gradient norm');
    hold on;
    radius_change = [0 sign(diff([info.Delta]))];
    text(xdata(radius_change > 0), ydata(radius_change > 0)/2, '+');
    text(xdata(radius_change < 0), ydata(radius_change < 0)/2, '-');
    text(xdata, ydata*2, num2str([info.numinner]'));
    hold off;

    % info
    % keyboard;

end

author: NicolasBoumal

repo: manopt


text:
function pf = loadPoleFigure_siemens(fname,varargin)
% load D5000 data file
%
% Syntax
%   pf = loadPoleFigure_D5000(fname)
%
% Input
%  fname - file name
%
% Output
%  pf    - @PoleFigure
%
% See also
% loadPoleFigure ImportPoleFigureData

fid = fopen(fname,'r');
d = {};
allH = {};
p = 0;

try
  while ~feof(fid)
    line = fgetl(fid);

    % new pole figure section
    if strfind(line,'*Pole figure:')
      p = p+1;
      allH{p} = string2Miller(line(14:end));
      d{p} = [];
      theta{p} = [];

      % new theta angle
    elseif strfind(line,'*Khi')
      theta{p} = [theta{p} sscanf(line(7:end),'%f')]; %#ok<*AGROW>
      bg = [];

      % new background line
    elseif strfind(line,'background')
      bg = [bg sscanf(line(end-8:end),'%8f')];
      if numel(bg)>1, bg = mean(bg); end

    elseif strfind(line,'*')
      %other information

    else % load intensities

      dd = cell2mat(textscan(line,'%n'));
      if ~isempty(bg), dd = dd - bg; end
      d{p} = [d{p},dd];
    end

  end
catch
  interfaceError(fname,fid);
end
fclose(fid);

% generate specimen directions
for p=1:length(allH)
  n = numel(d{p})/numel(theta{p});
  allR{p} = regularS2Grid('theta',theta{p}*degree,'rho',linspace(0,2*pi*(n-1/n),n));
end

% generate pole figure variable
pf = PoleFigure(allH,allR,d,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function [model] = svRegression(X,y, epsilon)
% Compute sizes
[n,d] = size(X);

% Add bias variable
Z = [ones(n,1) X];

one = zeros((d+1) + n, 1);
one(d+2:d+1+n) = 1;

A = [Z -eye(n, n);
    -Z -eye(n, n);
    zeros(n,d+1) -eye(n,n)];

b = [y+epsilon; -y+epsilon; zeros(n,1)];

x = linprog(one, A, b);
w = x(1:d+1);

model.w = w;
model.predict = @predict;

end

function [yhat] = predict(model,Xhat)
[t,d] = size(Xhat);
Zhat = [ones(t,1) Xhat];
yhat = Zhat*model.w;
end

author: gudbrandtandberg


text:
function checkMeshAdjacentFaces(vertices, edges, faces)
%CHECKMESHADJACENTFACES Check if adjacent faces of a mesh have similar orientation
%
%   checkMeshAdjacentFaces(VERTICES, EDGES, FACES)
%   The functions returns no output, but if two faces share a common edge
%   with the same direction (meaning that adjacent faces have normals in
%   opposite direction), a warning is displayed.
%
%   Example
%   [v e f] = createCube();
%   checkMeshAdjacentFaces(v, e, f);
%   % no output -> all faces have normal outwards of the cube
%
%   v = [0 0 0; 10 0 0; 0 10 0; 10 10 0];
%   e = [1 2;1 3;2 3;2 4;3 4];
%   f = [1 2 3; 2 3 4];
%   checkMeshAdjacentFaces(v, e, f);
%      Warning: Faces 1 and 2 run through the edge 3 (2-3) in the same direction
%
%   See also
%   meshes3d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-10-06,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

pattern = 'Faces %d and %d run through the edge %d (%d-%d) in the same direction';

edgeFaces = meshEdgeFaces(vertices, edges, faces);
Ne = size(edgeFaces, 1);

for i = 1:Ne
    % indices of extreimty vertices
    v1 = edges(i, 1);
    v2 = edges(i, 2);

    % index of adjacent faces
    indF1 = edgeFaces(i, 1);
    indF2 = edgeFaces(i, 2);

    % if one of the faces has index 0, then the edge is at the boundary
    if indF1 == 0 || indF2 == 0
        continue;
    end
    % vertices of adjacent faces
    face1 = meshFace(faces, indF1);
    face2 = meshFace(faces, indF2);

    % position of vertices in face vertex array
    ind11 = find(face1 == v1);
    ind12 = find(face1 == v2);
    ind21 = find(face2 == v1);
    ind22 = find(face2 == v2);

    % check if edge is traveled forward or backard
    direct1 = (ind12 == ind11+1) | (ind12 == 1 & ind11 == length(face1));
    direct2 = (ind22 == ind21+1) | (ind22 == 1 & ind21 == length(face2));

    % adjacent faces should travel the edge in opposite direction
    if direct1 == direct2
        warning(pattern, indF1, indF2, i, v1, v2); %#ok<WNTAG>
    end
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function lambda = rutis5_eigenvalues ( )

%*****************************************************************************80
%
%% RUTIS5_EIGENVALUES returns the eigenvalues of the RUTIS5 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real LAMBDA(4,1), the eigenvalues.
%
  n = 4;
  lambda = zeros ( n, 1 );

  lambda(1:4,1) = [ ...
    19.122479087555860; ...
    10.882816916492464; ...
     8.994169735037230; ...
     0.000534260914449 ];

  return
end

author: johannesgerer


text:
function [ abd, ipvt, info ] = sgbfa ( abd, lda, n, ml, mu )

%*****************************************************************************80
%
%% SGBFA factors a real band matrix by elimination.
%
%  Discussion:
%
%    SGBFA is usually called by SGBCO, but it can be called
%    directly with a saving in time if RCOND is not needed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real ABD(LDA,N), the matrix in band
%    storage.  The columns of the matrix are stored in the columns of ABD
%    and the diagonals of the matrix are stored in rows ML+1 through
%    2*ML+MU+1 of ABD.
%
%    Input, integer LDA, the leading dimension of the array ABD.
%    2*ML + MU + 1 <= LDA is required.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ML, MU, the number of diagonals below and above the
%    main diagonal.  0 <= ML < N, 0 <= MU < N.
%
%    Output, real ABD(LDA,N), an upper triangular matrix in band storage
%    and the multipliers which were used to obtain it.  The factorization
%    can be written A = L*U where L is a product of permutation and unit lower
%    triangular matrices and U is upper triangular.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, integer INFO, error flag.
%    0, normal value.
%    K, if U(K,K) == 0.0D+00.  This is not an error condition for this
%      subroutine, but it does indicate that SGBSL will divide by zero if
%      called.  Use RCOND in SGBCO for a reliable indication of singularity.
%
  m = ml + mu + 1;
  info = 0;
%
%  Zero initial fill-in columns.
%
  j0 = mu + 2;
  j1 = min ( n, m ) - 1;

  for jz = j0 : j1
    i0 = m + 1 - jz;
    abd(i0:ml,jz) = 0.0;
  end

  jz = j1;
  ju = 0;
%
%  Gaussian elimination with partial pivoting.
%
  for k = 1 : n-1
%
%  Zero out the next fill-in column.
%
    jz = jz + 1;
    if ( jz <= n )
      abd(1:ml,jz) = 0.0;
    end
%
%  Find L = pivot index.
%
    lm = min ( ml, n-k );
    l = isamax ( lm+1, abd(m:m+lm,k), 1 ) + m - 1;
    ipvt(k) = l + k - m;
%
%  Zero pivot implies this column already triangularized.
%
    if ( abd(l,k) == 0.0 )

      info = k;
%
%  Interchange if necessary.
%
    else

      if ( l ~= m )
        t = abd(l,k);
        abd(l,k) = abd(m,k);
        abd(m,k) = t;
      end
%
%  Compute multipliers.
%
      abd(m+1:m+lm,k) = - abd(m+1:m+lm,k) / abd(m,k);
%
%  Row elimination with column indexing.
%
      ju = min ( max ( ju, mu+ipvt(k) ), n );
      mm = m;

      for j = k+1 : ju
        l = l - 1;
        mm = mm - 1;
        t = abd(l,j);
        if ( l ~= mm )
          abd(l,j) = abd(mm,j);
          abd(mm,j) = t;
        end
        abd(mm+1:mm+lm,j) = ...
          saxpy ( lm, t, abd(m+1:m+lm,k), 1, abd(mm+1:mm+lm,j), 1 );
      end

    end

  end

  ipvt(n) = n;

  if ( abd(m,n) == 0.0 )
    info = n;
  end

  return
end


text:
function p = remove_redundant(p);

if isempty(p.F_struc)
    return
end

b = p.F_struc(1+p.K.f:p.K.l+p.K.f,1);
A = -p.F_struc(1+p.K.f:p.K.l+p.K.f,2:end);

redundant = find(((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <-1e-2));

if length(redundant)>1
    p.InequalityConstraintState(redundant) = inf;
end

if p.options.bmibnb.lpreduce
    b = p.lpcuts(:,1);
    A = -p.lpcuts(:,2:end);
    redundant = find(((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <-1e-2));
    if length(redundant)>1
        p.lpcuts(redundant,:) = [];
        p.cutState(redundant) = [];
    end
end

if any(p.K.f)
    b = p.F_struc(1:p.K.f,1);
    A = -p.F_struc(1:p.K.f,2:end);
    s1 = ((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <1e-6);
    s2 = ((-A>0).*(-A)*(p.ub-p.lb) - ((-b)-(-A)*p.lb) <1e-6);
    redundant = find(s1 & s2);
    if length(redundant)>1
        p.EqualityConstraintState(redundant) = inf;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
clear all; close all; clc

global VAR; %#ok<NUSED>
load_dir = GetOutputDataDir();
load(fullfile(load_dir,'VAR_new.mat'),'VAR'); % stores all clustering indices

%% Init load
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%% batch
totaltime = tic; %#ok<NASGU>
isFullData = 1;
data_masterdir = GetCurrentDataDir();

M_stimrange = GetStimRange();

range_fish = [12:15,17:18];

%% custom params here:
% numK1 = 20;
masterthres = 0.7;

%%
M_regthres = {0.7,0.5};
M_place = {1,2,3,1,4,5,6,7};
M_stimname = {'4x4','PT','OMR','defS','Spt','DF','Lm','Dot'};
%%
for i_count = 1,%%%%%%%%%
    masterthres = M_regthres{i_count};
%     clusParams = struct('merge',masterthres,'cap',masterthres,'reg1',masterthres,...
%         'reg2',masterthres,'minSize',10,'k1',numK1);

    for i_stimrange = 1:7
        if i_stimrange == 1,
            M_stimrange = GetStimRange('5');
        elseif i_stimrange == 2,
            M_stimrange = GetStimRange('P');
        elseif i_stimrange == 3,
            M_stimrange = GetStimRange('O');
        elseif i_stimrange == 4,
            M_stimrange = GetStimRange('M');
        elseif i_stimrange == 5,
            M_stimrange = GetStimRange('S');
        elseif i_stimrange == 6,
            M_stimrange = GetStimRange('D');
        elseif i_stimrange == 7,
            M_stimrange = GetStimRange('L');
%         elseif i_stimrange == 8,
%             M_stimrange = GetStimRange('Y');
        end

        for i = 1:length(range_fish),
            i_fish = range_fish(i);
            disp(i_fish);

            % check this loop
            stimrange = M_stimrange{i_fish};
            if isempty(stimrange),
                continue;
            end

            % Load fish
            LoadFullFish(hfig,i_fish,isFullData);

            %% 1.
            % setup
            absIX = getappdata(hfig,'absIX');

            i_ClusGroup = 2;
            i_Cluster = 1;

            % Load cluster data
            [cIX_load,gIX] = LoadCluster_Direct(i_fish,i_ClusGroup,i_Cluster,absIX);

            %% partitions for CV
            timelists = getappdata(hfig,'timelists');
            timelists_names = getappdata(hfig,'timelists_names');
            periods = getappdata(hfig,'periods');

            M_stim = M_stimrange{i_fish};

            timelistsCV_raw = cell(length(M_stim),2);
            timelistsCV = cell(1,2);

            for k_stim = 1:length(M_stim), % :3
                i_stim = M_stim(k_stim);
                TL = timelists{i_stim};
                seq = randperm(size(TL,2));
                halfpoint = floor(size(TL,2)/2);
                timelistsCV_raw{k_stim,1} = TL(seq(1:halfpoint));
                timelistsCV_raw{k_stim,2} = TL(seq(1+halfpoint:2*halfpoint));
%                 period = periods(i_stim);
%                 nrep = size(TL,2)/periods(i_stim); % integer
%                 n = floor(nrep/2);
%                 if n>0,
%                     timelistsCV_raw{k_stim,1} = TL(1:n*period);
%                     timelistsCV_raw{k_stim,2} = TL(1+n*period:2*n*period);% before 12/5/16: TL(1+n*period):TL(2*n*period);
%                 else % for spont, only one period
%                     halfperiod = floor(period/2);
%                     timelistsCV_raw{k_stim,1} = TL(1:halfperiod);
%                     timelistsCV_raw{k_stim,2} = TL(1+halfperiod:2*halfperiod);
%                 end
            end
            timelistsCV{1} = horzcat(timelistsCV_raw{:,1});
            timelistsCV{2} = horzcat(timelistsCV_raw{:,2});
            assert(length(timelistsCV{1})==length(timelistsCV{2}));

            %%
            for k = 1:2,% CV halves
                tIX = timelistsCV{k};
                M = GetTimeIndexedData_Default_Direct(hfig,cIX_load,tIX);
                M_0 = GetTimeIndexedData_Default_Direct(hfig,[],tIX,'isAllCells');

                % ------custom code here---------
                isWkmeans = true;
                isMakeFoxels = true;

                [cIX,gIX] = AutoClustering(cIX_load,gIX,M_0,cIX_load,isWkmeans,[],...
                    isMakeFoxels,masterthres);

                % save cluster
                name = ['Auto_',M_stimname{i_stimrange},'_M',num2str(masterthres),'_CV',num2str(k)];
                clusgroupID = 15+k; % 16 and 17 %13+k; % 14 and 15
                clusIDoverride = M_place{i_stimrange};
                SaveCluster_Direct(cIX,gIX,absIX,i_fish,name,clusgroupID,clusIDoverride);
            end
        end
    end
end
SaveVARwithBackup();
totaltime = toc;
disp(totaltime);

author: xiuyechen

repo: FishExplorer


text:
function out = ne(sR1,sR2)

out = ~eq(sR1,sR2);


text:
function test_suite=test_meeg_chan_neighborhood()
% tests for cosmo_meeg_chan_neighborhood
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_neighbors()
    % note: this tests assumes that meeg_chan_neighbors works properly
    if cosmo_skip_test_if_no_external('fieldtrip')
        return;
    end

    % only test every test_step-th channel
    test_step=10;

    ds=cosmo_synthetic_dataset('type','meeg','size','big');
    msk=ds.fa.chan<10 | mod(ds.fa.chan,5)==0 | ds.fa.chan>250;
    imsk=find(msk);
    rp=randperm(numel(imsk));
    ds=cosmo_slice(ds,imsk(rp),2);
    ds=cosmo_dim_prune(ds);
    n=numel(ds.a.fdim.values{1});
    ds.a.fdim.values{1}=ds.a.fdim.values{1}(randperm(n));

    nbrs=cosmo_meeg_chan_neighbors(ds,'chantype','meg_planar','radius',0);
    nh=cosmo_meeg_chan_neighborhood(ds,nbrs);

    assertEqual(nh.a.fdim.values,{{nbrs.label}});
    assertEqual(nh.a.fdim.labels,{'chan'});

    n=numel(nh.neighbors);
    for k=1:test_step:n
        idx=nh.neighbors{k};
        chan_label=ds.a.fdim.values{1}(ds.fa.chan(idx));
        assert(all(cosmo_match(chan_label,nh.a.fdim.values{1}{k})));
    end

    % test correspondence with neighbors
    args={'chantype','all','delaunay',true};
    nbrs=cosmo_meeg_chan_neighbors(ds,args{:});
    nh=cosmo_meeg_chan_neighborhood(ds,args{:});
    nh2=cosmo_meeg_chan_neighborhood(ds,nbrs);
    assertEqual(nh,nh2);

    ds_label=ds.a.fdim.values{1};

    n=numel(nh.neighbors);
    for k=1:test_step:n
        idx=nh.neighbors{k};

        chan_label=ds.a.fdim.values{1}(ds.fa.chan(idx));

        % ensure both arguments for intersect are column vectors,
        % because Octave behaves differently than Matlab
        assertEqual(intersect(ds_label(:),nbrs(k).neighblabel(:)),...
                        unique(chan_label'));

    end

    % try with dataset labels
    args={'chantype','all','count',5,'label','dataset'};
    nbrs=cosmo_meeg_chan_neighbors(ds,args{:});
    nh=cosmo_meeg_chan_neighborhood(ds,nbrs);
    n=numel(nh.neighbors);
    assertEqual(n, numel(ds.a.fdim.values{1}));
    assertEqual(1:n, nh.fa.chan);
    ds_label=ds.a.fdim.values{1};
    nh_label=nh.a.fdim.values{1};

    for k=1:test_step:n
        i=find(cosmo_match(nh_label,ds_label{k}));
        assertEqual(nbrs(i).label,ds_label{k});

        overlap=cosmo_overlap({ds_label(ds.fa.chan(nh.neighbors{i}))},...
                                {nbrs(i).neighblabel});
        assertEqual(overlap,1);
    end

    % test for number of channels
    nbrs=cosmo_meeg_chan_neighbors(ds,'count',5,...
                        'chantype','meg_combined_from_planar');
    nh=cosmo_meeg_chan_neighborhood(ds,nbrs);
    assertEqual(numel(nh.neighbors),102);
    h=cellfun(@numel,nh.neighbors)/7;
    assert(all(h>=5 & h<=10));
    mh=mean(h);
    assert(mh>6 && mh<9);

    nbrs=cosmo_meeg_chan_neighbors(ds,'count',5,'chantype','meg_planar');
    nh=cosmo_meeg_chan_neighborhood(ds,nbrs);
    assertEqual(numel(nh.neighbors),204);
    h=cellfun(@numel,nh.neighbors)/7;
    assert(all(h==5));

    % test tiny dataset
    ds=cosmo_synthetic_dataset('type','meeg');
    opt=struct();
    opt.delaunay=true;
    opt.label='dataset';
    opt.chantype='meg_axial';

    nh=cosmo_meeg_chan_neighborhood(ds,opt);
    assertEqual(nh.neighbors,{[1 4]});

    opt.chantype='meg_planar';
    nh=cosmo_meeg_chan_neighborhood(ds,opt);
    assertEqual(nh.neighbors,{[2 5 3 6]; [2 5 3 6]});

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [E,J]=SynthMeasWatsonSHCylSingleRadIsoV_GPD_B0(x, protocol, fibredir, roots)
% Substrate: Impermeable cylinders with one radius in a homogeneous background.
% Orientation distribution: Watson's distribution with SH approximation
% Signal approximation: Gaussian phase distribution.
% Notes: This version includes an isotropic diffusion compartment with its own
% diffusivity.
% Includes a free parameter for the measurement at b=0.
%
% [E,J]=SynthMeasWatsonSHCylSingleRadIsoV_GPD_B0(x, protocol, fibredir, roots)
% returns the measurements E according to the model and the Jacobian J of the
% measurements with respect to the parameters.  The Jacobian does not
% include derivates with respect to the fibre direction.
%
% x is the list of model parameters in SI units:
% x(1) is the volume fraction of the intracellular space.
% x(2) is the free diffusivity of the material inside and outside the cylinders.
% x(3) is the hindered diffusivity outside the cylinders in perpendicular directions.
% x(4) is the radius of the cylinders.
% x(5) is the concentration parameter of the Watson's distribution.
% x(6) is the volume fraction of the isotropic compartment.
% x(7) is the diffusivity of the isotropic compartment.
% x(8) is the measurement at b=0.
%
% protocol is the object containing the acquisition protocol.
%
% fibredir is a unit vector along the symmetry axis of the Watson's
% distribution.  It must be in Cartesian coordinates [x y z]' with size [3 1].
%
% roots contains solutions to the Bessel function equation from function
% BesselJ_RootsCyl.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

S0 = x(8);

% Call the other function to get normalized measurements.
if(nargout == 1)
    Enorm=SynthMeasWatsonSHCylSingleRadIsoV_GPD(x, protocol, fibredir, roots);
else
    [Enorm,Jnorm]=SynthMeasWatsonSHCylSingleRadIsoV_GPD(x, protocol, fibredir, roots);
end

E = Enorm*S0;

if(nargout>1)
    J = Jnorm*S0;
    [meas, params] = size(J);
    J(:,params+1) = Enorm;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function K = prtUtilQuadExpCovariance(x1, x2, theta0, theta1, theta2, theta3)
% xxx Need Help xxx
% K = prtUtilQuadExpCovariance(x1, x2, theta0, theta1, theta2, theta3)
%
% Common covariance function for Gaussian Processes
%
% Bishop, Pattern Recognition and Machine Learning Eq. 6.63

if nargin==3 && isvector(theta0) && numel(theta0)==4
    theta3 = theta0(4);
    theta2 = theta0(3);
    theta1 = theta0(2);
    theta0 = theta0(1);
else
    if nargin < 6 || isempty(theta3)
        theta3 = 0;
    end

    if nargin < 5 || isempty(theta2)
        theta2 = 0;
    end

    if nargin < 4 || isempty(theta1)
        theta1 = 1;
    end

    if nargin < 3 || isempty(theta0)
        theta0 = 1;
    end
end

K = theta0 * exp(-theta1/2 * prtDistanceEuclidean(x1,x2).^2) + theta2 + theta3*x1*x2';

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
classdef CorePair < handle
    % CorePair. A class for storing two values
    %
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        First
        Second
    end

    methods
        function obj = CorePair(first, second)
            obj.First = first;
            obj.Second = second;
        end

        function first = FirstList(obj)
            first = {obj.First};
        end

        function first = SecondList(obj)
            first = {obj.First};
        end

        function is_equal = eq(obj, other)
            % Equality operator

            if ~isa(other, 'CorePair')
                is_equal = false;
            else
                is_equal = isequal(obj.First, other.First) && isequal(obj.Second, other.Second);
            end
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function varargout=sprintcpx(varargin)
%
% For calling details please see v_sprintcpx.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_sprintcpx.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: sprintcpx.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_sprintcpx(varargin{:});
else
    v_sprintcpx(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%% simple example of use of renorm_sibling_layer
%% compute scattering
clear; close all;
Wop = wavelet_factory_2d_pyramid();
x = uiuc_sample;
Sx = scat(x, Wop);
%% extract the second layer
layer = Sx{3};
op = @(x)(sum(x,3));

%% to renormalized order 2, the sibling is all the node with same ancestor
sibling = @(p)(find(Sx{3}.meta.j(1,:) == Sx{3}.meta.j(1,p) & ...
    Sx{3}.meta.theta(1,:) == Sx{3}.meta.theta(1,p)));

%% renormalize
layer_renorm = renorm_sibling_layer(layer, op, sibling);

%% more sophistated example of use of renorm_sibling_layer
%% smooth a bit + L1 norm instead of just L1 norm
options.sigma_phi = 1;
options.P = 4;
filters = morlet_filter_bank_2d_pyramid(options);
h = filters.h.filter;
smooth = @(x)(conv_sub_2d(x, h, 0));
op = @(x)(smooth(sum(x,3)));

%% renormalize
layer_renorm = renorm_sibling_layer(layer, op, sibling);

%% display
figure(1);
image_scat_layer(layer, 0, 0);
title('second order of scattering')
figure(2);
image_scat_layer(layer_renorm, 0, 0);
title('normalized second order of scattering')

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [mm] = dm2mm(dm)
% Convert length from decimeters to millimeters.
% Chad A. Greene 2012
mm = dm*100;

author: Sable

repo: mcbench-benchmarks


text:
%% DEMO_febio_0029_contact_friction_benchmark_boxes
% Below is a demonstration for a FEBio demo presented at the FEBio workshop at
% WCB 2018. Four cubes are stacked on top of each other. Sliding-elastic
% contact with friction is defined between them. The top cube is forced
% downwards in a rotation motion thereby compressing the other cubes and
% also transferring some twist due to frictional forces.
%
% The demo includes
% * Building geometry for 4 cubes with hexahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * indentation
% * contact, sliding, sliding-elastic, friction
% * rigid body constraints
% * hexahedral elements, hex8
% * quadrilateral elements, quad4
% * static, solid
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=15;
faceAlpha1=0.8;
faceAlpha2=0.3;
markerSize=40;
markerSize2=20;
lineWidth=3;

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=fullfile(savePath,[febioFebFileNamePart,'.txt']); %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_strainEnergy=[febioFebFileNamePart,'_fsed_out.txt']; %Log file name for exporting strain energy density

%Specifying dimensions and number of elements for slab
sampleHeight=2; %Height
sampleWidth=sampleHeight; %Width
sampleThickness=sampleHeight; %Thickness

numElementsWidth=[4 3 5 4]; %Number of elemens in dir 1
numElementsThickness=numElementsWidth; %Number of elemens in dir 2
numElementsHeight=numElementsWidth; %Number of elemens in dir 3

%Material parameter set
youngsModuli   = [0.3 10  0.3 10 ];
poissonsRatios = [0.4 0.1 0.4 0.1];

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=25; %Optimum number of iterations
max_retries=10; %Maximum number of retires
dtmin=(1/numTimeSteps)/10; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size

%Contact parameters
contactInitialOffset=0;
boxOffsets=sampleHeight+contactInitialOffset;

%Prescribed displacement
prescribedDisplacement_Z=-2;
prescribedRotation_Z=pi/2;

%% Creating model geometry and mesh
% A box is created with tri-linear hexahedral (hex8) elements using the
% |hexMeshBox| function. The function offers the boundary faces with
% seperate labels for the top, bottom, left, right, front, and back sides.
% As such these can be used to define boundary conditions on the exterior.

E=[];
elementMaterialID=[];
V=[];
Fb=[];
Cb=[];
for q=1:1:4

    % Create a box with hexahedral elements
    beamDimensions=[sampleWidth sampleThickness sampleHeight]; %Dimensions
    beamElementNumbers=[numElementsWidth(q) numElementsThickness(q) numElementsHeight(q)]; %Number of elements
    outputStructType=2; %A structure compatible with mesh view
    [meshStruct]=hexMeshBox(beamDimensions,beamElementNumbers,outputStructType);

    %Access elements, nodes, and faces from the structure
    E1=meshStruct.elements; %The elements
    V1=meshStruct.nodes; %The nodes (vertices)
    Fb1=meshStruct.facesBoundary; %The boundary faces
    Cb1=meshStruct.boundaryMarker; %The "colors" or labels for the boundary faces
    elementMaterialIndices=ones(size(E1,1),1); %Element material indices

    V1(:,3)=V1(:,3)+(q-1)*boxOffsets;

    E=[E;E1+size(V,1)];
    Fb=[Fb;Fb1+size(V,1)];
    V=[V;V1];
    colorOffset=max(Cb(:));
    if isempty(colorOffset)
        colorOffset=0;
    end
    Cb=[Cb;Cb1+colorOffset];

    elementMaterialID=[elementMaterialID;q*ones(size(E1,1),1)];

end
V(:,3)=V(:,3)-min(V(:,3)); %Shift so bottom is at 0

%%
% Plotting model boundary surfaces and a cut view

hFig=cFigure;
hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
gpatch(Fb,V,Cb,'k',faceAlpha1);
colormap(gjet(250)); icolorbar;
axisGeom(gca,fontSize);
drawnow;

%%

logicTops=false(size(Cb,1),4);
logicBottoms=false(size(Cb,1),4);
for q=1:1:4
    logicTops(:,q)=Cb==6+(6*(q-1));
    logicBottoms(:,q)=Cb==5+(6*(q-1));
end

%%
% Plotting model boundary surfaces

hFig=cFigure;
hold on;
title('Contact faces','FontSize',fontSize);
gpatch(Fb,V,'kw','none',0.2);
for q=1:1:size(logicTops,2)-1
    gpatch(Fb(logicTops(:,q),:),V,q*ones(nnz(logicTops(:,q)),1),'g',1);
    gpatch(Fb(logicBottoms(:,q+1),:),V,q*ones(nnz(logicBottoms(:,q+1)),1),'y',1);
end

colormap(gjet(250)); icolorbar;
axisGeom(gca,fontSize);
drawnow;

%% Define boundary conditions

F_support=Fb(logicBottoms(:,1),:);

F_rigidBody=Fb(logicTops(:,end),:);
indNodesRigidBody=unique(F_rigidBody(:));
center_of_mass=mean(V(indNodesRigidBody,:),1);

%Supported nodes
bcSupportList=unique(F_support(:));

%Visualize BC's
hf=cFigure;
title('Boundary conditions model','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'kw','none',faceAlpha2);
hl(1)=gpatch(F_rigidBody,V,'rw','k',1);
hl(2)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);
hl(3)=plotV(center_of_mass,'r.','MarkerSize',50);

legend(hl,{'Rigid body','BC support','Rigid body center of mass'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.solver.symmetric_stiffness=0;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
for q=1:1:numel(youngsModuli)
    materialNameNow=['Material',num2str(q)];
    febio_spec.Material.material{q}.ATTR.name=materialNameNow;
    febio_spec.Material.material{q}.ATTR.type='neo-Hookean';
    febio_spec.Material.material{q}.ATTR.id=q;
    febio_spec.Material.material{q}.E=youngsModuli(q);
    febio_spec.Material.material{q}.v=poissonsRatios(q);
end

materialNameRigid='RigidMaterial';
matIdRigid=numel(youngsModuli)+1;
febio_spec.Material.material{matIdRigid}.ATTR.name=materialNameRigid;
febio_spec.Material.material{matIdRigid}.ATTR.type='rigid body';
febio_spec.Material.material{matIdRigid}.ATTR.id=numel(youngsModuli)+1;
febio_spec.Material.material{matIdRigid}.density=1;
febio_spec.Material.material{matIdRigid}.center_of_mass=center_of_mass;

%Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
n=1;
for q=1:1:numel(youngsModuli)
    logicMaterialNow=(elementMaterialID==q);
    partNameNow=['Part',num2str(q)];
    febio_spec.Mesh.Elements{q}.ATTR.name=partNameNow; %Name of this part
    febio_spec.Mesh.Elements{q}.ATTR.type='hex8'; %Element type of this set
    febio_spec.Mesh.Elements{q}.elem.ATTR.id=(n:1:(n-1+nnz(logicMaterialNow)))'; %Element id's
    febio_spec.Mesh.Elements{q}.elem.VAL=E(logicMaterialNow,:);
    n=n+nnz(logicMaterialNow);
end

partNameRigid='RigidPart';
n=max(febio_spec.Mesh.Elements{end}.elem.ATTR.id);
febio_spec.Mesh.Elements{numel(youngsModuli)+1}.ATTR.name=partNameRigid; %Name of this part
febio_spec.Mesh.Elements{numel(youngsModuli)+1}.ATTR.type='quad4'; %Element type of this set
febio_spec.Mesh.Elements{numel(youngsModuli)+1}.ATTR.mat=numel(youngsModuli)+1; %material index for this set
febio_spec.Mesh.Elements{numel(youngsModuli)+1}.elem.ATTR.id=(n:1:(n-1+size(F_rigidBody,1)))'; %Element id's
febio_spec.Mesh.Elements{numel(youngsModuli)+1}.elem.VAL=F_rigidBody;

% -> NodeSets
nodeSetName1='bcSupportList';
febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList(:);

%MeshDomains section
for q=1:1:numel(youngsModuli)
    partNameNow=['Part',num2str(q)];
    materialNameNow=['Material',num2str(q)];
    febio_spec.MeshDomains.SolidDomain{q}.ATTR.name=partNameNow;
    febio_spec.MeshDomains.SolidDomain{q}.ATTR.mat=materialNameNow;
end

febio_spec.MeshDomains.ShellDomain.ATTR.name=partNameRigid;
febio_spec.MeshDomains.ShellDomain.ATTR.mat=materialNameRigid;

% -> Surfaces
febio_spec.Mesh.Surface=[];
for q=1:1:numel(youngsModuli)-1
    F_contact_now=Fb(logicTops(:,q),:);
    c=numel(febio_spec.Mesh.Surface)+1;
    febio_spec.Mesh.Surface{c}.ATTR.name=['contact_',num2str(c)];
    febio_spec.Mesh.Surface{c}.quad4.ATTR.id=(1:1:size(F_contact_now,1))';
    febio_spec.Mesh.Surface{c}.quad4.VAL=F_contact_now;

    F_contact_now=Fb(logicBottoms(:,q+1),:);
    c=numel(febio_spec.Mesh.Surface)+1;
    febio_spec.Mesh.Surface{c}.ATTR.name=['contact_',num2str(c)];
    febio_spec.Mesh.Surface{c}.quad4.ATTR.id=(1:1:size(F_contact_now,1))';
    febio_spec.Mesh.Surface{c}.quad4.VAL=F_contact_now;
end

% -> Surface pairs
surfaceIndices=reshape(1:(2*(numel(youngsModuli)-1)),2,3)'; %Indices for surface pairs
for q=1:1:numel(youngsModuli)-1
    febio_spec.Mesh.SurfacePair{q}.ATTR.name=['Contact_',num2str(q)];

    numElements1=size(febio_spec.Mesh.Surface{surfaceIndices(q,1)}.quad4.VAL,1);
    numElements2=size(febio_spec.Mesh.Surface{surfaceIndices(q,2)}.quad4.VAL,1);
    if numElements1>numElements2
        febio_spec.Mesh.SurfacePair{q}.primary = febio_spec.Mesh.Surface{surfaceIndices(q,1)}.ATTR.name;
        febio_spec.Mesh.SurfacePair{q}.secondary = febio_spec.Mesh.Surface{surfaceIndices(q,2)}.ATTR.name;
    else
        febio_spec.Mesh.SurfacePair{q}.primary = febio_spec.Mesh.Surface{surfaceIndices(q,2)}.ATTR.name;
        febio_spec.Mesh.SurfacePair{q}.secondary = febio_spec.Mesh.Surface{surfaceIndices(q,1)}.ATTR.name;
    end
end

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x,y,z';

% %Rigid section
% -> Prescribed rigid body boundary conditions
febio_spec.Rigid.rigid_constraint{1}.ATTR.name='RigidFix_1';
febio_spec.Rigid.rigid_constraint{1}.ATTR.type='fix';
febio_spec.Rigid.rigid_constraint{1}.rb=matIdRigid;
febio_spec.Rigid.rigid_constraint{1}.dofs='Rx,Ry,Ru,Rv';

febio_spec.Rigid.rigid_constraint{2}.ATTR.name='RigidPrescribe';
febio_spec.Rigid.rigid_constraint{2}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{2}.rb=matIdRigid;
febio_spec.Rigid.rigid_constraint{2}.dof='Rz';
febio_spec.Rigid.rigid_constraint{2}.value.ATTR.lc=1;
febio_spec.Rigid.rigid_constraint{2}.value.VAL=prescribedDisplacement_Z;
febio_spec.Rigid.rigid_constraint{2}.relative=0;

febio_spec.Rigid.rigid_constraint{3}.ATTR.name='RigidPrescribe';
febio_spec.Rigid.rigid_constraint{3}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{3}.rb=matIdRigid;
febio_spec.Rigid.rigid_constraint{3}.dof='Rw';
febio_spec.Rigid.rigid_constraint{3}.value.ATTR.lc=1;
febio_spec.Rigid.rigid_constraint{3}.value.VAL=prescribedRotation_Z;
febio_spec.Rigid.rigid_constraint{3}.relative=0;

%Contact section
for q=1:1:numel(youngsModuli)-1
    febio_spec.Contact.contact{q}.ATTR.surface_pair=febio_spec.Mesh.SurfacePair{q}.ATTR.name;
    febio_spec.Contact.contact{q}.ATTR.type='sliding-elastic';
    febio_spec.Contact.contact{q}.two_pass=1;
    febio_spec.Contact.contact{q}.laugon=1;
    febio_spec.Contact.contact{q}.tolerance=0.2;
    febio_spec.Contact.contact{q}.gaptol=0;
    febio_spec.Contact.contact{q}.minaug=1;
    febio_spec.Contact.contact{q}.maxaug=10;
    febio_spec.Contact.contact{q}.search_tol=0.01;
    febio_spec.Contact.contact{q}.search_radius=0.1*sqrt(sum((max(V,[],1)-min(V,[],1)).^2,2));
    febio_spec.Contact.contact{q}.symmetric_stiffness=0;
    switch q
        case 1
            febio_spec.Contact.contact{q}.auto_penalty=1;
            febio_spec.Contact.contact{q}.penalty=5;
            febio_spec.Contact.contact{q}.fric_coeff=1;
        case 2
            febio_spec.Contact.contact{q}.auto_penalty=1;
            febio_spec.Contact.contact{q}.penalty=0.1;
            febio_spec.Contact.contact{q}.fric_coeff=0.08;
        case 3
            febio_spec.Contact.contact{q}.auto_penalty=1;
            febio_spec.Contact.contact{q}.penalty=0.1;
            febio_spec.Contact.contact{q}.fric_coeff=0.08;
    end
end

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_strainEnergy;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sed';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode='external';%'internal';

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

 %%
    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

     %%
    % Importing element stress from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_strainEnergy),1,1);

    %Access data
    E_energy=dataStruct.data;

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    [CV]=faceToVertexMeasure(E,V,E_energy(:,:,end));

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('$\sigma_{zz}$ [MPa]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),CV,'k',1); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    hp2=gpatch(F_rigidBody,V_DEF(:,:,end),'w','k',1); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([min(E_energy(:)) max(E_energy(:))]/2);
    axis(axisLim(V_DEF)); %Set axis limits statically
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments

        [CV]=faceToVertexMeasure(E,V,E_energy(:,:,qt));

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp hp2]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData','Vertices'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),CV,V_DEF(:,:,qt)}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON