text:
function [selected_file, selected_path] = rri_select_file(varargin)
%
%  USAGE: [selected_file, selected_path] = ...
%             rri_select_file(dir_name, fig_title)
%
%  Allow user to select a file from a list of Matlab competible
%	file format
%
%  Example:
%
%    [selected_file, selected_path] = ...
%		rri_select_file('/usr','Select Data File');
%
%  See Also RRI_GETFILES

%  -- Created June 2001 by Wilkin Chau, Rotman Research Institute
%
%  use rri_select_file to open & save Matlab recognized format
%  -- Modified Dec 2002 by Jimmy Shen, Rotman Research Institute
%

   if nargin == 0 | ischar(varargin{1}) 	% create rri_select_file figure

      dir_name = '';
      fig_title = 'Select a File';

      if nargin > 0
         dir_name = varargin{1};
      end

      if nargin > 1
         fig_title = varargin{2};
      end

      Init(fig_title,dir_name);
      uiwait;                           % wait for user finish

      selected_path = getappdata(gcf,'SelectedDirectory');
      selected_file = getappdata(gcf,'SelectedFile');

      cd (getappdata(gcf,'StartDirectory'));
      close(gcf);
      return;
   end;

   %  clear the message line,
   %
   h = findobj(gcf,'Tag','MessageLine');
   set(h,'String','');

   action = varargin{1}{1};

   %  change 'File format':
   %  update 'Files' & 'File selection' based on file pattern
   %
   if strcmp(action,'EditFilter'),
      EditFilter;

   %  run delete_fig when figure is closing
   %
   elseif strcmp(action,'delete_fig'),
      delete_fig;

   %  select 'Directories':
   %  go into the selected dir
   %  update 'Files' & 'File selection' based on file pattern
   %
   elseif strcmp(action,'select_dir'),
      select_dir;

   %  select 'Files':
   %  update 'File selection'
   %
   elseif strcmp(action,'select_file'),
      select_file;

   %  change 'File selection':
   %  if it is a file, select that,
   %  if it is more than a file (*), select those,
   %  if it is a directory, select based on file pattern
   %
   elseif strcmp(action,'EditSelection'),
      EditSelection;

   %  clicked 'Select'
   %
   elseif strcmp(action,'DONE_BUTTON_PRESSED'),
      h = findobj(gcf,'Tag','SelectionEdit');
      [filepath,filename,fileext] = fileparts(get(h,'String'));

      if isempty(filepath) | isempty(filename) | isempty(fileext)
         setappdata(gcf,'SelectedDirectory',[]);
         setappdata(gcf,'SelectedFile',[]);
      else
         if ~strcmp(filepath(end),filesep)		% not end with filesep
            filepath = [filepath filesep];		% add a filesep to filepath
         end

         setappdata(gcf,'SelectedDirectory',filepath);
         setappdata(gcf,'SelectedFile',[filename fileext]);
      end

      if getappdata(gcf,'ready')			% ready to exit
         uiresume;
      end

   %  clicked 'cancel'
   %
   elseif strcmp(action,'CANCEL_BUTTON_PRESSED'),
      setappdata(gcf,'SelectedDirectory',[]);
      setappdata(gcf,'SelectedFile',[]);
      set(findobj(gcf,'Tag','FileList'),'String','');
      uiresume;
   end;

   return;

% --------------------------------------------------------------------
function Init(fig_title,dir_name),

   StartDirectory = pwd;
   if isempty(StartDirectory),
       StartDirectory = filesep;
   end;

   filter_disp = {'JPEG image (*.jpg)', ...
	'TIFF image, compressed (*.tif)', ...
	'EPS Level 1 (*.eps)', ...
	'Adobe Illustrator 88 (*.ai)', ...
	'Enhanced metafile (*.emf)', ...
	'Matlab Figure (*.fig)', ...
	'Matlab M-file (*.m)', ...
	'Portable bitmap (*.pbm)', ...
	'Paintbrush 24-bit (*.pcx)', ...
	'Portable Graymap (*.pgm)', ...
	'Portable Network Graphics (*.png)', ...
	'Portable Pixmap (*.ppm)', ...
   };

   filter_string = {'*.jpg', ...
	'*.tif', ...
	'*.eps', ...
	'*.ai', ...
	'*.emf', ...
	'*.fig', ...
	'*.m', ...
	'*.pbm', ...
	'*.pcx', ...
	'*.pgm', ...
	'*.png', ...
	'*.ppm', ...
   };

%   filter_disp = char(filter_disp);
   filter_string = char(filter_string);

   margine = 0.05;
   line_height = 0.07;
   char_height = line_height*0.8;

   save_setting_status = 'on';
   rri_select_file_pos = [];

   try
      load('pls_profile');
   catch
   end

   if ~isempty(rri_select_file_pos) & strcmp(save_setting_status,'on')

      pos = rri_select_file_pos;

   else

      w = 0.4;
      h = 0.6;
      x = (1-w)/2;
      y = (1-h)/2;

      pos = [x y w h];

   end

   h0 = figure('parent',0, 'Color',[0.8 0.8 0.8], ...
        'Units','normal', ...
        'Name',fig_title, ...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'Position', pos, ...
        'deleteFcn','rri_select_file({''delete_fig''});', ...
        'WindowStyle', 'modal', ...
        'Tag','GetFilesFigure', ...
        'ToolBar','none');

   x = margine;
   y = 1 - 1*line_height - margine;
   w = 1-2*x;
   h = char_height;

   pos = [x y w h];

   h1 = uicontrol('Parent',h0, ...		% Filter Label
        'Style','text', ...
        'Units','normal', ...
        'BackgroundColor',[0.8 0.8 0.8], ...
	'fontunit','normal', ...
        'FontSize',0.5, ...
        'HorizontalAlignment','left', ...
        'Position', pos, ...
        'String','Choose one of the file format:', ...
        'Tag','FilterLabel');

   y = 1 - 2*line_height - margine + line_height*0.2;
   w = 1-2*x;

   pos = [x y w h];

   h_filter = uicontrol('Parent',h0, ...	% Filter list
        'Style','popupmenu', ...
        'Units','normal', ...
        'BackgroundColor',[1 1 1], ...
	'fontunit','normal', ...
        'FontSize',0.5, ...
        'HorizontalAlignment','left', ...
        'Position', pos, ...
        'String', filter_disp, ...
        'user', filter_string, ...
	'value', 1, ...
        'Callback','rri_select_file({''EditFilter''});', ...
        'Tag','FilterEdit');

   y = 1 - 3*line_height - margine;
   w = 0.5 - x - margine/2;

   pos = [x y w h];

   h1 = uicontrol('Parent',h0, ...            % Directory Label
        'Style','text', ...
        'Units','normal', ...
        'BackgroundColor',[0.8 0.8 0.8], ...
	'fontunit','normal', ...
        'FontSize',0.5, ...
        'HorizontalAlignment','left', ...
        'ListboxTop',0, ...
        'Position', pos, ...
        'String','Directories', ...
        'Tag','DirectoryLabel');

   x = 0.5;
   y = 1 - 3*line_height - margine;
   w = 0.5 - margine;

   pos = [x y w h];

   h1 = uicontrol('Parent',h0, ...            % File Label
        'Style','text', ...
        'Units','normal', ...
        'BackgroundColor',[0.8 0.8 0.8], ...
	'fontunit','normal', ...
        'FontSize',0.5, ...
        'HorizontalAlignment','left', ...
        'ListboxTop',0, ...
        'Position', pos, ...
        'String','Files', ...
        'Tag','FileLabel');

   x = margine;
   y = 4*line_height + margine;
   w = 0.5 - x - margine/2;
   h = 1 - 7*line_height - 2*margine;

   pos = [x y w h];

   h_dir = uicontrol('Parent',h0, ...            % Directory Listbox
        'Style','listbox', ...
        'Units','normal', ...
	'fontunit','normal', ...
        'FontSize',0.08, ...
        'HorizontalAlignment','left', ...
        'Interruptible', 'off', ...
        'ListboxTop',1, ...
        'Position', pos, ...
        'String', '', ...
        'Callback','rri_select_file({''select_dir''});', ...
        'Tag','DirectoryList');

   x = 0.5;
   y = 4*line_height + margine;
   w = 0.5 - margine;
   h = 1 - 7*line_height - 2*margine;

   pos = [x y w h];

   h_file = uicontrol('Parent',h0, ...            % File Listbox
        'Style','listbox', ...
        'Units','normal', ...
	'fontunit','normal', ...
        'FontSize',0.08, ...
        'HorizontalAlignment','left', ...
        'ListboxTop',1, ...
        'Position', pos, ...
        'String', '', ...
        'Callback','rri_select_file({''select_file''});', ...
        'Tag','FileList');

   x = margine;
   y = 3*line_height + margine - line_height*0.2;
   w = 1-2*x;
   h = char_height;

   pos = [x y w h];

   h1 = uicontrol('Parent',h0, ...            % Selection Label
        'Style','text', ...
        'Units','normal', ...
        'BackgroundColor',[0.8 0.8 0.8], ...
        'fontunit','normal', ...
        'FontSize',0.5, ...
        'HorizontalAlignment','left', ...
        'Position', pos, ...
        'String','File you selected:', ...
        'Tag','SelectionLabel');

   y = 2*line_height + margine;
   w = 1-2*x;

   pos = [x y w h];

   h_select = uicontrol('Parent',h0, ...            % Selection Edit
        'Style','edit', ...
        'Units','normal', ...
        'BackgroundColor',[1 1 1], ...
        'fontunit','normal', ...
        'FontSize',0.5, ...
        'HorizontalAlignment','left', ...
        'Position', pos, ...
        'String', '', ...
        'Callback','rri_select_file({''EditSelection''});', ...
        'Tag','SelectionEdit');

   x = 2*margine;
   y = line_height/2 + margine;
   w = 0.2;
   h = line_height;

   pos = [x y w h];

   h_done = uicontrol('Parent',h0, ...                      % DONE
        'Units','normal', ...
        'fontunit','normal', ...
        'FontSize',0.5, ...
        'ListboxTop',0, ...
        'Position', pos, ...
        'HorizontalAlignment','center', ...
        'String','Save', ...			% 'Select', ...
        'Callback','rri_select_file({''DONE_BUTTON_PRESSED''});', ...
        'Tag','DONEButton');

   x = 1 - x - w;

   pos = [x y w h];

   h_cancel = uicontrol('Parent',h0, ...                      % CANCEL
        'Units','normal', ...
        'fontunit','normal', ...
        'FontSize',0.5, ...
        'ListboxTop',0, ...
        'Position', pos, ...
        'HorizontalAlignment','center', ...
        'String','Cancel', ...
        'Callback','rri_select_file({''CANCEL_BUTTON_PRESSED''});', ...
        'Tag','CANCELButton');

   if isempty(dir_name)
      dir_name = StartDirectory;
   end

   set(h_select,'string',dir_name);

   filter_select = get(h_filter,'value');
   filter_pattern = filter_string(filter_select,:);

   setappdata(gcf,'FilterPattern',deblank(filter_pattern));
   setappdata(gcf,'filter_string',filter_string);

   setappdata(gcf,'h_filter', h_filter);
   setappdata(gcf,'h_dir', h_dir);
   setappdata(gcf,'h_file', h_file);
   setappdata(gcf,'h_select', h_select);
   setappdata(gcf,'h_done', h_done);
   setappdata(gcf,'h_cancel', h_cancel);
   setappdata(gcf,'StartDirectory',StartDirectory);

   EditSelection;

   h_file = getappdata(gcf,'h_file');
   if isempty(get(h_file,'string'))
      setappdata(gcf,'ready',0);
   else
      setappdata(gcf,'ready',1);
   end

   return;					% Init

%  called by all the actions, to update 'Directories' or 'Files'
%  based on filter_pattern. Select first file in filelist.
%
% --------------------------------------------------------------------

function update_dirlist;

   filter_path = getappdata(gcf,'curr_dir');
   filter_pattern = getappdata(gcf,'FilterPattern');

   if exist(filter_pattern) == 2	% user input specific filename
      is_single_file = 1;		% need manually take path out later
   else
      is_single_file = 0;
   end

   % take the file path out from filter_pattern
   %
   [fpath fname fext] = fileparts(filter_pattern);
   filter_pattern = [fname fext];

   dir_struct = dir(filter_path);
   if isempty(dir_struct)
      msg = 'ERROR: Directory not found!';
      uiwait(msgbox(msg,'File Selection Error','modal'));
      return;
   end;

   old_pointer = get(gcf,'Pointer');
   set(gcf,'Pointer','watch');

   dir_list = dir_struct(find([dir_struct.isdir] == 1));
   [sorted_dir_names,sorted_dir_index] = sortrows({dir_list.name}');

   dir_struct = dir([filter_path filesep filter_pattern]);
   if isempty(dir_struct)
      sorted_file_names = [];
   else
      file_list = dir_struct(find([dir_struct.isdir] == 0));

      if is_single_file			% take out path
         tmp = file_list.name;
         [fpath fname fext] = fileparts(tmp);
         file_list.name = [fname fext];
      end

      [sorted_file_names,sorted_file_index] = sortrows({file_list.name}');
   end;

   disp_dir_names = [];			% if need full path, use this
					% instead of sorted_dir_names
   for i=1:length(sorted_dir_names)
      tmp = [filter_path filesep sorted_dir_names{i}];
      disp_dir_names = [disp_dir_names {tmp}];
   end

   h = findobj(gcf,'Tag','DirectoryList');
   set(h,'String',sorted_dir_names,'Value',1);

   h = findobj(gcf,'Tag','FileList');
   set(h,'String',sorted_file_names,'value',1);

   h_select = getappdata(gcf,'h_select');
   if strcmp(filter_path(end),filesep)		% filepath end with filesep
      filter_path = filter_path(1:end-1);	% take filesep out
   end

   if isempty(sorted_file_names)
      set(h_select,'string',[filter_path filesep]);
   else
      set(h_select,'string',[filter_path filesep sorted_file_names{1}]);
   end

   set(gcf,'Pointer',old_pointer);

   return; 					% update_dirlist

%  change 'File format':
%  update 'Files' & 'File selection' based on file pattern
%
% --------------------------------------------------------------------

function EditFilter()

   filter_select = get(gcbo,'value');
   filter_string = getappdata(gcf,'filter_string');
   filter_pattern = filter_string(filter_select,:);
   filter_path = getappdata(gcf,'curr_dir');

   % update filter_pattern
   setappdata(gcf,'FilterPattern',deblank(filter_pattern));

   if isempty(filter_path),
       filter_path = filesep;
   end;

   update_dirlist;

   h_file = getappdata(gcf,'h_file');
   if isempty(get(h_file,'string'))
      setappdata(gcf,'ready',0);
   else
      setappdata(gcf,'ready',1);
   end

   return;					% EditFilter

%  select 'Directories':
%  go into the selected dir
%  update 'Files' & 'File selection' based on file pattern
%
% --------------------------------------------------------------------

function select_dir()

   listed_dir = get(gcbo,'String');
   selected_dir_idx = get(gcbo,'Value');
   selected_dir = listed_dir{selected_dir_idx};
   curr_dir = getappdata(gcf,'curr_dir');

   %  update the selection box
   %
   try
      cd ([curr_dir filesep selected_dir]);
   catch
      msg = 'ERROR: Cannot access directory';
      uiwait(msgbox(msg,'File Selection Error','modal'));
      return;
   end;

   if isempty(pwd)
      curr_dir = filesep;
   else
      curr_dir = pwd;
   end;

   setappdata(gcf,'curr_dir',curr_dir);
   update_dirlist;

   h_file = getappdata(gcf,'h_file');
   if isempty(get(h_file,'string'))
      setappdata(gcf,'ready',0);
   else
      setappdata(gcf,'ready',1);
   end

   return;					% select_dir

%  select 'Files':
%  update 'File selection'
%
% --------------------------------------------------------------------

function select_file()

   setappdata(gcf,'ready',1);
   listed_file = get(gcbo,'String');
   selected_file_idx = get(gcbo,'Value');
   selected_file = listed_file{selected_file_idx};
   curr_dir = getappdata(gcf,'curr_dir');

   if strcmp(curr_dir(end),filesep)		% filepath end with filesep
      curr_dir = curr_dir(1:end-1);	% take filesep out
   end

   h_select = getappdata(gcf,'h_select');
   set(h_select,'string',[curr_dir filesep selected_file]);

   return;					% select_file

%  change 'File selection':
%  if it is a file, select that,
%  if it is more than a file (*), select those,
%  if it is a directory, select based on file pattern
%
% --------------------------------------------------------------------

function EditSelection()

   filter_string = getappdata(gcf,'filter_string');
   h_select = getappdata(gcf,'h_select');
   selected_file = get(h_select,'string');

   if exist(selected_file) == 7			% if user enter a dir
      setappdata(gcf,'ready',0);
      setappdata(gcf,'curr_dir',selected_file);		% get new dir
      update_dirlist;
   else

      setappdata(gcf,'ready',1);

      [fpath fname fext]= fileparts(selected_file);
      if exist(fpath) ~=7			% fpath is not a dir
         setappdata(gcf,'ready',0);
         msg = 'ERROR: Cannot access directory';
         uiwait(msgbox(msg,'File Selection Error','modal'));
      end

      %  if the file format user entered is not supported by matlab
      if isempty(strmatch(['*',fext],filter_string,'exact'))
         setappdata(gcf,'ready',0);
         msg = 'ERROR: File format is not supported by Matlab.';
         uiwait(msgbox(msg,'File Selection Error','modal'));
      end

   end

   return;					% EditSelection

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

function delete_fig()

   try
      load('pls_profile');
      pls_profile = which('pls_profile.mat');

      rri_select_file_pos = get(gcbf,'position');

      save(pls_profile, '-append', 'rri_select_file_pos');
   catch
   end

   return;


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   example 2 DOF planar robot.
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   03/03/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

robot.name='Example 2DOF planar arm';

%kinematic data DH parameters
robot.DH.theta='[q(1) q(2)]';
robot.DH.d='[0  0]';
robot.DH.a='[1  1]';
robot.DH.alpha='[0  0]';

%number of degrees of freedom
robot.DOF = 2;

%rotational: R, translational: T
robot.kind=['R' 'R'];

%Jacobian matrix. It is easy to obtain the Jacobian matrix of this robot
% J is defined such that:
% v=[vn wn]', where vn is the linear speed of the end effector and wn in
% the angular speed. vn=[vx vy vz] and wn = [wx wy wz]. Being this a planar
% robot, vz=0 and wx=wy=0
robot.J=['[-a(1)*sin(q(1))-a(2)*sin(q(1)+q(2))  -a(2)*sin(q(1)+q(2));' ...
          'a(1)*cos(q(1))+a(2)*cos(q(1)+q(2))   a(2)*cos(q(1)+q(2));' ...
          '               0                                  0;' ...
          '               0                                  0;' ...
          '               0                                  0;' ...
          '               1                                  1]'];
%Function name to compute inverse kinematic
robot.inversekinematic_fn = 'inversekinematic_2dofplanar(robot, T)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-180) deg2rad(180); %Axis 1, minimum, maximum
                deg2rad(-180) deg2rad(180)]; %Axis 2, minimum, maximum

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(100);
                deg2rad(100)];

robot.accelmax=robot.velmax/3; % 0.1 is here an acceleration time

% end effectors maximum velocity
robot.linear_velmax = 0.5; %m/s,

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GRAPHICS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read graphics files
robot.graphical.has_graphics=1;
robot.graphical.color = [25 20 40];
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis = [-2.2 2.2 -2.2 2.2 0 2.2];
robot = read_graphics(robot);

% INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
% position, velocity and acceleration
robot.q=[0 0]';
robot.qd=[0 0]';
robot.qdd=[0 0]';
robot.time = [];

robot.q_vector=[];
robot.qd_vector=[];
robot.qdd_vector=[];

robot.last_target=directkinematic(robot, robot.q);
robot.last_zone_data = 'fine';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DYNAMIC PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robot.has_dynamics=1;

%consider friction or not
robot.dynamics.friction=0;

%link masses (kg)
robot.dynamics.masses=[2 2];

%COM of each link with respect to own reference system
robot.dynamics.r_com=[-0.5      0         0; %(rx, ry, rz) link 1, w/r to reference system 1
                      -0.5      0         0];%(rx, ry, rz) link 2

%link masses
m1 = robot.dynamics.masses(1);
m2 = robot.dynamics.masses(2);

%eval a to obtain parameters.
a=eval(robot.DH.a);
L1 = a(1);
L2 = a(2);

%Inertia matrices of each link
% Ixx	Iyy	Izz	Ixy	Iyz	Ixz, or each row
robot.dynamics.Inertia=[0   m1*L1^2/12   m1*L1^2/12    0	0	0;
                        0   m2*L2^2/12   m2*L2^2/12    0	0	0];

%Inertia of the rotor
robot.motors.Inertia=[0 0];
%Reduction ratio: motor_speed/joint speed
robot.motors.G=[10  10];

%Viscous friction factor of the motor
robot.motors.Viscous = [0.9  0.9];
%Coulomb friction of the motor
%Tc+, Tc-
robot.motors.Coulomb = [0	0;
                        0	0];

%Obtained from motor catalog under practicals/inverse_dynamics
%                        R(Ohm)  L(H)      Kv (V/rad/s):speed constant     Kp (Nm/A):torque constant        Max_current (A)
robot.motors.constants=[0.345  0.273e-3       2.3474e-05               84.9e-3                 139];%these correspond to Maxon, 167132;

repo: ARTE


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2017 Lime Microsystems Ltd.
%
% Licensed under the Apache License, Version 2.0 (the 'License');
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%    http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an 'AS IS' BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function NumptyGraphs(fpltname,iqDataTx,iqDataRxRaw,iqDataRx,binData,din,dout,FSR,osr,msgLen,saveGraphs)
  %
  %   D I S P L A Y   G R A P H I C A L   S U M M A R Y
  %
  subplot(2,3,1);
  plot(iqDataTx(1:((msgLen+1)*osr)),'b');
  xlabel('Time');
  ylabel('Level');
  title('Filtered BPSK Message');
  grid on;
  subplot(2,3,2);
  plot(binData,'go',dout,'rx-');
  legend('transmitted','received');
  xlabel('Time');
  ylabel('Level');
  title('Binary Message');
  grid on;
  subplot(2,3,3);
  fsc=FSR*((-(length(iqDataRx)-1)/2):((length(iqDataRx)-1)/2))/length(iqDataRx);
  fftr=fftshift(mag2db(abs(fft(iqDataTx)/length(iqDataTx))+1e-5));
  plot(fsc,fftr);
  grid on;
	title('Spectrum of RRC Filtered Tx Signal');
	xlabel('Frequency MHz');
	ylabel('Level (dB)');
	axis('tight');
  subplot(2,3,4);
  fftr=fftshift(mag2db(abs(fft(iqDataRxRaw)/length(iqDataRxRaw))+1e-5));
  plot(fsc,fftr);
  grid on;
	title('Spectrum of Raw Rx Signal');
	xlabel('Frequency MHz');
	ylabel('Level (dB)');
	axis('tight');
  subplot(2,3,5);
  fftr=fftshift(mag2db(abs(fft(iqDataRx)/length(iqDataRx))+1e-5));
  plot(fsc,fftr);
  grid on;
	title('Spectrum of RRC Filtered Rx Signal');
	xlabel('Frequency MHz');
	ylabel('Level (dB)');
	axis('tight');
  subplot(2,3,6);
 	plot(iqDataRx,'b-',din,'rx');
	grid on;
	title('Trajectory of RRC Filtered RX signal');
	xlabel('I level');
	ylabel('Q level');
  legend('Traj','Locked');
  if saveGraphs
    print( fpltname, '-dpng' );
  end
end

author: myriadrf


text:
function [data] = ft_rejectcomponent(cfg, comp, data)

% FT_REJECTCOMPONENT backprojects an ICA (or similar) decomposition to the
% channel level after removing the independent components that contain
% the artifacts. This function does not automatically detect the artifact
% components, you will have to do that yourself.
%
% Use as
%    [data] = ft_rejectcomponent(cfg, comp)
% or as
%    [data] = ft_rejectcomponent(cfg, comp, data)
%
% where the input comp is the result of FT_COMPONENTANALYSIS. The output
% data will have the same format as the output of FT_PREPROCESSING.
%
% An optional input argument data can be provided. In that case
% componentanalysis will do a subspace projection of the input data
% onto the space which is spanned by the topographies in the unmixing
% matrix in comp, after removal of the artifact components.  Please use
% this option of including data as input, if you wish to use the output
% data.grad in further computation, for example for leadfield computation.
%
% The configuration structure can contain
%   cfg.component  = list of components to remove, e.g. [1 4 7] or see FT_CHANNELSELECTION
%   cfg.demean     = 'no' or 'yes', whether to demean the input data (default = 'yes')
%   cfg.updatesens = 'no' or 'yes' (default = 'yes')
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_COMPONENTANALYSIS, FT_PREPROCESSING

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

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

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

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

% set the defaults
cfg.component       = ft_getopt(cfg, 'component',  []);
cfg.demean          = ft_getopt(cfg, 'demean',     'yes');
cfg.feedback        = ft_getopt(cfg, 'feedback',   'text');
cfg.updatesens      = ft_getopt(cfg, 'updatesens', 'yes');

% the data can be passed as input arguments or can be read from disk
hascomp = exist('comp', 'var');
hasdata = exist('data', 'var');

if hascomp && hasdata
  % check if the input data is valid for this function
  istlck  = ft_datatype(data, 'timelock');  % this will be temporary converted into raw
  data    = ft_checkdata(data, 'datatype', 'raw');
  comp    = ft_checkdata(comp, 'datatype', 'comp');
  label   = data.label;
  nchans  = length(data.label);
  ncomps  = length(comp.label);
elseif hascomp
  % check if the input data is valid for this function
  istlck  = ft_datatype(comp, 'timelock');  % this will be temporary converted into raw
  comp    = ft_checkdata(comp, 'datatype', 'raw+comp');
  label   = comp.topolabel;
  ncomps  = length(comp.label);
else
  ft_error('incorrect number of input arguments');
end

% cfg.component can be indicated by number or by label
cfg.component = ft_channelselection(cfg.component, comp.label);
reject = match_str(comp.label, cfg.component);

if isempty(reject)
  ft_warning('no components were selected for rejection');
end

if min(reject)<1
  ft_error('you cannot remove components that are not present in the data');
end

if max(reject)>ncomps
  ft_error('you cannot remove components that are not present in the data');
end

if hasdata && strcmp(cfg.demean, 'yes')
  % optionally perform baseline correction on each trial
  ft_info('baseline correcting data \n');
  for trial=1:numel(data.trial)
    data.trial{trial} = ft_preproc_baselinecorrect(data.trial{trial});
  end
end

% set the rejected component amplitudes to zero
ft_info('removing %d components\n', length(reject));
ft_info('keeping %d components\n',  ncomps-length(reject));

% create a projection matrix by subtracting the subspace spanned by the
% topographies of the to-be-removed components from identity
[seldat, selcomp] = match_str(label, comp.topolabel);

if hasdata && length(seldat)~=length(label)
  ft_warning('the subspace projection is not guaranteed to be correct for non-orthogonal components');
end

if hasdata
  mixing   = comp.topo(selcomp,:);
  unmixing = comp.unmixing(:,selcomp);

  % I am not sure about this, but it gives comparable results to the ~hasdata case
  % when comp contains non-orthogonal (=ica) topographies, and contains a complete decomposition

  montage     = [];
  montage.tra = eye(length(selcomp)) - mixing(:, reject)*unmixing(reject, :);
  % we are going from data to components, and back again
  montage.labelold = comp.topolabel(selcomp);
  montage.labelnew = comp.topolabel(selcomp);

  keepunused = 'yes'; % keep the original data which are not present in the mixing provided
  bname = 'reject';

else
  mixing = comp.topo(selcomp, :);
  mixing(:, reject) = 0;

  montage     = [];
  montage.tra = mixing;
  % we are going from components to data
  montage.labelold = comp.label;
  montage.labelnew = comp.topolabel(selcomp);

  keepunused = 'no'; % don't need to keep the original rejected components
  bname = 'invcomp';

  % create the initial data structure, remove all component details
  data = keepfields(comp, {'trial', 'time', 'label', 'fsample', 'grad', 'elec', 'opto', 'trialinfo', 'sampleinfo'});
end % if hasdata

% apply the linear projection to the data
data = ft_apply_montage(data, montage, 'keepunused', keepunused, 'feedback', cfg.feedback, 'showcallinfo', cfg.showcallinfo);

sensfield = cell(0,1);
if isfield(data, 'grad')
  sensfield{end+1} = 'grad';
end
if isfield(data, 'elec')
  sensfield{end+1} = 'elec';
end
if isfield(data, 'opto')
  sensfield{end+1} = 'opto';
end

% apply the linear projection also to the sensor description
if ~isempty(sensfield)
  if  strcmp(cfg.updatesens, 'yes')

    for m = 1:numel(sensfield)
      ft_info('also applying the backprojection matrix to the %s structure\n', sensfield{m});

      % the balance field is needed to keep the sequence of linear projections
      if ~isfield(data.(sensfield{m}), 'balance')
        data.(sensfield{m}).balance.current = 'none';
      end

      % keepunused = 'yes' is required to get back e.g. reference or otherwise
      % unused sensors in the sensor description. the unused components need to
      % be removed in a second step
      sens = ft_apply_montage(data.(sensfield{m}), montage, 'keepunused', 'yes', 'balancename', bname, 'feedback', cfg.feedback);

      % remove the unused channels from the grad/elec/opto
      [junk, remove]    = match_str(comp.label, sens.label);
      sens.tra(remove,:) = [];
      sens.label(remove) = [];
      sens.chanpos(remove,:) = [];
      if isfield(sens, 'chanori')
        sens.chanori(remove,:) = [];
      end

      % there could have been sequential subspace projections, so the
      % invcomp-field may have been renamed into invcompX. If this it the case,
      % take the one with the highest suffix
      invcompfield = bname;
      if  ~isfield(sens.balance, invcompfield)
        for k = 10:-1:1
          if isfield(sens.balance, [bname num2str(k)])
            invcompfield = [invcompfield num2str(k)];
            break;
          end

        end
      end

      % remove the unused components from the balancing
      [junk, remove]    = match_str(comp.label, sens.balance.(invcompfield).labelnew);
      sens.balance.(invcompfield).tra(remove, :)   = [];
      sens.balance.(invcompfield).labelnew(remove) = [];
      data.(sensfield{m})  = sens;
    end

  else
    for m = 1:numel(sensfield)
      ft_info('not applying the backprojection matrix to the %s structure\n', sensfield{m});
      % simply copy it over
      comp.(sensfield{m}) = data.(sensfield{m});
    end
  end
end % if sensfield

if istlck
  % convert the raw structure back into a timelock structure
  data = ft_checkdata(data, 'datatype', 'timelock');
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   comp data
ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function img = fdct_usfft_dispcoef(C)

% fdct_usfft_dispcoef - Returns an image containing all the curvelet coefficients
%
% Inputs
%     C         Curvelet coefficients
%
% Outputs
%     img       Image containing all the curvelet coefficents. The coefficents are rescaled so that
%               the largest coefficent in each subband has unit norm.
%

  [m,n] = size(C{end}{1});
  nbscales = floor(log2(min(m,n)))-3;

  img = C{1}{1};  img = img/max(max(abs(img))); %normalize
  for sc=2:nbscales-1
    nd = length(C{sc})/4;
    wcnt = 0;

    ONE = [];
    for w=1:nd
      ONE = [ONE, C{sc}{wcnt+w}];
    end
    wcnt = wcnt+nd;

    TWO = [];
    for w=1:nd
      TWO = [TWO; C{sc}{wcnt+w}];
    end
    wcnt = wcnt+nd;

    THREE = [];
    for w=1:nd
      THREE = [C{sc}{wcnt+w}, THREE];
    end
    wcnt = wcnt+nd;

    FOUR = [];
    for w=1:nd
      FOUR = [C{sc}{wcnt+w}; FOUR];
    end
    wcnt = wcnt+nd;

    [p,q] = size(img);
    [a,b] = size(ONE);
    [g,h] = size(TWO);
    m = 2*a+g;    n = 2*h+b; %size of new image
    scale = max(max( max(max(abs(ONE))),max(max(abs(TWO))) ), max(max(max(abs(THREE))), max(max(abs(FOUR))) )); %scaling factor

    new = 0.5 * ones(m,n); %background value
    new(a+1:a+g,1:h) = FOUR /scale;
    new(a+g+1:2*a+g,h+1:h+b) = THREE /scale;
    new(a+1:a+g,h+b+1:2*h+b) = TWO /scale;
    new(1:a,h+1:h+b) = ONE /scale; %normalize

    dx = floor((g-p)/2);    dy = floor((b-q)/2);

    new(a+1+dx:a+p+dx,h+1+dy:h+q+dy) = img;

    img = new;
  end

author: thomaskuestner


text:
function [xopt,fopt]=simann(func, x, LB, UB, sa_t, sa_rt, sa_nt, sa_ns,rseed)

% Simulated Annealing programmed for minimization problem
% INPUTS
% func, string variable containing name of function file to be optimized
% x, starting values
% LB, lower bound on optimization parameters
% UB, upper bound on optimization parameters
% sa_t, initial temperature
% sa_rt, temperature reduction factor, 0 < sa_rt < 1, try .85
% sa_nt, number of times through ns loop before temperature reduction (recommended value: 5)
% sa_ns, number of times through function before stepsize adjustment (recommended value: 20)
%
% OUTPUTS
% xopt, the optimal solution
%
%

LB=LB(:)';
UB=UB(:)';

rand('state',rseed);                      %sets seed for random number generator
sa_neps=4;                                %number of times eps
                                          %tolerance is achieved before termination
sa_eps=eps;                              %convergence criteria
sa_maxeval=60;%12000000;                      %maximum number of function evaluations

sa_nargs=length(LB);                      %number of parameters
sa_nobds=0;
sa_nacc=0;                                %number of acceptions
sa_nevals=0;                              %number of evaluations
sa_opteval=0;                             %optimum number of
                                          %function evaluations

fstar=Inf*ones(sa_neps,1);

%x=LB+(UB-LB).*rand(1, sa_nargs);         %starting values for model parameters
f=feval(func,x);                          %function evaluation with parameters x
%disp('initial loss function value:');disp(f);
sa_nevals=sa_nevals+1;
xopt=x;
fopt=f;
xtot=x;
fstar(1)=f;

VM=(UB-LB);%/2;                      %maximum step size

%LOOP
while 1

  nup=0;                                   %number of uphill movements
  nrej=0;                                  %number of rejections
  nnew=0;
  ndown=0;                                 %number of downhill movements
  lnobds=0;
  nacp=zeros(sa_nargs,1);
  C = progress('init','Determine initial hyperparameters for simplex...');
  for m=1:sa_nt
    for j=1:sa_ns
      for h=1:sa_nargs
        if sa_nevals>=sa_maxeval
          %disp('too many function evaluations')
          return
        end
        C = progress(C,sa_nevals/sa_maxeval);
        %workbar(sa_nevals/sa_maxeval,'Determine initial hyperparameters to build grid...','Progress')
        % generate xp, trial value of x
        xp=x;
        xp(h)=x(h)+VM(h)*(2*rand(1,1)-1.0);               %calculate new value for x (xp)
        if (xp(h)<LB(h)) | (xp(h)>UB(h))
          xp(h)=LB(h)+(UB(h)-LB(h))*rand(1,1);
          lnobds=lnobds+1;
          sa_nobds=sa_nobds+1;
        end
        % evaluate at xp and return as fp
        %disp ('current parameter vector:');disp(xp);
        fp=feval(func,xp);                                 %function evaluation with parameters xp
        %disp ('function value');disp(fp);
        sa_nevals=sa_nevals+1;

        % we minimize! accept if the function value decreases
        if fp<=f
          x=xp;
          f=fp;
          sa_nacc=sa_nacc+1;
          nacp(h)=nacp(h)+1;
          nup=nup+1;
          % if smaller than any previous point, record as new optimum
          if fp<fopt
            xopt=xp;
            fopt=fp;
            sa_opteval=sa_nevals;
            nnew=nnew+1;
          end
        else % function value increases
          p=exp((f-fp)/sa_t);                              %random number
          pp=rand(1,1);
          if pp<p
            x=xp;
            f=fp;
            sa_nacc=sa_nacc+1;
            nacp(h)=nacp(h)+1;
            ndown=ndown+1;
          else
            nrej=nrej+1;
          end
        end
      end
    end

    % adjust maximal step size vm
    c=ones(sa_nargs,1)*2;%??
    for i=1:sa_nargs
      ratio=nacp(i)/sa_ns;
      if ratio>0.6
        VM(i)=VM(i) * (1+c(i)*(ratio-0.6)/0.4);
      elseif ratio <0.4
        VM(i)=VM(i)/(1+c(i)*((0.4-ratio)/0.4));
      end
      if VM(i)>(UB(i)-LB(i))
        VM(i)=UB(i)-LB(i);
      end
    end

    % provide statistics about current state of optimization

%     disp('No. of evaluations');disp(sa_nevals);disp('  current temperature');disp(sa_t);
%     disp('current optimum function value');disp(fopt);
%     disp('No. of downhill steps');disp(nup);  % note misnomer in variable declaration!
%     disp('No. of accepted uphill steps');disp(ndown); % we minimize, thus downhill is always accepted!
%     disp('No. of rejections');disp(nrej);
%     disp('current parameter values');disp(xp);
%     disp('current optimum vector');disp(xopt);
%     disp('current step size');disp(VM);
    %disp('Variables used:');whos;

  for i=1:sa_nargs
     nacp(i) = 0;
  end
  end

  % check termination criteria
  fstar(1)=f;
  quit = ((fstar(1)-fopt) <= sa_eps);
  if any(abs(fstar-f)>sa_eps)
    quit=0;
  end

  if quit
    disp(['simulated annealing achieved termination after ', num2str(sa_nevals),' evals']);
    return
  end

  % reduce temperature
  sa_t=sa_t*sa_rt;
  fstar(2:4)=fstar(1:3);
  % continue from current optimum
  x=xopt;
  f=fopt;
end %while

author: peterhcharlton

repo: RRest


text:
function ehmm = hsupdate_ehmm(Xi,Gamma,T,ehmm)
%
% updates hidden state parameters of an HMM
%
% INPUT:
%
% Xi     probability of past and future state cond. on data
% Gamma  probability of current state cond. on data
% T      length of observation sequences
% hmm    single hmm data structure
%
% OUTPUT
% hmm    single hmm data structure with updated state model probs.
%
% Author: Diego Vidaurre, OHBA, University of Oxford

K = ehmm.K;

if isempty(Xi)  % non-exact estimation
    order = (sum(T)-size(Gamma,1)) / length(T);
    Xi = zeros(sum(T-1-order),K,2,2);
    for j = 1:length(T)
        indG = (1:(T(j)-order)) + sum(T(1:j-1)) - (j-1)*order;
        indXi =  (1:(T(j)-order-1)) + sum(T(1:j-1)) - (j-1)*(order+1);
        for k = 1:K
            g = [Gamma(indG,k) (1-Gamma(indG,k))];
            for t = 1:length(indXi)
                xi = g(t,:)' * g(t+1,:);
                xi = xi / sum(xi(:));
                Xi(indXi(t),k,:,:) = xi;
            end
        end
    end
end
if length(size(Xi))==4
    Xi = permute(sum(Xi),[2 3 4 1]);
end

% transitions
for k = 1:K
    ehmm.state(k).Dir2d_alpha = permute(Xi(k,:,:),[2 3 1]) ...
        + ehmm.state(k).prior.Dir2d_alpha;
    ehmm.state(k).P = zeros(2);
    for j = 1:2
        PsiSum = psi(sum(ehmm.state(k).Dir2d_alpha(j,:)));
        for j2 = 1:2
            ehmm.state(k).P(j,j2) = ...
                exp(psi(ehmm.state(k).Dir2d_alpha(j,j2))-PsiSum);
        end
        ehmm.state(k).P(j,:) = ehmm.state(k).P(j,:) ./ sum(ehmm.state(k).P(j,:));
    end
end

% initial state is always OFF for all chains
end

author: OHBA-analysis

repo: HMM-MAR


text:
function y = isnan(X)
    % True for Not-a-Number for symbolic arrays.
    %   ISNAN(X) returns an array that contains 1's where
    %   the elements of X are symbolic NaN's and 0's where they are not.
    %   For example, ISNAN(sym([pi NaN Inf -Inf])) is [0 1 0 0].
    %
    %   For any X, exactly one of ISFINITE(X), ISINF(X), or ISNAN(X)
    %   is 1 for each element.
    %
    %   See also SYM/ISFINITE, SYM/ISINF.

    % Convert inputs to SymExpression
    % X = SymExpression(X);

    % evaluate the operation in Mathematica and return the
    % expression string
    ret = eval_math(['Boole/@Map[SameQ, ' X.s ', NaN]'], 'math2matlab');

    y = logical(ret);
end

author: ayonga

repo: frost-dev

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


text:
function [jaccards,inters,false_pos,false_neg,true_areas] = eval_one(proposals, ground_truth)

    n_objs = length(ground_truth.masks);

    % Store true_areas
    true_areas = zeros(n_objs,1);
    for kk=1:n_objs
        true_areas(kk) = sum(ground_truth.masks{kk}(:));
    end

    % Case when proposals is a subfield
    if isfield(proposals,'proposals')
        proposals = proposals.proposals;
    end

    % Which type of result are we evaluating?
    % - Labels: Superpixel matrix + labels for each proposal
    % - Masks : 3D matrix of boolean masks
    % - Blobs : Boolean masks in a bounding box
    % You can add a case in this file to adapt the evaluation to your type
    % of result
    if isfield(proposals,'labels') || isfield(proposals,'indicator_matrix')
        if ~isfield(proposals,'labels')
            if isa(proposals.superpixels,'uint16')
                superpixels = uint32(proposals.superpixels); % That's RIGOR
            else
                superpixels = uint32(proposals.superpixels'); % That's GOP, transposed
            end

            assert(length(unique(superpixels))==size(proposals.indicator_matrix,1))
            n_cands = size(proposals.indicator_matrix,2);
            labels = cell(n_cands,1);
            for jj=1:n_cands
               labels{jj} = uint32(find(proposals.indicator_matrix(:,jj)'));
            end
        else
            labels = proposals.labels;
            superpixels = proposals.superpixels;
        end

        % Handle all as the multiple-superpixels case
        % [jaccards,inters,false_pos,false_neg,true_areas]
        if ~iscell(superpixels)
            superpixels = {superpixels};
            labels = {labels};
        end

        jaccards  = [];
        inters    = [];
        false_pos = [];
        false_neg = [];

        % Sweep all 'models' (Done for LPO)
        for mm=1:length(superpixels)
            % Overlap superpixels with GT objects
            n_leaves = length(unique(superpixels{mm}));
            superpixels{mm}(~ground_truth.valid_pixels) = 0;
            leave_fp   = zeros(n_objs,n_leaves);  % False positives
            leave_int  = zeros(n_objs,n_leaves);  % Intersection with GT
            for kk=1:n_objs
                tmp = hist(double(superpixels{mm}(:)).*(ground_truth.masks{kk}(:)),(0:n_leaves));
                leave_int(kk,:) = tmp(2:end);
                tmp = hist(double(superpixels{mm}(:)).*(~ground_truth.masks{kk}(:)),(0:n_leaves));
                leave_fp(kk,:) = tmp(2:end);
            end

            % Create matrix padded with zeros to be compatible with mex_eval_labels
            n_proposals  = length(labels{mm});
            n_max_labels = length(unique(superpixels{mm}));
            label_matrix = zeros(n_proposals,n_max_labels);
            for jj=1:n_proposals
                label_matrix(jj,1:length(labels{mm}{jj})) = labels{mm}{jj};
            end

            % Compute fp, fn etc. from these values on the superpixels
            this_inters     = zeros(n_objs,n_proposals);
            this_false_pos  = zeros(n_objs,n_proposals);
            this_false_neg  = zeros(n_objs,n_proposals);
            this_jaccards   = zeros(n_objs,n_proposals);
            if n_proposals>0
                for kk=1:n_objs
                    [this_inters(kk,:),this_false_pos(kk,:)] = mex_eval_labels(leave_int(kk,:),leave_fp(kk,:),label_matrix);
                    this_false_neg(kk,:) = true_areas(kk)-this_inters(kk,:);
                    this_jaccards(kk,:) = this_inters(kk,:)./(this_false_pos(kk,:)+true_areas(kk));
                end
            end
            assert(sum(isnan(this_jaccards(:)))==0)

            % Accumulate
            jaccards  = [jaccards  this_jaccards]; %#ok<AGROW>
            inters    = [inters    this_inters];   %#ok<AGROW>
            false_pos = [false_pos this_false_pos];%#ok<AGROW>
            false_neg = [false_neg this_false_neg];%#ok<AGROW>
        end
    elseif isfield(proposals,'masks')
        % Sweep and compare all masks
        [areas,inters,false_neg] = mex_eval_masks(proposals.masks,ground_truth.masks,ground_truth.valid_pixels);
        false_pos = repmat(areas,n_objs,1)-inters;
        jaccards = inters./(inters+false_pos+false_neg);
        assert(sum(isnan(jaccards(:)))==0)
    elseif isfield(proposals,'hBlobs')
        % Sweep and compare all blobs
        [areas,inters,false_neg] = mex_eval_blobs(proposals.hBlobs,ground_truth.masks,ground_truth.valid_pixels);
        false_pos = repmat(areas,n_objs,1)-inters;
        jaccards = inters./(inters+false_pos+false_neg);
        assert(sum(isnan(jaccards(:)))==0)
    end
end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function x = solve_triu(T,b)
% SOLVE_TRIU      Left division by upper triangular matrix.
% SOLVE_TRIU(T,b) is the same as T\b but requires T to be upper triangular
% and runs faster.

error('You must run install_lightspeed in order to use this function.')

author: andrewssobral

repo: lrslibrary


text:
function i4vec_sort_heap_a_test ( )

%*****************************************************************************80
%
%% I4VEC_SORT_HEAP_A_TEST tests I4VEC_SORT_HEAP_A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 20;
  b = 0;
  c = 3 * n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_SORT_HEAP_A_TEST\n' );
  fprintf ( 1, '  For a vector of integers,\n' );
  fprintf ( 1, '  I4VEC_SORT_HEAP_A ascending sorts,\n' );

  seed = 123456789;

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  Unsorted:' );

  a = i4vec_sort_heap_a ( n, a );

  i4vec_print ( n, a, '  Ascending sorted:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function shCoeff = ShearTransform3D(X,shearingFilter)
level=size(shearingFilter,2);
shCoeff.D=cell(3,level);

%Do the Band Pass of noisy Data
BP=DoPyrDec(X,level);
shCoeff.A=BP{1};
for pyrCone=1:3
  shCoeff.D(pyrCone,:) = ShDec(pyrCone,shearingFilter,BP,level,'single');
end


text:
function line_cvt_lloyd_test ( )

%*****************************************************************************80
%
%% LINE_CVT_LLOYD_TEST tests the line_cvt_lloyd library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_CVT_LLOYD_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the LINE_CVT_LLOYD library.\n' );

  line_cvt_lloyd_test01 ( );
  line_cvt_lloyd_test02 ( );
%
%  Repeat, using sorted initial points.
%
  line_cvt_lloyd_test03 ( );
  line_cvt_lloyd_test04 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_CVT_LLOYD_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
%% Experiment with the cnn_mnist_fc_bnorm

[net_bn, info_bn] = cnn_mnist(...
  'expDir', 'data/mnist-bnorm', 'useBnorm', true);

[net_fc, info_fc] = cnn_mnist(...
  'expDir', 'data/mnist-baseline', 'useBnorm', false);

figure(1) ; clf ;
subplot(1,2,1) ;
semilogy(info_fc.val.objective', 'o-') ; hold all ;
semilogy(info_bn.val.objective', '+--') ;
xlabel('Training samples [x 10^3]'); ylabel('energy') ;
grid on ;
h=legend('BSLN', 'BNORM') ;
set(h,'color','none');
title('objective') ;
subplot(1,2,2) ;
plot(info_fc.val.error', 'o-') ; hold all ;
plot(info_bn.val.error', '+--') ;
h=legend('BSLN-val','BSLN-val-5','BNORM-val','BNORM-val-5') ;
grid on ;
xlabel('Training samples [x 10^3]'); ylabel('error') ;
set(h,'color','none') ;
title('error') ;
drawnow ;

author: willard-yuan

repo: cnn-for-image-retrieval


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

function s = python_header()
%private

  persistent PyStrCache

  if (isempty(PyStrCache))
    thisdir = fileparts (mfilename ('fullpath'));
    pyfile = fullfile (thisdir, 'python_header.py');
    PyStrCache = fileread (pyfile);
  end

  s = PyStrCache;

end

repo: octsympy


text:
function f = define(f, s ,v) %#ok<INUSD>
%DEFINE   Deprecated function.
%   DEFINE(...) is deprecated and has been replaced by DEFINEPOINT() and
%   DEFINEINTERVAL(). However, most users will interact with both of these
%   methods via CHEBFUN/SUBSREF.
%
% See also DEFINEPOINT, DEFINEINTERVAL, CHEBFUN/SUBSREF.

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

error('CHEBFUN:CHEBFUN:define:deprecated', ...
    'CHEBFUN/DEFINE is deprecated. Use DEFINEPOINT or DEFINEINTERVAL.');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function b=num2ibm(x)
% num2ibm : convert IEEE 754 doubles to IBM 32 bit floating point format
%    b=num2ibm(x)
% x is a matrix of doubles
% b is a corresponding matrix of uint32
%
% The representations for NaN and inf are arbitrary
%
% See also ibm2num

%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
%
% (C) Brian Farrelly, 22 October 2001
%  mailto:Brian.Farrelly@nho.hydro.com          Norsk Hydro Research Centre
%  phone +47 55 99 68 74                 (((                  Postboks 7190
%  fax   +47 55 99 69 70                2oooS                 N-5020 Bergen
%  home  +47 55 13 78 49                HYDRO                        Norway
%

b=repmat(uint32(0),size(x));
err=zeros(size(x));

%format long

x(x> 7.236998675585915e+75)= inf;        % change big numbers to infinity
x(x<-7.236998675585915e+75)=-inf;        % 7.236998675585915e+75 is
                                         %    ibm2num(uint32(hex2dec('7fffffff')) or
                                         %    ibm2num(num2ibm(inf))

[F E]=log2(abs(x));

e=E/4;                         % exponent of base 16
ec=ceil(e);                    % adjust upwards to integer
p=ec+64;                       % offset exponent

f=F.*2.^(-4*(ec-e));           % correct mantissa for fractional part of exponent
f=round(f*2^24);               % convert to integer. Roundoff here can be as large as
                               % 0.5/2^20 when mantissa is close to 1/16 so that
                               % 3 bits of signifance are lost.

p(f==2^24)=p(f==2^24)+1;       % Roundoff can cause f to be 2^24 for numbers just under a
f(f==2^24)=2^20;               % power of 16, so correct for this

%format hex
psi=uint32(p*2^24);            % put exponent in first byte of psi.
phi=uint32(f);                 % put mantissa into last 3 bytes of phi

% make bit representation

b=bitor(psi,phi);                        % exponent and mantissa
b(x<0)=bitset(b(x<0),32);                % sign bit
%format long

% special cases

b(x==0)          =uint32(0)                  ;         %  bias is incorrect for zero
b(isnan(x))      =uint32(hex2dec('7fffffff'));         %  7.237005145973116e+75 in IBM format
b(isinf(x) & x>0)=uint32(hex2dec('7ffffff0'));         %  7.236998675585915e+75    ,,
b(isinf(x) & x<0)=uint32(hex2dec('fffffff0'));         % -7.236998675585915e+75    ,,
                                                       % Note that NaN > inf in IBM format

% check bit representation for normal cases

checkx=ibm2num(b);                      % note that use of base 16 in IBM format
z=x==0;                                 % can lead to a loss of 3 bits of precision
err(z)=0;                               % compared with an IEEE single.
q=(checkx(~z)-x(~z))./x(~z);
err(~z) = abs(q) > 5e-7;                % this is almost reached with numbers
                                        % of the form 16^n + 0.5*16^(n-5) where
                                        % the mantissa is 100001 hex. Roundoff
                                        % error is then 0.5/16^5=0.5/2^20=4.7684e-7

if any(err)
   warning('Conversion error in num2ibm for the following:')
   disp(x(logical(err)))
end

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function value = r8_beta ( a, b )

%*****************************************************************************80
%
%% R8_BETA evaluates the beta function of R8 arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real A, B, the arguments.
%
%    Output, real VALUE, the beta function of A and B.
%
  persistent alnsml
  persistent xmax

  if ( isempty ( xmax ) )
    [ xmin, xmax ] = r8_gaml ( );
    alnsml = log ( r8_mach ( 1 ) );
  end

  if ( a <= 0.0 || b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_BETA - Fatal error!\n' );
    fprintf ( 1, '  A and B must be greater than 0.\n' );
    error ( 'R8_BETA - Fatal error!' )
  end

  if ( a + b < xmax )
    value = r8_gamma ( a ) * r8_gamma ( b ) / r8_gamma ( a + b );
    return
  end

  value = r8_lbeta ( a, b );

  value = exp ( value );

  return
end


text:
%% Make plots
n_plot = floor(20/dt);
subplot(2,3,1);
hold on;
plot(vTime(1:n_plot),monetary_shock(1:n_plot),'linewidth',1.5,'linestyle',line_style,'color',color);
set(gcf,'color','w');
title('Monetary Policy Shock','interpreter','latex','fontsize',14);
ylabel('$\%$ deviation','interpreter','latex');
xlim([1 vTime(n_plot)]);
grid on;
hold off;

subplot(2,3,2);
hold on;
plot(vTime(1:n_plot),inflation(1:n_plot),'linewidth',1.5,'linestyle',line_style,'color',color);
set(gcf,'color','w');
title('Inflation','interpreter','latex','fontsize',14);
ylabel('$\%$ deviation','interpreter','latex');
xlim([1 vTime(n_plot)]);
grid on;
hold off;

subplot(2,3,3);
hold on;
plot(vTime(1:n_plot),consumption(1:n_plot),'linewidth',1.5,'linestyle',line_style,'color',color);
set(gcf,'color','w');
title('Consumption','interpreter','latex','fontsize',14);
ylabel('$\%$ deviation','interpreter','latex');
xlim([1 vTime(n_plot)]);
grid on;
hold off;

subplot(2,3,4);
hold on;
plot(vTime(1:n_plot),Y(1:n_plot),'linewidth',1.5,'linestyle',line_style,'color',color);
set(gcf,'color','w');
title('GDP','interpreter','latex','fontsize',14);
ylabel('$\%$ deviation','interpreter','latex');
xlim([1 vTime(n_plot)]);
grid on;
hold off;

subplot(2,3,5);
hold on;
plot(vTime(1:n_plot),lab_sup(1:n_plot),'linewidth',1.5,'linestyle',line_style,'color',color);
set(gcf,'color','w');
title('Labor Supply','interpreter','latex','fontsize',14);
ylabel('$\%$ deviation','interpreter','latex');
xlim([1 vTime(n_plot)]);
grid on;
hold off;

subplot(2,3,6);
hold on;
plot(vTime(1:n_plot),wage(1:n_plot),'linewidth',1.5,'linestyle',line_style,'color',color);
set(gcf,'color','w');
title('Wage','interpreter','latex','fontsize',14);
ylabel('$\%$ deviation','interpreter','latex');
xlim([1 vTime(n_plot)]);
grid on;
hold off;

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
function fcn_decodeCommandMoveRobot(str)
%% About:
% This function is used to decode the command string, and then it is used
% to move the robot along the decoded instruction.

%% Areguments:
% str: command string

% Copyright: Mohammad SAFEEA, 13th-Jul-2018

global read_state_var;
% check connection
if fcn_isConnected()
else
    message='Error in function (fcn_decodeCommandMoveRobot), not connected to the robot !!!';
    fcn_errorMessage(message);
    return;
end

% Turn off state read
read_state_var=false;

% decode the motion command
if strncmpi('moveLine',str,max(size('moveLine')))
    moveLinePTP(str);
elseif strncmpi('moveJoints',str,max(size('moveJoints')))
    moveJointsPTP(str);
elseif startsWith(str,'OuputPin')
    setOutputsOfMediaFlange(str);
else
end

% Turn on state read
read_state_var=true;
end

function moveLinePTP(str)
    global iiwa;
    delimiter='_';
    C = fcn_tokenize(str,delimiter);
    for i=2:7
        pos{i-1}=str2num(C{i});
    end
    relVel=10; % velocity mm/sec
    iiwa.movePTPLineEEF(pos, relVel);
end

function moveJointsPTP(str)
    global iiwa;
    delimiter='_';
    C = fcn_tokenize(str,delimiter);
    for i=2:8
        jPos{i-1}=str2num(C{i});
    end
    relVel=0.2; % overrised velocity
    iiwa.movePTPJointSpace(jPos, relVel);
end

%% Control touch-pneumatic media flange output
function setOutputsOfMediaFlange(str)
    global iiwa;
    off={'OuputPin_11_OFF','OuputPin_12_OFF','OuputPin_1_OFF','OuputPin_2_OFF'};
    on={'OuputPin_11_ON','OuputPin_12_ON','OuputPin_1_ON','OuputPin_2_ON'};
    % pin 11
    if strncmpi(off{1},str,max(size(str)))
        iiwa.setPin11Off();
        return;
    end
    if strncmpi(on{1},str,max(size(str)))
        iiwa.setPin11On();
        return;
    end
    % pin 12
	if strncmpi(off{2},str,max(size(str)))
        iiwa.setPin12Off();
        return;
    end
    if strncmpi(on{2},str,max(size(str)))
        iiwa.setPin12On();
        return;
    end
    % pin 1
    if strncmpi(off{3},str,max(size(str)))
        iiwa.setPin1Off();
        return;
    end
    if strncmpi(on{3},str,max(size(str)))
        iiwa.setPin1On();
        return;
    end
    % pin 2
    if strncmpi(off{4},str,max(size(str)))
        iiwa.setPin2Off();
        return;
    end
    if strncmpi(on{4},str,max(size(str)))
        iiwa.setPin2On();
        return;
    end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function val=intPowSin(u,n)
%%INTPOWSIN Evaluate the integral of u^n*sin(u) du.  A definite integral
%           can be evaluated, or an indefinite integral (with a particular
%           additive constant).
%
%INPUTS: u A 2XN (for definite integral) or a 1XN (for indefinite
%          integrals) set of N points. For definite integrals, u(1,:) are
%          the real lower bounds and u(2,:) are the real upper bounds. For
%          indefinite integrals, the integral is evaluated at the points in
%          u. The values in u should be real.
%        n The positive integer exponent of coefficient multiplying the
%          sine term.
%
%OUTPUTS: val The 1XN set of values of the integral of u^n*sin(u).
%
%This function simply implements the recursion that arises from integration
%by parts from basic calculus, as are given in the tables in the back of
%[1].
%
%REFERENCES:
%[1] J. Stewart, Calculus, 7th ed. Belmont, CA: Brooks/Cole, 2012.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(u,1);

if(isempty(u))
   val=[];
   return;
end

if(numDim==1)%An indefinite integral
    val=indefIntPowSin(u,n);
else%A definite integral
    val=indefIntPowSin(u(2,:),n)-indefIntPowSin(u(1,:),n);
end

end

function val=indefIntPowSin(u,n)

if(n==0)
    val=-cos(u);
    return;
elseif(n==1)
    val=sin(u)-u*cos(u);
    return;
end

%If here, n>=1
cosVal=cos(u);
sinVal=sin(u);

val=0;
if(mod(n,2)==0)%If n is even
    endVal=2;
else
    endVal=3;
end

coeffProd=1;
for k=n:-2:endVal
    %Do a sine step and a cosine step together!!!
    val=val-coeffProd*u.^k.*cosVal;
    coeffProd=coeffProd*k;
    val=val+coeffProd*u.^(k-1).*sinVal;
    coeffProd=-coeffProd*(k-1);
end

if(endVal==2)
    %The final integral is over u^0*sin(u)
    val=val-coeffProd*cosVal;
else%The final integral is over u^1*sin(u)
    val=val+coeffProd*(sinVal-u.*cosVal);
end
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 [pnt, tri, srf] = read_bv_srf(filename);

% READ_BV_SRF reads a triangulated surface from a BrainVoyager *.srf file
%
% Use as
%   [pnt, tri] = read_bv_srf(filename) or
%   [pnt, tri, srf] = read_bv_srf(filename)

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

% This documentation originates from
%   http://www.brainvoyager.com/BV2000OnlineHelp/BrainVoyagerWebHelp/mergedProjects/FileFormats/BrainVoyager_File_Formats.htm
%
% BYTES DATA TYPE DEFAULT DESCRIPTION
% 4 float 3 version number
% 4 int 0 reserved, must be '0'
% 4 int  NrOfVertices (number of vertices)
% 4 int  NrOfTriangles (number of triangles)
% 4 float 128.0 MeshCenterX
% 4 float 128.0 MeshCenterY
% 4 float 128.0 MeshCenterZ
% NrOfVertices*4 float  VertexX, sequence of X coordinates of all vertices
% NrOfVertices*4 float  VertexY, sequence of Y coordinates of all vertices
% NrOfVertices*4 float  VertexZ, sequence of Z coordinates of all vertices
% NrOfVertices*4 float  NormalX, sequence of X components of all vertex normals
% NrOfVertices*4 float  NormalY, sequence of Y components of all vertex normals
% NrOfVertices*4 float  NormalZ, sequence of Z components of all vertex normals
% 4 float 0.322 R component of convex curvature color (range: 0.0 - 1-0)
% 4 float 0.733 G component of convex curvature color (range: 0.0 - 1-0)
% 4 float 0.980 B component of convex curvature color (range: 0.0 - 1-0)
% 4 float 0.500 Alpha component of convex curvature color (range: 0.0 - 1-0)
% 4 float 0.100 R component of concave curvature color (range: 0.0 - 1-0)
% 4 float 0.240 G component of concave curvature color (range: 0.0 - 1-0)
% 4 float 0.320 B component of concave curvature color (range: 0.0 - 1-0)
% 4 float 0.500 Alpha component of concave curvature color (range: 0.0 - 1-0)
% NrOfVertices*4 int  MeshColor, sequence of color indices of all vertices (*1)
% 4 int  N, number of nearest neighbors of vertex 1
% 4 int  Nearest neighbor 1 of vertex 1
% : :  :
% 4 int  Nearest neighbor N of vertex 1
% :: ::  ::
% 4 int  N, number of nearest neighbors of vertex 'NrOfVertices'
% 4 int  Nearest neighbor 1 of vertex 'NrOfVertices'
% : :  :
% 4 int  Nearest neighbor N of vertex 'NrOfVertices'
% NrOfTriangels*3*4 int  Sequence of three indices to constituting vertices of each triangle
% 4 int  NrOfTriangleStripElements
% NrOfStripElements*4 int  Sequence of strip elements (if NrOfStripElements > 0)

fid = fopen(filename, 'rb', 'ieee-le');

srf.version_number                                        = fread(fid, 1, 'float');
srf.reserved                                              = fread(fid, 1, 'int'  );
srf.NrOfVertices                                          = fread(fid, 1, 'int'  );
srf.NrOfTriangles                                         = fread(fid, 1, 'int'  );
NrOfVertices  = srf.NrOfVertices;
NrOfTriangles = srf.NrOfTriangles;
srf.MeshCenterX                                           = fread(fid, 1, 'float');
srf.MeshCenterY                                           = fread(fid, 1, 'float');
srf.MeshCenterZ                                           = fread(fid, 1, 'float');
srf.X_coordinates                                         = fread(fid, NrOfVertices, 'float');
srf.Y_coordinates                                         = fread(fid, NrOfVertices, 'float');
srf.Z_coordinates                                         = fread(fid, NrOfVertices, 'float');
srf.X_components                                          = fread(fid, NrOfVertices, 'float');
srf.Y_components                                          = fread(fid, NrOfVertices, 'float');
srf.Z_components                                          = fread(fid, NrOfVertices, 'float');
srf.R_component_of_convex_curvature_color                 = fread(fid, 1, 'float');
srf.G_component_of_convex_curvature_color                 = fread(fid, 1, 'float');
srf.B_component_of_convex_curvature_color                 = fread(fid, 1, 'float');
srf.Alpha_component_of_convex_curvature_color             = fread(fid, 1, 'float');
srf.R_component_of_concave_curvature_color                = fread(fid, 1, 'float');
srf.G_component_of_concave_curvature_color                = fread(fid, 1, 'float');
srf.B_component_of_concave_curvature_color                = fread(fid, 1, 'float');
srf.Alpha_component_of_concave_curvature_color            = fread(fid, 1, 'float');
srf.MeshColor                                             = fread(fid, NrOfVertices, 'int'  );
for i=1:NrOfVertices
  number           = fread(fid, 1, 'int'  );
  srf.neighbour{i} = fread(fid, number, 'int'  );
end
srf.Triangles                                             = fread(fid, [3 NrOfTriangles], 'int'  );
srf.NrOfTriangleStripElements                             = fread(fid, 1, 'int'  );
srf.sequence_of_strip_elements                            = fread(fid, srf.NrOfTriangleStripElements, 'int'  );

fclose(fid);

pnt = [srf.X_coordinates(:) srf.Y_coordinates(:) srf.Z_coordinates(:)];
tri = srf.Triangles' + 1;

author: PatternRecognition

repo: OpenBMI


text:
function printPlot(fileName, directory, directoryHtml)

% PRINTPLOT Print a plot to eps and png files.
% FORMAT
% DESC prints a plot to the specified file name and directory.
% ARG fileName : the base name of the file to use.
% ARG directory : the directory to place the eps files in.
% ARG directoryPng : the directory to place png the file in.
%
% SEEALSO : preparePlot
%
% COPYRIGHT : Neil D. Lawrence, 2008

% NDLUTIL

  global printDiagram

  if isempty(printDiagram)
    printTrue = true; % User hasn't been explicit, assume print.
  else
    printTrue = printDiagram; % Use user preference.
  end
  if printTrue
    if nargin < 2
      directory = '.';
    end
    if nargin < 3
      png = false;
    else
      png = true;
    end
    fprintf('Printing eps plot ...\n');
    print('-depsc', [directory filesep fileName '.eps']) % OCTAVE doesn't
                                                       % include .eps automatically
    if isoctave
      % make a plot with text separated.
      fprintf('Printing epslatex plot ...\n');
      print('-depslatex', [directory filesep fileName '_latex.tex'])
    end
    cmap = colormap;
    if png
      fprintf('Printing png plot ...\n');
      % make smaller for PNG plot.
      pos = get(gcf, 'paperposition');
      origpos = pos;
      pos(3) = pos(3)/2;
      pos(4) = pos(4)/2;
      set(gcf, 'paperposition', pos);
      fontsize = get(gca, 'fontsize');
      set(gca, 'fontsize', fontsize/2);
      lineWidth = get(gca, 'lineWidth');
      set(gca, 'lineWidth', lineWidth*2);
      print('-dpng', [directoryHtml filesep fileName '.png']);
      set(gcf, 'paperposition', origpos);
      set(gca, 'fontsize', fontsize);
      set(gca, 'lineWidth', lineWidth);
    end

    fprintf('Printing black and white eps plot ...\n');
    colormap gray
    print('-deps', [directory filesep fileName 'NoColour.eps']) % Octave
                                                                % doesn't
                                                                % include
                                                                % .eps automatically
    colormap(cmap);
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [Classification_accuracy,p1,m1,classes]=simclass(datalearn, datatest,v,c, measure, p, m,pl)
%Inputs:
% datalearn: put your data file in matrix form, rows indicates samples and
% columns features of the data.
% datatest: put your data file in matrix form, rows indicates samples and
% columns features of the data.
% v:   how many columns of features you are using e.g. [1:4] means you will
% use first four columns as features of your data.
% c: column where you class labels are. They should be in numerical form.
% measure: Which quantifier in OWA operators you want to use.
% 1=Basic RIM quantifier
% 2=Quadratic quantifier
% 3=Exponential quantifier
% 4=Trigonometric quantifier
% 5=O'Hagans method.
% p: parameter in generalized Lukasiewics similarity, can be studied as a range of
% parameter values and then given as a vector i.e. p=[0.1:0.1:5].
% alpha: alpha value in OWA operators.
% Can be given as vector i.e. alpha=[0.1:0.1:5].
% pl: do you want to plot how the parameter changes in p and alpha changes the
% mean classification accuracies and variances.
%
%OUTPUTS:
% Mean_accuracy: Mean classification accuracy with best parameter values in
% p and m:
% p1 and m1: best parameter values w.r.t. mean classification accuracy.
% classes: Class information. To which class the sample was classified (in testing set data, datatest)

sv_name = 'results1'; % mat file where results are stored ('' = "no storing")

[datatest, lc1, cs1] = init_data(datatest,v,c);
[datalearn, lc2, cs2] = init_data(datalearn,v,c);

w_opt = 0; % Do we use weight optimization. Not implemented in this version.

%Initializations
N=1;
rn=1;
fitness = zeros(1,N);
fitness_id = zeros(1,N);
fitness_dif = zeros(1,N);
Means = zeros(length(p),length(m),length(rn));
Vars = zeros(length(p),length(m),length(rn));
Maxsf = zeros(length(p),length(m),length(rn));
Minsf = zeros(length(p),length(m),length(rn));

Means_fit_dif = zeros(length(p),length(m),length(rn));
Vars_fit_dir = zeros(length(p),length(m),length(rn));
Ideal_var = zeros(length(p),length(m), length(lc1),length(rn));

for n = 1 : 1
    rn_ideal = ones(1,length(lc1));
    for j = 1:length(m)
        for i = 1:length(p)
            y = [p(i), m(j),measure]; % p and m values and similarity measure

            for k = 1 : 1
                ideals(:,:,k) = idealvectors(datalearn, y); % idealvectors
                if w_opt == 0
                  [fitness(k), class, Simil] = calcfit(datatest, ideals(:,:,k), y);
                end
            end
            classinfo(i,j,1:length(class))=class;
            Means(i,j,n) = mean(fitness);
            Vars(i,j,n) = var(fitness);
            Maxsf(i,j,n) = max(fitness);
            Minsf(i,j,n) = min(fitness);
            fitness=[];
        end
    end
    tmp=max(max(Means));
    [p1,m1]=find(tmp==Means);
    Classification_accuracy=Means(p1(1),m1(1));
    classes=classinfo(p1(1),m1(1),:);
    if pl==1
        [X,Y] = meshgrid(m,p);
        figure
        surfc(X,Y,Means(:,:,n))
        title('Classification accuracies')
        xlabel('alpha-values')
        ylabel('p-values')
        zlabel('Classification accuracy')
    end
    clear Y
end
if length(sv_name) ~= 0
    save(sv_name)
end

author: Sable

repo: mcbench-benchmarks


text:
function d = empirical_dist(norm_type,x,y)

n = size(x,2);

% wasserstein integer conversion
Sc = 1e+6;

switch norm_type
    case {'wass1' 'wass2'}
        p = str2num(norm_type(5));
        C0 = distmat(x,y).^p;
        C = int32( Sc*C0 );
        [rho,varrho,u,v] = hungarianLSAP(C);
        d = sum( C0( (1:n) + (double(rho)'-1)*n ) ) / n;
        d = d^(1/p);
    case 'energy'
        k = @(a)-abs(a);
        mu = ones(n,1)/n;
        d = sqrt( rkhs_norm_invariant(k,x,mu,y,mu) );
    case 'gaussian'
        sigma = .3;
        k = @(a)exp( -(a).^2 / (2*sigma^2) );
        mu = ones(n,1)/n;
        d = sqrt( rkhs_norm_invariant(k,x,mu,y,mu) );
    otherwise
        error('Unknown');
end

end

author: optimaltransport


text:
function idxPrm = addPrm(obj, varargin)
% ADDPRM Add parameter(s) to problem

% Fields of prm structure
% - x0           -> mandatory
% - xhat         -> NaN
% - sig          -> NaN
% - const        -> can be omitted
% - idxAdj       -> NaN
% - lowerBound   -> can be omitted
% - upperBound   -> can be omitted
% - scale4report -> can be omitted
% - report       -> can be omitted
% - label        -> can be omitted

prm2add = paramvalue2struct(varargin{:}); % create struct

noPrm2add = numel(prm2add.x0); % no. of parameters to add

% Add default values if not specified
if ~isfield(prm2add, 'const'       ), prm2add.const        = false(       noPrm2add, 1); end
if ~isfield(prm2add, 'lowerBound'  ), prm2add.lowerBound   = repmat(-Inf, noPrm2add, 1); end
if ~isfield(prm2add, 'upperBound'  ), prm2add.upperBound   = repmat(+Inf, noPrm2add, 1); end
if ~isfield(prm2add, 'scale4report'), prm2add.scale4report = ones(        noPrm2add, 1); end
if ~isfield(prm2add, 'report'      ), prm2add.report       = true(        noPrm2add, 1); end
if ~isfield(prm2add, 'label'       ), prm2add.label        = repmat({''}, noPrm2add, 1); end

% Default values
prm2add.xhat   = NaN(noPrm2add,1);
prm2add.sig    = NaN(noPrm2add,1);
prm2add.idxAdj = NaN(noPrm2add,1);

% Add observations with 'direct indexing' (faster than other methods)
% Indices
if isempty(obj.prm)
    idxPrm2add = 1:noPrm2add;
else
    noPrm = numel(obj.prm.x0); % actual no. of parameters
    idxPrm2add = noPrm+1:noPrm+noPrm2add;
end

% Add
obj.prm.x0(          idxPrm2add,1) = prm2add.x0;
obj.prm.xhat(        idxPrm2add,1) = prm2add.xhat;
obj.prm.sig(         idxPrm2add,1) = prm2add.sig;
obj.prm.const(       idxPrm2add,1) = prm2add.const;
obj.prm.idxAdj(      idxPrm2add,1) = prm2add.idxAdj;
obj.prm.lowerBound(  idxPrm2add,1) = prm2add.lowerBound;
obj.prm.upperBound(  idxPrm2add,1) = prm2add.upperBound;
obj.prm.scale4report(idxPrm2add,1) = prm2add.scale4report;
obj.prm.report(      idxPrm2add,1) = prm2add.report;
obj.prm.label(       idxPrm2add,1) = prm2add.label;

% Indices of observations
idxPrm = uint32(idxPrm2add)';

% % Input parsing ----------------------------------------------------------------
%
% p = inputParser;
% p.addParamValue('x0'          , NaN  , @(x) size(x,2)==1);
% p.addParamValue('const'       , false, @(x) islogical(x) && size(x,2)==1);
% p.addParamValue('lowerBound'  , -Inf , @(x) isnumeric(x) && size(x,2)==1);
% p.addParamValue('upperBound'  , +Inf , @(x) isnumeric(x) && size(x,2)==1);
% p.addParamValue('scale4report', 1    , @(x) isnumeric(x) && size(x,2)==1);
% p.addParamValue('label'       , {''} , @iscell);
% p.parse(varargin{:});
% p = p.Results;
%
% % Add parameter(s) to problem --------------------------------------------------
%
% if numel(p.const)        == 1, p.const        = repmat(p.const       , numel(p.x0), 1); end
% if numel(p.lowerBound)   == 1, p.lowerBound   = repmat(p.lowerBound  , numel(p.x0), 1); end
% if numel(p.upperBound)   == 1, p.upperBound   = repmat(p.upperBound  , numel(p.x0), 1); end
% if numel(p.scale4report) == 1, p.scale4report = repmat(p.scale4report, numel(p.x0), 1); end
% if numel(p.label)        == 1, p.label        = repmat(p.label       , numel(p.x0), 1); end
%
% xhat   = NaN(numel(p.x0), 1);
% sig    = NaN(numel(p.x0), 1);
% idxAdj = NaN(numel(p.x0), 1);
%
% prm2add = table(                  p.x0, xhat  , sig  , p.const, idxAdj  , p.lowerBound, p.upperBound, p.scale4report, p.label, ...
%                 'VariableNames', {'x0', 'xhat', 'sig', 'const', 'idxAdj', 'lowerBound', 'upperBound', 'scale4report', 'label'});
%
% obj.prm = [obj.prm; prm2add];
%
% idxPrm = uint32([height(obj.prm)-height(prm2add)+1 : height(obj.prm)]');

end

author: pglira


text:
function flag = is_writable(name)
  % IS_WRITABLE Determine if a path is writable
  %
  % flag = is_writable(name)
  %
  % Inputs:
  %   name  path to check
  % Outputs:
  %   flag  whether writable
  if exist(name,'file')
    [~,attribs] = fileattrib(name);
    flag = attribs.UserWrite;
    return;
  end
  dirname = fileparts(name);
  if exist(dirname,'file')
    [~,attribs] = fileattrib(dirname);
    flag = attribs.UserWrite;
    return;
  end
  flag = false;
end

author: alecjacobson

repo: gptoolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 National Renewable Energy Laboratory and National
% Technology & Engineering Solutions of Sandia, LLC (NTESS).
% Under the terms of Contract DE-NA0003525 with NTESS,
% the U.S. Government retains certain rights in this software.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef constraintClass<handle
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The ``constraintClass`` creates a ``constraint`` object saved to the MATLAB
    % workspace. The ``constraintClass`` includes properties and methods used
    % to define constraints between the body motion relative to the global reference
    % frame or relative to other bodies.
    %
    %.. autoattribute:: objects.constraintClass.constraintClass
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (SetAccess = 'public', GetAccess = 'public') %input file
        hardStops                   = struct(...                    % (`structure`) Defines the constraint hardstop
          'upperLimitSpecify',          'off',...                   % (`string`) Initialize upper stroke limit. ``  'on' or 'off' `` Default = ``off``.
          'upperLimitBound',            1, ...                      % (`float`) Define upper stroke limit in m or deg. Only active if `lowerLimitSpecify` is `on` `` Default = ``1``.
          'upperLimitStiffness',        1e6, ...                    % (`float`) Define upper limit spring stiffness, N/m or N-m/deg. `` Default = ``1e6``.
          'upperLimitDamping',          1e3, ...                    % (`float`) Define upper limit damping, N/m/s or N-m/deg/s.  `` Default = ``1e3``.
          'upperLimitTransitionRegionWidth', 1e-4, ...              % (`float`) Define upper limit transition region, over which spring and damping values ramp up to full values. Increase for stability. m or deg. ``Default = 1e-4``
          'lowerLimitSpecify',          'off',...                   % Initialize lower stroke limit. ``  `on` or `off` `` Default = ``off``.
          'lowerLimitBound',            -1, ...                     % (`float`) Define lower stroke limit in m or deg. Only active if `lowerLimitSpecify` is `on` ``   `` Default = ``-1``.
          'lowerLimitStiffness',        1e6, ...                    % (`float`) Define lower limit spring stiffness, N/m or N-m/deg.  `` Default = ``1e6``.
          'lowerLimitDamping',          1e3, ...                    % (`float`) Define lower limit damping, N/m/s or N-m/deg/s.  `` Default = ``1e3``.
          'lowerLimitTransitionRegionWidth', 1e-4)                  % (`float`) Define lower limit transition region, over which spring and damping values ramp up to full values. Increase for stability. m or deg. ``Default = 1e-4``
        initial                     = struct(...                    %
            'displacement',             [0 0 0])                    % (`structure`) Defines the initial displacement of the constraint. ``displacement`` (`3x1 float vector`) is defined as the initial displacement of the constraint [m] in the following format [x y z], Default = [``0 0 0``].
        name                        = 'NOT DEFINED'                 % (`string`) Specifies the constraint name. For constraints this is defined by the user, Default = ``NOT DEFINED``.
        location                    = [999 999 999]                 % (`3x1 float vector`) Constraint location [m]. Defined in the following format [x y z]. Default = ``[999 999 999]``.
        orientation                 = struct(...                    % (`structure`) Defines the orientation axis of the constraint.
            'z',                        [0, 0, 1], ...              %
            'y',                        [0, 1, 0], ...              %
            'x',                        [], ...                     %
            'rotationMatrix',           [])                         % (`structure`) Defines the orientation axis of the constraint. ``z`` (`3x1 float vector`) defines the direciton of the Z-coordinate of the constraint, Default = [``0 0 1``]. ``y`` (`3x1 float vector`) defines the direciton of the Y-coordinate of the constraint, Default = [``0 1 0``]. ``x`` (`3x1 float vector`) internally calculated vector defining the direction of the X-coordinate for the constraint, Default = ``[]``. ``rotationMatrix`` (`3x3 float matrix`) internally calculated rotation matrix to go from standard coordinate orientation to the constraint coordinate orientation, Default = ``[]``.
    end

    properties (SetAccess = 'private', GetAccess = 'public') %internal
        number                      = []                            % Constraint number
    end

    methods (Access = 'public')
        function obj = constraintClass(name)
            % This method initilizes the ``constraintClass`` and creates a
            % ``constraint`` object.
            %
            % Parameters
            % ------------
            %     filename : string
            %         String specifying the name of the constraint
            %
            % Returns
            % ------------
            %     constraint : obj
            %         contraintClass object
            %
            if exist('name','var')
                obj.name = name;
            else
                error('The constraint class number(s) in the wecSimInputFile must be specified in ascending order starting from 1. The constraintClass() function should be called first to initialize each constraint with a name.')
            end
        end

        function obj = checkLoc(obj,action)
            % This method checks WEC-Sim user inputs and generate an error message if the constraint location is not defined in constraintClass.

            % Checks if location is set and outputs a warning or error. Used in mask Initialization.
            switch action
              case 'W'
                if obj.location == 999 % Because "Allow library block to modify its content" is selected in block's mask initialization, this command runs twice, but warnings cannot be displayed during the first initialization.
                    obj.location = [888 888 888];
                elseif obj.location == 888
                    obj.location = [0 0 0];
                    s1= ['For ' obj.name ': constraint.location was changed from [9999 9999 9999] to [0 0 0]'];
                    warning(s1)
                end
              case 'E'
                try
                    if obj.location == 999
                      s1 = ['For ' obj.name ': constraint.location needs to be specified in the WEC-Sim input file.' ...
                        ' constraint.location is the [x y z] location, in meters, for the pitch constraint.'];
                      error(s1)
                    end
                catch exception
                  throwAsCaller(exception)
                end
            end
        end

        function obj = setOrientation(obj)
            % This method calculates the constraint ``x`` vector and ``rotationMatrix`` matrix in the ``orientation`` structure based on user input.
            obj.orientation.z = obj.orientation.z / norm(obj.orientation.z);
            obj.orientation.y = obj.orientation.y / norm(obj.orientation.y);
            z = obj.orientation.z;
            y = obj.orientation.y;
            if abs(dot(y,z))>0.001
                error('The Y and Z vectors defining the constraint''s orientation must be orthogonal.')
            end
            x = cross(y,z)/norm(cross(y,z));
            x = x(:)';
            obj.orientation.x = x;
            obj.orientation.rotationMatrix  = [x',y',z'];
        end

        function setInitDisp(obj, relCoord, axisAngleList, addLinDisp)
            % Function to set a constraints's initial displacement
            %
            % This function assumes that all rotations are about the same relative coordinate.
            % If not, the user should input a relative coordinate of 0,0,0 and
            % use the additional linear displacement parameter to set the cg or location
            % correctly.
            %
            % Parameters
            % ------------
            %    relCoord : [1 3] float vector
            %        Distance from x_rot to the body center of gravity or the constraint
            %        or pto location as defined by: relCoord = cg - x_rot. [m]
            %
            %    axisAngleList : [nAngle 4] float vector
            %        List of axes and angles of the rotations with the
            %        format: [n_x n_y n_z angle] (angle in rad)
            %        Rotations applied consecutively in order of dimension 1
            %
            %    addLinDisp : [1 3] float vector
            %        Initial linear displacement (in addition to the
            %        displacement caused by rotation) [m]
            %

            % initialize quantities before for loop
            axisList = axisAngleList(:,1:3);
            angleList = axisAngleList(:,4);
            nAngle = size(axisList,1);
            rotMat = eye(3);

            % Loop through all axes and angles.
            for i=1:nAngle
                rotMat = axisAngle2RotMat(axisList(i,:),angleList(i))*rotMat;
            end

            % calculate net axis-angle rotation
%             [netAxis, netAngle] = rotMat2AxisAngle(rotMat);

            % calculate net displacement due to rotation
            rotatedRelCoord = relCoord*(rotMat');
            linDisp = rotatedRelCoord - relCoord;

            % apply rotation and displacement to object
            obj.initial.displacement = linDisp + addLinDisp;

        end

        function listInfo(obj)
            % This method prints constraint information to the MATLAB Command Window.
            fprintf('\n\t***** Constraint Name: %s *****\n',obj.name)
        end

        function setNumber(obj,number)
            % Method to set the private number property
            obj.number = number;
        end
    end
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
port = int16(2000);
client = py.carla.Client('localhost', port);
client.set_timeout(10.0);
world = client.get_world();

% Spawn Vehicle
blueprint_library = world.get_blueprint_library();
car_list = py.list(blueprint_library.filter("model3"));
car_bp = car_list{1};
spawn_point = py.random.choice(world.get_map().get_spawn_points());
tesla = world.spawn_actor(car_bp, spawn_point);
tesla.set_autopilot(true);

carla_is_running = true;

% Sensor 1
blueprint = world.get_blueprint_library().find('sensor.camera.rgb');
blueprint.set_attribute('image_size_x', '960')
blueprint.set_attribute('image_size_y', '540')
% blueprint.set_attribute('sensor_tick', '.1');
transform = py.carla.Transform(py.carla.Location(pyargs('x',-7.5, 'z',2.5)));
sensor = world.spawn_actor(blueprint, transform, pyargs('attach_to',tesla));

pyModule = sensorBind(sensor, "rgb", "rgb", "array");
currentImage = uint8(py.getattr(pyModule, 'array'));
imageHandle = imshow(currentImage);
set(gca,'units','pixels'); % set the axes units to pixels
x = get(gca,'position'); % get the position of the axes
set(gcf,'units','pixels'); % set the figure units to pixels
y = get(gcf,'position'); % get the figure position
set(gcf,'position',[y(1) y(2) x(3) x(4)]);% set the position of the figure to the length and width of the axes
set(gca,'units','normalized','position',[0 0 1 1]); % set the axes units to pixels
set(gcf,'menubar','none');
% Lidar
blueprint = world.get_blueprint_library().find('sensor.lidar.ray_cast');
blueprint.set_attribute('points_per_second', '56000');
blueprint.set_attribute('range', '5000');
blueprint.set_attribute('sensor_tick', '0.1');
transform = py.carla.Transform(py.carla.Location(pyargs('x',0.8, 'z',1.7)));
lidar = world.spawn_actor(blueprint, transform, pyargs('attach_to',tesla));

moduleLidar = sensorBind(lidar, 'lidar_file', 'lidar', 'array');

% Spawn NPC's
npc_bps = blueprint_library.filter("vehicle");
npc_to_spawn = 75;

% Preallocate memory
npc_list = cell(1, npc_to_spawn);

i = 1;
while i <= npc_to_spawn
    try
        npc_bp = py.random.choice(npc_bps);
        spawn_point = py.random.choice(world.get_map().get_spawn_points());
        npc_list{i} = world.spawn_actor(npc_bp, spawn_point);
        npc_list{i}.set_autopilot(true);
    catch
        % In case spawing fails due to collision, try again
        i = i - 1;
    end
    i = i + 1;
end

% A bounding box detector model.
detectorModel = HelperBoundingBoxDetector(...
    'XLimits',[-50 75],...              % min-max
    'YLimits',[-5 5],...                % min-max
    'ZLimits',[-2 5],...                % min-max
    'SegmentationMinDistance',1.6,...   % minimum Euclidian distance
    'MinDetectionsPerCluster',1,...     % minimum points per cluster
    'MeasurementNoise',eye(6),...       % measurement noise in detection report
    'GroundMaxDistance',0.3);           % maximum distance of ground points from ground plane

assignmentGate = [10 100]; % Assignment threshold;
confThreshold = [7 10];    % Confirmation threshold for history logic
delThreshold = [8 10];     % Deletion threshold for history logic
Kc = 1e-5;                 % False-alarm rate per unit volume

% IMM filter initialization function
filterInitFcn = @helperInitIMMFilter;

% A joint probabilistic data association tracker with IMM filter
tracker = trackerJPDA('FilterInitializationFcn',filterInitFcn,...
    'TrackLogic','History',...
    'AssignmentThreshold',assignmentGate,...
    'ClutterDensity',Kc,...
    'ConfirmationThreshold',confThreshold,...
    'DeletionThreshold',delThreshold,...
    'HasDetectableTrackIDsInput',true,...
    'InitializationThreshold',0);

% Create display
displayObject = HelperLidarExampleDisplay(uint8(py.getattr(pyModule, 'array')),...
    'PositionIndex',[1 3 6],...
    'VelocityIndex',[2 4 7],...
    'DimensionIndex',[9 10 11],...
    'YawIndex',8,...
    'MovieName','',...  % Specify a movie name to record a movie.
    'RecordGIF',false); % Specify true to record new GIFs

%% Loop Through Data
% Loop through the recorded lidar data, generate detections from the
% current point cloud using the detector model and then process the
% detections using the tracker.
start_time = cputime;

% Initiate all tracks.
allTracks = struct([]);

% Rotate point cloud
A = [0 -1 0 0; ...
     1  0 0 0; ...
     0  0 1 0; ...
     0  0 0 1];
tform = affine3d(A);

while carla_is_running
    try
        % Update time
        time = cputime - start_time;

        % Get current lidar scan
        XYZI = single(py.getattr(moduleLidar, 'array'));
        XYZ = lidarData(:, 1:3);

        % Flip the axis
        XYZ(:, 2) = -1 * XYZ(:, 2);
        currentLidar = pctransform(pointCloud(XYZ), tform);

        % Generator detections from lidar scan.
        [detections,obstacleIndices,groundIndices,croppedIndices] = detectorModel(currentLidar,time);

        % Calculate detectability of each track.
        detectableTracksInput = helperCalcDetectability(allTracks,[1 3 6]);

        % Pass detections to track.
        [confirmedTracks,tentativeTracks,allTracks] = tracker(detections,time,detectableTracksInput);

        % Get model probabilities from IMM filter of each track using
        % getTrackFilterProperties function of the tracker.
        modelProbs = zeros(2,numel(confirmedTracks));
        for k = 1:numel(confirmedTracks)
            c1 = getTrackFilterProperties(tracker,confirmedTracks(k).TrackID,'ModelProbabilities');
            modelProbs(:,k) = c1{1};
        end

        % Update display
        if isvalid(displayObject.PointCloudProcessingDisplay.ObstaclePlotter)
            % Get current image scan for reference image
            currentImage = uint8(py.getattr(pyModule, 'array'));

            set(imageHandle,'Cdata',currentImage);

            % Update display object
            displayObject(detections,confirmedTracks,currentLidar,obstacleIndices,...
                groundIndices,croppedIndices,[],modelProbs);
        end
    catch
        carla_is_running = false;
        close all;
    end
end

tesla.destroy();
sensor.destroy();
lidar.destroy();

for i=1:npc_to_spawn
   npc_list{i}.destroy();
end

author: darkscyla

repo: MATLAB-Carla-Interface