text:
classdef ClassifierStuff < handle
  % Class that includes the classifier plus classifier metadata, like the
  % classifier training parameters

  properties
    type
    params
    timeStamp
    confThresholds % 1x2 vec
    scoreNorm % AL does windowdata.scoreNorm dup .scoreNorm?
    postProcessParams
    trainingParams
    %featureNames
    windowdata
    savewindowdata
    predictOnlyCurrentFly
    selFeatures
  end

  % store featureCodes instead
  properties (Dependent)

    featureNames

  end

  % use featureCode instead of featureNames to save space
  properties (Access = public)

    featureCodes
    perframeFeatureNames
    statNames
    transNames

  end

  methods

    function self = ClassifierStuff(varargin)

      self.type = 'boosting';
      self.params = struct('dim',cell(0,1), ...
                           'error',[], ...
                           'dir',[], ....
                           'tr',[], ...
                           'alpha',[]);
      self.timeStamp = 0;  % default time stamp (bad idea?)
      self.confThresholds = [0 0];
      self.scoreNorm = nan;  % default time stamp, number of days since Jan 0, 0000 (typically non-integer)
      self.postProcessParams.method = 'Hysteresis';
      self.postProcessParams.hystopts(1) = struct('name','High Threshold','tag','hthres','value',0);
      self.postProcessParams.hystopts(2) = struct('name','Low Threshold','tag','lthres','value',0);
      self.postProcessParams.filtopts(1) = struct('name','Size','tag','size','value',1);
      self.postProcessParams.blen = 1;
      self.trainingParams = ...
        struct('iter',100, ...
               'iter_updates',10, ...
               'numSample',2500, ...
               'numBins',30, ...
               'CVfolds',7, ...
               'baseClassifierTypes','Decision Stumps', ...
               'baseClassifierSelected',1,...
               'nselfeatures',1000);

      % KB 20201208
      % replace featureNames with a code, make featureNames be dependent
      % variable
      self.featureCodes = nan(0,5);
      self.perframeFeatureNames = {};
      self.statNames = {};
      self.transNames = {};
      self.windowdata = WindowData.windowdata(0);
      self.savewindowdata = false;
      self.predictOnlyCurrentFly = false;
      self.selFeatures =  SelFeatures.createEmpty();

      % varargin should be key-value pairs
      % Keys should be property names, values should be initialization
      % values for the corresponding property
      keys = varargin(1:2:end);
      values = varargin(2:2:end);

      % Process the keys, assigning to any that match properties, and
      % throwing an error if a key matches no property
      propertyNames = properties(self);
      for iKey=1:length(keys)
        key = keys{iKey};
        iProperty = whichstr(key,propertyNames);
        if isempty(iProperty)
          error('Classifier:noPropertyMatchingKey', ...
                'Classifier: No property matching key %s',key);
        else
          propertyName = propertyNames{iProperty};
          value = values{iKey};
          self.(propertyName) = value;
        end
      end
    end

    function tfmodified = modernize(self)
      tfmodified = false;
      for i = 1:numel(self)
        if isempty(self(i).savewindowdata)
          self(i).savewindowdata = false;
          tfmodified = true;
        end
      end

      % MK: May 09 2016
      for i = 1:numel(self)
        % AL: prob don't need to check isprop(...,'selFeatures') b/c the
        % property is in the class now -- any/all ClassifierStuffs running
        % with this version should get the prop
        if ~isprop(self(i),'selFeatures') || isempty(self(i).selFeatures),
          tfmodified = true;
          self(i).selFeatures = SelFeatures.createEmpty();
        else
          [self(i).selFeatures,sfmod] = SelFeatures.modernize(self(i).selFeatures);
          tfmodified = tfmodified || sfmod;
        end
      end

      % MK: May 09 2016
      for i = 1:numel(self)
        if ~isfield(self(i).trainingParams,'nselfeatures'),
          tfmodified = true;
          self(i).trainingParams.nselfeatures = 1000;
        end
      end

      % MK: May 09 2016
      for i = 1:numel(self)
        % AL: isprop etc
        if ~isprop(self(i),'predictOnlyCurrentFly') || isempty(self(i).predictOnlyCurrentFly),
          tfmodified = true;
          self(i).predictOnlyCurrentFly = false;
        end
      end

%       % MK: modernize window data.
%       wdmodified = false;
%       for i = 1:numel(self)
%           [self(i).windowdata,wdm] =  WindowData.modernize(self(i).windowdata);
%           if wdm,
%               wdmodified = true;
%           end
%       end
%       if wdmodified,
%           tfmodified = true;
%       end

      % MK: modernize window data.
      % AL: fixed for classifiers without windowdata
      wdmodified = false;
      for i = 1:numel(self)
          if ~isequal(self(i).windowdata,[])
              [self(i).windowdata,wdm] =  WindowData.modernize(self(i).windowdata);
              if wdm,
                  wdmodified = true;
              end
          end
      end
      if wdmodified,
          tfmodified = true;
      end

    end

    function set.featureNames(obj,featureNames)

      [obj.featureCodes,obj.perframeFeatureNames,obj.statNames,obj.transNames] = ...
        windowFeatureNames2Codes(featureNames);

    end

    function featureNames = get.featureNames(obj)
      featureNames = windowFeatureCodes2Names(obj.featureCodes,obj.perframeFeatureNames,obj.statNames,obj.transNames);
    end

  end

end

author: kristinbranson

repo: JAABA


text:
function hplstat(HPL,HPE,HAL,src_name)
% HPL(1xN) - Horizontal protection level
% HPE(Nx1) - Horizontal position error = (easting.^2, northing.^2).^0.5
% HAL - Horizontal alert limit (define how large is the area for Normal operation
%       [the line between white and yellow])
% src_name = label of vertical axis which will go as subscript after 'HPL'

%*     Copyright c 1998 The board of trustees of the Leland Stanford     *
%*                      Junior University. All rights reserved.          *
%*     This script file may be distributed and used freely, provided     *
%*     this copyright notice is always kept with it.                     *
%*                                                                       *
%*     Questions and comments should be directed to Todd Walter at:      *
%*     walter@relgyro.stanford.edu                                       *
%HPLSTAT
%   HPLSTAT(FNAME, HAL, PATH) reads the given file of TMS Horizontal
%   Protection Limit statistics and plots the color coded histogram.
%   In addition the regions of integrity failures and unavailability
%   are shown.  The optional HAL argument is given in meters (default 30m).
%   See also : TMSSTAT, TRSNAMES, VERTSTAT, VPLSTAT
%   Example code for generating the file is included at the bottom of hplstat.m
%
%   NOTE: To plot on a non-color printer type:
%                 colormap('gray')
% % new colormap for colorbar (HPL and HPE plot)
 color_bar =[
     0.6980    0.6980    0.6980;
    0.6905    0.6905    0.6905;
     0.6830    0.6830    0.6830;
     0.6754    0.6754    0.6754;
     0.6679    0.6679    0.6679;
     0.6604    0.6604    0.6604;
     0.6528    0.6528    0.6528;
     0.6453    0.6453    0.6453;
     0.6378    0.6378    0.6378;
     0.6302    0.6302    0.6302;
     0.6227    0.6227    0.6227;
     0.6152    0.6152    0.6152;
     0.6076    0.6076    0.6076;
     0.6001    0.6001    0.6001;
     0.5926    0.5926    0.5926;
     0.5850    0.5850    0.5850;
     0.5775    0.5775    0.5775;
     0.5700    0.5700    0.5700;
     0.5624    0.5624    0.5624;
     0.5549    0.5549    0.5549;
     0.5474    0.5474    0.5474;
     0.5398    0.5398    0.5398;
     0.5323    0.5323    0.5323;
     0.5247    0.5247    0.5247;
     0.5172    0.5172    0.5172;
     0.5097    0.5097    0.5097;
     0.5021    0.5021    0.5021;
     0.4946    0.4946    0.4946;
     0.4871    0.4871    0.4871;
    0.4795    0.4795    0.4795;
     0.4720    0.4720    0.4720;
     0.4645    0.4645    0.4645;
     0.4569    0.4569    0.4569;
     0.4494    0.4494    0.4494;
     0.4419    0.4419    0.4419;
     0.4343    0.4343    0.4343;
     0.4268    0.4268    0.4268;
     0.4193    0.4193    0.4193;
     0.4117    0.4117    0.4117;
     0.4042    0.4042    0.4042;
     0.3967    0.3967    0.3967;
     0.3891    0.3891    0.3891;
     0.3816    0.3816    0.3816;
     0.3741    0.3741    0.3741;
     0.3665    0.3665    0.3665;
     0.3590    0.3590    0.3590;
     0.3515    0.3515    0.3515;
     0.3439    0.3439    0.3439;
     0.3364    0.3364    0.3364;
     0.3289    0.3289    0.3289;
     0.3213    0.3213    0.3213;
     0.2966    0.2966    0.2966;
     0.2719    0.2719    0.2719;
    0.2472    0.2472    0.2472;
     0.2225    0.2225    0.2225;
     0.1977    0.1977    0.1977;
     0.1730    0.1730    0.1730;
     0.1483    0.1483    0.1483;
     0.1236    0.1236    0.1236;
     0.0989    0.0989    0.0989;
     0.0742    0.0742    0.0742;
     0.0494    0.0494    0.0494;
     0.0247    0.0247    0.0247;
          0         0         0];
 colormap(color_bar)

% other colors used in the plot
%color1 = [134/255 134/255 134/255]; % - dark gray
color1 = [1 0.1 0.1];              % - red (original)

%color2 = [174/255 174/255 174/255]; % - middle gray
color2 = [1 0.55 0.55];            % - pink (original)

%color3 = [222/255 222/255 222/255]; % - light gray
color3 = [1 1 0.5];                % - yellow (original)

%color4 = color2;                    % - middle gray
color4 = [1 .55 0.3];              % - orange (original)

set(0,'DefaultTextFontName','Times');
set(0,'DefaultAxesFontName','Times');
set(0,'DefaultTextFontSize',18);
set(0,'DefaultAxesFontSize',18);

if nargin < 2
    error('Must input HPL and HPE!');
end
if nargin < 3
    HAL = 30;
end

% size of VPL, which should be the same for VPE as well
n = size(HPL,2);

% HPE
j = floor(2.0*abs(HPE))+1;
j(find(j>100)) = 100;
% HPL
k = floor(2.0*abs(HPL))+1;
k(find(k>100)) = 100;

% initialize hpl_stat
data     = zeros(100,100);
diagonal = zeros(100,1);
for i = 1:n
    % statistics
    data(k(i),j(i)) = data(k(i),j(i))+1;
    % diagonal
    if((k(i) == j(i)) && (abs(HPE(i,1)) < abs(HPL(1,i))))
        diagonal(k(i),1) = diagonal(k(i),1)+1;
    end
end

c = 1;
colors = 'brygcmw';
clf

err_bin  = 0.25:0.5:49.75;
sig_bin  = 0.25:0.5:49.75;

seconds  = n;
sec_available = n;
diag_cnt = sum(diagonal);

% determine the number of points and axis ranges
n_pts = sum(sum(data));
if sec_available == 1
    epochs = seconds;
else
    epochs = n_pts;
end

d_err_bin = mean(diff(err_bin));
x_lo_bnd  = min(err_bin) - d_err_bin/2;
x_up_bnd  = max(err_bin) + d_err_bin/2;

d_sig_bin = mean(diff(sig_bin));
y_lo_bnd  = min(sig_bin) - d_sig_bin/2;
y_up_bnd  = max(sig_bin) + d_sig_bin/2;

z_lo_bnd  = 1;
z_up_bnd  = max(max(data));

% clear plot
clf;

% plot each data point as a pixel
[i,j]=find(data);
face_mat=[[1 2 6 5]' [2 3 7 6]' [3 4 8 7]' ...
          [4 1 5 8]' [1 2 3 4]' [5 6 7 8]']';
colors=colormap;
for idx = 1:length(i)
  z = log10(data(i(idx),j(idx)));
  vtx_mat = [err_bin(j(idx))+[-0.5 0.5 0.5 -0.5 -0.5 0.5 0.5 -0.5]'*d_err_bin ...
             sig_bin(i(idx))+[-0.5 -0.5 0.5 0.5 -0.5 -0.5 0.5 0.5]'*d_sig_bin ...
             [0 0 0 0 z z z z]'];
  c_idx = ceil(63*(log10(data(i(idx),j(idx)))/log10(z_up_bnd))) + 1;
  patch('Vertices',  vtx_mat, ...
        'Faces',     face_mat, ...
        'FaceColor', colors(c_idx,:), ...
        'EdgeColor', 'none');
end

% determine availability and # of integrity failures
i_diag1 = find(err_bin == HAL);
i_diag2 = find(err_bin < HAL);
i_diag3 = find(err_bin > HAL);

i_fail1 = find(err_bin(j) >= HAL & sig_bin(i) < HAL);
n_fail1 = sum(sum(diag(data(i(i_fail1),j(i_fail1)))))...
         - sum(diagonal(i_diag1));
i_fail2 = find(err_bin(j)./sig_bin(i) >=1.0 & err_bin(j) < HAL);
n_fail2 = sum(sum(diag(data(i(i_fail2),j(i_fail2)))))...
         - sum(diagonal(i_diag2));
i_fail3 = find(err_bin(j)./sig_bin(i) >=1.0 & sig_bin(i) >= HAL);
n_fail3 = sum(sum(diag(data(i(i_fail3),j(i_fail3)))))...
         - sum(diagonal(i_diag3));
i_cont  = find(sig_bin(i) >= HAL);
n_cont  = sum(sum(diag(data(i(i_cont),j(i_cont)))));
%i_avail = find(err_bin(j) < HAL & sig_bin(i) < HAL);
i_avail = find(err_bin(j)./sig_bin(i) < 1.0 & sig_bin(i) < HAL);
n_avail = sum(sum(diag(data(i(i_avail),j(i_avail)))))...
         + sum(diagonal(i_diag2));

% set the axes limits and color values
set(gca,'XLim',[x_lo_bnd x_up_bnd]);
set(gca,'YLim',[y_lo_bnd y_up_bnd]);
set(gca,'CLim',[z_lo_bnd z_up_bnd]);

% show the region of normal operation
HT=text(0.37*(HAL - x_lo_bnd) + x_lo_bnd, ...
     0.93*(HAL - y_lo_bnd) + y_lo_bnd, ...
     0.25*sqrt(z_up_bnd*z_lo_bnd), ...
     '{\fontsize{14pt}Normal Operation}');
pos = get(HT, 'Position');
pos = pos + [.05 -0.05 0];
set(HT, 'Position', pos);
set(HT, 'Color', [1 1 1]); % white
text(0.37*(HAL - x_lo_bnd) + x_lo_bnd, ...
     0.93*(HAL - y_lo_bnd) + y_lo_bnd, ...
     0.25*sqrt(z_up_bnd*z_lo_bnd), ...
     '{\fontsize{14pt}Normal Operation}');
if n_avail/epochs >= .999995
  HT = text(0.33*(HAL - x_lo_bnd) + x_lo_bnd, ...
            0.86*(HAL - y_lo_bnd) + y_lo_bnd, ...
            0.25*sqrt(z_up_bnd*z_lo_bnd), ...
            '> 99.999%');
  pos = get(HT, 'Position');
  pos = pos + [.05 -0.05 0];
  set(HT, 'Position', pos);
  set(HT, 'Color', [1 1 1]); % white
  set(HT, 'FontSize', 14);
  HT = text(0.33*(HAL - x_lo_bnd) + x_lo_bnd, ...
            0.86*(HAL - y_lo_bnd) + y_lo_bnd, ...
            0.25*sqrt(z_up_bnd*z_lo_bnd), ...
            '> 99.999%');
  set(HT, 'FontSize', 14);
else
  HT = text(0.37*(HAL - x_lo_bnd) + x_lo_bnd, ...
            0.86*(HAL - y_lo_bnd) + y_lo_bnd, ...
            0.25*sqrt(z_up_bnd*z_lo_bnd), ...
            [num2str(100.0*n_avail/epochs,'%6.3f'), '%']);
  pos = get(HT, 'Position');
  pos = pos + [.05 -0.05 0];
  set(HT, 'Position', pos);
  set(HT, 'Color', [1 1 1]); % white
  set(HT, 'FontSize', 14);
  HT = text(0.37*(HAL - x_lo_bnd) + x_lo_bnd, ...
            0.86*(HAL - y_lo_bnd) + y_lo_bnd, ...
            0.25*sqrt(z_up_bnd*z_lo_bnd), ...
            [num2str(100.0*n_avail/epochs,'%6.3f'), '%']);
  set(HT, 'FontSize', 14);
end

% outline the region of integrity failures
patch([HAL HAL x_up_bnd x_up_bnd], ...
      [y_lo_bnd HAL HAL y_lo_bnd], ...
      -[0.5 0.5 0.5 0.5], ...
      color1);
HT = text(0.50*(x_up_bnd - HAL) + HAL, ...
          0.55*(HAL - y_lo_bnd) + y_lo_bnd, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          'HMI');
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
pos = get(HT, 'Position');
pos = pos + [.05 -0.05 0];
set(HT, 'Position', pos);
set(HT, 'Color', color1);
HT = text(0.50*(x_up_bnd - HAL) + HAL, ...
          0.55*(HAL - y_lo_bnd) + y_lo_bnd, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          'HMI');
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
HT = text(0.50*(x_up_bnd - HAL) + HAL, ...
          0.45*(HAL - y_lo_bnd) + y_lo_bnd, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          ['epochs: ', int2str(n_fail1)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
pos = get(HT, 'Position');
pos = pos + [.05 -0.05 0];
set(HT, 'Position', pos);
set(HT, 'Color', color1);
HT = text(0.50*(x_up_bnd - HAL) + HAL, ...
          0.45*(HAL - y_lo_bnd) + y_lo_bnd, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          ['epochs: ', int2str(n_fail1)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');

% outline the region of HPL failures
patch([x_lo_bnd HAL HAL], ...
      [y_lo_bnd HAL y_lo_bnd], ...
      -[0.5 0.5 0.5], ...
      color2);
HT = text(0.67*(HAL - x_lo_bnd) + x_lo_bnd, ...
          0.35*(HAL - y_lo_bnd) + y_lo_bnd, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          'MI');
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
pos = get(HT, 'Position');
pos = pos + [.05 -0.05 0];
set(HT, 'Position', pos);
set(HT, 'Color', color2);
HT = text(0.67*(HAL - x_lo_bnd) + x_lo_bnd, ...
          0.35*(HAL - y_lo_bnd) + y_lo_bnd, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          'MI');
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
HT = text(0.67*(HAL - x_lo_bnd) + x_lo_bnd, ...
          0.25*(HAL - y_lo_bnd) + y_lo_bnd, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          ['epochs: ', int2str(n_fail2)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
pos = get(HT, 'Position');
pos = pos + [.05 -0.05 0];
set(HT, 'Position', pos);
set(HT, 'Color', color2);
HT = text(0.67*(HAL - x_lo_bnd) + x_lo_bnd, ...
          0.25*(HAL - y_lo_bnd) + y_lo_bnd, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          ['epochs: ', int2str(n_fail2)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');

% outline the region of unavailability
patch([x_lo_bnd x_up_bnd x_up_bnd x_lo_bnd], ...
      [HAL HAL y_up_bnd y_up_bnd], ...
      -[0.5 0.5 0.5 0.5], ...
      color3);
HT = text(0.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, ...
          0.70*(y_up_bnd - HAL) + HAL, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          '{\fontsize{14pt}System Unavailable}');
set(HT,'HorizontalAlignment','Center');
pos = get(HT, 'Position');
pos = pos + [.05 -0.05 0];
set(HT, 'Position', pos);
set(HT, 'Color', color3);
HT = text(0.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, ...
          0.70*(y_up_bnd - HAL) + HAL, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          '{\fontsize{14pt}System Unavailable}');
set(HT,'HorizontalAlignment','Center');
HT = text(0.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, ...
          0.55*(y_up_bnd - HAL) + HAL, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          ['Alarm Epochs: ', int2str(n_cont)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
pos = get(HT, 'Position');
pos = pos + [.05 -0.05 0];
set(HT, 'Position', pos);
set(HT, 'Color', color3);
HT = text(0.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, ...
          0.55*(y_up_bnd - HAL) + HAL, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          ['Alarm Epochs: ', int2str(n_cont)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');

% outline the region where integrity failures and unavailability overlap
patch([HAL x_up_bnd x_up_bnd], ...
      [HAL y_up_bnd HAL], ...
      -z_lo_bnd*[0.45 0.45 0.45], ...
      color4);
HT = text(0.70*(x_up_bnd - HAL) + HAL, ...
          0.325*(y_up_bnd - HAL) + HAL, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          'MI');
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
pos = get(HT, 'Position');
pos = pos + [.05 -0.05 0];
set(HT, 'Position', pos);
set(HT, 'Color', color4);
HT = text(0.70*(x_up_bnd - HAL) + HAL, ...
          0.325*(y_up_bnd - HAL) + HAL, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          'MI');
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
HT = text(0.70*(x_up_bnd - HAL) + HAL, ...
          0.175*(y_up_bnd - HAL) + HAL, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          ['epochs: ', int2str(n_fail3)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');
pos = get(HT, 'Position');
pos = pos + [.05 -0.05 0];
set(HT, 'Position', pos);
set(HT, 'Color', color4);
HT = text(0.70*(x_up_bnd - HAL) + HAL, ...
          0.175*(y_up_bnd - HAL) + HAL, ...
          0.25*sqrt(z_up_bnd*z_lo_bnd), ...
          ['epochs: ', int2str(n_fail3)]);
set(HT, 'FontSize', 14);
set(HT,'HorizontalAlignment','Center');

hold on;
grid off;
% make the grid visible over the patch regions (integrity, unavailability)
ytick  = get(gca,'YTick');
xtick  = get(gca,'XTick');
nytick = length(ytick);
nxtick = length(xtick);
for i = 1:nytick
    plot3([x_lo_bnd x_up_bnd], ytick(i)*[1 1], [0.6 0.6], 'k:');
end
for i = 1:nxtick
    plot3(xtick(i)*[1 1], [y_lo_bnd y_up_bnd], [0.6 0.6], 'k:');
end

% label the axes and add a title
xlabel('{\fontsize{16pt}Error [m]}');
set(get(gca,'XLabel'),'Position', ...
    [ 0.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, ...
     -0.05*(y_up_bnd - y_lo_bnd) + y_lo_bnd, ...
     z_lo_bnd]);

ylabel(['{\fontsize{16pt}HPL_{',src_name,'} [m]}']);
set(get(gca,'YLabel'),'Position', ...
    [-0.06*(x_up_bnd - x_lo_bnd) + x_lo_bnd, ...
      0.50*(y_up_bnd - y_lo_bnd) + y_lo_bnd z_lo_bnd]);

title (['Horizontal Performance [',int2str(seconds),' seconds]']);

set(get(gca,'Title'),'FontSize',14);
set(get(gca,'Title'),'Position', ...
    [0.50*(x_up_bnd - x_lo_bnd) + x_lo_bnd, ...
     1.02*(y_up_bnd - y_lo_bnd) + y_lo_bnd, ...
     z_lo_bnd]);

% put the sigma_H_major scale on the right hand y-axis
% K_h_pa = 6.18;
% for i = ceil(y_lo_bnd/K_h_pa):floor(y_up_bnd/K_h_pa)
%   if abs(i*K_h_pa - (0.5*(y_up_bnd - y_lo_bnd) + y_lo_bnd)) > ...
%      0.05*(y_up_bnd - y_lo_bnd)
%     plot3([.99*(x_up_bnd - x_lo_bnd) + x_lo_bnd x_up_bnd], ...
%           i*K_h_pa*[1 1], ...
%           [0.65 0.65], ...
%           'k');
%     text(1.02*(x_up_bnd - x_lo_bnd) + x_lo_bnd, i*K_h_pa, 0.65, int2str(i));
%   end
% end
% HT = text(1.04*(x_up_bnd - x_lo_bnd) + x_lo_bnd, ...
%           0.50*(y_up_bnd - y_lo_bnd) + y_lo_bnd, ...
%           0.65, ...
%           '{\fontsize{14pt}\sigma}_{H_{major}} (m)');
% set(HT,'HorizontalAlignment','Center');
% set(HT,'Rotation',90);

% put the color scale up on the right hand side
H = colorbar('vert');
set(get(H,'Ylabel'),'String','{\fontsize{14pt}Number of Points per Pixel}');
%set(H,'YScale','log');
set(H,'YLim',[z_lo_bnd z_up_bnd]);
%set(H,'CLim', [z_lo_bnd z_up_bnd]);

if sec_available == 1
    display (['EGNOS differential for ', ...
              num2str(n_pts), ...
              ' out of ', ...
              num2str(seconds), ...
              ' seconds']);
end

set(H, 'linewidth',.8);
hold off
set(gca,'Fontsize',18);

%print -deps env_HPL

%%%%%%%%%%%% end hplstat.m  %%%%%%%%%%%%%%%%%%%%

author: yandld


text:
function [airway_mapped_image, airway_tree_root] = PTKMapAirwayCentrelineToImage(centreline_results, airway_image)
    % PTKMapAirwayCentrelineToImage.
    %
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2014.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    airway_mapped_image_raw = zeros(airway_image.ImageSize, 'uint16');
    airway_mapped_image = airway_image.BlankCopy;

    airway_tree_root = centreline_results.AirwayCentrelineTree;
    centreline_bronchi_to_do = CoreStack(airway_tree_root);

    bronchus_index = uint16(1);

    number_of_branches = airway_tree_root.CountBranches;

    parent_map = cell(number_of_branches, 1);
    child_map = cell(number_of_branches, 1);

    % Assign a label to each centreline bronchus, and mark the label
    % image with that index at each centreline voxel
    while ~centreline_bronchi_to_do.IsEmpty
        next_centreline_bronchus = centreline_bronchi_to_do.Pop;
        voxels = PTKTreeUtilities.GetCentrelineVoxelsForTheseBranches(next_centreline_bronchus, airway_image);
        airway_mapped_image_raw(voxels) = bronchus_index;
        next_centreline_bronchus.BronchusIndex = bronchus_index;

        % Add parent index to this branch, and add this branch index to
        % parent's child indices
        if ~isempty(next_centreline_bronchus.Parent)
            parent = next_centreline_bronchus.Parent;
            parent_index = parent.BronchusIndex;
            parent_map{bronchus_index} = parent_index;
            child_map{parent_index} = [child_map{parent_index}, bronchus_index];
        end

        centreline_bronchi_to_do.Push(next_centreline_bronchus.Children);
        bronchus_index = bronchus_index + 1;
    end

    % Find the nearest centreline point for every voxel in the airway
    % segmentation, and assign every voxel to that label
    [~, nearest_centreline_index] = bwdist(airway_mapped_image_raw > 0);
    airway_mapped_image_raw(:) = airway_mapped_image_raw(nearest_centreline_index(:));
    airway_mapped_image_raw(airway_image.RawImage ~= 1) = 0;
    airway_mapped_image.ChangeRawImage(airway_mapped_image_raw);

    airway_mapped_image.ChangeColorLabelParentChildMap(parent_map, child_map)
end

author: tomdoel

repo: pulmonarytoolkit


text:
function RIR = genMultiChannelRIR(t60, useGPU, filename)
if nargin<2
    useGPU = 0;
end

% the following geometry parameters can be randomly sampled from a
% distribution. For example, you want to cover different room sizes, just
% sample the length, width, and height of the room with a reasonable
% distribution.
nCh = 8;    % number of microphones in the array
room_size = [10 5 3];   % length, width, and height of the room
ctr_pos = [2 2 1];      % center position of the array
src_pos = [5 3 2.5];    % coordinate of the sound source
[mic_pos] = genMicPositionsWrtCenter(0);        % get the microphones' positions with respect to the center of the array
mic_pos = bsxfun(@plus, mic_pos, ctr_pos(:));   % get microphones' positions

% set up the settings
SetupStruc.Fs = 16000;                  % sampling frequency in Hz
SetupStruc.c = 343;                     % (optional) propagation speed of acoustic waves in m/s
SetupStruc.T60 = t60;                   % reverberation time T60, or define a T20 field instead!
SetupStruc.room = room_size;
SetupStruc.reflect_weights = ones(1,6) * exp(-13.82/sum(1./SetupStruc.room(1,:)*SetupStruc.c*t60));     %calculate the reflection coefficient.
                                        % It's good to first decide the T60 and room size, and then compute the corresponding to reflection rate.
SetupStruc.mic_pos = mic_pos';
SetupStruc.src_traj = src_pos;

% we are using a modified version of the ISM_RIR_bank function from Eric A. Lehmann
% The main change is that it now supports the use of GPU.
[RIR_cell] = ISM_RIR_bank_GPU(SetupStruc,'', 'useGPU', useGPU, 'SilentFlag', 1);     % generating RIR can be very slow for long T60. Using GPU makes it much faster
RIR = cell2mat(RIR_cell')';

if nargin==3
    wavReader = BinaryReader(nCh, 'int16');                 % a binary read/write object
    wavReader.write(filename, StoreWavInt16(RIR));          % we can write the RIR into binary format for later access
end

end

author: singaxiong

repo: SignalGraph


text:
function [data_train, labels_train, vid_ids_train_string, data_devel, labels_devel, vid_ids_devel_string, raw_devel, PC, means_norm, stds_norm, success_devel] = ...
    Prepare_HOG_AU_data_generic_intensity(train_users, devel_users, au_train, bp4d_dir, features_dir)

%%
addpath(genpath('../data extraction/'));

% First extracting the labels
[ labels_train, valid_ids_train, vid_ids_train ] = extract_BP4D_labels_intensity(bp4d_dir, train_users, au_train);
au_other = setdiff([6, 10, 12, 14, 17], au_train);
[ labels_other, ~, ~ ] = extract_BP4D_labels_intensity(bp4d_dir, train_users, au_other);
labels_other = cat(1, labels_other{:});

train_geom_data = Read_geom_files(train_users, features_dir);

% Reading in the HOG data (of only relevant frames)
[train_appearance_data, valid_ids_train_hog, vid_ids_train_string] = Read_HOG_files(train_users, features_dir);
train_appearance_data = cat(2, train_appearance_data, train_geom_data);

% Subsample the data to make training quicker
labels_train = cat(1, labels_train{:});
valid_ids_train = logical(cat(1, valid_ids_train{:}));

reduced_inds = false(size(labels_train,1),1);
reduced_inds(labels_train > 0) = true;

% make sure the same number of positive and negative samples is taken
pos_count = sum(labels_train > 0);
neg_count = sum(labels_train == 0);

num_other = floor(pos_count / (size(labels_other, 2)));

inds_all = 1:size(labels_train,1);

if(numel(train_users) > 0)
    if(numel(au_train) == 1)
        for i=1:size(labels_other, 2)+1

            if(i > size(labels_other, 2))
                % fill the rest with a proportion of neutral
                inds_other = inds_all(sum(labels_other,2)==0 & ~labels_train );
                num_other_i = min(numel(inds_other), pos_count - sum(labels_train(reduced_inds,:)==0));
            else
                % take a proportion of each other AU
                inds_other = inds_all(labels_other(:, i) & ~labels_train );
                num_other_i = min(numel(inds_other), num_other);
            end
            inds_other_to_keep = inds_other(round(linspace(1, numel(inds_other), num_other_i)));
            reduced_inds(inds_other_to_keep) = true;

        end
    end
    % Remove invalid ids based on CLM failing or AU not being labelled
    reduced_inds(~valid_ids_train) = false;
    reduced_inds(~valid_ids_train_hog) = false;

    labels_other = labels_other(reduced_inds, :);
    labels_train = labels_train(reduced_inds,:);
    train_appearance_data = train_appearance_data(reduced_inds,:);
    vid_ids_train_string = vid_ids_train_string(reduced_inds,:);

end

%% Extract devel data

% First extracting the labels
[ labels_devel, valid_ids_devel, vid_ids_devel ] = extract_BP4D_labels_intensity(bp4d_dir, devel_users, au_train);

devel_geom_data = Read_geom_files(devel_users, features_dir);
% Reading in the HOG data (of only relevant frames)
[devel_appearance_data, valid_ids_devel_hog, vid_ids_devel_string] = Read_HOG_files(devel_users, features_dir);
devel_appearance_data = cat(2, devel_appearance_data, devel_geom_data);

valid_ids_devel = logical(cat(1, valid_ids_devel{:}));

labels_devel = cat(1, labels_devel{:});

success_devel = valid_ids_devel;

% Peforming zone specific masking
if(au_train < 8 || au_train == 43 || au_train == 45) % upper face AUs ignore bottom face
    % normalise the data
    pca_file = '../../pca_generation/generic_face_upper.mat';
    load(pca_file);
elseif(au_train > 9) % lower face AUs ignore upper face and the sides
    % normalise the data
    pca_file = '../../pca_generation/generic_face_lower.mat';
    load(pca_file);
elseif(au_train == 9) % Central face model
    % normalise the data
    pca_file = '../../pca_generation/generic_face_rigid.mat';
    load(pca_file);
end

PC_n = zeros(size(PC)+size(devel_geom_data, 2));
PC_n(1:size(PC,1), 1:size(PC,2)) = PC;
PC_n(size(PC,1)+1:end, size(PC,2)+1:end) = eye(size(devel_geom_data, 2));
PC = PC_n;

means_norm = cat(2, means_norm, zeros(1, size(devel_geom_data,2)));
stds_norm = cat(2, stds_norm, ones(1, size(devel_geom_data,2)));

% Grab all data for validation as want good params for all the data
raw_devel = devel_appearance_data;

devel_appearance_data = bsxfun(@times, bsxfun(@plus, devel_appearance_data, -means_norm), 1./stds_norm);
data_devel = devel_appearance_data * PC;

if(numel(train_users) > 0)
    train_appearance_data = bsxfun(@times, bsxfun(@plus, train_appearance_data, -means_norm), 1./stds_norm);
    data_train = train_appearance_data * PC;
else
    data_train = [];
end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function pts = rotate(pts, phi, theta, psi)
% ROTATE - Apply rotation to a vector of points

% Define rotation matrix (right handed)
R_roll = [...
    1, 0, 0;...
    0, cos(phi), sin(phi);...
    0, -sin(phi), cos(phi)];
R_pitch = [...
    cos(theta), 0, -sin(theta);...
    0, 1, 0;...
    sin(theta), 0, cos(theta)];
R_yaw = [...
    cos(psi), sin(psi), 0;...
    -sin(psi), cos(psi), 0;...
    0, 0, 1];
R = R_roll * R_pitch * R_yaw;   % inertial to body
R = R';  % body to inertial

% Rotate vertices
pts = R*pts;

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
classdef AbstractGeometricVector < matlab.mixin.SetGet & matlab.mixin.Heterogeneous
    %AbstractGeometricVector Summary of this class goes here
    %   Detailed explanation goes here

    properties

    end

    methods
        vect = getVectorAtTime(obj, time, vehElemSet, inFrame)

        name = getName(obj)

        setName(obj, name)

        listboxStr = getListboxStr(obj)

        useTf = openEditDialog(obj)

        tf = isVehDependent(obj)

        origin = getOriginPointInViewFrame(obj, time, vehElemSet, viewFrame)

        tf = usesGroundObj(obj, groundObj)

        tf = usesGeometricPoint(obj, point)

        tf = usesGeometricVector(obj, vector)

        tf = usesGeometricCoordSys(obj, coordSys)

        tf = usesGeometricRefFrame(obj, refFrame)

        tf = usesGeometricAngle(~, ~)

        tf = usesGeometricPlane(~, ~)

        tf = isInUse(obj, lvdData)
    end

    methods(Sealed)
        function tf = eq(a,b)
            tf = eq@handle(a,b);
        end

        function tf = ne(a,b)
            tf = ne@handle(a,b);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function Lms=normapx(Ns, delta, epsilon)
K=sqrt(delta*(1-delta))*norminv(epsilon, 0, 1);
Lms = Ns.*(1-delta) + K.*sqrt(Ns);

author: yp-mit

repo: spectre

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


text:
function [ h, m, s ] = frac_to_hms ( f )

%*****************************************************************************80
%
%% FRAC_TO_HMS converts a fractional day into hours, minutes, seconds.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real F, a day fraction between 0.0 and 1.0.
%
%    Output, integer H, M, S, the equivalent hours, minutes
%    and seconds.
%
  f2 = f;

  f2 = 24.0 * f2;
  h = floor ( f2 );
  f2 = f2 - h;

  f2 = 60.0 * f2;
  m = floor ( f2 );
  f2 = f2 - m;

  f2 = 60.0 * f2;
  s = floor ( f2 );
  f2 = f2 - s;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [item, inputs] = fillvals(item, inputs, infcn)

% function [item, inputs] = fillvals(item, inputs, infcn)
% If ~all_set_item, try to set item.val{1} to inputs{1}. Validity checks
% are performed through subsasgn. If inputs{1} is not suitable for this
% item, it is discarded. If infcn is a function handle,
% [val sts] = infcn(item)
% will be called to obtain a value for this item. This call will be
% repeated until either val can be assigned to item or sts is true. val
% should be a cell array with 1 item and val{1} the value that is to be
% assigned to item.val{1}.
% This function is suitable for all cfg_leaf items.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: fillvals.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $'; %#ok

if ~all_set_item(item)
    if ~isempty(inputs)
        item = setval(item, inputs{1}, false);
        inputs = inputs(2:end);
    end;
    if ~all_set_item(item) && ~isempty(infcn) && subsasgn_check_funhandle(infcn)
        sts = false;
        while ~sts && ~all_set_item(item)
            [val, sts] = feval(infcn, item);
            if sts
                item = setval(item, val, false);
            end;
        end;
    end;
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [kW] = MW2kW(MW)
% Convert power from megawatts to kilowatts.
% Chad A. Greene 2012
kW = MW*1000;

author: Sable

repo: mcbench-benchmarks


text:
% pixCon:  Return a matrix of pixel connections in an image or video
%
% Note:  The mex-file version runs much faster than the m-file.
%
% Usage:
%   function pc = pixCon(sz,mode)

function pc = pixCon(sz,mode)

npix = prod(sz);
ndim = length(sz);
rd = cumprod([1 sz(1:end-1)]);
d = [-rd rd];
%pc = spdiags(ones(npix,2*ndim),d,npix,npix);
nbr = zeros(npix,2*ndim);
dz = cumprod(sz(1:end-1));
for i = 1:ndim
    dnbr = ones(sz);
    [subs{1:ndim}] = deal(':');
    subs{i} = sz(i);
    subsasgn(dnbr,struct('type','()','subs',{subs}),0);
    nbr(:,i) = dnbr(:);
    dnbr = ones(sz);
    [subs{1:ndim}] = deal(':');
    subs{i} = 1;
    subsasgn(dnbr,struct('type','()','subs',{subs}),0);
    nbr(:,ndim+i) = dnbr(:);
end;
pc = spdiags(nbr,d,npix,npix);


function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
%function [ydum,xdum,breaks]=varprior(nv,nx,lags,mnprior,vprior)
% ydum, xdum:   dummy observation data that implement the prior
% breaks:       vector of points in the dummy data after which new dummy obs's start
%                   Set breaks=T+[0;breaks], ydata=[ydata;ydum], xdum=[xdata;xdum], where
%                   actual data matrix has T rows, in preparing input for rfvar3
% nv,nx,lags: VAR dimensions
% mnprior.tight:Overall tightness of Minnesota prior
% mnprior.decay:Standard deviations of lags shrink as lag^(-decay)
% vprior.sig:   Vector of prior modes for diagonal elements of r.f. covariance matrix
% vprior.w:     Weight on prior on vcv.  1 corresponds to "one dummy observation" weight
%                   Should be an integer, and will be rounded if not.  vprior.sig is needed
%                   to scale the Minnesota prior, even if the prior on sigma is not used itself.
%                   Set vprior.w=0 to achieve this.
% Note:         The original Minnesota prior treats own lags asymmetrically, and therefore
%                   cannot be implemented entirely with dummy observations.  It is also usually
%                   taken to include the sum-of-coefficients and co-persistence components
%                   that are implemented directly in rfvar3.m.  The diagonal prior on v, combined
%                   with sum-of-coefficients and co-persistence components and with the unit own-first-lag
%                   prior mean generates larger prior variances for own than for cross-effects even in
%                   this formulation, but here there is no way to shrink toward a set of unconstrained
%                   univariate AR's.

% Original file downloaded from:
% http://sims.princeton.edu/yftp/VARtools/matlab/varprior.m

if ~isempty(mnprior)
    xdum = zeros(lags+1,nx,lags,nv);
    ydum = zeros(lags+1,nv,lags,nv);
    for il = 1:lags
        ydum(il+1,:,il,:) = il^mnprior.decay*diag(vprior.sig);
    end
    ydum(1,:,1,:) = diag(vprior.sig);
    ydum = mnprior.tight*reshape(ydum,[lags+1,nv,lags*nv]);
    ydum = flipdim(ydum,1);
    xdum = mnprior.tight*reshape(xdum,[lags+1,nx,lags*nv]);
    xdum = flipdim(xdum,1);
    breaks = (lags+1)*[1:(nv*lags)]';
    lbreak = breaks(end);
else
    ydum = [];
    xdum = [];
    breaks = [];
    lbreak = 0;
end
% elle = 0;
% for j1 = 1 : nv
%     for j2  =  1 : lags
%         elle = 1 + elle;
%         ydum(:,:,elle) = ydum(:,:,elle) * mnprior.unit_root_(j1);
%     end
% end
if ~isempty(vprior) && vprior.w>0
    ydum2 = zeros(lags+1,nv,nv);
    xdum2 = zeros(lags+1,nx,nv);
    ydum2(end,:,:) = diag(vprior.sig);
    for i = 1:vprior.w
        ydum = cat(3,ydum,ydum2);
        xdum = cat(3,xdum,xdum2);
        breaks = [breaks;(lags+1)*[1:nv]'+lbreak];
        lbreak = breaks(end);
    end
end
dimy = size(ydum);
ydum = reshape(permute(ydum,[1 3 2]),dimy(1)*dimy(3),nv);
xdum = reshape(permute(xdum,[1 3 2]),dimy(1)*dimy(3),nx);
breaks = breaks(1:(end-1));


text:
function [ y2, m2, d2, f2 ] = ymdf_next_islamic ( y1, m1, d1, f1 )

%*****************************************************************************80
%
%% YMDF_NEXT_ISLAMIC returns the Islamic YMDF date of the next day.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Output, integer Y2, M2, D2, real F2,
%    tomorrow's YMDF date.
%
  y2 = y1;
  m2 = m1;
  d2 = d1 + 1;
  f2 = f1;

  [ y2, m2, d2 ] = day_carry_islamic ( y2, m2, d2 );

  return
end


text:
function interleaved_vec = interleave_vectors(vectors)
    % Interleave multiple vectors element-by-element
    num_vectors = length(vectors);
    element_count = zeros(1, num_vectors);
    for i = 1:num_vectors
        element_count(i) = numel(vectors{i});
    end

    temp = cell(num_vectors, max(element_count));
    for i = 1:num_vectors
        temp(i, 1:element_count(i)) = num2cell(vectors{i});
    end

    interleaved_vec = [temp{:}];
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function F=PolFreq(n, X, Y,Z)
F=zeros(n,n,n);
for j =1:n
  for i=1:n
    for k=1:n
    F(X(i,k,j),Y(i,k,j),Z(i,k,j))=F(X(i,k,j),Y(i,k,j),Z(i,k,j))+1;
    end
  end
end

author: thomaskuestner


text:
function [varargout] = file_readBV(file, varargin)
% FILE_READBV - load EEG data which is stored in BrainVision format.
%                  C-functions are used for better performance. Use the
%                  slower FILE_LOADBV if this function does not work.
%
% Synopsis:
%   [CNT, MRK, HDR]= file_readBV(FILE, 'Property1',Value1, ...)
%
% Arguments:
%   FILE: file name (no extension),
%         relative to BTB.RawDir unless beginning with '/' (resp '\').
%         FILE may also contain the wildcard symbol '*'. In this case
%         make sure that the order of the files (printed to the terminal)
%         is appropriate.
%         FILE may also be a cell array of file names.
%
% Properties:
%   'CLab': Channels to load (labels or indices). Default all
%           (which can be explicitly specified by [])
%   'Fs': Sampling interval, must be an integer divisor of the
%         sampling interval of the raw data. fs may also be 'raw' which means
%         sampling rate of raw signals. Default: 'raw'.
%   'Ival': Interval to read, [start end] in msec. It is not checked
%           whether the whole interval could be loaded, or the file is shorter.
%   'IvalSa': Same as 'ival' but [start end] in samples of the downsampled data.
%   'Start': Start [msec] reading.
%   'MaxLen': Maximum length [msec] to be read.
%   'Filt': Filter to be applied to raw data before subsampling.
%           opt.Filt must be a struct with fields 'b' and 'a' (as used for
%           the Matlab function filter).
%           Note that using opt.Filt may slow down loading considerably.
%   'SubsamplePolicy': Function that is used for subsampling after filtering,
%           specified as as string or a vector.
%           Default 'subsampleByMean'. Other 'subsampleByLag'
%           If you specify a vector it has to be the same size as lag.
%   'LinearDerivation' : for creating bipolar channels (see
%   procutil_biplist2projection for details)
%
% Remark:
%   Properties 'Ival' and 'Start'+'MaxLen' are exclusive, i.e., you may only
%   specify one of them.
%
% Returns:
%   CNT: struct for contiuous signals
%        .x: EEG signals (time x channels)
%        .clab: channel labels
%        .fs: sampling interval
%        .scale: if this field is given, the real data are .x*.scale,
%           and .x is int
%   MRK: struct of marker information
%   HDR: struct of header information
%
% TODO: The function so far can only read specific formats, e.g.,
%       multiplexed, INT_16, ... The function does not even check, whether
%       the file is in this format!
%
% See also: file_* procutil_biplist2projection
%
%Hints:
% A low pass filter to get rid of the line noise can be designed as follows:
%  hdr= file_readBVheader(file);
%  Wps= [40 49]/hdr.fs*2;
%  [n, Ws]= cheb2ord(Wps(1), Wps(2), 3, 50);
%  [filt.b, filt.a]= cheby2(n, 50, Ws);
% You can also use filtdemo to design your own filters.

% Benjamin Blankertz

%   2008/06/20/ - Max Sagebaum
%               - refactored file_loadBV to read the data with read_bv.c
%               - removed code fragments marked as obsolete
%               - you can now use the ival option when concatinating
%                 multiple eeg files
%   2008/06/17  - Max Sagebaum
%               - the iir filter was not properly send to read_bv
%   2010/09/09  - Max Sagebaum
%               - There was an bug in the check for the lag

%% check if the mex file is present
readBV_status = exist('read_bv','file');
if not(readBV_status == 3)
    warning('Could not detect mex files for read_bv! Using file_loadBV instead.')
    varargout= cell(1, nargout);
    [varargout{:}]= file_loadBV(file, varargin{:});
    return;
end

%% start file_readBV
global BTB

props= {'CLab'               ''       'CHAR|CELL{CHAR}'
        'Fs'                 'raw'    'CHAR|DOUBLE'
        'Start'              0        'DOUBLE[1]'
        'MaxLen'             inf      'DOUBLE[1]'
        'Prec'               0        'DOUBLE[1]'
        'Ival'               []       'DOUBLE[2]'
        'IvalSa'             []       'DOUBLE[2]'
        'SubsamplePolicy'    'mean'   'CHAR(mean lag)|DOUBLE'
        'Filt'               []       'STRUCT(a b)'
        'LinearDerivation'   []       'STRUCT'
        'TargetFormat'       'bbci'   'CHAR'
        'Verbose'            1        'BOOL'
       };

props_readBVmarkers= file_readBVmarkers;
props_readBVheader= file_readBVheader;
all_props= opt_catProps(props, props_readBVmarkers, props_readBVheader);

if nargin==0,
  varargout= {all_props};
  return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, all_props);

misc_checkType(file, 'CHAR|CELL{CHAR}');

if ~isempty(opt.Ival),
  if ~isdefault.Start || ~isdefault.MaxLen,
    error('specify either ''Ival'' or ''Start''+''MaxLen'' but not both');
  end
  if ~isempty(opt.IvalSa),
    error('specify either ''Ival'' or ''IvalSa'' but not both');
  end
  opt.Start= opt.Ival(1);
  opt.MaxLen= diff(opt.Ival)+1;
end

%% read the headers an prepare the clab

if ~iscell(file)
  file = {file};
end

fileNames = cell(1,length(file));
fileTitle = cell(1,length(file));
% use BTB.RawDir as default dir
for filePos = 1:length(file)
  if fileutil_isAbsolutePath(file{filePos}),
    fileNames{filePos}= file{filePos};
    [dmy, fileTitle{filePos}]= fileparts(file{filePos});
  else
    fileNames{filePos} = fullfile(BTB.RawDir, file{filePos});
    fileTitle{filePos}= file{filePos};
  end
end

% get all files specified with the file object
fileNamesTemp = {};
for filePos = 1:length(file)

  if ischar(fileNames{filePos}) && ismember('*', fileNames{filePos},'legacy'),
    dd= dir([fileNames{filePos} '.eeg']);
    if isempty(dd),
      error('\nFile not found: %s\n', fileNames{filePos});
    end
    fc= cellfun(@(x)(x(1:end-4)), {dd.name}, 'UniformOutput',0);

    fileNamesTemp = cat(2,fileNamesTemp,strcat(fileparts(fileNames{filePos}), '/', fc));
  else
    fileNamesTemp = cat(2,fileNamesTemp,{fileNames{filePos}});
  end
end
fileNames = fileNamesTemp;
if length(fileNames)>1,
  if opt.Verbose,
    fprintf('concatenating files in the following order:\n');
    fprintf('  %s\n', fileNames{:});
  end
end

% now we read all headers and make some consistent checks if neeeded
hdr = cell(1,length(fileNames));
opt_readBVheader= opt_substruct(opt, props_readBVheader(:,1));
for filePos = 1:length(fileNames)
  hdr{filePos} = file_readBVheader(fileNames{filePos}, opt_readBVheader);

  % set the clabs and the raw_fs if we are in the loop for the first time
  if(filePos == 1)
    cnt.clab= hdr{filePos}.clab;
    raw_fs= hdr{filePos}.fs;
  end

  if ~isequal(cnt.clab, hdr{filePos}.clab),
    warning(['inconsistent clab structure will be repaired ' ...
             'by using the intersection']);
    cnt.clab = intersect(cnt.clab, hdr{filePos}.clab,'legacy');
  end
  if isequal(opt.Fs, 'raw')
    % if we want to read the raw data check if for each file the raw data
    % is the same
    if~isequal(raw_fs, hdr{filePos}.fs)
      error('inconsistent sampling rate');
    end
  else
    % if we have a specific fs check if for each file we have a positive
    % lag
    lag = hdr{filePos}.fs/opt.Fs;
    if lag~=round(lag) || lag<1,
      error('fs must be a positive integer divisor of every file''s fs');
    end
  end
end
clab_in_file= cnt.clab;

% select specified channels
if ~isempty(opt.CLab) && strcmp(opt.TargetFormat,'bbci'),
  cnt.clab= cnt.clab(util_chanind(cnt, opt.CLab));
end

% sort channels for memory efficient application of linear derivation:
% temporary channels are moved to the end
if ~isempty(opt.LinearDerivation),
  rm_clab= cell_flaten({opt.LinearDerivation.rm_clab});
  rmidx= util_chanind(cnt.clab, rm_clab);
  cnt.clab(rmidx)= [];
  cnt.clab= cat(2, cnt.clab, rm_clab);
end

%% prepare the output samples
firstFileToRead = 1;
firstFileSkip = 0;
lastFileToRead = length(fileNames);
lastFileLength = inf;

% check if we want to load the raw data
if isequal(opt.Fs, 'raw'),
  opt.Fs= raw_fs;
end
cnt.fs= opt.Fs;
cnt.title = str_vec2str(fileTitle);
cnt.file = str_vec2str(fileNames);

nChans= length(cnt.clab);

% get the skip and maxlen values for the data in samples for the new
% sampling rate
if ~isempty(opt.IvalSa),
  if ~isdefault.Start || ~isdefault.MaxLen,
    error('specify either <IvalSa> or <Start/MaxLen> but not both');
  end
  skip= opt.IvalSa(1)-1;
  maxlen = diff(opt.IvalSa)+1;
else
  skip= max(0, floor(opt.Start/1000*opt.Fs));
  maxlen = ceil(opt.MaxLen/1000*opt.Fs);
end

%get the number of samples for every file and check from which file we have
%to read
nSamples = 0;
dataSamples = 0;
dataSize = zeros(1,length(fileNames));
for filePos = 1:length(fileNames)
  % check if we can read the data with read_bv
  % currently only 16Bit Integers are supported
  switch hdr{filePos}.BinaryFormat
   case 'INT_16',
    cellSize= 2;
    readbv_binformat(filePos)=1;
   case 'INT_32',
    cellSize= 4;
    readbv_binformat(filePos)=2;
   case {'IEEE_FLOAT_32', 'FLOAT_32'},
    cellSize= 4;
    readbv_binformat(filePos)=3;
   case {'IEEE_FLOAT_64', 'FLOAT_64', 'DOUBLE'},
    cellSize= 8;
    readbv_binformat(filePos)=4;
   otherwise
    error('Precision %s not known.', hdr.BinaryFormat);
  end

  % open the file to get the size
  fid= fopen([fileNames{filePos} '.eeg'], 'r', hdr{filePos}.endian);
  if fid==-1, error('%s.eeg not found', fileNames{filePos}); end
  fseek(fid, 0, 'eof');
  fileLen= ftell(fid);
  fclose(fid);

  curChannels = length(hdr{filePos}.clab);
  curLag = hdr{filePos}.fs/opt.Fs;
  samples_in_file = floor(fileLen/(cellSize*curChannels));
  samples_after_subsample = floor(samples_in_file / curLag);
  dataSize(filePos) = samples_after_subsample;

  % set the new first file and the first data in this file
  if nSamples <= skip
    firstFileToRead = filePos;
    firstFileSkip = skip - nSamples;
    dataSamples = samples_after_subsample - firstFileSkip;
  else
    dataSamples = dataSamples + samples_after_subsample;
  end
  % advance to the end of the cur file
  nSamples = nSamples + samples_after_subsample;

  % if we reach the end set the last file and stop reading
  if nSamples >= (skip + maxlen)
    lastFileToRead = filePos;
    lastFileLength = samples_after_subsample - (nSamples - (skip + maxlen));
    dataSamples = dataSamples - (samples_after_subsample - lastFileLength);
     break;
  else
    % only if we have no maxlen
    lastFileLength = samples_after_subsample;
  end

end

%% reading the data
%create the data block for all samples
chosen_clab = cnt.clab;
cnt.x = zeros(dataSamples,nChans);
cnt.T = dataSize;

dataOffset = 0; % the offset for the current file
for filePos = firstFileToRead:lastFileToRead
  % get the channel id for this file
  chanids = util_chanind(clab_in_file,chosen_clab); % the -1 is for read_bv

  read_opt = struct('fs',cnt.fs, 'chanidx',chanids);

  if ~isempty(opt.Filt)
    read_opt.filt_b = opt.Filt.b;
    read_opt.filt_a = opt.Filt.a;
  end
  % set the subsample filter
  lag = hdr{filePos}.fs/opt.Fs;
  switch opt.SubsamplePolicy,
    case 'mean',
      read_opt.filt_subsample = ones(1,lag)/lag;
    case 'lag',
      read_opt.filt_subsample = [zeros(1,lag-1) 1];
    otherwise,
      read_opt.filt_subsample = opt.SubsamplePolicy;
  end

  read_hdr = struct('fs',hdr{filePos}.fs, ...
                    'nChans',hdr{filePos}.NumberOfChannels, ...
                    'scale',hdr{filePos}.scale, ...
                    'endian',hdr{filePos}.endian, ...
                    'BinaryFormat',readbv_binformat(filePos));

  % get the position for the data in the whole data set
  if firstFileToRead == filePos
    firstX = 1;
    firstData = firstFileSkip + 1;
  else
    firstX = lastX + 1;
    firstData = 1;
  end

  if lastFileToRead == filePos
    lastX = nSamples;
    lastData = lastFileLength;
  else
    lastX = firstX + dataSize(filePos) - 1 - firstFileSkip;
    lastData = dataSize(filePos);
  end

  read_opt.data = cnt.x;
  read_opt.dataPos = [firstX lastX firstData lastData] - 1;

  % read the data, read_bv will set the data in cnt.x because of the
  % read_opt.data options
  read_bv([fileNames{filePos} '.eeg'], read_hdr, read_opt);
  cnt.yUnit= hdr{filePos}.unit;

  %% Markers
  if nargout>1,
    opt_mrk= opt_substruct(opt, props_readBVmarkers(:,1));
    curmrk= file_readBVmarkers(fileNames{filePos}, opt_mrk);
    curmrk.time= curmrk.time + dataOffset*1000/cnt.fs;
    % find markers in the loaded interval
    inival= find(curmrk.time > skip*1000/cnt.fs & ...
                 curmrk.time <= (skip+maxlen)*1000/cnt.fs);
    % add special case: don't loose t=0 markers
    % NO: markers with time=0 make problems!
    %if skip==0,
    %  idxzero= find(curmrk.time==0);
    %  inival= [idxzero, inival];
    %end
    curmrk= mrk_selectEvents(curmrk, inival);
    %let the markers start at zero
    curmrk.time= curmrk.time - skip*1000/cnt.fs;

    if firstFileToRead == filePos
      mrk = curmrk;
    else
      mrk = mrk_mergeMarkers(mrk, curmrk);
    end
    dataOffset = dataOffset + dataSize(filePos);
  end
end
clear read_opt;

if ~isempty(opt.LinearDerivation),
  ld= opt.LinearDerivation;
  for cc= 1:length(ld),
    ci= util_chanind(cnt.clab, ld(cc).chan);
    support= find(ld(cc).filter);
    s2= util_chanind(cnt.clab, ld(cc).clab(support));
    cnt.x(:,ci)= cnt.x(:,s2) * ld(cc).filter(support);
    cnt.clab{ci}= ld(cc).new_clab;
  end
  % delete temporary channels: TODO in a memory efficient way
  idx= util_chanind(cnt, rm_clab);
  cnt.x(:,idx)= [];
  cnt.clab(idx)= [];
end

varargout= cell(1, nargout);

varargout{1}= cnt;
if nargout > 1,
  varargout{2} = mrk;
end
if nargout>2,
  if(1 == length(hdr))
    varargout{3} = hdr{1};
  else
      varargout{3}= hdr;
  end
end

author: bbci


text:
function OUT = MovAvg(DATA,window)
% =======================================================================
% Moving average of the matrix DATA (TxN). The moving average is
% computed down each column.
% =======================================================================
% OUT = MovAvg(DATA,window)
% -----------------------------------------------------------------------
% INPUT
%   DATA: T observations x N variables [double]
%   window: window of the moving average [double]
% -----------------------------------------------------------------------
% OUPUT
%   OUT: T x N matrix (the first windows-1 observations are NaN) [double]
% -----------------------------------------------------------------------
% EXAMPLE
%   X = rand(20,2);
%   OUT = MovAvg(X,2)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

if nargin<2,                error('Not enough input.'),          end
if window<=0,               error('Window must be positive.'),   end
if (window~=floor(window)), error('Window must be an integer.'), end

if min(size(DATA))==1
    DATA = DATA(:); % forces DATA to be a column vector
end

[nobs,nvar] = size(DATA);
if window>nobs
    error('window must not be greater than the length of DATA.')
end

temp=[];
for row=1:(nobs-window+1),
    temp = [temp; mean(DATA(row:(row+window-1),:))];
end

OUT = temp;
OUT = [nan(window-1,nvar); OUT]; % add nans to make conformable to original

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function dcmPatientCoordM = convertCerrToDcmPoints(cerrPointsM,scanNum,planC)
%
% function dcmPatientCoordM = convertCerrToDcmPoints(cerrPointsM,scanNum,planC)
%
% Example:
% global planC
% indexS = planC{end};
% scanNum = 1;
% cerrPointsM = planC{indexS.structures}(1).contour(35).segments(1).points;
% scanS = planC{indexS.scan}(scanNum);
% dcmPatientCoordM = convertCerrToDcmPoints(cerrPointsM,scanS);
%
% APA, 11/16/2021

if exist('planC','var')
    indexS = planC{end};
end
if isstruct(scanNum)
    scanS = scanNum;
else
    scanS = planC{indexS.scan}(scanNum);
end

% Convert contour points to DICOM r,c,s
[xs,ys,zs] = getScanXYZVals(scanS);
dx = xs(2)-xs(1);
dy = ys(2)-ys(1);
slice_distance = zs(2)-zs(1);
virPosMtx = [dx 0 0 xs(1);0 dy 0 ys(1); 0 0 slice_distance -zs(end); 0 0 0 1]; % (-)ve zs since CERR z is opposite of DICOM
cerrPointsM(:,3) = -cerrPointsM(:,3); % (-)ve since cerr z is opposite of DICOM
cerrPointsM = [cerrPointsM, ones(size(cerrPointsM,1),1)];
rcsM = virPosMtx \ cerrPointsM';

% Convert r,c,s to DICOM Patient coordinates (mm)
positionMatrix = getDicomImageToPatientAffineMat(scanS);
dcmPatientCoordM = positionMatrix * rcsM;
dcmPatientCoordM = dcmPatientCoordM(1:3,:)';

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%LAPLACIAN  Calculates the Laplacian of an image
%
%     dst = cv.Laplacian(src)
%     dst = cv.Laplacian(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ Source image.
%
% ## Output
% * __dst__ Destination image of the same size and the same number of channels
%   as `src`.
%
% ## Options
% * __DDepth__ Desired depth of the destination image. default -1, which uses
%   the same type as the input `src`.
% * __KSize__ Aperture size used to compute the second-derivative filters. See
%   cv.getDerivKernels for details. The size must be positive and odd.
%   default 1
% * __Scale__ Optional scale factor for the computed Laplacian values. By
%   default, no scaling is applied. See cv.getDerivKernels for details.
%   default 1
% * __Delta__ Optional delta value that is added to the results prior to
%   storing them in `dst`. default 0
% * __BorderType__ Pixel extrapolation method, see cv.copyMakeBorder.
%   Default 'Default'
%
% The function calculates the Laplacian of the source image by adding up the
% second x and y derivatives calculated using the Sobel operator:
%
%     dst = \Delta src = d^2(src)/dx^2 + d^2(src)/dy^2
%
% This is done when `KSize > 1`. When `KSize = 1`, the Laplacian is computed
% by filtering the image with the following 3x3 aperture:
%
%     [0, 1, 0; 1,-4, 1; 0, 1, 0]
%
% See also: cv.Sobel, cv.Scharr, imfilter, fspecial
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function p = prob_node(CPD, self_ev, pev)
% PROB_NODE Compute P(y|pa(y), theta) (tabular)
% p = prob_node(CPD, self_ev, pev)
%
% self_ev{m} is the evidence on this node in case m
% pev{i,m} is the evidence on the i'th parent in case m
% If there is a single case, self_ev can be a scalar instead of a cell array

ncases = size(pev, 2);

%assert(~any(isemptycell(pev))); % slow
%assert(~any(isemptycell(self_ev))); % slow

CPT = CPD_to_CPT(CPD);
sz = mysize(CPT);
nparents = length(sz)-1;
assert(nparents == size(pev, 1));

if ncases==1
  x = cat(1, pev{:});
  if iscell(y)
    y = self_ev{1};
  else
    y = self_ev;
  end
  switch nparents
   case 0, p = CPT(y);
   case 1, p = CPT(x(1), y);
   case 2, p = CPT(x(1), x(2), y);
   case 3, p = CPT(x(1), x(2), x(3), y);
   otherwise,
    ind = subv2ind(CPD.sizes, [x y]);
    p = CPT(ind);
  end
else
  x = num2cell(pev)'; % each row is a case
  y = cat(1, self_ev{:})';
  ind = subv2ind(CPD.sizes, [x y]);
  p = CPT(ind);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function B = chb(A,a,b)

% B = chb(A,a,b) change numerical base:
% a,b scalars from 2 to 64 , A,B strings.
% Number A in base a , expressed in base b by B.
% Symbols are 0 .. 9 A .. Z a .. z @ &.
%
% Example:
% '-54.13'-chb(chb('-54.13',6,27),27,6)
%
% Giampiero Campa 21-2-94

B=[];
l=length(A);
v=0;
c=0;
d=0;
y=0;

if l==0, error('A is unexisting')
end

if      ~isstr(A)|l==0|A=='-', error('A must be a string')
elseif   (a > 64)|(b > 64)|(a < 2)|(b < 2), error('invalid base')
end

if A(1)=='-',

s='-';
l=l-1;
A=A(2:l+1);

else s='';
end

for i = 1:l,

if      (A(i)=='&' & a<64)|(A(i)=='@' & a<63), error('invalid base or symbol')
elseif  (A(i)~='&' & A(i)~='@'),
	if  A(i)<'.'|A(i)=='/'|A(i)>'z'|...
	    (A(i)>'9' & A(i)<'A')|(A(i)>'Z' & A(i)<'a'), error('invalid base or symbol')
	end
end

if   (abs(A(i)) > a+60)|(abs(A(i))>a+54 & a<37)|(abs(A(i))>a+47 & a<11)|(A(i)=='.' & v>0),
     error('invalid base or symbol')

elseif  A(i)=='.', v=i;
end

end

if v==0, v=l+1;
end

for i=1:v-1;

if       A(i)<':', y=A(i)-'0';
elseif   A(i)<'[', y=A(i)-55;
else     y=A(i)-61;
end

if       A(i)=='@', y=62;
elseif   A(i)=='&', y=63;
end

c=c+y*a^(v-1-i);
end

for i=v+1:l;

if       A(i)<':', y=A(i)-'0';
elseif   A(i)<'[', y=A(i)-55;
else     y=A(i)-61;
end

if       A(i)=='@', y=62;
elseif   A(i)=='&', y=63;
end

d=d+y*a^(v-i);
end

while c>0,

z=round(b*(c/b-fix(c/b)));
c=fix(c/b);

if      z<10, q=setstr(48+z);
elseif  z<36, q=setstr(55+z);
else          q=setstr(61+z);
end

if z=='62', q='@';
elseif z>='63', q='&';
end

B=[q,B];
end

B=[s,B];

if d>0,

B=[B,'.'];
for i=1:38;

z=fix(d*b);
d=d*b-z;

if z==b,
d=1-1e-10;
z=b-1;
end

if      z<10, q=setstr(48+z);
elseif  z<36, q=setstr(55+z);
else          q=setstr(61+z);
end

if z=='62', q='@';
elseif z>='63', q='&';
end

B=[B,q];
end

while B(length(B))=='0',
B=B(1:length(B)-1);
end

end


text:
function [xnew,ynew]=splineValue(x,y)
%+++ spline for interpolation.
%+++
%+++ Aug., 2010.

minx=min(x);
maxx=max(x);
d=(maxx-minx)/128;
xnew=linspace(minx-d,maxx+d,128)';
ynew=spline(x,y,xnew);

author: viggin

repo: domain-adaptation-toolbox


text:
function precision = precision_table ( rule )

%*****************************************************************************80
%
%% PRECISION_TABLE returns the precision of a Lebedev rule.
%
%  Modified:
%
%    15 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule, between 1 and 65.
%
%    Output, integer PRECISION, the precision of the rule.
%
  rule_max = 65;

  table = [ ...
      3,   5,   7,   9,  11,  13,  15,  17,  19,  21, ...
     23,  25,  27,  29,  31,  33,  35,  37,  39,  41, ...
     43,  45,  47,  49,  51,  53,  55,  57,  59,  61, ...
     63,  65,  67,  69,  71,  73,  75,  77,  79,  81, ...
     83,  85,  87,  89,  91,  93,  95,  97,  99, 101, ...
    103, 105, 107, 109, 111, 113, 115, 117, 119, 121, ...
    123, 125, 127, 129, 131 ]';

  if ( rule < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PRECISION_TABLE - Fatal error!\n' );
    fprintf ( 1, '  RULE < 1.\n' );
    error ( 'PRECISION_TABLE - Fatal error!' );
  elseif ( rule_max < rule )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PRECISION_TABLE - Fatal error!\n' );
    fprintf ( 1, '  RULE_MAX < RULE.\n' );
    error ( 'PRECISION_TABLE - Fatal error!' );
  end

  precision = table(rule);

  return
end


text:
function [ o, x, w ] = en_r2_05_6 ( n )

%*****************************************************************************80
%
%% EN_R2_05_6 implements the Stroud rule 5.6 for region EN_R2.
%
%  Discussion:
%
%    The rule has order O = ( N + 1 ) * 2^N.
%
%    The rule has precision P = 5.
%
%    EN_R2 is the entire N-dimensional space with weight function
%
%      w(x) = exp ( - x1^2 - x2^2 ... - xn^2 )
%
%    The rule requires 5 <= N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%    5 <= N.
%
%    Output, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  if ( n < 5 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EN_R2_05_6 - Fatal error!\n' );
    fprintf ( 1, '  5 <= N is required.\n' );
    error ( 'EN_R2_05_6 - Fatal error!' );
  end

  o = 2^n * ( n + 1 );
  volume = sqrt ( pi^n );

  a = volume / 2^n / ( n + 1 );

  r = sqrt ( ( n - sqrt ( 2 ) + ( n - 1 ) * sqrt ( 2 * ( n + 1 ) ) ) / 2 / n );
  s = sqrt ( ( n - sqrt ( 2 ) -             sqrt ( 2 * ( n + 1 ) ) ) / 2 / n );
  t = sqrt ( ( 1 + sqrt ( 2 ) ) / 2 );

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  k = 0;
%
%  N * 2^N points.
%
  for i = 1 : n

    k = k + 1;
    x(1:n,k) = - s;
    x(i,k)   = - r;
    w(k) = a;

    more = 1;

    while ( more )
      more = 0;
      for j = n : -1 : 1
        if ( x(j,k) < 0.0 )
          k = k + 1;
          x(1:n,k) = x(1:n,k-1);
          x(j,k)     =   abs ( x(j,k) );
          x(j+1:n,k) = - abs ( x(j+1:n,k) );
          w(k) = a;
          more = 1;
          break;
        end
      end
    end

  end
%
%  2^N points.
%
  k = k + 1;
  x(1:n,k) = - t;
  w(k) = a;
  more = 1;
  while ( more )
    more = 0;
    for j = n : -1 : 1
      if ( x(j,k) < 0.0 )
        k = k + 1;
        x(1:n,k) = x(1:n,k-1);
        x(j,k)     =   abs ( x(j,k) );
        x(j+1:n,k) = - abs ( x(j+1:n,k) );
        w(k) = a;
        more = 1;
        break;
      end
    end
  end

  return
end


text:
load -ascii letter_A
load -ascii letter_T

[N, m] = size(letter_A)

class = N

app  = letter_A;size(app)
test = letter_T;size(test)

Napp = size(app,2);
Ntest = size(test,2);

unique(app(class,:))
unique(test(class,:))

ns1 = max(letter_A');
ns2 = max(letter_T');
ns = max(ns1,ns2)
clear letter_A letter_T ns1 ns2

% N, ns(class), Napp, Ntest, mean(ns),


text:
function [px, py, gx, gy] = position_5_link(z,P)
%[px,py,gx,gy] = POSITION_5_LINK(Z,P)
%
%FUNCTION:  This function computes the cartesian positions
%    of the CoM and tip of each link
%INPUTS:
%
%
%OUTPUTS:
%    px = [nLink X nTime] position of the end of each link
%    py = [nLink X nTime] position of the end of each link
%    gx = [nLink X nTime] position of the CoM of each link
%    gy = [nLink X nTime] position of the CoM of each link
%
%NOTES:
%    This file was automatically generated by writePosition.m

g  = P.g ; %gravity
m1 = P.m(1); % Link 1 mass
m2 = P.m(2); % Link 2 mass
m3 = P.m(3); % Link 3 mass
m4 = P.m(4); % Link 4 mass
m5 = P.m(5); % Link 5 mass
l1 = P.l(1); % Link 1 length
l2 = P.l(2); % Link 2 length
l3 = P.l(3); % Link 3 length
l4 = P.l(4); % Link 4 length
l5 = P.l(5); % Link 5 length
I1 = P.I(1); % Link 1 moment of inertia about its center of mass
I2 = P.I(2); % Link 2 moment of inertia about its center of mass
I3 = P.I(3); % Link 3 moment of inertia about its center of mass
I4 = P.I(4); % Link 4 moment of inertia about its center of mass
I5 = P.I(5); % Link 5 moment of inertia about its center of mass
d1 = P.d(1); % Link 1 distance between center of mass and parent joint
d2 = P.d(2); % Link 2 distance between center of mass and parent joint
d3 = P.d(3); % Link 3 distance between center of mass and parent joint
d4 = P.d(4); % Link 4 distance between center of mass and parent joint
d5 = P.d(5); % Link 5 distance between center of mass and parent joint

th1 = z(1,:);
th2 = z(2,:);
th3 = z(3,:);
th4 = z(4,:);
th5 = z(5,:);

nTime = length(th1);
px = zeros(5,nTime);
py = zeros(5,nTime);
gx = zeros(5,nTime);
gy = zeros(5,nTime);

px(1,:) = l1.*cos(th1);
py(1,:) = l1.*sin(th1);
gx(1,:) = d1.*cos(th1);
gy(1,:) = d1.*sin(th1);

px(2,:) = l1.*cos(th1) + l2.*cos(th2);
py(2,:) = l1.*sin(th1) + l2.*sin(th2);
gx(2,:) = d2.*cos(th2) + l1.*cos(th1);
gy(2,:) = d2.*sin(th2) + l1.*sin(th1);

px(3,:) = l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3);
py(3,:) = l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3);
gx(3,:) = d3.*cos(th3) + l1.*cos(th1) + l2.*cos(th2);
gy(3,:) = d3.*sin(th3) + l1.*sin(th1) + l2.*sin(th2);

px(4,:) = l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4);
py(4,:) = l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4);
gx(4,:) = d4.*cos(th4) + l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3);
gy(4,:) = d4.*sin(th4) + l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3);

px(5,:) = l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4) + l5.*cos(th5);
py(5,:) = l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4) + l5.*sin(th5);
gx(5,:) = d5.*cos(th5) + l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3) + l4.*cos(th4);
gy(5,:) = d5.*sin(th5) + l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3) + l4.*sin(th4);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [C,S,a,c]=getEMMCoeffs(M,year,fullyNormalize)
%%GETEMMCOEFFS Obtain spherical harmonic coefficients for the 2017
%              version of the National Oceanic and Atmospheric
%              Administration's (NOAA's) Enchaned Magnetic Model (EMM) at a
%              particular time or at the reference epoch (2017). The model
%              is considered valid down to 10km underground.
%
%INPUTS: M The integer maximum order of the spherical harmonic coefficients
%          obtained. This is a value between 1 and 740. If this parameter
%          is omitted, the default value is 740. If one wishes to load all
%          coefficients, one can also just pass Inf.
%     year A decimal number indicating a year in the Gregorian calendar as
%          specified by universal coordinated time (UTC). For example,
%          halfway through the year 2017 would be represented as 2017.5.
%          The precision of the model is not sufficiently fine that leap
%          seconds matter. If this parameter is omitted, then the
%          reference year of the EMM model is used. In this instance,
%          2017.0.
% fullyNormalize Geomagnetic models are normally given in terms of Schmidt
%          semi-normalized Legendre functions. If fullyNormalize=true, then
%          the coefficients are converted for use with fully normalized
%          associated Legendre functions, as are commonly used with
%          spherical harmonic algorithms for gravitational models. If this
%          parameter is omitted, the default value is true.
%
%OUTPUTS: C An array holding the coefficient terms that are multiplied by
%           cosines in the harmonic expansion. This can be given to a
%           CountingClusterSet so that C(n+1,m+1) is the coefficient of
%           degree n and order m. The coefficients have units of Tesla. The
%           coefficients are normalized according to the fullyNormalize
%           term.
%         S An array holding the coefficient terms that are multiplied by
%           sines in the harmonic expansion. The format of S is the same as
%           that of C.
%         a The numerator in the (a/r)^n term in the spherical harmonic sum
%           having units of meters.
%         c The constant value by which the spherical harmonic series is
%           multiplied, having units of squared meters.
%
%Details on the normalization of the coefficients is given in the comments
%to the function spherHarmonicEval.
%
%This function first checks for a .mat file with the coefficients in it.
%The .mat file is small and can be read quite quickly. However, if one does
%not exist, then it tries to read the EMM2015.COF and EMM2015SV.COF text
%files that one can obtain directly from the NOAA. Reading from the text
%files is very slow.
%
%Documentation is in [1] and [2]. Documentation as well as some code
%containing containing the coefficients is available at
%http://www.ngdc.noaa.gov/geomag/EMM/
%Being created by the U.S. government, the model is not subject to
%copyright.
%
%The data is kept in zipped files in the ./data folder. If all of the data
%is being loaded for a particular year (M=Inf), then a .mat file will be
%created in the ./data folder with the data for that year so that it can be
%loaded more quickly in the future.
%
%REFERENCES:
%[1] S. Maus, "An ellipsoidal harmonic representation of Earth's
%    lithospheric magnetic field to degree and order 720," Geochemistry,
%    Geophysics, Geosystems, vol. 11, no. 6, Jun. 2010.
%[2] Information on Schmidt semi-normalized Legendre functions is given in
%    D. E. Winch, D. J. Ivers, J. P. R. Turner, and R. J. Stening,
%    "Geomagnetism and Schmidt quasi-normalization," Geophysical Journal
%    International, vol. 160, no. 2, pp. 487-504, Feb. 2005.
%
%June 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<1||isempty(M))
    M=Inf;%Load all of the coefficients.
end

yearRef=2015.0;
if(nargin<2||isempty(year))
    year=yearRef;
end

if(nargin<3||isempty(fullyNormalize))
    fullyNormalize=true;
end

saveMatFile=(M==Inf);

%The EMM2015 magnetic coefficient data file, should be located in a data
%folder that is in the same folder as this file. This find the path to this
%file.
ScriptPath=mfilename('fullpath');
ScriptFolder=fileparts(ScriptPath);

%We have to load the version of the EMM that predicts for the year
%given. To know which one to load, we will see which years are in the
%zip archive. The years are given in the filenames.
filePathAndName=fileNamesInZipArchive([ScriptFolder,'/data/EMM_Coefficients.zip']);
secVarFilePathAndName=fileNamesInZipArchive([ScriptFolder,'/data/EMM_Secular_Variations.zip']);

numFiles=length(filePathAndName);
EMMYears=zeros(numFiles,1);
for curFile=1:numFiles
    [~,fileName] = fileparts(filePathAndName{curFile});
    %The filenames all begin with EMM, so skip the first three
    %characters to get the years.
    EMMYears(curFile)=sscanf(fileName(4:end),'%f',1);
end

%We have to find the EMM model year that is equal to the integer part
%of year. If there is none, then we have to find the closest year.
yearIdx2Load=find(EMMYears==fix(year),1);
if(isempty(yearIdx2Load)||yearIdx2Load~=year)
    [~,yearIdx2Load]=min(abs(EMMYears-fix(year)));
end
yearRef=EMMYears(yearIdx2Load);

%First, see if a .mat file with all of the data for the appropriate year
%exists. If so, then use that to load the coefficients from which
%interpolation must be performed.
matFile=[ScriptFolder,'/data/EMM',int2str(yearRef),'.mat'];

if(exist(matFile,'file'))
    load(matFile,'CCoeffs','SCoeffs','C1Coeffs','S1Coeffs');

    %Keep only as many coefficients as the maximum order provided.
    totalNumCoeffs=(M+1)*(M+2)/2;
    if(length(CCoeffs)<totalNumCoeffs)
        totalNumCoeffs=length(CCoeffs);
    end
    C=CountingClusterSet(CCoeffs(1:totalNumCoeffs));
    S=CountingClusterSet(SCoeffs(1:totalNumCoeffs));
    M=C.numClust-1;

    %If the drift coefficients need to be reduced in size.
    totalNumDriftCoeffs=min(length(C1Coeffs),totalNumCoeffs);
    if(length(C1Coeffs)<totalNumCoeffs)
        totalNumDriftCoeffs=length(C1Coeffs);
    end

    C1=CountingClusterSet(C1Coeffs(1:totalNumDriftCoeffs));
    S1=CountingClusterSet(S1Coeffs(1:totalNumDriftCoeffs));
else
    %Otherwise, just read the data from the text files.

    %We now know which coefficients and secular variation terms to load.
    temp=readZipArchive([ScriptFolder,'/data/EMM_Coefficients.zip'],filePathAndName{yearIdx2Load});
    coefficientFile=char(temp{2});
    temp=readZipArchive([ScriptFolder,'/data/EMM_Secular_Variations.zip'],secVarFilePathAndName{yearIdx2Load});
    secVarTerms=char(temp{2});

    %The first line of the coefficient file is just a description, so find
    %the first occurence of a return character and only take the file after
    %that character
    startIdx=1;
    %char(10) is the newline character used in the file.
    while(strcmp(coefficientFile(startIdx),newline)~=1)
        startIdx=startIdx+1;
    end
    startIdx=startIdx+1;
    %Load all of the coefficients.
    dataCoeffs=sscanf(coefficientFile(startIdx:end),'%d %d %f %f',Inf);
    numCoeffEntries=length(dataCoeffs)/4;
    dataCoeffs=reshape(dataCoeffs,4,numCoeffEntries);

    %The number of coefficients.
    M=min(M,dataCoeffs(1,end));
    totalNumCoeffs=(M+1)*(M+2)/2;

    %Load all of the drift coefficients.
    driftCoeffs=sscanf(secVarTerms,'%d %d %f %f',Inf);
    numDriftEntries=length(driftCoeffs)/4;
    driftCoeffs=reshape(driftCoeffs,4,numDriftEntries);
    MDrift=driftCoeffs(1,end);
    totalNumDriftCoeffs=(MDrift+1)*(MDrift+2)/2;

    %Allocate space for the coefficients.
    emptyData=zeros(totalNumCoeffs,1);
    C=CountingClusterSet(emptyData);
    S=CountingClusterSet(emptyData);
    %These will hold the time-varying terms.
    emptyData=zeros(totalNumDriftCoeffs,1);
    C1=CountingClusterSet(emptyData);
    S1=CountingClusterSet(emptyData);

    %Put the data into the CountingClusterSet class instances.
    for curRow=1:numCoeffEntries
        n=dataCoeffs(1,curRow);%The degree of the coefficient.
        m=dataCoeffs(2,curRow);%The order of the coefficient.
        C(n+1,m+1)=dataCoeffs(3,curRow);
        S(n+1,m+1)=dataCoeffs(4,curRow);
    end

    for curRow=1:numDriftEntries
        n=driftCoeffs(1,curRow);%The degree of the coefficient.
        m=driftCoeffs(2,curRow);%The order of the coefficient.
        C1(n+1,m+1)=driftCoeffs(3,curRow);
        S1(n+1,m+1)=driftCoeffs(4,curRow);
    end

    %Change the units fron Nanotesla to Tesla.
    C(:)=10^(-9)*C(:);
    S(:)=10^(-9)*S(:);
    C1(:)=10^(-9)*C1(:);
    S1(:)=10^(-9)*S1(:);

    %If the coefficients should be saved into a .mat file for more
    %efficient loading in the future.
    if(saveMatFile)
        CCoeffs=C.clusterEls;
        SCoeffs=S.clusterEls;

        C1Coeffs=C1.clusterEls;
        S1Coeffs=S1.clusterEls;

        save(matFile,'CCoeffs','SCoeffs','C1Coeffs','S1Coeffs');
    end
end

%If interpolation to other dates must be performed.
if(year~=yearRef)
    yearDiff=year-yearRef;
    %Perform linear interpolation.
    idx=1:totalNumDriftCoeffs;
    C.clusterEls(idx)=C.clusterEls(idx)+yearDiff*C1.clusterEls(idx);
    S.clusterEls(idx)=S.clusterEls(idx)+yearDiff*S1.clusterEls(idx);
end

%If the coefficients should be fully normalized.
if(fullyNormalize~=false)
     for n=0:M
        k=1/sqrt(1+2*n);
        C(n+1,:)=k*C(n+1,:);
        S(n+1,:)=k*S(n+1,:);
     end
end

%Return S and C as arrays, not as a CountingClusterSet classes.
C=C.clusterEls;
S=S.clusterEls;

%The EMM2015 model uses the same reference ellipse as the WMM2010.
a=Constants.WMM2010SphereRad;%meters
c=a^2;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [propulsionModel] = definePropulsionModel(quadRotorParams,plotflag)
% [propulsionModel] = definePropulsionModel(quadrotorParams)
%
% Defines a propulsion parameter plant model for use with dynBodyFrame.
% Supports any number of thrusters, n.
%
% First row of 'thrustLocations','thrustAxes','isSpinDirectionCCW' are
% associated with the first motor.
% Second row of 'thrustLocations','thrustAxes','isSpinDirectionCCW', are associate with the second
% motor, and so forth.
%
% Following parameters are the same for each thruster
%   (i.e.assumes same properties of prop-motor combination at all locations):
%       maxThrust, maxRPM, maxTorque, d_prop
%
% Inputs:
%   quadRotorParams [struct] parameter struct with the following fields
%       .thrustLocations = [n x 3] [m] location of motors
%       .thrustAxes = [n x 3] [unitvectors] direction of thrust axis
%       .isSpinDirectionCCW = [n x 1] [bool] if set
%       .maxThrust [n x 1] thrust at 100% throttle (N)
%       .maxRPM [n x 1] RPM at 100% throttle (RPM)
%       .maxTorque [n x 1] torque at 100% throttle (Nm)
%       .d_prop [n x 1] propeller diameter (m)
%
% Output:
%   quadRotorModel = [struct] with n elements, one for each motor.
%
% Written by Conrad McGreal 2020

% give shorter name
p = quadRotorParams ;

%% Compute propulsion system coefficients
% See: https://web.mit.edu/16.unified/www/FALL/thermodynamics/notes/node86.html
% air density during propulsion data collection (kg/m^3)
% (used to determine propulsion system coefficients).
rho = 1.225 ;

%% Assign to struct
propulsionModel = struct() ; % initialize output struct

for i=1:size(p.thrustLocations,1)
    propulsionModel(i).thrustAxis = p.thrustAxes(i,:) ;      % [port, nose, top] % nose faces north when body and world axis are aligned (world is "East North Up")
    propulsionModel(i).thrustLocation = p.thrustLocations(i,:) ;
    propulsionModel(i).isSpinDirectionCCW = p.isSpinDirectionCCW(i,:) ;
    propulsionModel(i).maxRPM = p.maxRPM(i) ;
    propulsionModel(i).maxTorque = p.maxTorque(i) ;
    propulsionModel(i).d_prop = p.d_prop(i) ;

    % compute coefficients.
    C_t = p.maxThrust(i) / (rho * (p.maxRPM(i)/60)^2 * p.d_prop(i)^4) ; % thrust coefficient
    propulsionModel(i).C_t = C_t ;

    C_q = p.maxTorque(i) / (rho * (p.maxRPM(i)/60)^2 * p.d_prop(i)^5) ; % torque coefficient
    propulsionModel(i).C_q = C_q ;

end

% plot propulsion model, if desired
if ~exist('plotflag','var')
else
    if plotflag
        showPropulsionModel(propulsionModel) ;
    end
end

repo: OptimTraj


text:
function hypercube_integrals_test01 ( )

%*****************************************************************************80
%
%% HYPERCUBE_INTEGRALS_TEST01: estimate integrals over the unit hypercube in 3D
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2014
%
%  Author:
%
%    John Burkardt
%
  m = 3;
  n = 4192;
  test_num = 20;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  Compare exact and estimated integrals\n' );
  fprintf ( 1, '  over the interior of the unit hypercube in 3D.\n' );
%
%  Get sample points.
%
  seed = 123456789;
  [ x, seed ] = hypercube01_sample ( m, n, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of sample points used is %d\n', n );
%
%  Randomly choose exponents.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ex  Ey  Ez     MC-Estimate           Exact      Error\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    [ e, seed ] = i4vec_uniform_ab ( m, 0, 4, seed );

    value = monomial_value ( m, n, e, x );

    result = hypercube01_volume ( m ) * sum ( value(1:n) ) / n;
    exact = hypercube01_monomial_integral ( m, e );
    error = abs ( result - exact );

    fprintf ( 1, '  %2d  %2d  %2d  %14.6g  %14.6g  %10.2e\n', ...
      e(1:m), result, exact, error );

  end

  return
end


text:
function J = vl_imup(I)
% VL_IMUP Upsample an image by two
%   J=VL_IMUP(I) doubles the resolution of the image I by using
%   bilinear interpolation.
%
%   See also: VL_IMDOWN(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

[M,N,K] = size(I) ;

J = zeros(2*M,2*N,K) ;

J(1:2:end,1:2:end,:) = I ;

J(2:2:end,1:2:end,:) = 0.5*(I+[I(2:end,:,:);I(end,:,:)]) ;
J(1:2:end,2:2:end,:) = 0.5*(I+[I(:,2:end,:),I(:,end,:)]) ;
J(2:2:end,2:2:end,:) = ...
  0.25*(...
  J(2:2:end,1:2:end-1,:)+...
  J(1:2:end-1,2:2:end,:)+...
  [J(2:2:end,3:2:end,:),J(2:2:end,end-1,:)]+...
  [J(3:2:end,2:2:end,:);J(end-1,2:2:end,:)]) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [res,ires]=dydt1(t,y,yp)

res(1)=yp(1)+4*y(1);
res(2)=yp(2)+6*y(2);
res(3)=y(3)-y(1)-2.0;
ires=0;


text:
 function xpad = zpad0(x, nxp, nyp, nzp)
%function xpad = zpad0(x, nxp, nyp, nzp)
% zero pad an input signal x symmetrically around "0" (image center)
% by A. Yendiki

if nargin < 1, ir_usage, end
if nargin == 1 && streq(x, 'test'), zpad0_test, return, end

if nargin == 2
	xpad = ir_pad_into_center(x, nxp);
elseif nargin == 3
	xpad = ir_pad_into_center(x, [nxp nyp]);
elseif nargin == 4
	xpad = ir_pad_into_center(x, [nxp nyp nzp]);
else
	error 'Too many arguments'
end

function zpad0_test
zpad0([1 2 1], 7)
zpad0([1 2 1], 7, 3)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function hout=suptitle(str,varargin)
    % biafra ahanonu
    % started updating: 2016.03.10
    % this is an overloaded suptitle to remove focus stealing

    %SUPTITLE puts a title above all subplots.
    %
    %	SUPTITLE('text') adds text to the top of the figure
    %	above all subplots (a "super title"). Use this function
    %	after all subplot commands.
    %
    %   SUPTITLE is a helper function for yeastdemo.

    %   Copyright 2003-2010 The MathWorks, Inc.

    % Warning: If the figure or axis units are non-default, this
    % will break.

    %========================
    % Amount of the figure window devoted to subplots
    options.plotregion = 0.94;
    % Y position of title in normalized coordinates
    options.titleypos = 0.96;
    % font size
    options.fontSize = [];
    % Str or rgb vector: color of suptitle
    options.Color = 'k';
	% Str: interpreter, latex or tex
	options.Interpreter = 'latex';
    % get options
    options = ciapkg.io.getOptions(options,varargin);
    % display(options)
    % unpack options into current workspace
    % fn=fieldnames(options);
    % for i=1:length(fn)
    %   eval([fn{i} '=options.' fn{i} ';']);
    % end
    %========================

    % Parameters used to position the supertitle.

    % Amount of the figure window devoted to subplots
    plotregion = options.plotregion;

    % Y position of title in normalized coordinates
    titleypos  = options.titleypos;

    % Fontsize for supertitle
    if isempty(options.fontSize)
        fs = get(gcf,'defaultaxesfontsize')+4;
    else
        fs = options.fontSize;
    end

    % Fudge factor to adjust y spacing between subplots
    fudge=1;

    haold = gca;
    figunits = get(gcf,'units');

    % Get the (approximate) difference between full height (plot + title
    % + xlabel) and bounding rectangle.

    if (~strcmp(figunits,'pixels')),
        set(gcf,'units','pixels');
        pos = get(gcf,'position');
        set(gcf,'units',figunits);
    else
        pos = get(gcf,'position');
    end
    ff = (fs-4)*1.27*5/pos(4)*fudge;

    % The 5 here reflects about 3 characters of height below
    % an axis and 2 above. 1.27 is pixels per point.

    % Determine the bounding rectangle for all the plots

    % h = findobj('Type','axes');

    % findobj is a 4.2 thing.. if you don't have 4.2 comment out
    % the next line and uncomment the following block.

    h = findobj(gcf,'Type','axes');  % Change suggested by Stacy J. Hills

    max_y=0;
    min_y=1;
    oldtitle = NaN;
    numAxes = length(h);
    thePositions = zeros(numAxes,4);
    for i=1:numAxes
        pos=get(h(i),'pos');
        thePositions(i,:) = pos;
        if (~strcmp(get(h(i),'Tag'),'suptitle')),
            if (pos(2) < min_y)
                min_y=pos(2)-ff/5*3;
            end;
            if (pos(4)+pos(2) > max_y)
                max_y=pos(4)+pos(2)+ff/5*2;
            end;
        else
            oldtitle = h(i);
        end
    end
    if max_y > plotregion,
        scale = (plotregion-min_y)/(max_y-min_y);
        for i=1:numAxes
            pos = thePositions(i,:);
            pos(2) = (pos(2)-min_y)*scale+min_y;
            pos(4) = pos(4)*scale-(1-scale)*ff/5*3;
            set(h(i),'position',pos);
        end
    end

    np = get(gcf,'nextplot');
    set(gcf,'nextplot','add');
    if ishghandle(oldtitle)
        delete(oldtitle);
    end
    warning off
    axes('pos',[0 1 1 1],'visible','off','Tag','suptitle','SortMethod','depth');
    warning on
    ht = text('position',[.5 titleypos-1],'Interpreter',options.Interpreter,'String',str);
	set(ht,'horizontalalignment','center','fontsize',fs,'Color',options.Color);
	set(ht,'FontName',get(0,'DefaultAxesFontName'));
    set(gcf,'nextplot',np);
    % focus stealing occurs here, turn off
    set(gcf, 'CurrentAxes', haold)
    % axes(haold); %#ok<MAXES>
    if nargout,
        hout=ht;
    end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
 function ob = Ginterp1(mask, varargin)
%function ob = Ginterp1(mask, options)
%|
%| Construct Ginterp1 object for 1D image registration.
%| This method does not have an adjoint.
%|
%| See Ginterp1_test() below for example usage.
%|
%| in
%|	mask	size(image)	logical array of object support.
%|
%| options
%|	'interp1_arg'	{}	arguments to interp1 for 'forward' operation
%|
%| out
%|	ob [nd np]	np = sum(mask(:)), so it is already "masked"
%|			nd = numel(interp1_arg{2})
%|
%| Copyright 2013-12-05, Jeff Fessler, University of Michigan

if nargin == 1 && streq(mask, 'test'), Ginterp1_test, return, end
if nargin < 1, ir_usage, end

arg.mask = mask;

% option defaults
arg.interp1_x = [];
arg.interp1_xq = [];
arg.interp1_arg = {'linear', 0};

% options specified by name/value pairs
arg = vararg_pair(arg, varargin);

if isempty(arg.interp1_x) || isempty(arg.interp1_xq)
	fail '"interp1_arg" and "interp1_x" required'
end

arg.ndim = ndims(mask);
if arg.ndim == 2 && size(mask,2) == 1
	arg.ndim = 1;
else
	fail 'only 1D mask allowed'
end

if streq(arg.interp1_arg{1}, 'linear')
	abs_arg = {'abs', @(ob) ob}; % linear interpolation uses nonnegative coefficients
else
	abs_arg = {}; % unknown
end

arg.fun_forw = @(arg, x) ...
	interp1(arg.interp1_x, x, arg.interp1_xq, arg.interp1_arg{:});
arg.fun_back = @(arg, y) fail('interp adjoint not done');

% build object
idim = size(mask);
if numel(idim) == 2 && idim(2) == 1
	idim = idim(1); % 1d
end
ob = fatrix2('mask', mask, 'arg', arg, ...
	'idim', idim, 'odim', numel(arg.interp1_xq), ...
	abs_arg{:}, 'forw', arg.fun_forw, 'back', arg.fun_back);

% Ginterp1_test()
function Ginterp1_test

nx = 8;
mask = true(nx,1);
mask(1) = false;
xq = linspace(0,nx+1, 101);

A = Ginterp1(mask, 'interp1_x', 1:nx, 'interp1_xq', xq);

if 0 % todo
	fatrix2_tests(A, 'complex', 0, 'halt', 0, ...
		'check1', false, 'full', false) % because of bad adjoint
%	test_adjoint(A, 'complex', 1);
else
	warn 'adjoint not tested'
end

im plc 1 2
Af = full(A);
im(1, Af'), axis xy, axis normal

im subplot 2
x = [1:nx]' .* mask;
y = A * x;

plot(1:nx, x, 'o', xq, y, '.-')


text:
function [M,P] = ghkf_predict(M,P,f,Q,f_param,p)
% GHKF_PREDICT - Gauss-Hermite Kalman filter prediction step
%
% Syntax:
%   [M,P] = GHKF_PREDICT(M,P,[f,Q,f_param,p])
%
% In:
%   M - Nx1 mean state estimate of previous step
%   P - NxN state covariance of previous step
%   f - Dynamic model function as a matrix A defining
%       linear function f(x) = A*x, inline function,
%       function handle or name of function in
%       form f(x,param)                   (optional, default eye())
%   Q - Process noise of discrete model   (optional, default zero)
%   f_param - Parameters of f               (optional, default empty)
%   p - Degree of approximation (number of quadrature points)
%
% Out:
%   M - Updated state mean
%   P - Updated state covariance
%
% Description:
%   Perform additive form Gauss-Hermite Kalman Filter prediction step.
%
%   Function f(.) should be such that it can be given a
%   DxN matrix of N sigma Dx1 points and it returns
%   the corresponding predictions for each sigma
%   point.
%
% See also:
%   GHKF_UPDATE, GHRTS_SMOOTH, GH_TRANSFORM

% History:
%   Aug 5,  2010 - Renamed from 'gh_predict' to 'ghkf_predict' (asolin)

% Copyright (C) 2009 Hartikainen, Särkkä, Solin
%
% $Id: gh_predict.m,v 1.2 2009/07/01 06:34:40 ssarkka Exp $
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.
%%

  %
  % Check which arguments are there
  %
  if nargin < 2
    error('Too few arguments');
  end
  if nargin < 3
    f = [];
  end
  if nargin < 4
    Q = [];
  end
  if nargin < 5
    f_param = [];
  end
  if nargin < 6
     p = [];
  end

  %
  % Apply defaults
  %
  if isempty(f)
    f = eye(size(M,1));
  end
  if isempty(Q)
    Q = zeros(size(M,1));
  end
  if isempty(p)
    p = 10;
  end

  %
  % Do transform and add process noise
  %
  tr_param = {p};
  [M,P] = gh_transform(M,P,f,f_param,tr_param);
  P = P + Q;

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
%% <../index.html MATLAB xUnit Test Framework>: How to Test Using a Floating-Point Tolerance
% MATLAB performs arithmetic operations using the floating-point
% hardware instructions on your processor. Because
% almost all floating-point operations are subject to round-off
% error, arithmetic operations can sometimes produce surprising
% results. Here's an example.

a = 1 + 0.1 + 0.1 + 0.1

%%
a == 1.3

%%
% So why doesn't |a| equal 1.3? Because 0.1, 1.3, and most other
% decimal fractions do not have exact representations in the binary
% floating-point number representation your computer uses.  The
% first line above is doing an approximate addition of 1 plus an
% approximation of 0.1, plus an approximation of 0.1, plus an
% approximation of 0.1.  The second line compares the result of all
% that with an approximation of 1.3.
%
% If you subtract 1.3 from |a|, you can see that the computed result
% for |a| is _extremely close_ to the floating-point approximation
% of 1.3, but it is not exactly the same.

a - 1.3

%%
% As a general rule, when comparing the results of floating-point
% calculations for equality, it is necessary to use a tolerance
% value.  Two types of tolerance comparisons are commonly used: absolute
% tolerance and relative tolerance.  An absolute tolerance comparison of _a_ and _b_
% looks like:
%
% $$|a-b| \leq T$$
%
% A relative tolerance comparison looks like:
%
% $$|a-b| \leq T\max(|a|,|b|) + T_f$$
%
% where _Tf_ is called the _floor tolerance_. It acts as an absolute tolerance
% when _a_ and _b_ are very close to 0.
%
% For example, suppose that _a_ is 100, _b_ is 101, and T is 0.1.  Then _a_ and
% _b_ would not be considered equal using an absolute tolerance, because 1 >
% 0.1.  However, _a_ and _b_ would be considered equal using a relative
% tolerance, because they differ by only 1 part in 100.
%
% MATLAB xUnit provides the utility assertion functions called
% |assertElementsAlmostEqual| and |assertVectorAlmostEqual|. These functions
% make it easy to write tests involving floating-point tolerances.
%
% |assertElementsAlmostEqual(A,B)| applies the tolerance test independently to
% every element of |A| and |B|.  The function uses a relative tolerance test by
% default, but you make it use an absolute tolerance test, or change the
% tolerance values used, by passing additional arguments to it.
%
% |assertVectorsAlmostEqual(A,B)| applies the tolerance test to the vectors |A|
% and |B| in the L2-norm sense.  For example, suppose |A| is |[1 1e10|], |B|
% is |[2 1e10]|, and the tolerance is 1e-8.  Then |A| and |B| would fail an
% elementwise relative tolerance comparison, because the relative difference
% between the first elements is 0.5.  However, they would pass a vector relative
% tolerance comparison, because the relative vector difference between |A| and
% |B| is only about 1 part in 1e10.
%
% The |examples_general| directory contains a portion of a unit test for the
% |sin| function.  The output of |sin| can sometimes be a bit surprising because
% of floating-point issues.  For example:

sin(pi)

%%
% That's very close but not exactly equal to 0.  Here's how the
% |sin| unit test uses |assertElementsAlmostEqual| to write the |sin(pi)|
% test with a minimum of fuss.

cd examples_general
type testSin

%%
% Run the test using |runtests|.

runtests testSin

%%
% <../index.html Back to MATLAB xUnit Test Framework>

%%
% Copyright 2008-2010 The MathWorks, Inc.

author: ricedsp