text:
%% Example of using KST class for interfacing with KUKA iiwa robots

% Collaborative pick and place application, the robot detects collision
% with its surrounding while performing the pick and place application,
% when a collision is detected, the robots stops its motion, when the
% operator gives the robot a double touch on the end-effector in the Z
% direcdtion, the robot returns back to doing its operation.

% First start the server on the KUKA iiwa controller
% Then run this script using Matlab

% This example works with Sunrise application version KST_1.7  and higher.
% Copyright Mohammad SAFEEA, 26th-June-2018

close all;clear;clc;
warning('off')
%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
    return;
end
pause(1);
disp('An example on using the interruptible PTP motion functions for Human robot collaboration')

%% Move PTP in joint space to some inital position
jPos={0,0,0,-pi/2,0,pi/2,0};
relVel=0.15;
iiwa.movePTPJointSpace(jPos, relVel); % move to initial configuration
%% Move PTP in joint space with torque limits
[ Pos ] = iiwa.getEEFPos();
deltaY=250;
deltaZ=332-150;

joints_indices=[1,2,3,4,5,6,7]; % joint index start from one
max_torque=ones(7,1)*8;
max_torque(1)=3;
min_torque=-max_torque;
VEL=175;
w=17; % weight of tool newton
%% Go down near the table
Pos{3}=Pos{3}-150;
iiwa.moveWaitForDTWhenInterrupted(  ...
     Pos, VEL,joints_indices,max_torque,min_torque,w);
%% Open gripper
iiwa.setPin1Off();
pause(1);
iiwa.setPin11On();
pause(1);
%% Motion twards the box
Pos{2}=Pos{2}-deltaY;
iiwa.moveWaitForDTWhenInterrupted( ...
     Pos, VEL,joints_indices,max_torque,min_torque,w);

Pos{3}=Pos{3}-deltaZ;
iiwa.moveWaitForDTWhenInterrupted(  ...
     Pos, VEL,joints_indices,max_torque,min_torque,w);
%% Close gripper
iiwa.setPin1On();
pause(1);
iiwa.setPin11Off();
pause(1);
%% Picking the box twards the goal
Pos{3}=Pos{3}+deltaZ;
iiwa.moveWaitForDTWhenInterrupted(  ...
     Pos, VEL,joints_indices,max_torque,min_torque,w);

Pos{2}=Pos{2}+2*deltaY;
iiwa.moveWaitForDTWhenInterrupted(  ...
     Pos, VEL,joints_indices,max_torque,min_torque,w);

Pos{3}=Pos{3}-deltaZ;
iiwa.moveWaitForDTWhenInterrupted(  ...
     Pos, VEL,joints_indices,max_torque,min_torque,w);
%% Release the box
% Open gripper
iiwa.setPin1Off();
pause(1);
iiwa.setPin11On();
pause(1);

Pos{3}=Pos{3}+deltaZ;
iiwa.moveWaitForDTWhenInterrupted(  ...
     Pos, VEL,joints_indices,max_torque,min_torque,w);

%% Turn off the server
iiwa.net_turnOffServer();

repo: KST-Kuka-Sunrise-Toolbox


text:
% JM: 2006/07/19
% viewer for data Tc on domain Omega

function varargout = viewImage(varargin);
% this is just output and hence not commented; fell free to improve (-;
% function varargout = viewImage(Tc,Omega,m,varargin);

persistent PARA

if nargin == 0,
  varargout{1} = PARA;
  return;
end;

% set persistent parameter
if nargin >= 1 & strcmp(varargin{1},'set'),
  PARA = setPARA(PARA,varargin{2:end});
  varargout{1} = PARA;
  return;
end;

if isstr(varargin{1}) & strcmp(varargin{1},'clear'),
  PARA = [];
  return;
end;

if isempty(PARA),
  PARA = setPARA(PARA,'viewer','image2D','colmap',gray(256));
  varargout{1} = PARA;
  warning('!!! no PARA given yet, set viewer to image2D !!!')
end;

% -----------------------------------------------------------------------------
% start to work

Tc    = varargin{1};
Omega = varargin{2};
m     = varargin{3};
varargin(1:3) = [];

% start with global parameters and update
para = PARA;
for j=1:length(varargin)/2,
  field = varargin{2*j-1};
  value = varargin{2*j};
  para  = setfield(para,field,value);
end;

if ~para.plots, return;  end;          % do nothing, just return;

% build parameter list
fields = fieldnames(para);
for j=1:length(fields);
  optn{2*j-1} = fields{j};
  optn{2*j}   = getfield(para,fields{j});
end;

% scale image to [0,255]
if isfield(para,'scale')  & para.scale,
  maxTc = max(Tc);
  Tc = 255/(maxTc+(maxTc==0))*Tc;
end;

% invert image
if isfield(para,'invert') & para.invert,
  Tc = max(Tc)-Tc;
end;

% setup figure
if ~isfield(para,'fig'), para.fig = [];   end;
if isempty(para.fig),    para.fig = gcf;  end;

fig = figure(para.fig);
if any(para.sub ~= 1),
  subplot(para.sub(1),para.sub(2),para.sub(3)); cla;
end;

if ~isempty(para.figname),
  set(fig,'numbertitle','off','name',sprintf('[JM-%d]: %s',fig,para.figname));
end;

% call the viewer
fh = feval(para.viewer,Tc,Omega,m,optn{:});
str = sprintf('%s',para.name);
title(str);

if nargout>0,
  varargout{1} = fh;
end;

return;
%==============================================================================
function PARA = setPARA(PARA,varargin);

% set default parameter
p.viewer   = [];
p.figname  = [];
p.fig      = [];
p.sub      = [1,1,1];
p.name     = '';
p.plots    = 1;
p.scale    = 0;
p.invert   = 0;

% get default values
fn = fieldnames(p);
for j=1:length(fn),
  default{2*j-1} = fn{j};
  default{2*j}   = getfield(p,fn{j});
end;

% get old values
if isempty(PARA),
  current = {};
else
  fn = fieldnames(PARA);
  for j=1:length(fn),
    current{2*j-1} = fn{j};
    current{2*j}   = getfield(PARA,fn{j});
  end;
end;

% get new values
if strcmp(varargin{1},'para') | strcmp(varargin{1},'PARA'),
  PARA = varargin{2};
  fn = fieldnames(PARA);
  for j=1:length(fn),
    new{2*j-1} = fn{j};
    new{2*j}   = getfield(PARA,fn{j});
  end;
else
  new = varargin;
end;

optn = {default{:},current{:},new{:}};

% rewrite default parameter
PARA = [];
for j=1:length(optn)/2,
  field = optn{2*j-1};
  value = optn{2*j};
  PARA  = setfield(PARA,field,value);
end;

fprintf('%s: set view parameter, viewer=%s\n',mfilename,PARA.viewer);
return;
% -----------------------------------------------------------------------------

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function Control = setupController(P)

if P.Sim.springDoublePendulum

    Control.legOne.nomPos = 1; %(m)
    Control.legOne.Kp = 35;

    Control.legTwo.nomPos = 1; %(m)
    Control.legTwo.Kp = 35;

    %Then set everything to zero except for the proportional gains on the
    %leg length actuators. This turns them both into springs, which makes
    %the system conserve energy. This is a handle error check.

    %The Foot One actuator operates on the absolute angle of leg one
    Control.footOne.nomPos = 0;
    Control.footOne.nomVel = 0;
    Control.footOne.Kp = 0;
    Control.footOne.Kd = 0;

    %The Foot Two actuator operates on the absolute angle of leg two
    Control.footTwo.nomPos = 0;
    Control.footTwo.nomVel = 0;
    Control.footTwo.Kp = 0;
    Control.footTwo.Kd = 0;

    %The Leg One actuator operates on the length of leg one
    Control.legOne.nomVel = 0;
     Control.legOne.Kd = 0;

    %The Leg Two actuator operates on the length of leg two
    Control.legTwo.nomVel = 0;
    Control.legTwo.Kd = 0;

    %The Hip actuator operates on the angle between the legs: th2-th1
    Control.hip.nomPos = 0;
    Control.hip.nomVel = 0;
    Control.hip.Kp = 0;
    Control.hip.Kd = 0;
else
    %Allow system to be fully controlled. This will not necessarily
    %conserve energy.

    %The Foot One actuator operates on the absolute angle of leg one
    Control.footOne.nomPos = 0;
    Control.footOne.nomVel = 0;
    Control.footOne.Kp = 0;
    Control.footOne.Kd = 0;

    %The Foot Two actuator operates on the absolute angle of leg two
    Control.footTwo.nomPos = 0;
    Control.footTwo.nomVel = 0;
    Control.footTwo.Kp = 0;
    Control.footTwo.Kd = 0;

    %The Leg One actuator operates on the length of leg one
    Control.legOne.nomPos = 1; %(m)
    Control.legOne.nomVel = 0;
    Control.legOne.Kp = 100;
    Control.legOne.Kd = 50;

    %The Leg Two actuator operates on the length of leg two
    Control.legTwo.nomPos = 1; %(m)
    Control.legTwo.nomVel = 0;
    Control.legTwo.Kp = 100;
    Control.legTwo.Kd = 500;

    %The Hip actuator operates on the angle between the legs: th2-th1
    Control.hip.nomPos = 0;
    Control.hip.nomVel = 0;
    Control.hip.Kp = 0;
    Control.hip.Kd = 0;
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function x = r8pbu_sl ( n, mu, a_lu, b )

%*****************************************************************************80
%
%% R8PBU_SL solves a R8PBU system factored by R8PBU_FA.
%
%  Discussion:
%
%    The R8PBU storage format is for a symmetric positive definite band matrix.
%
%    To save storage, only the diagonal and upper triangle of A is stored,
%    in a compact diagonal format that preserves columns.
%
%    The diagonal is stored in row MU+1 of the array.
%    The first superdiagonal in row MU, columns 2 through N.
%    The second superdiagonal in row MU-1, columns 3 through N.
%    The MU-th superdiagonal in row 1, columns MU+1 through N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, Philadelphia, 1979.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer MU, the number of superdiagonals of the matrix.
%    MU must be at least 0 and no more than N-1.
%
%    Input, real A_LU(MU+1,N), the matrix, as factored by R8PBU_FA.
%
%    Input, real B(N), the right hand side of the linear system.
%
%    Output, real X(N), the solution vector.
%
  x(1:n) = b(1:n);
%
%  Solve L * Y = B.
%
  for k = 1 : n
    ilo = max ( 1, k - mu );
    x(k) = ( x(k) - x(ilo:k-1) * a_lu(mu+1+ilo-k:mu,k)  ) ...
      / a_lu(mu+1,k);
  end
%
%  Solve U * X = Y.
%
  for k = n : -1 : 1

    x(k) = x(k) / a_lu(mu+1,k);

    ilo = max ( 1, k - mu );
    for i = ilo : k-1
      x(i) = x(i) - x(k) * a_lu(mu+1+i-k,k);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function bool=sde_issquare(V)
%SDE_ISSQUARE  True for square matrices.
%
%   See also: SDE_ISMATRIX, SDE_ISDIAG

%   Andrew D. Horchler, horchler @ gmail . com, Created 4-30-13
%   Revision: 1.2, 5-4-13

bool = (ndims(V) == 2 && size(V,1) == size(V,2));	%#ok<ISMAT>

author: horchler

repo: SDETools

save_path: github-repos/MATLAB/horchler-SDETools


text:
%MAPM Mapping to execute arbitray Matlab function on dataset
%
%    [B,OUT] = MAPM(A,COMMAND,PAR1,PAR2,....)
%    [B,OUT] = A*MAPM([],COMMAND,PAR1,PAR2,....)
%    [B,OUT] = A*MAPM(COMMAND,PAR1,PAR2,....)
%
% INPUT
%    A               Dataset or datafile or double
%    COMMAND         String with function name
%    {PAR1,...}      Cell array with optional parameters for COMMAND
%
% OUTPUT
%    B               Resulting dataset, datafile or double array
%    OUT             (only for datasets or doubles) possible additional
%                    outputs
%
% DESCRIPTION
% On the data in A, the double array A the following command is executed
%
%    [B,OUT] = COMMAND(DATA,PAR1,PAR2,....)
%
% If A is a dataset then B is converted to a similar dataset. If A is a
% cell array, then results are combined into a cell array as well.
%
% This command differs from the similar command FILTM which operates object
% by object (row by row). The COMMAND in MAPM operates on the entire data
% matrix.
%
% EXAMPLES
% argmin = mapm('min')*out2;
% % A*argmin returns for every column the row indices of the minimum
% % A'*argmin returns for every row the column indices of its minimum
%
% nexpm = mapm('uminus')*mapm('exp');
% % A*nexpm returns exp(-A), useful incombination with other commands
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, IM2OBJ, DATA2IM, FILTM, FILTIM, OUT2
% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function varargout = mapm(varargin)

  if nargin == 0
    error('No command found')
  end
  varargin = shiftargin(varargin,'char');
  if isempty(varargin{1})
    varargout{1} = prmapping(mfilename,'fixed_cell',varargin(2:end));
    varargout{1} = setname(varargout{1},'AnyMap');
    return
  end

  % now we have: mapm(data,command,pars)
  % data might be dataset, datafile, cell array, double, structure array
  pars = cell(1,nargin-2);     % input pars
  [a,command,pars{:}] = deal(varargin{:}); % use nice names
  if numel(pars) > 1 % allow argument list as cell array or not.
   pars = {pars}
  end
  varargout = cell(1,nargout);

  if isdatafile(a)               % for datafiles processing is stored

    if isempty(getpostproc(a)) & ~ismapping(command)
                                 % as preprocessing (if no postproc defined)
      b = addpreproc(a,command,pars);
    else                         % or as mapping as postprocessing
      v = prmapping(mfilename,'fixed',{command,pars});
      v = setname(v,mapname);
      b = addpostproc(a,v);      % store mapping
    end
    varargout{1} = b;
    return

  elseif isdataset(a) % are executed here

    try
      [varargout{:}] = my_feval(command,a,pars);
    catch me
      [varargout{:}] = my_feval(command,+a,pars);
      if nargout > 0
        if size(varargout{1},1) == size(a,1)
          varargout{1} = setdata(a,varargout{1});
        else
          varargout{1} = prdataset(varargout{1});
          varargout{1} = setuser(varargout{1},getuser(a));
        end
      end
    end

  elseif iscell(a)

    b = cell(size(a));
    varargout = cell(numel(a),numel(varargout));
    n = numel(b);
    s = sprintf('Processing %i objects: ',n);
    prwaitbar(n,s);
    for i=1:n
      prwaitbar(n,i,[s int2str(i)]);
      [varargout{i,:}] = my_feval(command,a{i},pars);
    end
    prwaitbar(0);

  else % doubles

    [varargout{:}] = my_feval(command,a,pars{:});

  end

return

function varargout = my_feval(command,a,pars)
  if nargin < 3, pars = []; end
  varargout = cell(1,nargout);
  if isempty(pars)
    [varargout{:}] = feval(command,a);
  elseif iscell(pars)
    [varargout{:}] = feval(command,a,pars{:});
  else
    [varargout{:}] = feval(command,a,pars);
  end
return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function defaultOptions = prtOptionsGetFactory()
% prtOptionsGetFactory  Get the factory options for the PRT
%   PRT factory options are hard coded in MATLAB class m files stored
%   within the prtOptions package. prtOptionsGetFactory returns a structure
%   with each of these classes loaded into a structure with field names
%   equal to the class names
%
%   defaultOptions = prtOptionsGetFactory()
%
% See also. prtOptionsGet, prtOptionsSet, prtOptionsSetDefault

optionsPackageInfo = meta.package.fromName('prtOptions');

if isempty(optionsPackageInfo)
    error('Package prtOptions was not found on path');
end

for iClass = 1:length(optionsPackageInfo.Classes)
    defaultOptions.(strrep(optionsPackageInfo.Classes{iClass}.Name,'prtOptions.','')) = feval(optionsPackageInfo.Classes{iClass}.Name);
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function value = angle_deg_2d ( p1, p2, p3 )

%*****************************************************************************80
%
%% ANGLE_DEG_2D returns the angle swept out between two rays in 2D.
%
%  Discussion:
%
%    Except for the zero angle case, it should be true that
%
%      ANGLE_DEG_2D(P1,P2,P3) + ANGLE_DEG_2D(P3,P2,P1) = 360.0
%
%        P1
%        /
%       /
%      /
%     /
%    P2--------->P3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), P3(2,1), define the rays
%    P1 - P2 and P3 - P2 which in turn define the angle.
%
%    Output, real VALUE, the angle swept out by the rays, measured
%    in degrees.  0 <= VALUE < 360.  If either ray has zero length,
%    then VALUE is set to 0.
%
  p(1,1) = ( p3(1,1) - p2(1,1) ) * ( p1(1,1) - p2(1,1) ) ...
         + ( p3(2,1) - p2(2,1) ) * ( p1(2,1) - p2(2,1) );

  p(2,1) = ( p3(1,1) - p2(1,1) ) * ( p1(2,1) - p2(2,1) ) ...
         - ( p3(2,1) - p2(2,1) ) * ( p1(1,1) - p2(1,1) );

  if ( p(1,1) == 0.0 & p(2,1) == 0.0 )
    value = 0.0;
    return
  end

  value = atan2 ( p(2,1), p(1,1) );

  if ( value < 0.0 )
    value = value + 2.0 * pi;
  end

  value = radians_to_degrees ( value );

  return
end


text:
function obj = viewMovieFiltering(obj)
	% Allows testing of spatial filtering.
	% Biafra Ahanonu
	% branched from controllerAnalysis: 2014.08.01 [16:09:16]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	[fileIdxArray idNumIdxArray nFilesToAnalyze nFiles] = obj.getAnalysisSubsetsToAnalyze();

	movieSettings = inputdlg({...
			'start:end frames (leave blank for all)',...
			'FFT filter type (bandpass, lowpass, highpass)',...
			'FFT filter mask type (gaussian or binary)',...
			'movie regular expression',...
			'HDF5 dataset name',...
			'Show division and dfof of movie'...
		},...
		'Settings for movie spatial filtering check',[1 100],...
		{...
			'1:25',...
			'bandpass',...
			'gaussian',...
			'concat',...
			'/1',...
			'1'...
		}...
	);
	frameList = str2num(movieSettings{1});
	% options.normalizeFreqLow = str2num(movieSettings{1});
	% options.normalizeFreqHigh = str2num(movieSettings{1});
	options.normalizeBandpassType = movieSettings{2};
	options.normalizeBandpassMask = movieSettings{3};
	fileFilterRegexp = movieSettings{4};
	inputDatasetName = movieSettings{5};
	testDuplicateDfof = str2num(movieSettings{6});

	% loop over all directories and copy over a specified chunk of the movie file
	for thisFileNumIdx = 1:nFilesToAnalyze
			thisFileNum = fileIdxArray(thisFileNumIdx);
			obj.fileNum = thisFileNum;
			display(repmat('=',1,21))
			display([num2str(thisFileNumIdx) '/' num2str(nFilesToAnalyze) ' (' num2str(thisFileNum) '/' num2str(nFiles) '): ' obj.fileIDNameArray{obj.fileNum} 10 obj.inputFolders{obj.fileNum}]);
			movieList = getFileList(obj.inputFolders{obj.fileNum}, fileFilterRegexp);
			movieList = movieList{1};

			% get list of frames and correct for user input
			if isempty(frameList)
				frameListTmp = frameList;
			else
				movieDims = loadMovieList(movieList,'convertToDouble',0,'frameList',[],'inputDatasetName',inputDatasetName,'getMovieDims',1);
				nMovieFrames = sum(movieDims.z);
				display(['movie frames: ' num2str(nMovieFrames)]);
				if nMovieFrames<nanmax(frameList)
					frameListTmp = frameList;
					frameListTmp(frameListTmp>nMovieFrames) = [];
				else
					frameListTmp = frameList;
				end
			end

			% get movie, normalize, and display
			[primaryMovie] = loadMovieList(movieList,'convertToDouble',0,'frameList',frameListTmp(:),'inputDatasetName',inputDatasetName);
			primaryMovie = single(primaryMovie);

			stopTesting = 1;
			while ~isempty(stopTesting)
				[inputMovieNormalized] = normalizeMovie(primaryMovie,'normalizationType','matlabFFT_test','bandpassType',options.normalizeBandpassType,'bandpassMask',options.normalizeBandpassMask,'testDuplicateDfof',testDuplicateDfof);

				playMovie(inputMovieNormalized,'colormapColor','gray');

				stopTesting = strmatch(questdlg('Repeat?', ...
					'Run filtering for this movie again?', ...
					'yes','no','yes'),'yes');
			end
	end

end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
% Uncomment individual Example lines.

%---- Example 1: E5 (nyquist frequency)
fs = 120*1.023e6;
n_E5_periods = 1;

%---- Example 2: E5 (nyquist frequency, several periods)
% fs = 120*1.023e6;
% n_E5_periods = 8.2;

%---- Example 3: E5 (sub-nyquist frequency, several periods)
% fs = 10*1.023e6;
% n_E5_periods = 2.2;

%---- Example 4: E5 (higher frequency, several periods)
% fs = 130*1.023e6;
% n_E5_periods = 1.2;

%---- Example 4: E5 (higher frequency non multiple, several periods)
% fs = 130*1e6;
% n_E5_periods = 1.2;

%---- Example 4: E5 (lower frequency non multiple, several periods)
% fs = 10*1e6;
% n_E5_periods = 1.2;

coh_time = n_E5_periods*1e-3;
fh = BOCgen;
[sd1, sd2, sp1 sp2]  = fh.BOC_E5(coh_time,fs);

%==========================================================================
%% Plots
%==========================================================================
t_boc_chips = linspace(0,1.023e6*120*coh_time,length(sd1));

figure;
        hold on;
        stairs(t_boc_chips,sd1,'b');
        stairs(t_boc_chips,sd2,'r');
        grid on;
        title('Data spread codes');
        xlabel('BOC symbol period')
        ylabel('amplitude');
        legend('Data spread code', 'Data spread code delayed');

figure; hold on;
        stairs(t_boc_chips,sp1,'b');
        stairs(t_boc_chips,sp2,'r');
        grid on;
        title('Pilot spread codes');
        xlabel('BOC symbol period')
        ylabel('amplitude');
        legend('Pilot spread code', 'Pilot spread code delayed');

% Next plots just shows that consecutive BOC symbols are identical
sub_period = round(fs/1.023e6*120);
n_periods = floor(length(sd1)/sub_period);
figure, hold on;
for i=1:n_periods-1
     plot(sd2(sub_period*(i-1)+1:sub_period*i)-sd2(sub_period*i+1:sub_period*(i+1)));
end

figure, hold on;
for i=1:n_periods-1
     plot(sd1(sub_period*(i-1)+1:sub_period*i)-sd1(sub_period*i+1:sub_period*(i+1)));
end

figure, hold on;
for i=1:n_periods-1
     plot(sp1(sub_period*(i-1)+1:sub_period*i)-sp1(sub_period*i+1:sub_period*(i+1)));
end

figure, hold on;
for i=1:n_periods-1
     plot(sp2(sub_period*(i-1)+1:sub_period*i)-sp2(sub_period*i+1:sub_period*(i+1)));
end

author: danipascual

repo: GNSS-matlab


text:
function [g, c, tmp, m]=kMeansCluster(m,k,isRand)
% kMeansCluster - Simple k means clustering algorithm
% Author: Kardi Teknomo, Ph.D.
%
% Purpose: classify the objects in data matrix based on the attributes
% Criteria: minimize Euclidean distance between centroids and object points
% For more explanation of the algorithm, see http://people.revoledu.com/kardi/tutorial/kMean/index.html
% Output: matrix data plus an additional column represent the group of each object
%
% Example: m = [ 1 1; 2 1; 4 3; 5 4]  or in a nice form
%          m = [ 1 1;
%                2 1;
%                4 3;
%                5 4]
%          k = 2
% kMeansCluster(m,k) produces m = [ 1 1 1;
%                                   2 1 1;
%                                   4 3 2;
%                                   5 4 2]
% Input:
%   m      - required, matrix data: objects in rows and attributes in columns
%   k      - optional, number of groups (default = 1)
%   isRand - optional, if using random initialization isRand=1, otherwise input any number (default)
%            it will assign the first k data as initial centroids
%
% Local Variables
%   f      - row number of data that belong to group i
%   c      - centroid coordinate size (1:k, 1:maxCol)
%   g      - current iteration group matrix size (1:maxRow)
%   i      - scalar iterator
%   maxCol - scalar number of rows in the data matrix m = number of attributes
%   maxRow - scalar number of columns in the data matrix m = number of objects
%   temp   - previous iteration group matrix size (1:maxRow)
%   z      - minimum value (not needed)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3,        isRand=0;   end
if nargin<2,        k=1;        end
tmp = [];

[maxRow, maxCol]=size(m);
if maxRow<=k,
    y=[m, 1:maxRow];
else

	% initial value of centroid
    if isRand,
        p = randperm(size(m,1));      % random initialization
        for i=1:k
            c(i,:)=m(p(i),:)
    	end
    else
        for i=1:k
           c(i,:)=m(i,:);        % sequential initialization
    	end
    end

	temp=zeros(maxRow,1);   % initialize as zero vector

	while 1,
        d=DistMatrix(m,c);  % calculate objcets-centroid distances
        [z,g]=min(d,[],2);  % find group matrix g
        if g==temp,
            break;          % stop the iteration
        else
            temp=g;         % copy group matrix to temporary variable
        end
        for i=1:k
            f=find(g==i);
            if f            % only compute centroid if f is not empty
                c(i,:)=mean(m(find(g==i),:),1);
            end
        end
	end

	y=[m,g];

end

%The Matlab function kMeansCluster above call function DistMatrix as shown in the code below. It works for multi-dimensional Euclidean distance. Learn about other type of distance here.

    function d=DistMatrix(A,B)
                 %%%%%%%%%%%%%%%%%%%%%%%%%
                 % DISTMATRIX return distance matrix between points in A=[x1 y1 ... w1] and in B=[x2 y2 ... w2]
                 % Copyright (c) 2005 by Kardi Teknomo,  http://people.revoledu.com/kardi/
                 %
                 % Numbers of rows (represent points) in A and B are not necessarily the same.
                 % It can be use for distance-in-a-slice (Spacing) or distance-between-slice (Headway),
                 %
                 % A and B must contain the same number of columns (represent variables of n dimensions),
                 % first column is the X coordinates, second column is the Y coordinates, and so on.
                 % The distance matrix is distance between points in A as rows
                 % and points in B as columns.
                 % example: Spacing= dist(A,A)
                 % Headway = dist(A,B), with hA ~= hB or hA=hB
                 %          A=[1 2 3; 4 5 6; 2 4 6; 1 2 3]; B=[4 5 1; 6 2 0]
                 %          dist(A,B)= [ 4.69   5.83;
                 %                       5.00   7.00;
                 %                       5.48   7.48;
                 %                       4.69   5.83]
                 %
                 %          dist(B,A)= [ 4.69   5.00     5.48    4.69;
                 %                       5.83   7.00     7.48    5.83]
                 %%%%%%%%%%%%%%%%%%%%%%%%%%%
                 [hA,wA]=size(A);
                 [hB,wB]=size(B);
                 if wA ~= wB,  error(' second dimension of A and B must be the same'); end
                 for k=1:wA
                      C{k}= repmat(A(:,k),1,hB);
                      D{k}= repmat(B(:,k),1,hA);
                 end
                 S=zeros(hA,hB);
                 for k=1:wA
                      S=S+(C{k}-D{k}').^2;
                 end
                 d=sqrt(S);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function triangulation_test26 ( )

%*****************************************************************************80
%
%% TEST26 tests TRIANGULATION_ORDER6_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST26\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER6_PRINT prints out\n' );
  fprintf ( 1, '  an order6 triangulation.\n' );

  [ node_num, triangle_num, hole_num ] = ...
    triangulation_order6_example1_size ( );

  [ node_xy, triangle_node, triangle_neighbor ] = ...
    triangulation_order6_example1 ( node_num, triangle_num );

  triangulation_order6_print ( node_num, triangle_num, node_xy, ...
    triangle_node, triangle_neighbor );

  return
end


text:
%% Copyright (C) 2016 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/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym bessely (@var{alpha}, @var{x})
%% Symbolic Bessel function of the second kind.
%%
%% Example:
%% @example
%% @group
%% syms n x
%% A = 2*bessely(n, x)
%%   @result{} A = (sym) 2⋅bessely(n, x)
%% diff(A)
%%   @result{} (sym) bessely(n - 1, x) - bessely(n + 1, x)
%% @end group
%% @end example
%%
%% @seealso{@@sym/besselj, @@sym/besseli, @@sym/besselk}
%% @end defmethod

function Y = bessely(n, x)

  if (nargin ~= 2)
    print_usage ();
  end

  Y = elementwise_op ('bessely', sym(n), sym(x));

end

%!test
%! X = [1 2 3; 4 5 6];
%! ns = [sym(0) 1 -2; sym(1)/2 -sym(3)/2 pi];
%! n = double(ns);
%! A = double(bessely(ns, X));
%! B = bessely(n, X);
%! assert (all (all (abs (A - B) < 50*eps*abs(A))))

%!test
%! % roundtrip
%! syms x
%! A = bessely(2, 10);
%! q = bessely(2, x);
%! h = function_handle(q);
%! B = h(10);
%! assert (abs (A - B) <= eps*abs(A))

%!error bessely(sym('x'))

repo: octsympy


text:
function [thrs, fars, frrs] = slroc(scores, signs, thres, op)
%SLROC Computes the ROC
%
% $ Syntax $
%   - [thrs, fars, frrs] = slroc(scores, signs, thres, op)
%
% $ Arguments $
%   - scores:           the scores representing the signal intensities
%   - signs:            the signs representing the groundtruth (0 or 1)
%   - thres:            the descriptor indicating how to sample the
%                       thresholds at which the rate is computed
%   - op:               the option stating the attributes of the scores.
%   - thrs:             the sampled threshold values
%   - fars:             the false accept rates at the sampled thresholds
%   - frrs:             the false reject rates at the sampled thresholds
%
% $ Description $
%   - [thrs, fars, frrs] = slroc(scores, signs, thres, op) Computes the
%     ROC of a receiver from the scores and groudtruth signs. The argument
%     thres specifies the sampled thresholds at which the false accept
%     rate and false reject rate is evaluated. If thres is an integer,
%     say n, then n equal-interval integers from lowest to highest scores
%     are taken as samples. The thres can also be a vector containing the
%     sampled thresholds, which should be arranged in ascending order.
%     op states the attributes of the scores, which
%     takes either of the two values: 'high' or 'low'. If op is 'high', a
%     higher score indicates a better match; if op is 'low', a lower score
%     indicates a better match. The cumulative score will be computed
%     up to the number of all classes.
%     For output, if n thresholds are sampled, then thrs, fars, and frrs
%     for all n x 1 vectors, containing the sampled threshold values,
%     and the corresponding false accept rates, and false reject rates.
%
% $ History $
%   - Created by Dahua Lin on Jun 9th, 2005
%   - Modified by Dahua Lin on May 1st, 2006
%     - Base on the sltoolbox v4
%   - Modified by Dahua Lin on Aug 8th, 2006
%     - Base on slhistroc
%

%% parse and verify the input arguments
if nargin < 4
    raise_lackinput('slroc', 4);
end
if ~isequal(size(scores), size(signs))
    error('sltoolbox:sizmismatch', ...
        'The sizes of scores and signs are not match');
end

% the following two statements are disabled in the 2006-08-08 modification
% since scores(signs) will automatic serialize the values
% there is no need to write this statement, just a waste of time and mem
% scores = scores(:);
% signs = logical(signs(:));

if numel(thres) == 1
    n = thres;
    highscore = max(max(scores));
    lowscore = min(min(scores));
    thrs = linspace(lowscore, highscore, n)';
else
    thrs = thres(:);
end

%% compute
scores_a = scores(signs);
scores_r = scores(~signs);

hist_a = histc(scores_a, thrs);
hist_r = histc(scores_r, thrs);

[thrs, fars, frrs] = slhistroc(hist_a, hist_r, thrs, op);

author: lmthang


text:
function [img_pca,img_pca_null]=connectopic_laplacian(s,ind_ins,N,Vn,Null)

% Input:
% s: similarity matrx
% ind_ins: index of subcortical voxels
% N: size of image (MNI152:91x109x91)
% K: degree, used when Disparity=1
% Vn: index of gradient to compute: Vn=2 -> Gradint I; Vn=3 -> Gradient II; Vn=4 -> Gradient III
% Null: property of null data,structural variable including the following:

% Null.NumNull: Number of randomizations if compute the null data
% Null.T: number of time points in the empirical data
% Null.FWHM: Gaussian smooth kernel used in the empirical data
% Null.voxelsize: voxel size in mm in the empirical data

% Output:
% img_pca: eigenmap (main output)
% img_pca_null: eigenmap computed from null data

%Global thresholding. Haak et al 2017
img_pca=zeros(N);
w=squareform(pdist(s));  %similarity to distance mapping

fprintf('Thresholding to minimum density needed for graph to remain connected\n');
ind_upper=find(triu(ones(length(w),length(w)),1));
[~,ind_srt]=sort(w(ind_upper));
w_thresh=zeros(length(w),length(w));
dns=linspace(0.001,1,1000);
for i=1:length(dns)
    ttl=ceil(length(ind_upper)*dns(i));
    w_thresh(ind_upper(ind_srt(1:ttl)))=s(ind_upper(ind_srt(1:ttl)));
    [~,comp_sizes]=get_components(~~w_thresh+~~w_thresh');
    if length(comp_sizes)==1
        break
    end
end

fprintf('Density=%0.2f%%\n',100*(length(find(~~w_thresh))/length(ind_upper)));
dns=dns(i);
w_thresh=w_thresh+w_thresh';

fprintf('Computing Laplacian\n');
L=diag(sum(w_thresh))-w_thresh;

fprintf('Finding eigenvectors\n');
[v,d]=eig(L);d=diag(d);

% Variance explained
per=1./d(2:end);
per=per/sum(per)*100;

%figure; plot(per(1:20));

if v(1,Vn)>v(end,Vn)
    y=v(:,Vn);
else
    y=-v(:,Vn);
end
min_val=min(y);
y=y-min_val;
img_pca(ind_ins)=y;

if nargin==5
    fprintf('Null Model: Synthetic data + MST + Geometry\n');
    NumNull=Null.NumNull;
    T=Null.T;
    FWHM=Null.FWHM;
    voxelsize=Null.voxelsize;

    % Reference matrix
    fprintf('Generating reference matrix\n')
    M=zeros(length(ind_ins),length(ind_ins));
    for i=1:length(ind_ins)
        for j=1:length(ind_ins)
            [xx1,yy1,zz1]=ind2sub(N,ind_ins(i));
            [xx2,yy2,zz2]=ind2sub(N,ind_ins(j));
            dd=sqrt((xx1-xx2)^2 + (yy1-yy2)^2 + (zz1-zz2)^2);
            if dd<=sqrt(2)
                M(i,j)=1;
            end
        end
    end
    ind_m_upper=find(triu(M,1)); % All the available locations (neighboring only)

    % Randomizations

    img_pca_null=zeros([N,NumNull]);
    for nn=1:NumNull
        fprintf('Simulating random data %d\n',nn)
        x=randn([N,T]);
        x_ins=zeros(T,length(ind_ins));
        frst=0;
        for i=1:T
            x(:,:,:,i)=imgaussfilt3(x(:,:,:,i),FWHM/voxelsize/2.355);
            tmp=x(:,:,:,i);
            x_ins(i,:)=tmp(ind_ins);
            show_progress(i,T,frst);frst=1;
        end
        clear x
        % Normalization
        x_ins=detrend(x_ins,'constant'); x_ins=x_ins./repmat(std(x_ins),T,1);

        % Correlation
        c=x_ins'*x_ins; c=c/T;z=atanh(c);

        [~,ind_srt_z]=sort(z(ind_upper),'descend');
        ind_z_upper=ind_upper(ind_srt_z(1:ttl)); % Available locations;ttl is computed from acutual data

        % MST ensures that the graph is fully connected
        % Only allow MST found in the neighboring locations
        ss=zeros(size(s));
        ss(ind_m_upper)=rand(length(ind_m_upper),1);
        ss=ss+ss';

        mst=adjacency(minspantree(graph(ss)));
        ind_mst_upper=find(triu(mst,1));

        % Randomise top weighted edges.
        ind_srt_ttl=ind_upper(ind_srt(1:ttl));
        ind_rand=randperm(length(ind_srt_ttl));
        ind_srt_ttl=ind_srt_ttl(ind_rand);

        % Add edges to mst locations first
        Nm=length(ind_mst_upper);
        w_thresh_null=zeros(size(s));
        w_thresh_null(ind_mst_upper)=s(ind_srt_ttl(1:Nm));

        % Then, add remaining edges to remaining desired locations
        ind_remain=setdiff(ind_z_upper,ind_mst_upper);
        ind_rand=randperm(length(ind_remain));
        ind_remain=ind_remain(ind_rand);
        w_thresh_null(ind_remain(1:ttl-Nm))=s(ind_srt_ttl(Nm+1:ttl));

        [~,comp_sizes]=get_components(~~w_thresh_null+~~w_thresh_null');
        if length(comp_sizes)~=1
            fprintf('Warning: Null model is not fully connected\n')
        end
        dns_null=length(find(~~w_thresh_null))/length(ind_upper);
        fprintf('Density null=%0.2f%%\n',dns_null*100);

        w_thresh_null=w_thresh_null + w_thresh_null';

        fprintf('Computing Laplacian\n');
        L=diag(sum(w_thresh_null))-w_thresh_null;

        fprintf('Finding eigenvectors\n');
        [v,d]=eig(L);d=diag(d);

        if v(1,Vn)>v(end,Vn)
            y=v(:,Vn);
        else
            y=-v(:,Vn);
        end
        min_val_null=min(y);
        y=y-min_val_null;
        tmp=zeros(N);
        tmp(ind_ins)=y;
        img_pca_null(:,:,:,nn)=tmp;

    end
else
    img_pca_null=[];
end

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function cfl = fd1d_heat_implicit_cfl ( k, t_num, t_min, t_max, x_num, ...
  x_min, x_max )

%*****************************************************************************80
%
%% FD1D_HEAT_IMPLICIT_CFL: compute the Courant-Friedrichs-Loewy coefficient.
%
%  Discussion:
%
%    The equation to be solved has the form:
%
%      dUdT - k * d2UdX2 = F(X,T)
%
%    over the interval [X_MIN,X_MAX] with boundary conditions
%
%      U(X_MIN,T) = U_X_MIN(T),
%      U(X_MIN,T) = U_X_MAX(T),
%
%    over the time interval [T_MIN,T_MAX] with initial conditions
%
%      U(X,T_MIN) = U_T_MIN(X)
%
%    The code uses the finite difference method to approximate the
%    second derivative in space, and an explicit forward Euler approximation
%    to the first derivative in time.
%
%    The finite difference form can be written as
%
%      U(X,T+dt) - U(X,T)                  ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) )
%      ------------------  = F(X,T) + k *  ------------------------------------
%               dt                                   dx * dx
%
%    or, assuming we have solved for all values of U at time T, we have
%
%      U(X,T+dt) = U(X,T) + cfl * ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) ) + dt * F(X,T)
%
%    Here "cfl" is the Courant-Friedrichs-Loewy coefficient:
%
%      cfl = k * dt / dx / dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Lindfield, John Penny,
%    Numerical Methods Using MATLAB,
%    Second Edition,
%    Prentice Hall, 1999,
%    ISBN: 0-13-012641-1,
%    LC: QA297.P45.
%
%  Parameters:
%
%    Input, real K, the heat conductivity coefficient.
%
%    Input, integer T_NUM, the number of time values, including the initial
%    value.
%
%    Input, real T_MIN, T_MAX, the minimum and maximum times.
%
%    Input, integer X_NUM, the number of nodes.
%
%    Input, real X_MIN, X_MAX, the minimum and maximum spatial coordinates.
%
%    Output, real CFL, the Courant-Friedrichs-Loewy coefficient.
%
  x_step = ( x_max - x_min ) / ( x_num - 1 );
  t_step = ( t_max - t_min ) / ( t_num - 1 );
%
%  Check the CFL condition, print out its value, and quit if it is too large.
%
  cfl = k * t_step / x_step / x_step;

  if ( 0.5 <= cfl )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FD1D_HEAT_EXPLICIT_CFL - Fatal error!\n' );
    fprintf ( 1, '  CFL condition failed.\n' );
    fprintf ( 1, '  0.5 <= K * dT / dX / dX = %f\n', cfl );
    error ( 'FD1D_HEAT_EXPLICIT_CFL - Fatal error!' );
  end

  return
end


text:
% interp_table_test.m

hr = @(t, J) (1-abs(t/(J/2))) .* (abs(t) <= J/2);
hi = @(t, J) (1-abs(t/(J/2))).^2 .* (abs(t) <= J/2);

% 1D
if 1, printm '1d'
	L = 10;
	J = 6;
	s = [-J/2*L:J/2*L]'/L;
	hs = hr(s, J) + 1i * hi(s, J);
	if length(hs) ~= J*L+1, error 'size', end

	if im
		clf, subplot(211)
		plot(s, real(hs), 'c.-', s, imag(hs), 'y.-')
	end

	K = 20;
	ck = zeros(K,1);
	ck(2+1) = 1;
	ck = complexify(ck);
	tm = linspace(-2*K, 2*K, 2001)';
	fm = interp1_table_mex(ck, hs, int32(J), int32(L), tm);

	if im
		subplot(212)
		plot(tm, real(fm), 'c.-', tm, imag(fm), 'y.-')
	end
prompt
end

% 2D
if 1, printm '2d'
	L = [2^5 2^4];
	J = [6 4];
	s1 = [-J(1)/2*L(1):J(1)/2*L(1)]'/L(1);
	s2 = [-J(2)/2*L(2):J(2)/2*L(2)]'/L(2);
	h1 = 0*hr(s1, J(1)) + 1i * hi(s1, J(1));
	h2 = 0*hr(s2, J(2)) + 1i * hi(s2, J(2));
%	h1 = complexify(h1);
%	h2 = complexify(h2);

	if im
		clf, subplot(211)
		plot(	s1, real(h1), 'c.-', s1, imag(h1), 'y.-', ...
			s2, real(h2), 'g.-', s2, imag(h2), 'm.-')
	end

	K = [8 12];
	ck = zeros(K);
	ck(0+1, 0+1) = -1i;
%	ck = complexify(ck);
	t1 = linspace(-2*K(1), 2*K(1), 201)';
	t2 = linspace(-2*K(2), 2*K(2), 199)';
	t1 = linspace(0, K(1), 201)';
	t2 = linspace(0, K(2), 199)';
	[tt1 tt2] = ndgrid(t1, t2);
	tm = [tt1(:) tt2(:)];
%	tic
	fm = interp2_table_mex(ck, h1, h2, int32(J), int32(L), tm);
%	toc
	fm = reshape(fm, size(tt1));

	im(121, t1, t2, real(fm), 'real'), cbar
	im(122, t1, t2, imag(fm), 'imag'), cbar
prompt
end

% 3D
if 1, printm '3d'
	L = [2^5 2^4 2^6];
	J = [5 3 4];
	s1 = [-J(1)/2*L(1):J(1)/2*L(1)]'/L(1);
	s2 = [-J(2)/2*L(2):J(2)/2*L(2)]'/L(2);
	s3 = [-J(3)/2*L(3):J(3)/2*L(3)]'/L(3);
	h1 = 1 * hr(s1, J(1)) + 0i * hi(s1, J(1));
	h2 = 1 * hr(s2, J(2)) + 0i * hi(s2, J(2));
	h3 = 1 * hr(s3, J(3)) + 0i * hi(s3, J(3));
	h1 = complexify(h1);
	h2 = complexify(h2);
	h3 = complexify(h3);

	if im
		clf, subplot(211)
		plot(	s1, real(h1), 'c.-', s1, imag(h1), 'y.-', ...
			s2, real(h2), 'g.-', s2, imag(h2), 'm.-', ...
			s3, real(h3), 'y.-', s3, imag(h3), 'w.-')
	end

	K = [8 12 4];
	ck = zeros(K);
	ck(0+1, 0+1) = 1i;
	ck = complexify(ck);
	t1 = linspace(-2*K(1), 2*K(1), 69)';
	t2 = linspace(-2*K(2), 2*K(2), 89)';
	t3 = linspace(-2*K(3), 2*K(3), 9)';
%	t1 = linspace(0, K(1), 201)';
%	t2 = linspace(0, K(2), 199)';
	t3 = [-1 0 1];
	[tt1 tt2 tt3] = ndgrid(t1, t2, t3);
	tm = [tt1(:) tt2(:) tt3(:)];
	tic
	fm = interp3_table_mex(ck, h1, h2, h3, int32(J), int32(L), tm);
	toc
	fm = reshape(fm, size(tt1));

	im(121, t1, t2, real(fm), 'real'), cbar
	im(122, t1, t2, imag(fm), 'imag'), cbar
	if length(t3) == 1 % compare to 2D
		f2 = interp2_table_mex(ck, h1, h2, ...
			int32(J(1:2)), int32(L(1:2)), tm(:,1:2));
		f2 = reshape(f2, size(tt1));
		max_percent_diff(f2, fm)
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [H,V] = lp_arn_m(name,Bf,Kf,k,r)
%
%  Arnoldi method w.r.t. inv(F), where F = A-Bf*Kf'.
%
%  Calling sequence:
%
%    [H,V] = lp_arn_m(name,Bf,Kf,k)
%    [H,V] = lp_arn_m(name,Bf,Kf,k,r)
%
%  Input:
%
%    name      basis name of the m-file which generates matrix
%              operations with A, e.g., 'as';
%    Bf        matrix Bf;
%              Set Bf = [] if not existing or zero!
%    Kf        matrix Kf;
%              Set Kf = [] if not existing or zero!
%    k         number of Arnoldi steps (usually k << n, where
%              n is the order of the system);
%    r         initial n-vector
%              (optional - chosen by random, if omitted).
%
%  Output:
%
%    H         matrix H ((k+1)-x-k matrix, upper Hessenberg);
%    V         matrix V (n-x-(k+1) matrix, orthogonal columns).
%
%  User-supplied functions called by this function:
%
%    '[name]_m', '[name]_i'
%
%  Method:
%
%    The ("inverse") Arnoldi method produces matrices V and H such that
%
%      V(:,1) in span{r},
%      V'*V = eye(k+1),
%      inv(F)*V(:,1:k) = V*H.
%
%  Remark:
%
%    This implementation does not check for (near-)breakdown!
%
%
%  LYAPACK 1.0 (Thilo Penzl, May 1999)

% Input data not completely checked!

na = nargin;

with_BK = length(Bf)>0;

eval(lp_e( 'n = ',name,'_m;' ));                    % Get system order.
if k >= n-1, error('k must be smaller than the order of A!'); end
if na<5, r = randn(n,1); end

V = zeros(n,k+1);
H = zeros(k+1,k);

V(:,1) = (1.0/norm(r))*r;

beta = 0;

if with_BK,      % SM = inv(F)*Bf*inv(I-Kf'*inv(F)*Bf)
                 % (This is the main part of the term needed for the low
                 % rank correction in the Sherman-Morrison formula.)
  Im = eye(size(Bf,2));
  TM = Bf;
  eval(lp_e( 'TM = ',name,'_l(''N'',TM);' ));
  SM = TM/(Im-Kf'*TM);
end

for j = 1:k

  if j > 1
    H(j,j-1) = beta;
    V(:,j) = (1.0/beta)*r;
  end

  w = V(:,j);
  eval(lp_e( 'w = ',name,'_l(''N'',w);' ));
  if with_BK, w = w + SM*(Kf'*w); end     % LR correction by SM formula
  r = w;

  for i = 1:j
    H(i,j) = V(:,i)'*w;
    r = r-H(i,j)*V(:,i);
  end

  beta = norm(r);
  H(j+1,j) = beta;

end

V(:,k+1) = (1.0/beta)*r;

author: Sable

repo: mcbench-benchmarks


text:
function [fileList] = getFileList(inputDir, filterExp,varargin)
	% [fileList] = getFileList(inputDir, filterExp,varargin)
	% Gathers a list of files or folders in a directory based on an input regular expression.
	% Biafra Ahanonu
	% started: 2013.10.08 [11:02:31]
	% inputs
	%	inputDir - directory to gather files from and regexp filter for files
	%	filterExp - regexp used to find files/folders.
	% outputs
	%	fileList - cell array of strings containing identified files or folders.

	% changelog
		% 2014.03.21 - added feature to input cell array of filters
		% 2016.03.20 - added exclusion filter to function
		% 2019.03.08 [13:12:59] - added support for natural sorting of files
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2021.09.10 [03:17:56] - Added support to exclude adding the input directory to each file path.
		% 2022.05.26 [22:24:42] - Improved support multiple directory input.
		% 2022.06.12 [16:42:24] - Added folder filter.
	% TODO
		% [DONE] Fix "recusive" (sp) option to recursive in a backwards compatible way.

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% Binary: 1 = recursively find files in all sub-directories. 0 = only find files in inputDir directory.
	options.recursive = '';
	% Append
	options.regexpWithFolder = 0;
	% String: filter to exclude.
	options.excludeFilter = '';
	% Binary: 1 = add inputDir to file paths, 0 = do not add input directory.
	options.addInputDirToPath = 1;
	% Char: lexicographic (e.g. 1 10 11 2 21 22 unless have 01 02 10 11 21 22) or numeric (e.g. 1 2 10 11 21 22) or natural (e.g. 1 2 10 11 21 22)
	% options.sortMethod = 'lexicographic';
	options.sortMethod = 'natural';
	% Binary: 1 = only include folders in the output. 0 = include folders and files.
	options.onlyFolders = 0;
	% DEPRECIATED 1 = recursively find files in all sub-directories. 0 = only find files in inputDir directory.
	options.recusive = 0;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	%     eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================
	if ~isempty(options.recursive)
		options.recusive = options.recursive;
	end
	% options.recursive = options.recusive;

	if ~iscell(filterExp)
		filterExp = {filterExp};
	end
	if ischar(inputDir)
		inputDir = {inputDir};
	end

	fileList = {};
	nDirs = length(inputDir);
	%for thisDir = inputDir
	for i = 1:nDirs
		thisDir = inputDir{i};
		% thisDirHere = thisDir{i};
		thisDirHere = thisDir;
		if options.recusive==0
			files = dir(thisDirHere);
		else
			files = dirrec(thisDirHere)';
		end
		for file=1:length(files)
			if options.recusive==0
				filename = files(file,1).name;
				if options.regexpWithFolder==1
					filename = [options.regexpWithFolder filesep filename];
				end
				% If option selected, remove non-folders.
				if options.onlyFolders==1&&isfolder(filename)==0
					continue
				end
				% Add found file/folder to the list.
				if(~isempty(cell2mat(regexpi(filename, filterExp))))
					if options.addInputDirToPath==1
						fileList{end+1} = [thisDirHere filesep filename];
					else
						[~,filename,filenameExt] = fileparts(filename);
						filename = [filename filenameExt];
						fileList{end+1} = [filename];
					end
				end
			else
				filename = files(file,:);
				filename = filename{1};
				% If option selected, remove non-folders.
				if options.onlyFolders==1&&isfolder(filename)==0
					continue
				end
				% Add found file/folder to the list.
				if(~isempty(cell2mat(regexpi(filename, filterExp))))
					if options.addInputDirToPath==1
						fileList{end+1} = [filename];
					else
						[~,filename,filenameExt] = fileparts(filename);
						filename = [filename filenameExt];
						fileList{end+1} = [filename];
					end
				end
			end
		end
	end
	if ~isempty(options.excludeFilter)
		% excludeIdx = find(~cellfun(@isempty,(regexp(fileList,options.excludeFilter))));
		excludeIdx = ~cellfun(@isempty,(regexp(fileList,options.excludeFilter)));
		fileList(excludeIdx) = [];
		% includeIdx = setdiff(1:length(fileList),excludeIdx)
	end
	if ~isempty(fileList)&&(strcmp(options.sortMethod,'natural')|strcmp(options.sortMethod,'numeric'))
		fileList = natsortfiles(fileList);
	end
end


text:
% At_f.m
%
% Adjoint for "scrambled Fourier" measurements.
%
% Usage: x = At_f(b, N, OMEGA, P)
%
% b - K vector = [real part; imag part]
%
% N - length of output x
%
% OMEGA - K/2 vector denoting which Fourier coefficients to use
%         (the real and imag parts of each freq are kept).
%
% P - Permutation to apply to the input vector.  Fourier coeffs of
%     x(P) are embedded.
%     Default = 1:N (no scrambling).
%
% Written by: Justin Romberg, Caltech
% Created: October 2005
% Email: jrom@acm.caltech.edu
%

function x = At_f(b, N, OMEGA, P)

if (nargin < 4),  P = 1:N;  end

K = length(b);
fx = zeros(N,1);
fx(OMEGA) = sqrt(2)*b(1:K/2) + i*sqrt(2)*b(K/2+1:K);
x = zeros(N,1);
x(P) = sqrt(N)*real(ifft(fx));

repo: REST


text:
classdef TestAddWeighted
    %TestAddWeighted

    methods (Static)
        function test_images
            img1 = cv.imread(fullfile(mexopencv.root(),'test','img001.jpg'), 'ReduceScale',2);
            img2 = randi([0 255], size(img1), 'uint8');

            out = cv.addWeighted(img1,1.0, img2,0.4, 0.2);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            if mexopencv.require('images')
                expected = my_add_weighted(img1,1.0, img2,0.4, 0.2);
                assert(isequal(out, expected));
            end
        end

        function test_2d_matrices
            A = uint8([100 150; 200 250]);
            B = uint8([200 100; 100 100]);

            out = cv.addWeighted(A,0.2, B,0.8, 0.0);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = A*0.2 + B*0.8 + 0;
            assert(isequal(out, expected));

            out = cv.addWeighted(A,0.2, 1,0.8, 0.0);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = A*0.2 + 1*0.8 + 0;
            assert(isequal(out, expected));

            out = cv.addWeighted(1,0.2, B,0.8, 0.0);
            validateattributes(out, {class(B)}, {'size',size(B)});
            expected = 1*0.2 + B*0.8 + 0;
            assert(isequal(out, expected));

            out = cv.addWeighted(1,0.2, 1,0.8, 0.0);
            validateattributes(out, {'double'}, {'scalar'});
            expected = 1*0.2 + 1*0.8 + 0;
            assert(isequal(out, expected));
        end

        function test_output_depth
            out = cv.addWeighted(uint16(1),0.5, int8(1),0.5, 0.0, 'DType','single');
            validateattributes(out, {'single'}, {'scalar'});
            expected = single(1*0.5 + 1*0.5 + 0);
            assert(isequal(out, expected));
        end

        function test_error_argnum
            try
                cv.addWeighted();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

function out = my_add_weighted(src1, a, src2, b, c, dtype)
    %MY_ADD_WEIGHTED  Similar to cv.addWeighted using imlincomb from IPT

    if nargin < 6, dtype = class(src1); end

    % add two weighted images with specified output class
    out = imlincomb(a,src1, b,src2, c, dtype);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function point = intersectLines(line1, line2, varargin)
%INTERSECTLINES Return all intersection points of N lines in 2D
%
%   PT = intersectLines(L1, L2);
%   returns the intersection point of lines L1 and L2. L1 and L2 are 1-by-4
%   row arrays, containing parametric representation of each line (in the
%   form [x0 y0 dx dy], see 'createLine' for details).
%
%   In case of colinear lines, returns [Inf Inf].
%   In case of parallel but not colinear lines, returns [NaN NaN].
%
%   If each input is [N*4] array, the result is a [N*2] array containing
%   intersections of each couple of lines.
%   If one of the input has N rows and the other 1 row, the result is a
%   [N*2] array.
%
%   PT = intersectLines(L1, L2, EPS);
%   Specifies the tolerance for detecting parallel lines. Default is 1e-14.
%
%   Example
%   line1 = createLine([0 0], [10 10]);
%   line2 = createLine([0 10], [10 0]);
%   point = intersectLines(line1, line2)
%   point =
%       5   5
%
%   See also
%   lines2d, edges2d, intersectEdges, intersectLineEdge
%   intersectLineCircle
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

%   HISTORY
%   2004-02-19 add support for multiple lines.
%   2007-03-08 update doc
%   2011-10-07 code cleanup

%% Process input arguments

% extract tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% check size of each input
N1 = size(line1, 1);
N2 = size(line2, 1);
N = max(N1, N2);
if N1 ~= N2 && N1*N2 ~= N
    error('matGeom:IntersectLines:IllegalArgument', ...
        'The two input arguments must have same number of lines');
end

%% Check parallel and colinear lines

% coordinate differences of origin points
dx = bsxfun(@minus, line2(:,1), line1(:,1));
dy = bsxfun(@minus, line2(:,2), line1(:,2));

% indices of parallel lines
denom = line1(:,3) .* line2(:,4) - line2(:,3) .* line1(:,4);
par = abs(denom) < tol;

% indices of colinear lines
col = abs(dx .* line1(:,4) - dy .* line1(:,3)) < tol & par ;

% initialize result array
x0 = zeros(N, 1);
y0 = zeros(N, 1);

% initialize result for parallel lines
x0(col) = Inf;
y0(col) = Inf;
x0(par & ~col) = NaN;
y0(par & ~col) = NaN;

% in case all line couples are parallel, return
if all(par)
    point = [x0 y0];
    return;
end

%% Extract coordinates of itnersecting lines

% indices of intersecting lines
inds = ~par;

% extract base coordinates of first lines
if N1 > 1
    line1 = line1(inds,:);
end
x1 =  line1(:,1);
y1 =  line1(:,2);
dx1 = line1(:,3);
dy1 = line1(:,4);

% extract base coordinates of second lines
if N2 > 1
    line2 = line2(inds,:);
end
x2 =  line2(:,1);
y2 =  line2(:,2);
dx2 = line2(:,3);
dy2 = line2(:,4);

% re-compute coordinate differences of origin points
dx = bsxfun(@minus, line2(:,1), line1(:,1));
dy = bsxfun(@minus, line2(:,2), line1(:,2));

%% Compute intersection points

denom = denom(inds);
x0(inds) = (x2 .* dy2 .* dx1 - dy .* dx1 .* dx2 - x1 .* dy1 .* dx2) ./ denom ;
y0(inds) = (dx .* dy1 .* dy2 + y1 .* dx1 .* dy2 - y2 .* dx2 .* dy1) ./ denom ;

% concatenate result
point = [x0 y0];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


%LEARNALGOMLEUCLIDEAN Simulates plain L2 distance as learning algorithm.
classdef LearnAlgoMLEuclidean < LearnAlgo

    properties
       p %parameters
       s %struct
    end

    properties (Constant)
        type = 'identity'
    end

    methods
        function obj = LearnAlgoMLEuclidean(p)
            if nargin < 1
               p = struct();
            end

            if ~isfield(p,'roccolor')
                p.roccolor = 'r';
            end

            obj.p = p;
        end

        function s = learnPairwise(obj,X,idxa,idxb,matches)
            s.M = eye(size(X,1));
            s.t = 0.0;
            s.learnAlgo = obj;
            s.roccolor = obj.p.roccolor;
        end

        function s = learn(obj,X,y)
            s.M = eye(size(X,1));
            s.t = 0.0;
            s.learnAlgo = obj;
            s.roccolor = obj.p.roccolor;
        end

        function d = dist(obj, s, X, idxa,idxb)
            d = cdistM(s.M,X,idxa,idxb);
        end
    end
end


text:
function [positions, time] = tracker(video_path, img_files, pos, target_sz, ...
	padding, kernel, lambda, output_sigma_factor, interp_factor, cell_size, ...
	features, show_visualization)
%TRACKER Kernelized/Dual Correlation Filter (KCF/DCF) tracking.
%   This function implements the pipeline for tracking with the KCF (by
%   choosing a non-linear kernel) and DCF (by choosing a linear kernel).
%
%   It is meant to be called by the interface function RUN_TRACKER, which
%   sets up the parameters and loads the video information.
%
%   Parameters:
%     VIDEO_PATH is the location of the image files (must end with a slash
%      '/' or '\').
%     IMG_FILES is a cell array of image file names.
%     POS and TARGET_SZ are the initial position and size of the target
%      (both in format [rows, columns]).
%     PADDING is the additional tracked region, for context, relative to
%      the target size.
%     KERNEL is a struct describing the kernel. The field TYPE must be one
%      of 'gaussian', 'polynomial' or 'linear'. The optional fields SIGMA,
%      POLY_A and POLY_B are the parameters for the Gaussian and Polynomial
%      kernels.
%     OUTPUT_SIGMA_FACTOR is the spatial bandwidth of the regression
%      target, relative to the target size.
%     INTERP_FACTOR is the adaptation rate of the tracker.
%     CELL_SIZE is the number of pixels per cell (must be 1 if using raw
%      pixels).
%     FEATURES is a struct describing the used features (see GET_FEATURES).
%     SHOW_VISUALIZATION will show an interactive video if set to true.
%
%   Outputs:
%    POSITIONS is an Nx2 matrix of target positions over time (in the
%     format [rows, columns]).
%    TIME is the tracker execution time, without video loading/rendering.
%
%   Joao F. Henriques, 2014

	%if the target is large, lower the resolution, we don't need that much
	%detail
	resize_image = (sqrt(prod(target_sz)) >= 100);  %diagonal size >= threshold
	if resize_image,
		pos = floor(pos / 2);
		target_sz = floor(target_sz / 2);
	end

	%window size, taking padding into account
	window_sz = floor(target_sz * (1 + padding));

% 	%we could choose a size that is a power of two, for better FFT
% 	%performance. in practice it is slower, due to the larger window size.
% 	window_sz = 2 .^ nextpow2(window_sz);

	%create regression labels, gaussian shaped, with a bandwidth
	%proportional to target size

	output_sigma = sqrt(prod(target_sz)) * output_sigma_factor / cell_size;
    if isfield(features, 'deep') && features.deep
        yf = fft2(gaussian_shaped_labels(output_sigma, ceil(window_sz / cell_size)));
%         sz = ceil(window_sz/cell_size)-1+4-4;
%         yf = fft2(gaussian_shaped_labels(output_sigma, sz));

    else
        yf = fft2(gaussian_shaped_labels(output_sigma, floor(window_sz / cell_size)));
    end

	%store pre-computed cosine window
	cos_window = hann(size(yf,1)) * hann(size(yf,2))';

	if show_visualization,  %create video interface
		update_visualization = show_video(img_files, video_path, resize_image);
	end

	%note: variables ending with 'f' are in the Fourier domain.

	time = 0;  %to calculate FPS
	positions = zeros(numel(img_files), 2);  %to calculate precision

	for frame = 1:numel(img_files),
		%load image
		im = imread([video_path img_files{frame}]);
        if ~isfield(features, 'deep')
            if size(im,3) > 1,
                im = rgb2gray(im);
            end
        end

		if resize_image,
			im = imresize(im, 0.5);
		end

		tic()

		if frame > 1,
			%obtain a subwindow for detection at the position from last
			%frame, and convert to Fourier domain (its size is unchanged)
			patch = get_subwindow(im, pos, window_sz);
			zf = fft2(get_features(patch, features, cell_size, cos_window));

			%calculate response of the classifier at all shifts
			switch kernel.type
			case 'gaussian',
				kzf = gaussian_correlation(zf, model_xf, kernel.sigma);
			case 'polynomial',
				kzf = polynomial_correlation(zf, model_xf, kernel.poly_a, kernel.poly_b);
			case 'linear',
				kzf = linear_correlation(zf, model_xf);
			end
			response = real(ifft2(model_alphaf .* kzf));  %equation for fast detection

			%target location is at the maximum response. we must take into
			%account the fact that, if the target doesn't move, the peak
			%will appear at the top-left corner, not at the center (this is
			%discussed in the paper). the responses wrap around cyclically.
			[vert_delta, horiz_delta] = find(response == max(response(:)), 1);
			if vert_delta > size(zf,1) / 2,  %wrap around to negative half-space of vertical axis
				vert_delta = vert_delta - size(zf,1);
			end
			if horiz_delta > size(zf,2) / 2,  %same for horizontal axis
				horiz_delta = horiz_delta - size(zf,2);
			end
			pos = pos + cell_size * [vert_delta - 1, horiz_delta - 1];
		end

		%obtain a subwindow for training at newly estimated target position
		patch = get_subwindow(im, pos, window_sz);
		xf = fft2(get_features(patch, features, cell_size, cos_window));

		%Kernel Ridge Regression, calculate alphas (in Fourier domain)
		switch kernel.type
		case 'gaussian',
			kf = gaussian_correlation(xf, xf, kernel.sigma);
		case 'polynomial',
			kf = polynomial_correlation(xf, xf, kernel.poly_a, kernel.poly_b);
		case 'linear',
			kf = linear_correlation(xf, xf);
		end
		alphaf = yf ./ (kf + lambda);   %equation for fast training

		if frame == 1,  %first frame, train with a single image
			model_alphaf = alphaf;
			model_xf = xf;
		else
			%subsequent frames, interpolate model
			model_alphaf = (1 - interp_factor) * model_alphaf + interp_factor * alphaf;
			model_xf = (1 - interp_factor) * model_xf + interp_factor * xf;
		end

		%save position and timing
		positions(frame,:) = pos;
		time = time + toc();

		%visualization
		if show_visualization,
			box = [pos([2,1]) - target_sz([2,1])/2, target_sz([2,1])];
			stop = update_visualization(frame, box);
			if stop, break, end  %user pressed Esc, stop early

			drawnow
% 			pause(0.05)  %uncomment to run slower
		end

	end

	if resize_image,
		positions = positions * 2;
	end
end

repo: KCF


text:
function [ y, m, d ] = day_carry_common ( y, m, d )

%*****************************************************************************80
%
%% DAY_CARRY_COMMON carries days to months in a Common date.
%
%  Discussion:
%
%    While ( number of days in M ) < D:
%      decrease the day D by the number of days in the month M;
%      increase M by 1;
%      if necessary, adjust Y.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, the YMD date.
%
%    Output, integer Y, integer M, integer D, the YMD date.
%    On output, D is between 1 and the number of days in M.
%

%
%  If the date is in the transition month, deflate it,
%  so we can perform ordinary arithmetic.
%
  [ y, m, d ] = deflate_common ( y, m, d );

  days = month_length_common ( y, m );

  while ( days < d )

    d = d - days;
    m = m + 1;
    days = month_length_common ( y, m );
%
%  Make sure the month isn't too big.
%
    [ y, m ] = month_carry_common ( y, m );

  end
%
%  If the date is in the transition month, inflate it.
%
  [ y, m, d ] = inflate_common ( y, m, d );

  return
end


text:
%{
Author: Deepti Ghadiyaram

Description: Given an input RGB image, this method extracts the 564 FRIQUEE
features proposed in [1] in three different color spaces.

1. We first convert the given input image into different color spaces -
LAB, LMS, and HSI.

2. We then extract the required components from each of these color spaces
eg: Hue and Saturation components, M and S components.

3. We also compute the chroma map from the A and B components of the image
in LAB color space.

Input: a MXNX3 array upon reading any image.

Output: Features extracted in luminance, chroma, LMS color components of the
image.

Dependencies: This method depends on the following methods:
computeSigmaMap.m
divisiveNormalization.m
%}
function feat = DoGFeat(imGray)

    %Get the sigma map of the image.
    sigmaMap = computeSigmaMap(imGray);

    % Construct two Gaussian windows where sig2 = 1.5*sig1;
    k=1.5;

    window1 = fspecial('gaussian',7,7/6);
    window1 = window1/sum(sum(window1));

    window2 = fspecial('gaussian',7,7*k/6);
    window2 = window2/sum(sum(window2));

    %Compute DOG of the sigmaMap
    DoGSigma= filter2((window1-window2), sigmaMap, 'same');

    %Compute the sigmaMap of DoGSigma
    DoGSigma1 = computeSigmaMap(DoGSigma);

    %DNT of DOG of sigma.
    divNormDoGSigma = divisiveNormalization(DoGSigma);

    %shape, variance, skewness.
    [alpha,lsigma] = estimateggdparam(divNormDoGSigma(:));

    % Apply DNT on DoGSigma1
    divNormDoGSigma1 = divisiveNormalization(DoGSigma1);

    % The model features are the GGD parameters of divNormDoGSigma (DNT of
    % dogSigma)
    modelFeats = [alpha lsigma];

    % The sample statistical features are extracted from both
    % divNormDoGSigma and divNormDoGSigma1
    sampleFeats = [skewness(divNormDoGSigma(:)) kurtosis(divNormDoGSigma(:)) skewness(divNormDoGSigma1(:)) kurtosis(divNormDoGSigma1(:))];

    feat = [modelFeats sampleFeats];
end

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
function [T, DISTR, df, Ybar, S] = symNormalTest2(TEST_TYPE, Ybar1, S1, n1, Ybar2, S2, n2, COV_TYPE, COV_DIFF)

% Two-sample tests for symmetric matrices.
%
%   [T, DISTR, df, Ybar, S] = symNormalTest2(TEST_TYPE, Ybar1, S1, n1, Ybar2, [S2], [n2], [COV_TYPE], [COV_DIFF])
%
% Input:
%   TEST_TYPE   Controls the type of test (assuming equal variances
%               between the two groups):
%                   'full': H0: both groups have the same mean tensors.
%                   'val' : H0: both groups have the same eigenvalues,
%                           with possibly different unknown eigenvectors.
%                   'vec' : H0: both groups have the same eigenvectors,
%                           with common unknown eigenvalues.
%   Ybar1, Ybar2  pxpx[sz] arrays of mean matrices for each group
%   S1, S2        qxqx[sz] arrays of covariance matrices for each group, q=p(p+1)/2
%                      (default S2 = 0)
%   n1, n2        Number of subjects in each group (default n2 = 1).
%   COV_TYPE      Type of covariance: 'spherical', 'orth-inv' or 'full' (default).
%   COV_DIFF      Assume different covariances: 1 (yes - default) or 0 (no)
% Output:
%   T           [sz] array of test statistics
%   DISTR       'f' or 'gamma'
%   df          2x1 degrees of freedom of F, or 2x[sz] parameters of Gamma (nu/2, 2a)
%   Ybar        pxpx[sz] array of pooled mean matrices
%   S           qxqx[sz] array of pooled covariance matrices
%
% E.g.:
%   [Ybar1, S1, n1] = symNormalStats(Y(:,:,g1));
%   [Ybar2, S2, n2] = symNormalStats(Y(:,:,g2));
%   [T, DISTR, df] = symNormalTest2('vec', Ybar1, S1, n1, Ybar2, S2, n2);
%
% Copyright by Armin Schwartzman, 2009

% HISTORY:
%   2008.12.30 ASH (armins@hsph.harvard.edu) wrote it.
%

% Check inputs
if ~exist('S2'),
    S2 = 0;
end
if ~exist('n2'),
    n2 = 1;
end
if (size(Ybar1,1) ~= size(Ybar1,2) | size(S1,1) ~= size(S1,2) | ...
    size(Ybar2,1) ~= size(Ybar2,2) | size(S2,1) ~= size(S2,2)),
    error('Wrong input format');
end
if (size(Ybar1,1) ~= size(Ybar2,1) | size(S1,1) ~= size(S2,1)),
    error('Wrong input format');
end
if ~exist('COV_TYPE'), COV_TYPE = 'full'; end
if (~strmatch(COV_TYPE,'spherical') & ~strmatch(COV_TYPE,'rot-inv') & ~strmatch(COV_TYPE,'full')),
    error('Only spherical, rot-inv and full covariance types supported.')
end
if ~exist('COV_DIFF'), COV_DIFF = 1; end

% Constants
n = n1 + n2;
p = size(Ybar1, 1);
q = size(S1, 1);
if (q ~= p*(p+1)/2),
    error('Wrong input format');
end

% Pooled mean and covariance
Ybar = (n1*Ybar1 + n2*Ybar2)/n;
S = ((n1-1)*S1 + (n2-1)*S2)/(n-2);

% Test type
switch TEST_TYPE,
    case 'full',
        d = permute(vecd(Ybar1 - Ybar2), [1 ndims(Ybar1)+1 2:ndims(Ybar1)]);
        switch COV_TYPE,
        case 'spherical',
            T = n1*n2/n * ndfun('mult', permute(d, [2 1 3:ndims(d)]), d);
            DISTR = 'f';
            if COV_DIFF,
                error('Full test, different spherical cov., not implemented.')
            end
            s2 = S(1,1,:);
            df = [q; q*(n-2)];
            T = df(2)/df(1) * T./(q*(n-2)*s2);
        case 'orth-inv',
            error('Full test, orth-inv cov., not implemented.')
        case 'full',
            if COV_DIFF,
                S = S1/n1 + S2/n2;
                Sinv = ndfun('inv', S);
                T = ndfun('mult', permute(d, [2 1 3:ndims(d)]), ndfun('mult', Sinv, d));
                Sinv1 = ndfun('mult', Sinv, ndfun('mult', S1/n1, Sinv));
                T1 = ndfun('mult', permute(d, [2 1 3:ndims(d)]), ndfun('mult', Sinv1, d));
                Sinv2 = ndfun('mult', Sinv, ndfun('mult', S2/n2, Sinv));
                T2 = ndfun('mult', permute(d, [2 1 3:ndims(d)]), ndfun('mult', Sinv2, d));
                m = 1./((T1./T).^2/(n1-1) + (T2./T).^2/(n2-1));
                df(2,:) = m-q+1; df(1,:) = q;
                DISTR = 'f';
                T = (m-q+1)./(q*m) .* T;
                df = shiftdim(df);
            else
                Sinv = ndfun('inv', S);
                T = n1*n2/n * ndfun('mult', permute(d, [2 1 3:ndims(d)]), ndfun('mult', Sinv, d));
                df = [q; n-q-1];
                DISTR = 'f';
                T = df(2)./df(1) .* T./(n-2);
            end
        end

    case 'val',
        [V1,L1] = ndSymEig(Ybar1); % [V1,L1] = ndfunm('eig', Ybar1);
        [V2,L2] = ndSymEig(Ybar2); % [V2,L2] = ndfunm('eig', Ybar2);
        T = n1*n2/n * sum(sum((L1 - L2).^2, 1), 2);
        switch COV_TYPE,
        case 'spherical',
            DISTR = 'f';
            if COV_DIFF,
                error('Eigval test, different spherical cov., not implemented.')
            end
            s2 = S(1,1,:);
            df = [p; q*(n-2)];
            T = df(2)/df(1) * T./(q*(n-2)*s2);
        case 'orth-inv',
            error('Full test, orth-inv cov., not implemented.')
        case 'full',
            Omega = 0;
            for i=1:p,
                W = cat(1, v(i,i,V1), -v(i,i,V2));
                Omega = Omega + ndfun('mult', W, permute(W, [2 1 3:ndims(W)]));
            end
            Omega = (n1*n2)/n * Omega;
            SS = zeros(size(Omega));  % (2q)x(2q)x[]
            if COV_DIFF,
                SS(1:q,1:q,:) = S1(1:q,1:q,:)/n1; SS(q+1:2*q,q+1:2*q,:) = S2(1:q,1:q,:)/n2;
            else
                SS(1:q,1:q,:) = S(1:q,1:q,:)/n1; SS(q+1:2*q,q+1:2*q,:) = S(1:q,1:q,:)/n2;
            end
            [a,nu] = chi2approx(SS, Omega);
            DISTR = 'gamma';
            df = cat(1, nu/2, 2*a);  % shape parameter, scale parameter
        end

    case 'vec',
        [V1,L1] = ndSymEig(Ybar1); % [V1,L1] = ndfunm('eig', Ybar1);
        [V2,L2] = ndSymEig(Ybar2); % [V2,L2] = ndfunm('eig', Ybar2);
        T = 2*n1*n2/n * (sum(sum(L1.*L2, 1), 2) - sum(sum(Ybar1.*Ybar2, 1), 2));
        switch COV_TYPE,
        case 'spherical',
            DISTR = 'f';
            if COV_DIFF,
                error('Eigvec test, different spherical cov., not implemented.')
            end
            s2 = S(1,1,:);
            df = [q-p; q*(n-2)];
            T = df(2)/df(1) * T./(q*(n-2)*s2);
        case 'orth-inv',
            error('Eigvec test, orth-inv cov., not implemented.')
        case 'full',
            Omega = 0;
            for i=1:p,
                for j=1:p,
                    W = omega(i,j,n1,n2,V1,V2);
                    Omega = Omega + ndfun('mult', W, permute(W, [2 1 3:ndims(W)]));
                end
            end
            Omega = (n1*n2)/n * Omega;
            SS = zeros(size(Omega));  % (2q)x(2q)x[]
            if COV_DIFF,
                SS(1:q,1:q,:) = S1(1:q,1:q,:)/n1; SS(q+1:2*q,q+1:2*q,:) = S2(1:q,1:q,:)/n2;
            else
                SS(1:q,1:q,:) = S(1:q,1:q,:)/n1; SS(q+1:2*q,q+1:2*q,:) = S(1:q,1:q,:)/n2;
            end
            [a,nu] = chi2approx(SS, Omega);
            DISTR = 'gamma';
            df = cat(1, nu/2, 2*a);  % shape parameter, scale parameter
        end
end

% Adjust output
T = shiftdim(T, 2);

end

%------------------------------------------------------------------------
% Auxiliary functions

function v = v(i,j,U)
    Ui = zeros(size(U));
    Ui(:,1,:) = U(:,i,:);
    Uj = zeros(size(U));
    Uj(1,:,:) = U(:,j,:);
    v = ndfun('mult', Ui, Uj);
    v = vecd((v + permute(v, [2 1 3:ndims(v)]))/2);
    v = permute(v, [1 ndims(v)+1 2:ndims(v)]);
end

function W = omega(i,j,n1,n2,U1,U2)
    sz = size(U1);
    p = sz(1);
    v1 = v(j, i, permute(U1, [2 1 3:ndims(U1)]))/2;
    v2 = v(j, i, permute(U2, [2 1 3:ndims(U1)]))/2;
    b = zeros([p 1 sz(3:end)]);
    b(1:p,:) = (n1*v2(1:p,:) + n2*v1(1:p,:))/(n1+n2);
    Eji = zeros(size(U1)); Eji(j,i,:) = 1/2;  Eji(i,j,:) = 1/2;
    Eji = permute(vecd(Eji), [1 ndims(Eji)+1 2:ndims(Eji)]);
    w1 = Eji - ndfun('mult', B(U1), b);
    w2 = Eji - ndfun('mult', B(U2), b);
    W = cat(1, w1, -w2);
end

function B = B(U)
    sz = size(U);
    p = sz(1);
    q = p*(p+1)/2;
    B = zeros([q p sz(3:end)]);
    for i=1:p,
        vv = v(i,i,U);
        B(:,i,:) = vv(1:q,1,:);
    end
end

function [a,nu] = chi2approx(S,Q)
    SQ = ndfun('mult',S,Q);
    k1 = ndfunm('trace',SQ);
    k2 = ndfunm('trace',ndfun('mult',SQ,SQ));
    a = k2./k1; nu = k1.^2./k2;
end

%------------------------------------------------------------------------
% Debugging
% M = zeros(3);
% S = eye(6);
% Y = symNormalRnd(M, S, [100 4]);
% [Ybar1, S1, n1] = symNormalStats(Y(:,:,1:40,:), 'full');
% [Ybar2, S2, n2] = symNormalStats(Y(:,:,41:100,:), 'full');
% [T, DISTR, df, Ybar, S] = symNormalTest2('val', Ybar1, S1, n1, Ybar2, S2, n2)

repo: vistasoft


text:
% READTXTFILE - Read text file into a Matlab variable
%
% Usage: >> str = readtxtfile( filename );
%
% Input:
%   filename - [string] name of the file.
%
% Output:
%   str - [string] content of the file.
%
% Author: Arnaud Delorme, SCCN / INC / UCSD, August 2009

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

function com = readtxtfile(filename)

    com = '';
    if ~exist(filename)
        disp([ 'Cannot find option file ' filename ]);
    else
        fid = fopen(filename, 'r');
        if fid == -1
            disp([ 'Cannot open option file ' filename ]);
        else
            com = '';
            while ~feof(fid)
                a = fgetl(fid);
                com = [ com 10 a ];
            end
            fclose(fid);
        end
    end


text:
% Line Data for B-Bus (Shunt Admittance)Formation.

function bbus = bbusppg()     % Returns B-bus..

linedata = linedata30();
fb = linedata(:,1);
tb = linedata(:,2);
b = linedata(:,5);
nbus = max(max(fb),max(tb));    % no. of buses...
nbranch = length(fb);           % no. of branches...
bbus = zeros(nbus,nbus);

 for k=1:nbranch
     bbus(fb(k),tb(k)) = b(k);
     bbus(tb(k),fb(k)) = bbus(fb(k),tb(k));
 end


text:
function [ss, dist_map] = makeSphericalSection(radius, height, width, plot_section, binary)
%MAKESPHERICALSECTION     Create a binary map of a sphere segment within a 3D grid.
%
% DESCRIPTION:
%       makeSphericalSection creates a binary map of a section of a
%       spherical surface within a three-dimensional matrix. The sphere is
%       created using an extension of the midpoint circle algorithm. A
%       single grid point is taken as the sphere center so the total
%       diameter will always be an odd number of grid points. The sphere is
%       then truncated based on the values for height and width (a diagram
%       of the input sizes is given below). The face of the spherical
%       section faces in the positive x-direction and the optional width
%       parameter truncates the size in the y-direction.
%
%       If the optional input parameter "binary" is set to false (the
%       default), the section map is returned as a double precision matrix.
%       If it is set to true, the map is returned as a logical matrix. The
%       average distance between each grid point in the spherical section
%       and its contiguous neighbours can also be returned. This is given
%       as a ratio compared to the average neighbour distance for a flat
%       surface.
%
% SYNTAX:
%       ss = makeSphericalSection(radius, height)
%       ss = makeSphericalSection(radius, height, width)
%       ss = makeSphericalSection(radius, height, width, plot_section)
%       ss = makeSphericalSection(radius, height, [], plot_section)
%       ss = makeSphericalSection(radius, height, width, plot_section, binary)
%       ss = makeSphericalSection(radius, height, [], [] binary)
%
%       [ss, dist_map] = makeSphericalSection(radius, height)
%       [ss, dist_map] = makeSphericalSection(radius, height, width)
%       [ss, dist_map] = makeSphericalSection(radius, height, [], plot_section)
%       [ss, dist_map] = makeSphericalSection(radius, height, width, plot_section)
%       [ss, dist_map] = makeSphericalSection(radius, height, width, plot_section, binary)
%       [ss, dist_map] = makeSphericalSection(radius, height, [], [] binary)
%
% INPUTS:
%       radius          - radius of curvature [grid points]
%       height          - transducer height [grid points]
%
% OPTIONAL INPUTS:
%       width           - section width (must be specified as an odd
%                         number) [grid points]
%       plot_section    - Boolean controlling whether the spherical section
%                         is plotted using voxelPlot (default = false)
%       binary          - Boolean controlling whether the spherical section
%                         is returned as a double precision matrix (false)
%                         or a logical matrix (true) (default = false)
%
% OUTPUTS:
%       ss              - binary matrix containing spherical section
%       dist_map        - ratio of average neighbour distance for each grid
%                         point within the spherical section compared to a
%                         flat surface
%
% ABOUT:
%       author          - Bradley Treeby
%       date            - 3rd February 2012
%       last update     - 20th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also makeBall, makeSphere

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

%#ok<*UNRCH>

% option to use the spherical sections or square sections of the sphere
% returned by makeSphere (true by default)
use_spherical_sections = true;

% force inputs to be integers
radius = round(radius);
height = round(height);

% check for width input
if nargin < 3 || (nargin == 4 && isempty(width))

    % set width truncation flag to false
    use_width = false;

else

    % set width truncation flag to true
    use_width = true;

    % force input to an integer
    width = round(width);

    % check that it's an odd number
    if ~rem(width, 2)
        error('width must be an odd number');
    end

end

% check for plot input
if nargin < 4 || isempty(plot_section)
    plot_section = false;
end

% check for data type input
if nargin < 5 || isempty(binary)
    binary = false;
end

% calculate minimum grid dimensions to fit entire sphere
Nx = 2*radius + 1;

% create sphere
ss = makeSphere(Nx, Nx, Nx, radius, [], binary);

% truncate to given height
if use_spherical_sections
    ss = ss(1:height, :, :);
else
    ss = permute(ss(:, 1:height, :), [2, 3, 1]);
end

% flatten transducer and store the maximum and indices
mx = squeeze(max(ss, [], 1));

% calculate the total length/width of the transducer
length = sum(mx(ceil(end/2), :), 2);

% truncate transducer grid based on length (removes empty rows and columns)
offset = (Nx - length)/2;
ss = ss(:, 1 + offset:end - offset, 1 + offset:end - offset);

% also truncate to given width if defined by user
if use_width

    % check the value is appropriate
    if width > length
        error('input for width must be less than or equal to transducer length');
    end

    % calculate offset
    offset = (length - width)/2;

    % truncate transducer grid
    ss = ss(:, 1 + offset:end - offset, :);

end

% compute average distance between each grid point and it's contiguous
% neighbours if dist_map output is defined
if nargout == 2

    % calculate x-index of each grid point in the spherical section, create
    % mask and remove singleton dimensions
    [mx, mx_ind] = max(ss, [], 1);
    mask = squeeze(mx ~= 0);
    mx_ind = squeeze(mx_ind).*mask;

    % double check there there is only one value of spherical section in
    % each matrix column
    if sum(mx(:)) ~= sum(ss(:))
        error('mean neighbour distance cannot be calculated uniquely due to overlapping points in the x-direction');
    end

    % calculate average distance to grid point neighbours in the flat case
    x_dist = repmat([1 0 1], 3, 1);
    y_dist = x_dist.';
    flat_dist = sqrt(x_dist.^2 + y_dist.^2);
    flat_dist = mean(flat_dist(:));

    % compute distance map
    dist_map = zeros(size(mx_ind));
    sz = size(mx_ind);
    for m = 1:sz(1)
        for n = 1:sz(2)

            % clear map
            local_heights = zeros(3, 3);

            % extract the height (x-distance) of the 8 neighbouring grid
            % points
            if m == 1 && n == 1
                local_heights(2:3, 2:3) = mx_ind(m:m + 1, n:n + 1);
            elseif m == sz(1) && n == sz(2)
                local_heights(1:2, 1:2) = mx_ind(m - 1:m, n - 1:n);
            elseif m == 1 && n == sz(2)
                local_heights(2:3, 1:2) = mx_ind(m:m + 1, n - 1:n);
            elseif m == sz(1) && n == 1
                local_heights(1:2, 2:3) = mx_ind(m - 1:m, n:n + 1);
            elseif m == 1
                local_heights(2:3, :) = mx_ind(m:m + 1, n - 1:n + 1);
            elseif m == sz(1)
                local_heights(1:2, :) = mx_ind(m - 1:m, n - 1:n + 1);
            elseif n == 1
                local_heights(:, 2:3) = mx_ind(m - 1:m + 1, n:n + 1);
            elseif n == sz(2)
                local_heights(:, 1:2) = mx_ind(m - 1:m + 1, n - 1:n);
            else
                local_heights = mx_ind(m - 1:m + 1, n - 1:n + 1);
            end

            % compute average variation from center
            local_heights_var = abs(local_heights - local_heights(2, 2));

            % threshold no neighbours
            local_heights_var(local_heights == 0) = 0;

            % calculate total distance from centre
            dist = sqrt(x_dist.^2 + y_dist.^2 + local_heights_var.^2);

            % average and store as a ratio
            dist_map(m, n) = 1 + (mean(dist(:)) - flat_dist)./flat_dist;

        end
    end

    % threshold out the non-transducer grid points
    dist_map(mask ~= 1) = 0;

    % flattened plot and distance values
    if plot_section
        figure;
        subplot(2, 1, 1);
        imagesc(mx_ind);
        axis image;
        colorbar;
        title('Height of Transducer Face');
        subplot(2, 1, 2);
        dist_map_plot = dist_map;
        dist_map_plot(dist_map_plot ~= 0) = dist_map_plot(dist_map_plot ~= 0) -1;
        imagesc(100*dist_map_plot);
        axis image;
        colorbar;
        title('Percentage Increase in Average Neighbour Distance Compared to a Flat Surface');
    end

end

% plot if required
if plot_section

    voxelPlot(double(ss));
    view(150, 20);

%     % create surface plot
%     figure;
%     sz = size(ss);
%     [x, y] = meshgrid(1:sz(3),1:sz(2));
%     surface_ind = find(mx ~= 0);
%     plot3(x(surface_ind), y(surface_ind), mx_ind(mx ~= 0), 'k.');
%     axis image;

end

repo: Aero-matlab


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%

% discrete time convolution of
% x[n]=[1,2],n=0,1 and h[n]=[2,1,1,1],n=0,1,2,3

%the signals
 kx=[0 1];
 x=[1 2];
 stem(kx,x);
 legend('x[k]')
 axis([-.1 3.1 -.1 2.1])

 figure
 kh=0:3;
 h=[2,1,1,1];
 stem(kh,h);
 axis([-.1 3.1 -.1 2.1])
 legend('h[k]')

 %reflection
 stem(-kh,h)
 axis([-3.1 .1 -.1 2.1])
 legend('h[-k]')

 %shifting
 %zero overlap
 figure
  stem(kx,x)
 axis([-5.1 3.1 -.1 2.1])
 legend('x[k]')
 figure
 n=-2;
 stem(-kh+n,h);
 axis([-5.1 3.1 -.1 2.1])
 legend('h[-2-k]')

  %overlap
 figure
 stem(kx,x)
 axis([-5.1 3.1 -.1 2.1])
 legend('x(k)')
 figure
 n=0;
 stem(-kh+n,h);
 axis([-5.1 3.1 -.1 2.1])
 legend('h(n-k)=h(0-k)')

 figure;
 n=1;
 stem(-kh+n,h);
 axis([-5.1 3.1 -.1 2.1])
 legend('h[n-k]=h[1-k]')

 figure;
 n=2;
 stem(-kh+n,h);
 axis([-5.1 3.1 -.1 2.1])
 legend('h[n-k]=h[2-k]')

 figure
 n=3;
 stem(-kh+n,h);
 axis([-2.1 6.1 -.1 2.1])
 legend('h[n-k]=h[3-k]')

 figure
 n=4;
 stem(-kh+n,h);
 axis([-2.1 6.1 -.1 2.1])
 legend('h[n-k]=h[4-k]')

%zero overlap
 n=5;
 stem(-kh+n,h);
 axis([-2.1 6.1 -.1 2.1])
 legend('h(n-k)=h(5-k)')

 %result
  n=0:4;
  y=[2,5,3,3,2];
  stem(n,y);
  axis([-.1 4.1 -.1 5.1])
  legend('y[n]')


text: plot(rho,h,'LineWidth',2,'Color','g'),grid,ylabel('Altitude: h [m]'),xlabel('Density: \rho [kg/m^3]')


text:
global CURRENT_PROJECT
projpath = 0;
while 1
    menulist = {'Keyboard input'
        'Quit'
        'Open help'
        'Choose/create project folder\n'};
    if any(projpath) && ~isempty(CURRENT_PROJECT)
        menulist(5) = {'Load specifications from Requirements.xls\n'};
        try %#ok<TRYNC>
            load ([projpath '\Requirements.mat'])
            menulist(6) = {'Iteration menu'};
        end
    end
    mainchoice = txtmenu(['Main Menu | Current Project: '....
        CURRENT_PROJECT],menulist);

    switch mainchoice
        case 0 % keyboard
            disp('Type ''return'' or hit F5 to exit keyboard mode and return to main menu')
            disp(' ')
            keyboard;

        case 1 %quit
            break

        case 2 %help
            open([CORE_PARENT_DIR '\CORE\core_help.html'])

        case 3 % Choose or create a project to work on
            [CURRENT_PROJECT projpath] = COREchooseproject();
            %sets current directory to project folder

        case 4 % Create mission parameters and objective functions
            try
                COREmkrequirements();
            catch ME
                disp('FAILURE: error creating requirements')
                disp(ME)
                beep;
            end

        case 5 % Iterate to find the best aircraft for the job
            CORE_ITERATE(Spec);

    end

end


text:
%% Test Problems for SEDUMI Interface
clc
clear

%Check SEDUMI is available
if(~optiSolver('sedumi',0))
    fprintf(2,'\n\nSEDUMI IS NOT AVAILABLE, NO SEDUMI TESTS CAN BE RUN\n\n');
    return;
end

files = {'sdp_truss1.mat','sdp_truss2.mat','sdp_truss3.mat','sdp_truss4.mat','sdp_hinf01.mat',...
         'sdp_hinf02.mat','sdp_hinf03.mat','sdp_hinf04.mat','sdp_arch0.mat','sdp_arch2.mat',...
         'sdp_arch4.mat','sdp_qap05.mat','sdp_qap06.mat','sdp_qap07.mat','sdp_qap08.mat',...
         'sdp_mcp100.mat','sdp_mcp124-1.mat','sdp_mcp124-2.mat','sdp_mcp124-3.mat','sdp_mcp124-4.mat',...
         'sdp_mcp250-1.mat','sdp_mcp250-2.mat','sdp_mcp250-3.mat','sdp_mcp250-4.mat'};

len = length(files);
cfval = zeros(len,1);
sfval = zeros(len,1);
optsC = optiset('solver','csdp','maxtime',15);
optsS = optiset('solver','sedumi','maxtime',15);

fprintf('SEDUMI SDP Checking\n');
%For each problem
for i = 1:len
    %Read In Problem
    prob = sdpRead(files{i});
    %Solve
    fprintf('Running Problem %s....',files{i});
    [~,cfval(i)] = solve(opti(prob,optsC));
    [~,sfval(i)] = solve(opti(prob,optsS));
    fprintf(' Done\n');
end

%Display fvals
display([cfval sfval])
if(norm(cfval-sfval) > 1e-1)
    error('Failed SeDuMi Checking');
else
    fprintf('Passed SeDuMi Checking\n');
end

% %% Complex Checking
% % Create a positive (semi)definite matrix:
% R=randn(5)+j*randn(5); R=R*R';
% % Make sure it's Hermitian:
% R=(R+R')/2;
%
% % Primal formulation
% clear At b c K
% c=vec(eye(5));
% At=vec(R);
% b=1;
% K.s=5;
% K.scomplex=1;
%
% [x,y,info]=sedumi(At,b,c,K);
%
% % Normed eigenvector:
% e=x(1:5)/norm(x(1:5));
% % Verify the eigenvalue
% [real(e'*R*e), max(eig(R))]

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function u = dirichlet_condition ( node_num, node_xy, time )

%*****************************************************************************80
%
%% DIRICHLET_CONDITION sets the value of a Dirichlet boundary condition.
%
%  Discussion:
%
%    For efficiency, this routine may be called once, with all the
%    nodes input.  Of course, only those nodes that lie on the boundary
%    will need a meaningful value of U to be assigned.  Other nodes
%    may be given a dummy value for U.
%
%    This routine is set for the unit square.
%
%    We assume that the equation to be solved is
%
%      dUdT - Laplacian U + K * U = F
%
%    with
%
%      K(X,Y,T) = 0
%      F(X,Y,T) = (2*pi*pi-1)*sin(pi*x)*sin(pi*y)*exp(-t).
%
%    The exact solution is:
%
%      U = sin(pi*x) * sin(pi*y) * exp(-t)
%
%  Modified:
%
%    07 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of the points.
%
%    Input, real TIME, the current time.
%
%    Output, real U(NODE_NUM), the value of the solution at the
%    the points.
%
  u(1:node_num) = sin ( pi * node_xy(1,1:node_num) ) ...
               .* sin ( pi * node_xy(2,1:node_num) ) * exp ( - time );

  return
end

repo: jburkardt-m