text:
function [positions, time] = tracker(video_path, img_files, pos, target_sz, ...
    padding, kernel, lambda1, lambda2, 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;
	yf = fft2(gaussian_shaped_labels(output_sigma, floor(window_sz / cell_size)));

	%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
    offset = [-target_sz(1) 0; 0 -target_sz(2); target_sz(1) 0; 0 target_sz(2)];

	for frame = 1:numel(img_files),
		%load image
		im = imread([video_path img_files{frame}]);
		if size(im,3) > 1,
			im = rgb2gray(im);
		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));

			response = real(ifft2(model_wf .* zf));  %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));

		kf = conj(xf) .* xf;

        kfn = zeros([size(xf) length(offset)]);
        for j=1:length(offset)
            %obtain a subwindow close to target for regression to 0
            patch = get_subwindow(im, pos+offset(j,:), window_sz);
            xfn = fft2(get_features(patch, features, cell_size, cos_window));
            kfn(:,:,j) = conj(xfn) .* xfn;
        end

        num = conj(xf) .* yf;
        den = kf + lambda1 + lambda2.*sum(kfn,3);
        wf = num ./ den;

		if frame == 1,  %first frame, train with a single image
			model_wf = wf;
		else
			%subsequent frames, interpolate model
			model_wf = (1 - interp_factor) * model_wf + interp_factor * wf;
		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: Context-Aware-CF-Tracking


text:
function varargout = drawLabels3d(varargin)
%DRAWLABELS3D Draw text labels at specified 3D positions.
%
%   drawLabels3d(X, Y, Z, LBL) draw labels LBL at position X and Y.
%   LBL can be either a string array, or a number array. In this case,
%   string are created by using sprintf function, with '%.2f' mask.
%
%   drawLabels3d(POS, LBL) draw labels LBL at position specified by POS,
%   where POS is a N-by-3 int array.
%
%   drawLabels3d(..., NUMBERS, FORMAT) create labels using sprintf function,
%   with the mask given by FORMAT (e. g. '%03d' or '5.3f'), and the
%   corresponding values.
%
%   See also
%     drawLabels
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/01/2019.
%

%   HISTORY
%   2019-01-31 write 3D version from drawLabels

%% Parse input arguments

% check if enough inputs are given
if isempty(varargin)
    error('wrong number of arguments in drawLabels3d');
end

% extract handle of axis to draw on
if isscalar(varargin{1}) && ishandle(varargin{1})
    ax = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

% process input parameters
var = varargin{1};
if size(var, 2) == 1
    % coordinates given as separate arguments
    if length(varargin) < 4
        error('wrong number of arguments in drawLabels');
    end
    px  = var;
    py  = varargin{2};
    pz  = varargin{3};
    lbl = varargin{4};
    varargin(1:4) = [];
else
    % parameters given as a packed array
    if length(varargin) < 2
        error('wrong number of arguments in drawLabels');
    end
    if size(var, 2) < 3
        error('Requires coordinates array to have at least three columns');
    end
    px  = var(:,1);
    py  = var(:,2);
    pz  = var(:,3);
    lbl = varargin{2};
    varargin(1:2) = [];
end

% parse format for displaying numeric values
format = '%.2f';
if ~isempty(varargin)
    format = varargin{1};
end
if size(format, 1) == 1 && size(px, 1) > 1
    format = repmat(format, size(px, 1), 1);
end

%% compute the strings that have to be displayed
labels = cell(length(px), 1);
if isnumeric(lbl)
    for i = 1:length(px)
        labels{i} = sprintf(format(i,:), lbl(i));
    end
elseif ischar(lbl)
    for i = 1:length(px)
        labels{i} = lbl(i,:);
    end
elseif iscell(lbl)
    labels = lbl;
end
labels = char(labels);

%% display the text
h = text(px, py, pz, labels, 'parent', ax);

%% format output
if nargout > 0
    varargout = {h};
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function jac = p12_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P12_JAC evaluates the jacobian for problem p12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac = zeros ( neqn, neqn );

  for i = 1 : neqn - 1
    jac(i,i) = - i;
  end

  for i = 2 : neqn
    jac(i,i-1) = i - 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function labels = loadMNISTLabels(filename)
%loadMNISTLabels returns a [number of MNIST images]x1 matrix containing
%the labels for the MNIST images

fp = fopen(filename, 'rb');
assert(fp ~= -1, ['Could not open ', filename, '']);

magic = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(magic == 2049, ['Bad magic number in ', filename, '']);

numLabels = fread(fp, 1, 'int32', 0, 'ieee-be');

labels = fread(fp, inf, 'unsigned char');

assert(size(labels,1) == numLabels, 'Mismatch in label count');

fclose(fp);

end

repo: MachineLearning


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.

function CstAlignWCSwithFace(mws)

WCS = invoke(mws,'WCS');
invoke(WCS,'AlignWCSWithSelected','Face');
end

repo: CST-MATLAB-API


text:
function [yd3] = ft32yd3(ft3)
% Convert volume from cubic feet to cubic yards.
% Chad Greene 2012
yd3 = ft3*0.037037037037;

author: Sable

repo: mcbench-benchmarks


text:
%Copyright (c) October,15 2008 by Varsha Hedau, UIUC.  All rights reserved.
function [vp p All_lines]=getVP(imdir,imagename,DO_DISPLAY,savedir)
% getVP Get a triplet of orthogonal vanishing points for an image.
%For details see [1] Varsha Hedau, Derek Hoiem, David Forsyth, Recovering the Spatial
%     Layout of Cluttered Rooms, in the Twelfth IEEE International Conference
%     on Computer Vision, 2009.

% INPUT:
%    image -imagename and image directory

%OUTPUT:
% vp - three orthognal vanishing points
% p - vote for each line and each vanishing point
% All_lines - detected line segments in the image [x1 x2 y1 y2 theta r]

vp=[];p=[];lines=[];
img=imread([imdir imagename]);
%img=imresize(img,500/size(img,1));
[h w k]=size(img);
grayIm=rgb2gray(img);
[All_lines] = getLargeConnectedEdges(grayIm,30);

% chucking out the lines near image boundaries imaging artifacts
inds = find(sum(double(All_lines(:,1:2)>10),2) & sum(double(All_lines(:,1:2)<w-10),2) & ...
    sum(double(All_lines(:,3:4)>10),2) & sum(double(All_lines(:,3:4)<h-10),2));
All_lines = All_lines(inds,:);
All_lines=[All_lines sqrt(((All_lines(:,1)-All_lines(:,2)).^2+(All_lines(:,3)-All_lines(:,4)).^2))];
maxl=max(All_lines(:,7));
imsize = size(grayIm);

%Computing intersections of all the lines
lines = All_lines;
Xpnts = ComputeIntersectionPoints(lines);
inds = find(~isnan(Xpnts(:,1)) & ~isnan(Xpnts(:,2)) & ...
    ~isinf(Xpnts(:,1)) & ~isinf(Xpnts(:,2)));
Xpnts = Xpnts(inds,:);

%Computing votes for every point from all lines
VoteArr = ComputeLinePtVote(lines,Xpnts);
Vote=sum(VoteArr,1);

%get the first point & remove the lines of this point
[vv ii]=sort(Vote,'descend');
vp(1:2)=Xpnts(ii(1),1:2);
Vote1 = VoteArr(:,ii(1));
active_lines = find((Vote1*maxl./All_lines(:,7))<0.8);
inactive_lines = find((Vote1*maxl./All_lines(:,7))>=0.8);
Vote1 = [Vote1(active_lines);Vote1(inactive_lines)];
lines = All_lines(active_lines,:);

%work with the remaining lines
Xpnts = ComputeIntersectionPoints(lines);
inds = find(~isnan(Xpnts(:,1)) & ~isnan(Xpnts(:,2)) & ...
    ~isinf(Xpnts(:,1)) & ~isinf(Xpnts(:,2)));
Xpnts = Xpnts(inds,:);
VoteArr = ComputeLinePtVote([lines;All_lines(inactive_lines,:)],Xpnts);
Vote=sum(VoteArr(1:size(lines,1),:),1);
[vv ii]=sort(Vote,'descend');
Vote = vv(:);
Xpnts=Xpnts(ii,:);
VoteArr = VoteArr(:,ii);
%Remove some of the points
[Xpnts,Vote,VoteArr] = RemoveRedundantPoints2(Xpnts,Vote,VoteArr,w,h);

% Vectorized orthogonality check
[pts2,pts1]=find(~triu(ones(length(Vote))));
npts=length(pts1);
orthochks=[];
for pt=1:100000:npts
    tempinds = [pt:min(pt+100000-1,npts)];
    temp_orthochks=chckothrogonalityvector(...
        ones(length(tempinds),1)*vp(1:2),...
        Xpnts(pts1(tempinds),:),...
        Xpnts(pts2(tempinds),:),w,h);
    orthochks = [orthochks;temp_orthochks(:)];
end
orthos = find(orthochks);
pts1 = pts1(orthos);
pts2 = pts2(orthos);
npts=length(pts1);

% Total vote computation for these points
totVote = zeros(npts,1);
for ln=1:length(Vote1)
    Votes = [Vote1(ln)*ones(npts,1) VoteArr(ln,pts1)' VoteArr(ln,pts2)'];
    Votes = max(Votes,[],2);
    totVote = totVote+Votes;
end
totVote = [pts1(:) pts2(:) totVote(:)];
%     lines = All_lines;

if size(totVote,1) > 0
    [vv ii]=sort(totVote(:,3),'descend');
    vp(3:4) = Xpnts(totVote(ii(1),1),:);
    vp(5:6) = Xpnts(totVote(ii(1),2),:);

    VoteArrTemp = ComputeLinePtVote(All_lines,[vp(1) vp(2);vp(3) vp(4);vp(5) vp(6)]);
    p=[VoteArrTemp.*maxl./repmat(All_lines(:,7),[1 3]) zeros(size(All_lines,1),1)];%4th vp is outliers
    ind=find(max(p(:,1:3),[],2)< 0.5);
    p(ind,4)=1;
    p=p./repmat(sum(p,2),[1 4]);
    %     [vv linemem] = max(VoteArrTemp,[],2);
    [vv linemem] = max(p,[],2);
    %Plot three vps
    if DO_DISPLAY
        figure(1000);plot(vp(1),vp(2),'r*');hold on;
        imagesc(img);hold on;
        plot(vp(1),vp(2),'r*');
        plot(vp(3),vp(4),'g*');
        plot(vp(5),vp(6),'b*');
        % linemem(vv==0) = 4;
        grp1=find(linemem==1);
        grp2=find(linemem==2);
        grp3=find(linemem==3);
        grp4=find(linemem==4);
        plot(All_lines(grp1, [1 2])', All_lines(grp1, [3 4])','r');
        plot(All_lines(grp2, [1 2])', All_lines(grp2, [3 4])','g');
        plot(All_lines(grp3, [1 2])', All_lines(grp3, [3 4])','b');
        plot(All_lines(grp4, [1 2])', All_lines(grp4, [3 4])','c');
        axis ij;axis equal;
        saveas(1000,[savedir imagename(1:end-3) 'fig']);
        close all

    end
    filename=fullfile(savedir,[imagename(1:end-4) '_vp.mat']);
    save(filename,'vp','p','VoteArrTemp','All_lines');
end
return

author: zouchuhang

repo: LayoutNet


text:
% run_simPressShaped.m
% Dana Goerzen and Jamie Near, McGill University 2021.
%
% USAGE:
% This script is run simply by editing the input parameters and then
% clicking "Run".
%
% DESCRIPTION:
% This script simulates a PRESS experiment with fully shaped refocusing
% pulses. Coherence order filtering is employed to only simulate desired signals
% This results in a 4x speed up compared to phase cycling (see deprecated run_simPressShaped_phCyc.m)
% Furthermore, simulations are run at various locations in space to account for the
% within-voxel spatial variation of the metabolite signal.  Summation
% across spatial positions is performed. The MATLAB parallel computing toolbox
% (parfor loop) was used to accelerate the simulations.  Acceleration
% is currently performed in the direction of the slice selective pulse along
% the x-direction, but this can be changed.  Up to a factor of 12 acceleration
% can be achieved using this approach.  To enable the use of the MATLAB
% parallel computing toolbox, initialize the multiple worked nodes using
% "matlabpool size X" where "X" is the number of available processing
% nodes.  If the parallel processing toolbox is not available, then replace
% the "parfor" loop with a "for" loop.
%
% INPUTS:
% To run this script, edit the parameters below as desired and then click
% "run":
% refocWaveform     = name of refocusing pulse waveform.
% refTp             = duration of refocusing pulses[ms]
% Bfield            = Magnetic field strength in [T]
% Npts              = number of spectral points
% sw                = spectral width [Hz]
% Bfield            = magnetic field strength [Tesla]
% lw                = linewidth of the output spectrum [Hz]
% thkX              = slice thickness of x refocusing pulse [cm]
% thkY              = slice thickness of y refocusing pulse [cm]
% fovX              = full simulation FOV in the x direction [cm]
% fovY              = full simulation FOV in the y direction [cm]
% nX                = number of spatial grid points to simulate in x-direction
% nY                = number of spatial grid points to simulate in y-direction
% taus              = vector of pulse sequence timings  [ms]
% spinSys           = spin system to simulate
%
% OUTPUTS:
% out_posxy         = Simulation results, spatially resolved.
% out               = Simulation results, summed over all space.

% ************INPUT PARAMETERS**********************************
refocWaveform='sampleRefocPulse.pta'; %name of refocusing pulse waveform.
refTp=5; %duration of refocusing pulses[ms]
Npts=2048; %number of spectral points
sw=2000; %spectral width [Hz]
Bfield=3; %magnetic field strength [Tesla]
lw=2; %linewidth of the output spectrum [Hz]
thkX=1.66; %slice thickness of x refocusing pulse [cm]
thkY=1.66; %slice thickness of y refocusing pulse [cm]
fovX=4; %size of the full simulation Field of View in the x-direction [cm]
fovY=4; %size of the full simulation Field of View in the y-direction [cm]
nX=16; %Number of grid points to simulate in the x-direction
nY=16; %Number of grid points to simulate in the y-direction
tau1=6; %TE1 for first spin echo [ms]
tau2=6; %TE2 for second spin echo [ms]
spinSys='Lac'; %spin system to simulate
flipAngle=180; %Refocusing pulse flip angle [degrees]
centreFreq=1.3; %Centre frequency of refocusing pulses [ppm]
% ************END OF INPUT PARAMETERS**********************************

%set up spatial grid
x=linspace(-fovX/2,fovX/2,nX); %X positions to simulate [cm]
y=linspace(-fovY/2,fovY/2,nY); %y positions to simulate [cm]

%Load RF waveform
refRF=io_loadRFwaveform(refocWaveform,'ref',0);

gamma=42577000; %gyromagnetic ratio

%Load spin systems
load spinSystems
sys=eval(['sys' spinSys]);

%Resample refocusing RF pulse from 400 pts to 100 pts to reduce
%computational workload
refRF=rf_resample(refRF,100);

Gx=(refRF.tbw/(refTp/1000))/(gamma*thkX/10000); %[G/cm]
Gy=(refRF.tbw/(refTp/1000))/(gamma*thkY/10000); %[G/cm]

%Initialize structures:
out_posxy=cell(length(x),length(y));
out=struct([]);

%loop through space: Don't forget to initialize the parallel processing
%toolbox workers using 'matlabpool open N' (for N workers, 12 max).

%for X=1:length(x);
parfor X=1:length(x)
    for Y=1:length(y)
                disp(['Executing X-position ' num2str(X) ' of ' num2str(length(x)) ', '...
                    'Y-position ' num2str(Y) ' of ' num2str(length(y))])
                out_posxy{X}{Y}=sim_press_shaped(Npts,sw,Bfield,lw,sys,tau1,tau2,...
                    refRF,refTp,x(X),y(Y),Gx,Gy,flipAngle,centreFreq);

        out=op_addScans(out,out_posxy{X}{Y});

    end %end of spatial loop (parfor) in y direction.
end %end of spatial loop (parfor) in x direction.

%For consistent scaling across different shaped simulations, we need to :
%1.  Scale down by the total number of simulations run (since these were
%    all added together.
numSims=(nX*nY);
out=op_ampScale(out,1/numSims);

%2.  Scale by the total size of the simulated region, relative to the size
%    of the voxel.
voxRatio=(thkX*thkY)/(fovX*fovY);
out=op_ampScale(out,1/voxRatio);

%Uncomment this part to view the results of the spatially resolved simulation:
 if length(x)>1 && length(y)>1
     ppmmin=1.0;
     ppmmax=1.6;
     sim_make2DSimPlot(out_posxy,ppmmin,ppmmax);
 else
     plot(out.ppm,out.specs);
 end
op_plotspec(out)

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function out = CO_trev(y,tau)
% CO_trev   Normalized nonlinear autocorrelation, trev function of a time series
%
% Calculates the trev function, a normalized nonlinear autocorrelation,
% mentioned in the documentation of the TSTOOL nonlinear time-series analysis
% package (available here: http://www.physik3.gwdg.de/tstool/).
%
% The quantity is often used as a nonlinearity statistic in surrogate data
% analysis, cf. "Surrogate time series", T. Schreiber and A. Schmitz, Physica D,
% 142(3-4) 346 (2000).
%
%---INPUTS:
%
% y, time series
%
% tau, time lag (can be 'ac' or 'mi' to set as the first zero-crossing of the
%       autocorrelation function, or the first minimum of the automutual
%       information function, respectively)
%
%---OUTPUTS:
% The raw trev expression, its magnitude, the numerator and its magnitude, and
% the denominator.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
%% Set defaults:
if nargin < 2 || isempty(tau)
    tau = 'ac';
end

%-------------------------------------------------------------------------------
% Can set the time lag, tau, to be 'ac' or 'mi':
if strcmp(tau,'ac')
    tau = CO_FirstCrossing(y,'ac',0,'discrete');
    % tau is first zero crossing of the autocorrelation function
elseif strcmp(tau,'mi')
    tau = CO_FirstMin(y,'mi');
    % tau is the first minimum of the automutual information function
end
if isnan(tau)
    error('No valid setting for time delay. (Is the time series too short?)');
end

%-------------------------------------------------------------------------------
% Compute trev quantities:

yn = y(1:end-tau);
yn1 = y(1+tau:end); % yn, tau steps ahead

% The trev expression used in TSTOOL:
out.raw = mean((yn1-yn).^3)/(mean((yn1-yn).^2))^(3/2);

% The magnitude
out.abs = abs(out.raw);

% The numerator
out.num = mean((yn1-yn).^3);
out.absnum = abs(out.num);

% The denominator
out.denom = (mean((yn1-yn).^2))^(3/2);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function L = gaussianNoiseExpectationLogLikelihood(noise, mu, varsigma, y)

% GAUSSIANNOISEEXPECTATIONLOGLIKELIHOOD Likelihood of the data under the GAUSSIAN noise model.
% FORMAT
% DESC returns the likelihoods for data points under the  Gaussian noise model.
% ARG noise : the noise structure for which the likelihood is required.
% ARG mu : input mean locations for the likelihood.
% ARG varSigma : input variance locations for the likelihood.
% ARG y : target locations for the likelihood.
%
% SEEALSO : gaussianNoiseParamInit, gaussianNoiseLogLikelihood, noiseLikelihood
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

N = size(y, 1);
D = size(y, 2);
varsigma = varsigma + noise.sigma2;
for i = 1:D
  mu(:, i) = mu(:, i) + noise.bias(i);
end
arg = (mu - y)./sqrt(varsigma);
L = (2*pi*varsigma).^(-1/2).*exp( - .5*arg.*arg);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function g = givens2(a,b)
%GIVENS2 find a Givens rotation.
% Example:
%   g = givens2(a,b)
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

if (b == 0)
    c = 1 ; s = 0 ;
elseif (abs (b) > abs (a))
    tau = -a/b ; s = 1 / sqrt (1+tau^2) ; c = s*tau ;
else
    tau = -b/a ; c = 1 / sqrt (1+tau^2) ; s = c*tau ;
end
g = [c -s ; s c] ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [err,Info, BRIKinfo] = BrikInfo (BrickName)
%
%   [err,Info] = BrikInfo (BrickName)
%
%Purpose:
%
%   returns some field values in the .HEAD files
%
%Input Parameters:
%   Brick filename
%  (also works for 1D files but Info has very limited info.
%   You should really avoid using BrikInfo on 1D file
%   because the 1D file is read in its entirety (and then cleared)
%   before Info can be returned. Better used BrikLoad directly or
%   Read_1D)
%
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Mucho Problems
%   Info is a structure with the following fields
%           Fieldnames in UPPER CASE correspond to the fields in AFNI
%           Field names in lower case correspond to some interpretation of the uppercase ones
%      .RootName : the brik name without .HEAD or .BRIK.
%      .TypeName: The interpretation of BRICK_TYPES (short, int, floats, etc...)
%      .TypeBytes : The corresponding byte size of BRICK_TYPES
%		 .ByteOrder: interpretation of BYTEORDER_STRING
%      .Orientation : interpretation of ORIENT_SPECIFIC, Orientation(:,1)' forms the three letter orientation code
%        The ith row of Orientation describes the orientation along that dimension.
%   see AFNI's README.attributes for a complete description of these fields
%   MANDATORY FIELDS AS SPECIFIED IN: ~cox/README.attributes (added April 6 2001)
%      .DATASET_RANK : ASK BOB. I think the first is always 3 (3D) and the other is the number of sub-bricks
%      .DATASET_DIMENSIONS : Number of voxels in each of the three dimensions
%      .TYPESTRING: Determines if the data set is of type Anat (underlay) or Func (overlay)
%      .SCENE_DATA : The three integer codes describing the dataset type
%		 .ORIENT_SPECIFIC : orintation code
%      .ORIGIN : The xyz coordinates of the center of the (0, 0, 0) voxel in the data set
%      .DELTA : the increment (in mm) to go from one voxel to the next (could be +ve or -ve depending on slices)
%      .TAXIS_NUMS: see readme file
%      .TAXIS_FLOATS
%      .TAXIS_OFFSETS
%   ALMOST MANDATORY FIELDS
%      .IDCODE_STRING
%      .IDCODE_DATE
%      .BYTEORDER_STRING : Byte order string
%      .BRICK_STATS : range of values in brick (min to max)
%         (Do not apply the scaling factor to these values)
%      .BRICK_TYPES : types of values in .BRIK
%      .BRICK_FLOAT_FACS : float factors to apply to bricks for recovering original values
%      .BRICK_LABS : The Sub-brick labels (~delimited)
%      .BRICK_STATAUX : Auxilliary Statistical Information
%      .STAT_AUX
%      .HISTORY_NOTE
%      .NOTES_COUNT
%      .NOTE_NUMBER_001
%      .TAGALIGN_MATVEC
%      .VOLREG_CENTER_OLD
%      .VOLREG_CENTER_BASE
%      .VOLREG_ROTPARENT_IDCODE
%      .VOLREG_ROTPARENT_NAME
%      .VOLREG_GRIDPARENT_IDCODE
%      .VOLREG_GRIDPARENT_NAME
%      .VOLREG_INPUT_IDCODE
%      .VOLREG_INPUT_NAME
%      .VOLREG_BASE_IDCODE
%      .VOLREG_BASE_NAME
%      .VOLREG_ROTCOM_NUM
%      .VOLREG_MATVEC_xxxxxx
%      .VOLREG_ROTCOM_xxxxxx
%      .IDCODE_ANAT_PARENT
%      .TO3D_ZPAD
%      .IDCODE_WARP_PARENT
%      .WARP_TYPE
%      .WARP_DATA
%      .MARKS_XYZ
%      .MARKS_LAB
%      .MARKS_HELP
%      .MARKS_FLAGS
%      .TAGSET_NUM
%      .TAGSET_FLOATS
%      .TAGSET_LABELS
%      .LABEL_1
%      .LABEL_2
%      .DATASET_NAME
%      .DATASET_KEYWORDS
%      .BRICK_KEYWORDS
%      .HISTORY_NOTE
%      .NOTES_COUNT
%      .NOTE_NUMBER_xxx
%      .TEMPLATE_SPACE
%      .ATLAS_LABEL_TABLE
%
%  To implement in the future,
%      VOLREG_MATVEC_xxxxxx , VOLREG_ROTCOM_xxxxxx
%
%  The following fields were added to support 1D file format.
%      .FileFormat: 'BRIK' ('1D' is allowed for 1D files but 1D files do not use BrikInfo)
%      .Extension_1D: The extension of the 1D filename
%  BrikLoad and WriteBrik now read and write 1D files
%
%  If the .HEAD contains a TEMPLATE_SPACE and ATLAS_LABEL_TABLE, these are added as well
%
%Key Terms:
%   The 1st, second and third dimensions refer to the dimensions the slices were entered into to3d
%More Info :
%    afni's README.attributes
%
%     Author : Ziad Saad
%     Date : Sun Oct 17 00:13:49 CDT 1999 , Major modification, April 06 2001

%Define the function name for easy referencing
FuncName = 'BrikInfo';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;
Info = [];

%check if this is a 1D file ..
is1D = 0;
[St, xtr] = Remove1DExtension(BrickName);
if (~isempty(xtr))
   is1D = 1;
end

if (is1D) % 1D land
   [err, V, Info] = Read_1D(BrickName, 1);
   if (err)
      ErrMessage = sprintf ('%s: Failed to read %s file', FuncName, BrickName);
      err = ErrEval(FuncName,'Err_1D file could not be read');
      return;
   end
   clear V; %biggest waste since 2002
   return;
end

	vtmp = findstr(BrickName,'.BRIK');
	if (~isempty(vtmp)) %remove .BRIK
		BrickName = BrickName(1:vtmp(1)-1);
	end

	vtmp = findstr(BrickName,'.HEAD');
	if (isempty(vtmp)) %add .HEAD
		BrickName = sprintf('%s.HEAD',BrickName);
	end

	if (exist(BrickName) ~= 2)
		err = ErrEval(FuncName,'Err_Could not find data set');
		return;
	end

	%store the name without the extension
		vtmp = findstr(BrickName,'.HEAD');
		Info.RootName = BrickName(1:vtmp-1);

	fidin=fopen(BrickName,'r');
	if (fidin < 0),	err = ErrEval(FuncName,'Err_Could not read .HEAD file'); return;	end
		BRIKinfo = fscanf(fidin,'%c');
	fclose(fidin);

%get subBrik info
	%Nx, Ny, Nz
		[err,Info.DATASET_DIMENSIONS] = BrikInfo_SectionValue (BRIKinfo, 'DATASET_DIMENSIONS');

	%DATASET_RANK
		[err, Info.DATASET_RANK] = BrikInfo_SectionValue(BRIKinfo, 'DATASET_RANK');

	%BRICK_TYPES
		[err, Info.BRICK_TYPES] = BrikInfo_SectionValue(BRIKinfo, 'BRICK_TYPES');

	%BRICK_STATS
		[err, Info.BRICK_STATS] = BrikInfo_SectionValue(BRIKinfo, 'BRICK_STATS');

	%BRICK_FLOAT_FACS
		[err, Info.BRICK_FLOAT_FACS] = BrikInfo_SectionValue(BRIKinfo, 'BRICK_FLOAT_FACS');

	%BYTEORDER_STRING
		[err, Info.BYTEORDER_STRING] = BrikInfo_SectionValue(BRIKinfo, 'BYTEORDER_STRING');

	%ORIENT_SPECIFIC
		[err, Info.ORIENT_SPECIFIC] = BrikInfo_SectionValue(BRIKinfo, 'ORIENT_SPECIFIC');

	%ORIGIN
		[err, Info.ORIGIN] = BrikInfo_SectionValue(BRIKinfo, 'ORIGIN');

	%DELTA
		[err, Info.DELTA] = BrikInfo_SectionValue(BRIKinfo, 'DELTA');

	%BRICK_LABS
		[err, Info.BRICK_LABS] = BrikInfo_SectionValue(BRIKinfo, 'BRICK_LABS');

	%BRICK_KEYWORDS
		[err, Info.BRICK_KEYWORDS] = BrikInfo_SectionValue(BRIKinfo, 'BRICK_KEYWORDS');

	%SCENE_DATA
		[err, Info.SCENE_DATA] = BrikInfo_SectionValue(BRIKinfo, 'SCENE_DATA');

	%TYPESTRING
		[err, Info.TYPESTRING] = BrikInfo_SectionValue(BRIKinfo, 'TYPESTRING');

	%TAXIS_NUMS
		[err, Info.TAXIS_NUMS] = BrikInfo_SectionValue(BRIKinfo, 'TAXIS_NUMS');

	%TAXIS_FLOATS
		[err, Info.TAXIS_FLOATS] = BrikInfo_SectionValue(BRIKinfo, 'TAXIS_FLOATS');

	%TAXIS_OFFSETS
		[err, Info.TAXIS_OFFSETS] = BrikInfo_SectionValue(BRIKinfo, 'TAXIS_OFFSETS');

	%IDCODE_STRING
		[err, Info.IDCODE_STRING] = BrikInfo_SectionValue(BRIKinfo, 'IDCODE_STRING');

	%IDCODE_DATE
		[err, Info.IDCODE_DATE] = BrikInfo_SectionValue(BRIKinfo, 'IDCODE_DATE');

	%BRICK_STATAUX
		[err, Info.BRICK_STATAUX] = BrikInfo_SectionValue(BRIKinfo, 'BRICK_STATAUX');

	%STAT_AUX
		[err, Info.STAT_AUX] = BrikInfo_SectionValue(BRIKinfo, 'STAT_AUX');

	%HISTORY_NOTE
		[err, Info.HISTORY_NOTE] = BrikInfo_SectionValue(BRIKinfo, 'HISTORY_NOTE');

	%NOTES_COUNT
		[err, Info.NOTES_COUNT] = BrikInfo_SectionValue(BRIKinfo, 'NOTES_COUNT');

	%NOTE_NUMBER_001
		[err, Info.NOTE_NUMBER_001] = BrikInfo_SectionValue(BRIKinfo, 'NOTE_NUMBER_001');

	%TAGALIGN_MATVEC
		[err, Info.TAGALIGN_MATVEC] = BrikInfo_SectionValue(BRIKinfo, 'TAGALIGN_MATVEC');

	%VOLREG_CENTER_OLD
		[err, Info.VOLREG_CENTER_OLD] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_CENTER_OLD');

	%VOLREG_CENTER_BASE
		[err, Info.VOLREG_CENTER_BASE] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_CENTER_BASE');

	%VOLREG_ROTPARENT_IDCODE
		[err, Info.VOLREG_ROTPARENT_IDCODE] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_ROTPARENT_IDCODE');

	%VOLREG_ROTPARENT_NAME
		[err, Info.VOLREG_ROTPARENT_NAME] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_ROTPARENT_NAME');

	%VOLREG_GRIDPARENT_IDCODE
		[err, Info.VOLREG_GRIDPARENT_IDCODE] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_GRIDPARENT_IDCODE');

	%VOLREG_GRIDPARENT_NAME
		[err, Info.VOLREG_GRIDPARENT_NAME] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_GRIDPARENT_NAME');

	%VOLREG_INPUT_IDCODE
		[err, Info.VOLREG_INPUT_IDCODE] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_INPUT_IDCODE');

	%VOLREG_INPUT_NAME
		[err, Info.VOLREG_INPUT_NAME] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_INPUT_NAME');

	%VOLREG_BASE_IDCODE
		[err, Info.VOLREG_BASE_IDCODE] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_BASE_IDCODE');

	%VOLREG_BASE_NAME
		[err, Info.VOLREG_BASE_NAME] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_BASE_NAME');

	%VOLREG_ROTCOM_NUM
		[err, Info.VOLREG_ROTCOM_NUM] = BrikInfo_SectionValue(BRIKinfo, 'VOLREG_ROTCOM_NUM');

   %Inbal's matrix,
   %TO get AFNI's RAI coords from AFNI i, j, k, indices (start at 0):
   % M = reshape(Info.IJK_TO_DICOM_REAL, 4, 3)'; I = [i j k 1]';
   % X = M*I
      [err, Info.IJK_TO_DICOM_REAL ] = BrikInfo_SectionValue(BRIKinfo, 'IJK_TO_DICOM_REAL');

	if (~err)
		for i=1:1:Info.VOLREG_ROTCOM_NUM
			sp = pad_strn(num2str(i-1), '0', 6, 1);
			spad = sprintf('[err, Info.VOLREG_MATVEC_%s] = BrikInfo_SectionValue(BRIKinfo, ''VOLREG_MATVEC_%s'');' , sp, sp);
			eval([ spad ]);
			spad = sprintf('[err, Info.VOLREG_ROTCOM_%s] = BrikInfo_SectionValue(BRIKinfo, ''VOLREG_ROTCOM_%s'');' , sp, sp);
			eval([ spad ]);
		end
	end

	%IDCODE_ANAT_PARENT
		[err, Info.IDCODE_ANAT_PARENT] = BrikInfo_SectionValue(BRIKinfo, 'IDCODE_ANAT_PARENT');

	%TO3D_ZPAD
		[err, Info.TO3D_ZPAD] = BrikInfo_SectionValue(BRIKinfo, 'TO3D_ZPAD');

	%IDCODE_WARP_PARENT
		[err, Info.IDCODE_WARP_PARENT] = BrikInfo_SectionValue(BRIKinfo, 'IDCODE_WARP_PARENT');

	%WARP_TYPE
		[err, Info.WARP_TYPE] = BrikInfo_SectionValue(BRIKinfo, 'WARP_TYPE');

	%WARP_DATA
		[err, Info.WARP_DATA] = BrikInfo_SectionValue(BRIKinfo, 'WARP_DATA');

	%MARKS_XYZ
		[err, Info.MARKS_XYZ] = BrikInfo_SectionValue(BRIKinfo, 'MARKS_XYZ');

	%MARKS_LAB
		[err, Info.MARKS_LAB] = BrikInfo_SectionValue(BRIKinfo, 'MARKS_LAB');

	%MARKS_HELP
		[err, Info.MARKS_HELP] = BrikInfo_SectionValue(BRIKinfo, 'MARKS_HELP');

	%MARKS_FLAGS
		[err, Info.MARKS_FLAGS] = BrikInfo_SectionValue(BRIKinfo, 'MARKS_FLAGS');

	%TAGSET_NUM
		[err, Info.TAGSET_NUM] = BrikInfo_SectionValue(BRIKinfo, 'TAGSET_NUM');

	%TAGSET_FLOATS
		[err, Info.TAGSET_FLOATS] = BrikInfo_SectionValue(BRIKinfo, 'TAGSET_FLOATS');

	%TAGSET_LABELS
		[err, Info.TAGSET_LABELS] = BrikInfo_SectionValue(BRIKinfo, 'TAGSET_LABELS');

	%LABEL_1
		[err, Info.LABEL_1] = BrikInfo_SectionValue(BRIKinfo, 'LABEL_1');

	%LABEL_2
		[err, Info.LABEL_2] = BrikInfo_SectionValue(BRIKinfo, 'LABEL_2');

	%DATASET_NAME
		[err, Info.DATASET_NAME] = BrikInfo_SectionValue(BRIKinfo, 'DATASET_NAME');

	%DATASET_KEYWORDS
		[err, Info.DATASET_KEYWORDS] = BrikInfo_SectionValue(BRIKinfo, 'DATASET_KEYWORDS');

	%BRICK_KEYWORDS
		[err, Info.BRICK_KEYWORDS] = BrikInfo_SectionValue(BRIKinfo, 'BRICK_KEYWORDS');

	%HISTORY_NOTE
		%[err, Info.HISTORY_NOTE] = BrikInfo_SectionValue(BRIKinfo, 'HISTORY_NOTE');

	%NOTES_COUNT
		%[err, Info.NOTES_COUNT] = BrikInfo_SectionValue(BRIKinfo, 'NOTES_COUNT');

	if (~err)
		for i=1:1:Info.NOTES_COUNT
			sp = pad_strn(num2str(i), '0', 3, 1);
			spad = sprintf('[err, Info.NOTE_NUMBER_%s] = BrikInfo_SectionValue(BRIKinfo, ''NOTE_NUMBER_%s'');' , sp, sp);
			eval([ spad ]);
		end
	end

	%SOMETHING
		%[err, Info.SOMETHING] = BrikInfo_SectionValue(BRIKinfo, 'SOMETHING');

	%DOF
		[err, Info.WORSLEY_DF] = BrikInfo_SectionValue(BRIKinfo, 'WORSLEY_DF');

	%NONJ
		[err, Info.WORSLEY_NCONJ] = BrikInfo_SectionValue(BRIKinfo, 'WORSLEY_NCONJ');

	%FWHM
		[err, Info.WORSLEY_FWHM] = BrikInfo_SectionValue(BRIKinfo, 'WORSLEY_FWHM');

        %TEMPLATE_SPACE
		[err, Info.TEMPLATE_SPACE] = BrikInfo_SectionValue(BRIKinfo, 'TEMPLATE_SPACE');

        %ATLAS_LABEL_TABLE (if present)
                [err, atlastable] = BrikInfo_SectionValue(BRIKinfo, 'ATLAS_LABEL_TABLE');

    try
      Info.ATLAS_LABEL_TABLE = parse_atlastable(atlastable);
    end

    itype = unique(Info.BRICK_TYPES);

 	if (length(itype) > 1)
		Info.TypeName = 'Mutliple Types';
	else
		switch itype
			case 0
				Info.TypeName = 'byte';
				Info.TypeBytes = 1;
			case 1
				Info.TypeName = 'short';
				Info.TypeBytes = 2; %Platform dependent
			case 2
				Info.TypeName = 'int';
				Info.TypeBytes = 4; %Platform dependent
			case 3
				Info.TypeName = 'float';
				Info.TypeBytes = 4; %Platform dependent
			otherwise
				Info.TypeName = 'Dunno';
				Info.TypeBytes = 0;
		end
	end

	 if (isempty(Info.BYTEORDER_STRING)) %field not found go native
	 	Info.ByteOrder = 'unspecified';
	 else
		 if (~isempty(strmatch('MSB_FIRST', Info.BYTEORDER_STRING)))
				Info.ByteOrder = 'ieee-be'; %Big Endian
		else
			if (~isempty(strmatch('LSB_FIRST', Info.BYTEORDER_STRING)))
					Info.ByteOrder = 'ieee-le'; %Little Endian
			else
				err = ErrEval(FuncName,'Err_Could not understand BYTEORDER_STRING');
				return;
			end
		end
	 end

	for i=1:1:3
		switch Info.ORIENT_SPECIFIC(i)
			case 0
				Info.Orientation(i,:) = 'RL';	%right to left
			case 1
				Info.Orientation(i,:) = 'LR';
			case 2
				Info.Orientation(i,:) = 'PA';
			case 3
				Info.Orientation(i,:) = 'AP';
			case 4
				Info.Orientation(i,:) = 'IS';
			case 5
				Info.Orientation(i,:) = 'SI';
      otherwise
				err = ErrEval(FuncName,'Err_Cannot understand Orientation code');
				return;
		end
	end

   Info.FileFormat = 'BRIK';
   Info.Extension_1D = '';
err = 0;
return;

function T = parse_atlastable(atlastable)

% subfunction to convert a string (read from an afni .HEAD file) with
% atlastable information into a struct

% assume that the string is well-behaved, and that it consists of a series
% of atlas points that has the attributes STRUCT, VAL, OKEY, GYoAR, COG
ap = regexp(atlastable, '<ATLAS_POINT');
st = regexp(atlastable, 'STRUCT');
v  = regexp(atlastable, 'VAL');
o  = regexp(atlastable, 'OKEY');
g  = regexp(atlastable, 'GYoAR');
c  = regexp(atlastable, 'COG');

% the sub-brick label is optional, only present for probabilistic atlases
s  = regexp(atlastable, 'SB_LABEL');
isprobabilistic = ~isempty(s);

assert(numel(ap)==numel(st));
assert(numel(ap)==numel(v));
assert(numel(ap)==numel(o));
assert(numel(ap)==numel(g));
assert(numel(ap)==numel(c));
if isprobabilistic
  assert(numel(ap)==numel(s));
end

ap(end+1) = numel(atlastable);

T = [];
for k = 1:numel(st)
  T(k,1).struct = getcontent(atlastable(st(k):(v(k)-1)));
  T(k,1).val    = str2double(getcontent(atlastable(v(k):(o(k)-1))));
  T(k,1).okey   = str2double(getcontent(atlastable(o(k):(g(k)-1))));
  T(k,1).gyoar  = str2double(getcontent(atlastable(g(k):(c(k)-1))));
  if isprobabilistic
    T(k,1).cog = str2num(getcontent(atlastable(c(k):(s(k)-1))));
    tmp        = getcontent(atlastable(s(k):(ap(k+1)-1)));
    % the sub-brick label can be numeric, or string
    if ~isempty(str2num(tmp))
      T(k,1).sb_label = str2num(tmp);
    else
      T(k,1).sb_label = tmp;
    end
  else
    T(k,1).cog = str2num(getcontent(atlastable(c(k):(ap(k+1)-1))));
  end
end

function out = getcontent(substring)

% anecdotally, I have come across 2 situations: one with the values between
% double quotation marks, or without, try to accommodate both situations
sel = regexp(substring, '"');
if isempty(sel)
  sel = regexp(substring, '=');
  out = deblank(substring((sel+1):end));
elseif numel(sel)>=2
  out = substring((sel(1)+1):(sel(2)-1));
else
  error('something unexpected happened');
end


text:
% test for vector field reorientation

n = 128/2;
h = compute_gaussian_filter([61 61],20/(2*n),[n n]);
Phi = randn(n);
Phi = perform_convolution(Phi, h);
v0 = grad(Phi);

% corrupt
v = v0 .* repmat( sign(randn(n)), [1 1 2] );

% reorient
options.method = 'randomized';
options.niter_reorient = 500;
options.method = 'propagation';
options.method = 'laplacian';
v1 = v;
v1 = perform_vf_reorientation(v1, options);

clf;
imageplot(v0, 'original', 1,3,1);
imageplot(v, 'noisy', 1,3,2);
imageplot(v1, 'recovered', 1,3,3);

author: gpeyre

repo: matlab-toolboxes


text:
function v = vectorize(m)

v = m(:);

author: kristinbranson

repo: JAABA


text:
function marginal = marginal_nodes(engine, nodes)
% MARGINAL_NODES Compute the marginal on the specified query nodes (likelihood_weighting)
% marginal = marginal_nodes(engine, nodes)

bnet = bnet_from_engine(engine);
ddom = myintersect(nodes, bnet.dnodes);
cdom = myintersect(nodes, bnet.cnodes);
nsamples = size(engine.samples, 1);
ns = bnet.node_sizes;

%w = normalise(engine.weights);
w = engine.weights;
if mysubset(nodes, ddom)
  T = 0*myones(ns(nodes));
  P = prod(ns(nodes));
  indices = ind2subv(ns(nodes), 1:P);
  samples = reshape(cat(1, engine.samples{:,nodes}), nsamples, length(nodes));
  for j = 1:P
    rows = find_rows(samples, indices(j,:));
    T(j) = sum(w(rows));
  end
  T = normalise(T);
  marginal.T = T;
elseif subset(nodes, cdom)
  samples = reshape(cat(1, engine.samples{:,nodes}), nsamples*sum(ns(nodes)), length(nodes));
  [marginal.mu, marginal.Sigma] =  wstats(samples', normalise(w));
else
  error('can''t handle mixed marginals yet');
end

marginal.domain = nodes;

%%%%%%%%%

function rows = find_rows(M, v)
% FINDROWS Find rows which are equal to a specified vector
% rows = findrows(M, v)
% Each row of M is a sample

temp = abs(M - repmat(v, size(M, 1), 1));
rows = find(sum(temp,2) == 0);

%%%%%%%%

function [mu, Sigma] = wstats(X, w)

% Computes the weighted mean and weighted covariance matrix for a given
% set of observations X(:,i), and a set of normalised weights w(i).
% Each column of X is a sample.

d = X - repmat(X * w', 1, size(X, 2));
mu = sum(X .* repmat(w, size(X, 1), 1), 2);
Sigma = d * diag(w) * d';

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [xvec, yvec]=generatepointsrandomboundCircle(npoints,diamcirc,origin, minsep)
% [xvec, yvec]=generatepointsrandomboundCircle(npoints,diamcirc,origin,
% minsep)
% generates points randomly distributed in the circle with 2r=diamcirc
% with center at origin. Minimum separation between points is minsep.
if ~exist('minsep','var')
    minsep = 0;
end

n=1;
xvec = 0;
yvec = 0;
while n<=npoints
    pos=origin+diamcirc*(rand(1,2)-.5);
    if sum((pos-origin).^2) <= (diamcirc/2)^2
        if or(n==1, min(pdist([[xvec'; pos(1)],[yvec'; pos(2)]], 'euclidean'))>minsep)
            xvec(n)=pos(1);
            yvec(n)=pos(2);
            n=n+1;
        end
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function bb = bb_click_move(img,bb)
% Change

[c,r,p] = impixel(bb_burn(img,bb));

bbW = bb_width(bb);
bbH = bb_height(bb);

bb = [c(1) - bbW/2; r(1) - bbH/2; c(1) + bbW/2; r(1) + bbH/2];

repo: OpenTLD


text:
function [B,Lambda,CLS,exitCode]=jointMatDiagLS(A,w,maxIter,CLSTol,RelTol,AbsTol)
%%JOINTMATDIAGLS Given K real or complex Hermitian or symmetric matrices in
%              A, try to find a common diagonalizing matrix B and sets of
%              values Lambda such that A(:,:,k)=B*diag(Lambda(:,k))*B' for
%              all k. An exact solution is only guaranteed for positive
%              (semi-)definite matrices with K=2. However, for general
%              problems the B and Lambda values are chosen trying to
%              minimize the least squares cost function
%              CLS=\sum_{k=1}^K w(k)*norm(A_k-B*diag(lambda_k)*B','fro')^2
%              when the matrices A are Hermitian, lambda_k is an NX1 vector
%              and
%              CLS=\sum_{k=1}^K w(k)*norm(A_k-B*diag(lambda_k)*B.','fro')^2
%              when the matrices A re symmetric. Real matrices count as
%              Hermitian; one cannot mix Hermitian and Symmetric complex
%              matrices in a single A.
%
%INPUTS: A An NXNXK set of K NXN Hermitian or symmetric matrices. This
%          cannot contain Inf or NaN values.
%        w An optional set of positive that go into the cost function.
%          These do not need to sum to any particular value. The default if
%          omitted or an empty matrix is passed is ones(K,1).
%  maxIter The maximum number of iterations to perform. The default if
%          omitted or an empty matrix is passed is 10000.
%   CLSTol If the least squares cost is <= this threshold, then the
%          algorithm will terminate early. This value is checked every 10
%          iterations. The default if omitted or an empty matrix is passed
%          is 1e-6.
% RelTol, AbsTol Absolute and relative tolerances on the parameters being
%          estimated. let theta=[B(:);Lambda(:)] and thetaPrev the estimate
%          from the previous step. Say that diff=abs(theta-thetaPrev).
%          Convergence is declared if
%          all((diff<=AbsTol)|diff<=RelTol*abs(theta))
%          The default values if omitted or empty matrices are passed are
%          1e-9 and 1e-12.
%
%OUTPUTS: B The NXN matrix that approximately turns the Lambda values into
%           A. These are the diagonalizer values. B is generally not
%           an orthonormal matrix.
%    Lambda The NXK set of the diagonals of the K diagonal matrices.
%       CLS The value of the cost function on termination.
%  exitCode A parameter indicating how the function terminated. Possible
%           values are:
%           -1 The RelTol and/or AbsTol criteria were satisfied.
%            0 The CLSTol criterion was satisfied.
%            1 The maximum number of iterations was reached.
%
%This function implements the algorithm of [1]. Though exact solutions are
%possible for pairs of positive semidefinite matrices, note that
%convergence is often faster for pairs of positive definite matrices than
%for pairs of positive semidefinite matrices, particularly when the
%matrices are complex.
%
%EXAMPLE 1:
%Here, we show that we can obtain an exact solution for a pair of positive
%definite real, symmetric (Hermitian) matrices where one is positive
%semi-definite.
% C1=[57,  7, 12, 17;
%      7, 47, 17, 22;
%     12, 17, 37, 27;
%     17, 22, 27, 27];
% [V,D,U] = svd(C1);
% D(2,2)=0;%Make uninformative.
% D(3,3)=0;%Make uninformative.
% C1=V*D*U';
% C2=[87, 25, 18, 31;
%     25, 63, 20, 17;
%     18, 20, 65, 29;
%     31, 17, 29, 65];
% A=zeros(4,4,2);
% A(:,:,1)=C1;
% A(:,:,2)=C2;
% [B,Lambda,CLS,exitCode]=jointMatDiagLS(A)
% %The CLS is verified as:
% CLSCheck=norm(A(:,:,1)-B*diag(Lambda(:,1))*B','fro')^2+norm(A(:,:,2)-B*diag(Lambda(:,2))*B','fro')^2
%One will see that convergence is by the CLS bound.
%
%EXAMPLE 2:
%Here, we have two positive definite complex Hermitian matrices. Again, an
%exact solution is found rapidly. If the additive identity matrices to C1
%and C2 are removed, then the matrices are no longer positive definite and
%only a least-squares approximation can be found. In such an instance, the
%algorithm is much slower.
% C1=[  9+  0*1i,   -65+  0*1i,  -11-153*1i,  -91-173*1i;
%     -65+  0*1i,    83+  0*1i,   54- 38*1i,   31+ 28*1i;
%     -11+153*1i,    54+ 38*1i,  130+  0*1i,   16- 47*1i;
%     -91+173*1i,    31- 28*1i,   16+ 47*1i,   22+  0*1i]+215*eye(4);
%
% C2=[-16+  0*1i,  -32- 56*1i,  -12-128*1i,   16+114*1i;
%     -32+ 56*1i,   79+  0*1i,  -87- 67*1i,  -48+ 51*1i;
%     -12+128*1i,  -87+ 67*1i,   76+  0*1i,   -7- 96*1i;
%      16-114*1i,  -48- 51*1i,   -7+ 96*1i, -147+  0*1i]+259*eye(4);
% A=zeros(4,4,2);
% A(:,:,1)=C1;
% A(:,:,2)=C2;
% [B,Lambda,CLS,exitCode]=jointMatDiagLS(A)
% %The CLS is verified as:
% CLSCheck=norm(A(:,:,1)-B*diag(Lambda(:,1))*B','fro')^2+norm(A(:,:,2)-B*diag(Lambda(:,2))*B','fro')^2
%Again, convergence is due to satisfying the cost function criterion.
%
%EXAMPLE 3:
%Though only approximations can be obtained for more than two matrices in
%general, here we show that if three matrices originate from a common set
%of basis vectors, then this function will find the correct
%diagonalization.
% [V,D]=eig(magic(4));
% C1=V*D*V';
% C2=V*diag([1;2;3;5])*V';
% C3=V*diag([12;1;0;8])*V';
% A=zeros(4,4,3);
% A(:,:,1)=C1;
% A(:,:,2)=C2;
% A(:,:,3)=C3;
% [B,Lambda,CLS,exitCode]=jointMatDiagLS(A)
%Convergence of the cost function is quickly achieved for the three
%matrices. Note, however, that B is not the same as V.
%
%REFERENCES:
%[1] A. Yeredor, "Non-orthogonal joint diagonalization in the least-squares
%    sense with application in blind source separation," IEEE Transactions
%    on Signal Processing, vol. 50, no. 7, pp. 1545-1553, Jul. 2002.
%
%July 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    K=size(A,3);
    N=size(A,1);

    if(nargin<6||isempty(AbsTol))
        AbsTol=1e-12;
    end

    if(nargin<5||isempty(RelTol))
        RelTol=1e-9;
    end

    if(nargin<4||isempty(CLSTol))
        CLSTol=1e-6;
    end

    if(nargin<3||isempty(maxIter))
        maxIter=10000;
    end

    if(nargin<2||isempty(w))
        w=ones(K,1);
    end

    %Check whether all (assume symmetric) matrices are Hermitian. If so,
    %then a different algorithm can be used.
    isHermitian=true;
    if(~isreal(A))
        for k=1:K
            if(any(any(A(:,:,k)~=A(:,:,k)')))
                isHermitian=false;
                break;
            end
        end
    end

    w=w(:);

    %Initial values. This doesn't necessarily help speed up the
    %convergence, but it will work even if the matrices in A are singular.
    B=eye(N,N);
    Lambda=ones(N,K);

    exitCode=1;
    if(isHermitian)
        for curIter=1:maxIter
            thetaPrev=[B(:);Lambda(:)];

            for iter=1:4
                for l=1:N
                    B=ACPhase(A,Lambda,w,B,l);
                end
            end
            Lambda=DCPhase(A,B);

            %Check for convergence.
            if(mod(curIter,10))
                %Cost criterion.
                CLS=0;
                for k=1:K
                    CLS=CLS+w(k)*norm(A(:,:,k)-B*diag(Lambda(:,k))*B','fro')^2;
                end
                if(CLS<=CLSTol)
                    exitCode=0;
                    return;
                end

                %Parameter change criterion.
                theta=[B(:);Lambda(:)];

                diff=abs(theta-thetaPrev);
                if(all((diff<=AbsTol)|diff<=RelTol*abs(theta)))
                    exitCode=-1;
                    break;
                end
            end
        end
    else%Assume it is symmetric, not Hermitian.
        for curIter=1:maxIter
            thetaPrev=[B(:);Lambda(:)];

            for iter=1:4
                for l=1:N
                    B=ACPhaseSym(A,Lambda,w,B,l);
                end
            end
            Lambda=DCPhaseSym(A,B);

            %Check for convergence.
            if(mod(curIter,10))
                %Cost criterion.
                CLS=0;
                for k=1:K
                    CLS=CLS+w(k)*norm(A(:,:,k)-B*diag(Lambda(:,k))*B.','fro')^2;
                end
                if(CLS<=CLSTol)
                    exitCode=0;
                    return;
                end

                %Parameter change criterion.
                theta=[B(:);Lambda(:)];

                diff=abs(theta-thetaPrev);
                if(all((diff<=AbsTol)|diff<=RelTol*abs(theta)))
                    exigtCode=-1;
                    break;
                end
            end
        end
    end

    %If the algorithm did not convergence, compute the current CLS if it should
    %be returned.
    if(nargout>2)
        CLS=0;
        if(isHermitian)
            for k=1:K
                CLS=CLS+w(k)*norm(A(:,:,k)-B*diag(Lambda(:,k))*B','fro')^2;
            end
        else
            for k=1:K
                CLS=CLS+w(k)*norm(A(:,:,k)-B*diag(Lambda(:,k))*B.','fro')^2;
            end
        end
    end
end

function B=ACPhase(A,Lambda,w,B,l)

    N=size(A,1);
    K=size(A,3);

    %Step 1
    P=zeros(N,N);
    for k=1:K
        innerSum=B*diag(Lambda(:,k))*B'-Lambda(l,k)*(B(:,l)*B(:,l)');

        P=P+w(k)*Lambda(l,k)*(A(:,:,k)-innerSum);
    end

    %Step 2
    [V,D]=eig(P);
    d=diag(D);
    [~,idx]=max(d);
    mu=d(idx);
    betaVec=V(:,idx);

    %The paper suggests using the first nonzero index of betaVec. We just
    %use the index having the largest magnitude value. in it.
    [~,nonZeroIdx]=max(abs(betaVec));

    %Force the first nonzero element of the unit eigenvector to be real and
    %positive.
    if(isreal(betaVec))
        if(betaVec(nonZeroIdx)<0)
            betaVec=-betaVec;
        end
    else
        betaVec=exp(-angle(betaVec(nonZeroIdx))*1i)*betaVec;
        %betaVec(nonZeroIdx) should be real and positive, but this line
        %deals avoids any possible finite precision issues.
        betaVec(nonZeroIdx)=abs(betaVec(nonZeroIdx));
    end

    if(mu<0)
        B(:,l)=zeros(N,1);
    else
        B(:,l)=betaVec*sqrt(mu)/sqrt(sum(w.*(Lambda(l,:).^2).'));
    end
end

function Lambda=DCPhase(A,B)

    N=size(A,1);
    K=size(A,3);

    %Step 1
    prodVal=B'*B;
    prodVal=conj(prodVal).*prodVal;
    %Use a pseudoinverse for poorly conditioned problems.
    if(rcond(prodVal)<=1e-14)
        G=pinv(prodVal);
    else
        G=inv(prodVal);
    end

    %Step 2
    Lambda=zeros(N,K);
    for k=1:K
        Lambda(:,k)=G*diag(B'*A(:,:,k)*B);
    end
end

function B=ACPhaseSym(A,Lambda,w,B,l)

    N=size(A,1);
    K=size(A,3);

    %Step 1
    P=zeros(N,N);
    for k=1:K
        innerSum=A(:,:,k);
        for n=1:N
            if(n~=l)
                innerSum=innerSum-Lambda(n,k)*(B(:,n)*B(:,n).');
            end
        end

        P=P+w(k)*Lambda(l,k)*conj(innerSum);
    end

    %Step 2
    [V,D]=eig([real(P), -imag(P);
              -imag(P), -real(P)]);

    %D should be real. This line is just in case some finite precision
    %issue in matlab makes it not real.
    D=real(D);

    d=diag(D);
    [~,idx]=max(d);
    mu=d(idx);
    xi=V(:,idx);

    %Step 3
    if(mu<0)
        B(:,l)=zeros(N,1);
    else
        gammaVal=xi(1:N);
        delta=xi((N+1):(2*N));

        B(:,l)=(gammaVal+1j*delta)*sqrt(mu)/sqrt(sum(w.*(abs(Lambda(l,:)).^2).'));
    end
end

function Lambda=DCPhaseSym(A,B)

    N=size(A,1);
    K=size(A,3);

    %Step 1
    prodVal=B'*B;
    prodVal=prodVal.*prodVal;
    %Use a pseudoinverse for poorly conditioned problems.
    if(rcond(prodVal)<=1e-14)
        G=pinv(prodVal);
    else
        G=inv(prodVal);
    end

    %Step 2
    Lambda=zeros(N,K);
    for k=1:K
        Lambda(:,k)=G*diag(B'*A(:,:,k)*conj(B));
    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:
classdef LinearSequence < IncrementalSequence

    methods (Access = public)

        function obj = LinearSequence(x0,x1,nSteps,initialValue,finalValue)
            obj.init(x0,x1,nSteps,initialValue,finalValue);
            obj.generateAlphaSequence();
        end

    end

    methods (Access = protected)

        function generateAlphaSequence(obj)
             obj.alpha = linspace(obj.x0,obj.x1,obj.nSteps);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% e = mohsst5_performance_rmsew(Yh, I)

function e = mohsst5_performance_rmsew(Yh, I)

load('mohsst5_data.mat', 'observations');

Y = observations;

if all(size(Yh) == [1727 1632])
  % The seas have been removed
  disp('Adding the land areas to the matrix.')
  sea = sum(~isnan(Y),2) > 0;
  F = nan(size(Y));
  F(sea,:) = Yh;
  Yh = F;
  Itmp = false(size(Y));
  Itmp(sea,:) = I;
  I = Itmp;
end

Y(~I) = nan;

e = mohsst5_rmsew(Y,Yh);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function sT = rotate(sT, rot)

sT.vertices = rotate(sT.vertices,rot);
sT.edges = rotate(sT.edges,rot);

end

author: mtex-toolbox

repo: mtex

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


text:
classdef TestInvert
    %TestInvert

    methods (Static)
        function test_nonsingular
            A = randn(20);
            [B,status] = cv.invert(A);
            assert(isscalar(status));
            assert(ismatrix(B) && isequal(size(B),size(A)));
            assert(norm(A*B - eye(20)) < 1e-9);
        end

        function test_singular
            A = [1 0 0; -2 0 0; 4 6 1];
            %B = inv(A);
            [B,status] = cv.invert(A);
        end

        function test_close_to_singular
            A = 5*eye(5) - ones(5);
            %B = inv(A);
            [B,status] = cv.invert(A);
        end

        function test_hilbert
            A = hilb(7);
            B = inv(A);
            B = cv.invert(A, 'Method','LU');
            B = cv.invert(A, 'Method','Cholesky');
        end

        function test_pseudoinv
            A = magic(8); A = A(:,1:6);
            B = pinv(A); c = 1./cond(A);
            [B,c] = cv.invert(A, 'Method','SVD');
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function legendre_2d_exactness ( a, b, n, x, y, w, t )

%*****************************************************************************80
%
%% LEGENDRE_2D_EXACTNESS: monomial exactness for the 2D Legendre integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(2), the lower limits of integration.
%
%    Input, real B(2), the upper limits of integration.
%
%    Input, integer N, the number of points in the rule.
%
%    Input, real X(N), Y(N), the quadrature points.
%
%    Input, real W(N), the quadrature weights.
%
%    Input, integer T, the maximum total degree.
%    0 <= T.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Quadrature rule for the 2D Legendre integral.\n' );
  fprintf ( 1, '  Number of points in rule is %d\n', n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   D   I       J          Relative Error\n' );

  for tt = 0 : t

    fprintf ( 1, '  %2d\n', tt );

    for j = 0 : tt

      i = tt - j;

      p(1) = i;
      p(2) = j;

      s = legendre_2d_monomial_integral ( a, b, p );

      v(1:n,1) = x(1:n,1) .^ p(1) .* y(1:n,1) .^ p(2);

      q = w' * v;

      if ( s == 0.0 )
        e = abs ( q );
      else
        e = abs ( q - s ) / abs ( s );
      end

      fprintf ( 1, '  %6d  %6d  %24.16g\n', p(1), p(2), e );

    end
  end

  return
end

repo: jburkardt-m


text:
function oM  = PatalaOrientationMapping(varargin)

warning('Starting with MTEX 5.1 "PatalaOrientationMapping" has been replaced by "PatalaColorKey".')
oM = PatalaColorKey(varargin{:});

end


text:
function p02_story ( )

%*****************************************************************************80
%
%% P02_STORY prints the "story" for problem p02.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    ETY Lee,
%    Choosing Nodes in Parametric Curve Interpolation,
%    Computer-Aided Design,
%    Volume 21, Number 6, July/August 1989, pages 363-370.
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This example is due to ETY Lee of Boeing.\n' );
  fprintf ( 1, '  Data near the corners is more dense than in regions of small curvature.\n' );
  fprintf ( 1, '  A local interpolation method will produce a more plausible\n' );
  fprintf ( 1, '  interpolant than a nonlocal interpolation method, such as\n' );
  fprintf ( 1, '  cubic splines.\n' );

  return
end


text:
function c = block_fwt( f, w, J)
%BLOCK_FWT FWT func. wrapper for a block processing
%   Usage: c = block_fwt( f, w, J);
%
%   Input parameters:
%         f     : Input data.
%         w     : Analysis Wavelet Filterbank.
%         J     : Number of filterbank iterations.
%
%   Output parameters:
%         c      : Coefficient vector.
%
%   `c = block_fwt(f,w,J)` accepts suitably extended block of data *f*
%   and produces correct coefficients using the SegDWT algorithm (based on
%   overlap-save block convolution) with wavelet filters defined by *w*
%   and *J* levels. *f* is expected to be a column vector or a matrix and
%   the processing is done column-wise.
%
%   Do not call this function directly. The function is called from
%   |blockana| when used with frame type 'fwt' and 'segola' block transform
%   handling see |blockframeaccel|.
%
%   Function should be independent of block_interface.
%
%   See also: block, block_ifwt
%
%   References: ltfatnote026

if nargin<3
  error('%s: Too few input parameters.',upper(mfilename));
end;

% Initialize the wavelet filters structure
%h = fwtinit(h,'ana');

if any(w.a~=w.a(1))
   error('%s: Non-equal subsampling factors are not supported.',upper(mfilename));
end

w = fwtinit(w);
% Extended block length
Ls = size(f,1);
% Low-pass filter length
m = numel(w.h{1}.h);
% Low-pass subsampling factor
a = w.a(1);
% Extension length
rred = (a^J-1)/(a-1)*(m-a);
% Block boundaries
blocksize=w.a(1)^J;
% Input signal samples to be processed

% This is effectivelly the "negative" right extension described in chapter
% 4.1.4 in the reference.
L=rred+floor((Ls-rred)/blocksize)*blocksize;

levelLen = L;
filtNo = length(w.h);
subbNo = (filtNo-1)*J+1;
Lc = zeros(subbNo,1);
runPtr = 0;
for jj=1:J
   for ff=filtNo:-1:2
      Lc(end-runPtr) = floor((levelLen-m-1)/w.a(ff));
      runPtr = runPtr + 1;
   end
   levelLen = floor((levelLen-m-1)/w.a(1));
end
Lc(1)=levelLen;

%
%[Lc, L] = fwtclength(Ls,h,J,'valid');

% Crop to the right length
if(Ls>L)
   f=postpad(f,L);
end

if Ls<rred+a^J
   error('%s: Insufficient input signal length for the %s flag. Minimum is %i.',upper(mfilename),'''valid''',rred+a^J);
end

c = comp_fwt(f,w.h,w.a,J,'valid');

% Do the cropping
runPtr = 0;
for jj=1:J-1
   for ff=filtNo:-1:2
      cstart = (a^(J-jj)-1)/(a-1)*(m-a);
      c{end-runPtr} = c{end-runPtr}(cstart+1:end,:);
      runPtr = runPtr + 1;
   end
end

% To the pack format
c = cell2mat(c);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function y=rhartley(x,n)
%RHARTLEY Calculate the Hartley transform of real data Y=(X,N)
% Data is truncated/padded to length N if specified.
% The inverse transformation is x=hartley(y,n)/n

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: rhartley.m,v 1.5 2008/05/22 20:00:28 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 2
  y=fft(real(x));
else
  y=fft(real(x),n);
end
y=real(y)-imag(y);

author: decouples


text:
function ax = axis_pct(pct)
% AXIS_PCT       Set reasonable axis limits.
% AXIS_PCT(pct) sets axis limits to extend pct% beyond limits of plotted
% objects.  Default is 5%.
% Works for linear or log scale.
% Unfortunately, the axes won't change when new points are plotted.

if nargin < 1
  pct = 0.05;
end
ax = [Inf -Inf Inf -Inf Inf -Inf];

% find bounding box of plotted objects
children = get(gca,'children');
for child = children'
  if strcmp(get(child,'type'),'text')
    xyz = get(child,'position');
    % need to determine bounding box of the text
    c([1 2]) = xyz(1);
    c([3 4]) = xyz(2);
    c([5 6]) = xyz(3);
  else
    x = get(child,'xdata');
    c(1) = min(x);
    c(2) = max(x);
    y = get(child,'ydata');
    c(3) = min(y);
    c(4) = max(y);
    z = get(child,'zdata');
    if isempty(z)
      c([5 6]) = 0;
    else
      c(5) = min(z);
      c(6) = max(z);
    end
  end
  ax([1 3 5]) = min(ax([1 3 5]), c([1 3 5]));
  ax([2 4 6]) = max(ax([2 4 6]), c([2 4 6]));
end
if strcmp(get(gca,'xscale'), 'log')
  ax([1 2]) = log(ax([1 2]));
end
if strcmp(get(gca,'yscale'), 'log')
  ax([3 4]) = log(ax([3 4]));
end
dx = ax(2)-ax(1);
if dx == 0
  dx = 1;
end
dy = ax(4)-ax(3);
if dy == 0
  dy = 1;
end
dz = ax(6)-ax(5);
if dz == 0
  dz = 1;
end
ax = ax + [-dx dx -dy dy -dz dz]*pct;
if strcmp(get(gca,'xscale'), 'log')
  ax([1 2]) = exp(ax([1 2]));
end
if strcmp(get(gca,'yscale'), 'log')
  ax([3 4]) = exp(ax([3 4]));
end
% clip for 2D
ax = ax(1:length(axis));
axis(ax);
if nargout < 1
  clear ax
end


text:
function vdisp(s,varargin)

global prevCharCnt;
if isempty(prevCharCnt), prevCharCnt = 0; end

if ~check_option(varargin,'silent') && ~getMTEXpref('generatingHelpMode')
  fprintf(repmat('\b',1, prevCharCnt));
  prevCharCnt = 0;
  disp(s);
end


text:
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% S_ESCAPE_TEX de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, the string to be de-escaped.
%
%    Output, string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' || ...
         s1(s1_pos) == '_' || ...
         s1(s1_pos) == '^' || ...
         s1(s1_pos) == '{' || ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  end

  return
end

author: johannesgerer


text:
function prtUtilSplitLegend(axisH,symbolH,labelH)
% prtUtilSplitLegend split the current legend into two parts on either
%   side of the axes
%
% prtUtilSplitLegend with no input arguments splits the current legend in half
%
% prtUtilSplitLegend(axisHandle,lineHandles,strings) makes a legend on the current
% axes for the specified lines, with the specified strings.
%
% % Examples:
% close all;
% h = plot(rand(5,40));
% strings = prtUtilCellPrintf('%d',num2cell(1:40));
% legend(h,strings);
% prtUtilSplitLegend
%
% % Or:
% h = plot(rand(5,40));
% strings = prtUtilCellPrintf('%d',num2cell(1:40));
% prtUtilSplitLegend(gca,h,strings);
%

if nargin == 0
    axisH = gca;
    [legendH,~,symbolH] = legend;
    labelH = get(legendH,'String');
end

border = 0.05;        % Normalized
innerBorder = 0;

% Test legend width
midPt = round(length(symbolH)/2);
lgEastH = legend(gca,symbolH(midPt+1:end),labelH(midPt+1:end),'Location','EastOutside');
lgEastPos = get(lgEastH,'Position');
legendWidth = lgEastPos(3);

% Set axis location
legend(axisH,'off')
currentAxisPos = get(axisH,'OuterPosition');
newAxisPos = [...
    (border+innerBorder+legendWidth),...
    currentAxisPos(2),...
    (1 -(2*border)- (2*innerBorder)-(2*legendWidth)),...
    currentAxisPos(4)];
set(axisH,'OuterPosition',newAxisPos);

% Set up east legend
eastLegendPos = [...
    ((2*innerBorder)+border+legendWidth+newAxisPos(3)),...
    lgEastPos(2),...
    legendWidth,...
    lgEastPos(4)];
lgEastH = legend(gca,symbolH(midPt+1:end),labelH(midPt+1:end),'Location','EastOutside');
set(lgEastH,'Position',eastLegendPos);

% Set up west legend
extraAxisH = axes('Position',get(gca,'Position'),'Visible','off');
lgWestH = legend(extraAxisH,symbolH(1:midPt),labelH(1:midPt),'Location','WestOutside');
lgWestPos = get(lgWestH,'Position');
westLegendPos = [...
    border,...
    lgWestPos(2),...
    legendWidth,...
    lgWestPos(4)];
set(lgWestH,'Position',westLegendPos)

axes(axisH);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function n = nfold(cs,axis)
% maximal n-fold of symmetry axes

if nargin == 1
  switch cs.LaueName
    case {'112/m','2/m11','12/m1','mmm'}
      n = 2;
    case {'m-3','-3','-3m1','-31m'}
      n = 3;
    case {'4/m','4/mmm','m-3m'}
      n = 4;
    case {'6/m','6/mmm'}
      n = 6;
    otherwise
      n = 1;
  end
else
  axis = vector3d(axis);
  n = ones(size(axis));
  for i = 1:length(axis)
    ind = isnull(angle(cs.rot.axis,axis(i))) & cs.rot.angle>0;
    if any(ind(:))
      n(i) = 2*pi / min(cs.rot(ind).angle);
    end
  end

end


text:
%% manage paths

path = strsplit(pwd,"sim_vehicle_dynamics");
path2module = strcat(path(1), "sim_vehicle_dynamics");
path2frictionmodel_folder = strcat(path2module, "\variable_friction\frictionmodel");
path2input = strcat(path2frictionmodel_folder, "\GUI\inputs");
path2output = strcat(path2frictionmodel_folder, "\GUI\outputs");

filepath2output_sldd=strcat(path2output, "\frictionmaps.sldd");
filepath2output_config=strcat(path2output, "\Gridmap_Config.txt");

%Handover Referenceline path and resolution from base workspace
Grid_mode=evalin('base', 'Grid_mode');
InterpTime=evalin('base','InterpTime');

%Set Variables

fid=fopen(filepath2output_config, 'wt');
Gridmap_Ident=zeros(3,3);

%% generate data dictionaries containing friction maps

switch Grid_mode

    case 1
        filename_1 = dir(strcat(path2input, '\1*.mat'));
    	filename_1 = filename_1.name;

        Gridmap_1 = importdata(strcat(path2input, '\', filename_1));
        Gridmap_2 = int8(zeros(size(Gridmap_1,1),size(Gridmap_1,2)));
        Gridmap_3 = Gridmap_2;
        Gridmap_4 = Gridmap_2;

        fprintf(fid, ['Gridmap_1: ', filename_1, ...
            '\nGridmap_2: 0 \nGridmap_3: 0 \nGridmap_4: 0 ']);

        fclose(fid);
        Gridmap_Ident(1:2,1:2)  = Gridmap_1(1:2,1:2);
        Gridmap_Ident(1,3)      = Gridmap_1(1,3);
        Gridmap_Ident(3,1)      = size(Gridmap_1,1);
        Gridmap_Ident(3,2)      = size(Gridmap_1,2);

    case 2
        filename_1 = dir(strcat(path2input, '\1*.mat'));
        filename_2 = dir(strcat(path2input, '\2*.mat'));
     	filename_1 = filename_1.name;
        filename_2 = filename_2.name;

        Gridmap_1 = importdata(strcat(path2input, '\', filename_1));
        Gridmap_2 = importdata(strcat(path2input, '\', filename_2));
        Gridmap_3 = int8(zeros(size(Gridmap_1,1),size(Gridmap_1,2)));
        Gridmap_4 = Gridmap_3;

        fprintf(fid, ['Gridmap_1: ', filename_1, ...
            '\nGridmap_2: ', filename_2, ...
            '\nGridmap_3: 0 \nGridmap_4: 0 ']);

        fclose(fid);
        Gridmap_Ident(1:2,1:2)  = Gridmap_1(1:2,1:2);
        Gridmap_Ident(1,3)      = Gridmap_1(1,3);
        Gridmap_Ident(3,1)      = size(Gridmap_1,1);
        Gridmap_Ident(3,2)      = size(Gridmap_1,2);

    case 3
        filename_1 = dir(strcat(path2input, '\1*.mat'));
        filename_2 = dir(strcat(path2input, '\2*.mat'));
     	filename_3 = dir(strcat(path2input, '\3*.mat'));
     	filename_1 = filename_1.name;
        filename_2 = filename_2.name;
    	filename_3 = filename_3.name;

        Gridmap_1 = importdata(strcat(path2input, '\', filename_1));
        Gridmap_2 = importdata(strcat(path2input, '\', filename_2));
        Gridmap_3 = importdata(strcat(path2input, '\', filename_3));
        Gridmap_4 = int8(zeros(size(Gridmap_1,1),size(Gridmap_1,2)));

        fprintf(fid, ['Gridmap_1: ', filename_1, ...
            '\nGridmap_2: ', filename_2, ...
            '\nGridmap_3: ', filename_3, ...
            '\nGridmap_4: 0 ']);

        fclose(fid);
        Gridmap_Ident(1:2,1:2)  = Gridmap_1(1:2,1:2);
        Gridmap_Ident(1,3)      = Gridmap_1(1,3);
        Gridmap_Ident(3,1)      = size(Gridmap_1,1);
        Gridmap_Ident(3,2)      = size(Gridmap_1,2);

    case 4
        filename_1 = dir(strcat(path2input, '\1*.mat'));
        filename_2 = dir(strcat(path2input, '\2*.mat'));
     	filename_3 = dir(strcat(path2input, '\3*.mat'));
        filename_4 = dir(strcat(path2input, '\4*.mat'));
     	filename_1 = filename_1.name;
        filename_2 = filename_2.name;
    	filename_3 = filename_3.name;
      	filename_4 = filename_4.name;

        Gridmap_1 = importdata(strcat(path2input, '\', filename_1));
        Gridmap_2 = importdata(strcat(path2input, '\', filename_2));
        Gridmap_3 = importdata(strcat(path2input, '\', filename_3));
        Gridmap_4 = importdata(strcat(path2input, '\', filename_4));

        fprintf(fid, ['Gridmap_1: ', filename_1, ...
            '\nGridmap_2: ',filename_2, ...
            '\nGridmap_3: ',filename_3, ...
            '\nGridmap_4: ',filename_4]);

        Gridmap_Ident(1:2,1:2)  = Gridmap_1(1:2,1:2);
        Gridmap_Ident(1,3)      = Gridmap_1(1,3);
        Gridmap_Ident(3,1)      = size(Gridmap_1,1);
        Gridmap_Ident(3,2)      = size(Gridmap_1,2);

end

%Create Datadict
if exist(filepath2output_sldd, 'file')==0
    myDictionaryObj = Simulink.data.dictionary.create(filepath2output_sldd);
    dDataSectObj = getSection(myDictionaryObj,'Design Data');
    addEntry(dDataSectObj,'Grid_map_1',Gridmap_1);
    addEntry(dDataSectObj,'Grid_map_2',Gridmap_2);
    addEntry(dDataSectObj,'Grid_map_3',Gridmap_3);
    addEntry(dDataSectObj,'Grid_map_4',Gridmap_4);
    addEntry(dDataSectObj,'InterpTime',InterpTime);
    addEntry(dDataSectObj,'Gridmap_Ident',Gridmap_Ident);

else
    myDictionaryObj = Simulink.data.dictionary.open(filepath2output_sldd);
    dDataSectObj = getSection(myDictionaryObj,'Design Data');
    assignin(dDataSectObj,'Grid_map_1',Gridmap_1);
    assignin(dDataSectObj,'Grid_map_2',Gridmap_2);
    assignin(dDataSectObj,'Grid_map_3',Gridmap_3);
    assignin(dDataSectObj,'Grid_map_4',Gridmap_4);
    assignin(dDataSectObj,'InterpTime',InterpTime);
    assignin(dDataSectObj,'Gridmap_Ident',Gridmap_Ident);

end

%Save Changes
saveChanges(myDictionaryObj);

%Finish
clc
"Datadict generated successfully! It is located in \frictionmodel\GUI\outputs" %#ok<NOPTS>

author: TUMFTM


text:
function [pulseCleanedTemplate, pulseTemplate] = tapas_physio_get_template_from_pulses(...
    c, cpulse, halfTemplateWidthInSamples, verbose, varargin)
% Computes mean pulse template given time stamp of detected pulses and raw
% time series; removes outlier pulses with correlation to initial guess of
% mean template
%
%   [pulseCleanedTemplate, pulseTemplate] = tapas_physio_get_template_from_pulses(...
%                c, cpulse2ndGuess, halfTemplateWidthInSamples, ...
%                doNormalizeTemplate, dt)
% IN
%   c                   cardiac time series (raw physiological recording)
%   cpulse              index vector (wrt c) of detected pulse events (1= event onset/max)
%   halfTemplateWidthInSamples
%                       half the width in samples of template to be
%                       generated, i.e. numel(pulseTemplate) = 2*this+1
%   verbose             verbose substructure of physio, holds all figure
%                       handles created and verbose.level to trigger visualization
%
%   optional as propertyName/Value pairs:
%
%   thresholdHighQualityCorrelation
%                           default: 0.95
%                           outliers pulses below that correlation will be
%                           removed for averaging when retrieving final template
%   minFractionIncludedPulses
%                           default: 0.1
%                           defines minimum fraction of total number of
%                           found pulses to be included in final template
%                           (after removing outliers). If there are not
%                           enough clean pulses, this value is enforced by
%                           lowering quality thresholds.
%   doNormalizeTemplate     default: true; if true, subtracts mean and
%                           scales max to 1 for each pulse, before
%                           averaging
%   dt                      default: 1; defines time vector only used for
%                           output plot
%
%
% OUT
%   pulseCleanedTemplate    pulseTemplate, mean after removing outliers for
%                           correlation quality
%   pulseTemplate           average of all detected pulse shapes
% EXAMPLE
%   tapas_physio_get_template_from_pulses
%
%   See also tapas_physio_get_cardiac_pulse_template tapas_physio_get_cardiac_pulses_auto_matched

% Author:   Lars Kasper
% Created:  2019-01-29
% Copyright (C) 2019 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under
% the terms of the GNU General Public License (GPL), version 3. You can
% redistribute it and/or modify it under the terms of the GPL (either
% version 3 or, at your option, any later version). For further details,
% see the file COPYING or <http://www.gnu.org/licenses/>.

% z-transform to have all templates (for averaging) have
% same norm & be mean-corrected
defaults.doNormalizeTemplate = true;
defaults.minFractionIncludedPulses = 0.1;
defaults.thresholdHighQualityCorrelation = 0.95;
defaults.dt = 1; % for visualization only

args = tapas_physio_propval(varargin, defaults);
tapas_physio_strip_fields(args);

doDebug = verbose.level >= 3;

nSamplesTemplate = halfTemplateWidthInSamples * 2 + 1;
nPulses = numel(cpulse);
template = zeros(nPulses-3,nSamplesTemplate);

for n=2:nPulses-2
    startTemplate = cpulse(n)-halfTemplateWidthInSamples;
    endTemplate = cpulse(n)+halfTemplateWidthInSamples;

    template(n,:) = c(startTemplate:endTemplate);

    if doNormalizeTemplate
        template(n,:) = template(n,:) - mean(template(n,:),2);

        % std-normalized...
        %template(n,:) = template(n,:)./std(template(n,:),0,2);
        % max-norm:
        template(n,:) = template(n,:)./max(abs(template(n,:)));
    end

end

%delete first zero-elements of the template
template(1,:) = [];

% template as average of the found representative waves
pulseTemplate = mean(template);

if doDebug
    fh = verbose.fig_handles(end);
    figure(fh);
    subplot(3,1,3);
    tTemplate = dt*(0:2*halfTemplateWidthInSamples);
    plot(tTemplate, template');
    hold all;
    hp(1) = plot(tTemplate, pulseTemplate', '.--g', 'LineWidth', 3, 'Marker', ...
        'o');
    xlabel('t (seconds)');
    title('Templates of cycle time course and mean template');
end

%% Final template for peak search from best-matching templates
% delete the peaks deviating from the mean too
% much before building the final template
[~, pulseTemplate] = tapas_physio_corrcoef12(pulseTemplate, pulseTemplate);
isZtransformed = [0 1];

nTemplates = size(template,1);
similarityToTemplate = zeros(nTemplates,1);
for n=1:nTemplates
    similarityToTemplate(n) = tapas_physio_corrcoef12(template(n,:),pulseTemplate, ...
        isZtransformed);
end

% minimal number of high quality templates to be achieved, otherwise
% enforced
nMinHighQualityTemplates = ceil(minFractionIncludedPulses * nPulses);
indHighQualityTemplates = find(similarityToTemplate > ...
    thresholdHighQualityCorrelation);

% if threshold to restrictive, try with new one:
% best nMinHighQualityTemplates / nPulses of all found templates used for
% averaging
if numel(indHighQualityTemplates) < nMinHighQualityTemplates
    thresholdHighQualityCorrelation = tapas_physio_prctile(similarityToTemplate, ...
        1 - nMinHighQualityTemplates/nPulses);
    indHighQualityTemplates = find(similarityToTemplate > ...
        thresholdHighQualityCorrelation);
end
pulseCleanedTemplate = mean(template(indHighQualityTemplates, :));

if doDebug
    stringTitle = 'Preproc: Iterative Template Creation Single Cycle';
    hp(2) = plot(tTemplate, pulseCleanedTemplate, '.-g', 'LineWidth', 4, ...
        'Marker', 'x');
    legend(hp, 'mean of templates', 'mean of most similar, chosen templates');
    tapas_physio_suptitle(stringTitle);
end

author: translationalneuromodeling

repo: tapas


text:
function [out] = prep_selectTime(dat, varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prep_selectTime (Pre-processing procedure):
%
% Synopsis:
%     [out] = prep_selectTime(DAT,<OPT>)
%
% Example :
%      out = prep_selectTime(dat, {'Time',[1000 3000]})
%
% Arguments:
%     dat - Structure. epoched data
%     varargin - struct or property/value list of optional properties:
%          : time - Time interval to be selected (ms)
%
% Returns:
%     out - Data structure which has selected time from epoched data
%
%
% Description:
%     This function selects the part of a specific time interval
%     from continuous or epoched data.
%     (i)  For continuous data, this function selects data in specifie time
%      interval from the whole data.
%     (ii) For epoched data, this function selects time interval in each trial.
%      If you want to select trials in specific time interval, you can use
%      a function 'prep_selectTrials'
%     continuous data should be [time * channels]
%     epoched data should be [time * channels * trials]
%
% See also 'https://github.com/PatternRecognition/OpenBMI'
%
% Min-ho Lee, 12-2017
% mh_lee@korea.ac.kr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(varargin)
    warning('OpenBMI: Time interval should be specified')
    out = dat;
    return
end
opt = opt_cellToStruct(varargin{:});
if ~isfield(opt,'Time')
    warning('OpenBMI: Time interval should be specified.')
    return
end
ival = opt.Time;

if ~isfield(dat,'x') || ~isfield(dat,'t') || ~isfield(dat,'fs')
    warning('OpenBMI: Data must have fields named ''x'',''t'',''fs''')
    return
end

if isfield(dat,'y_dec') && isfield(dat,'y_logic') && isfield(dat,'y_class')
    a=1;
else
    a=0;
    warning('OpenBMI: Data should have fields named ''y_dec'',''y_logic'',''y_class''')
end

d = ndims(dat.x);
is = ceil(ival(1)*dat.fs/1000);
ie = floor(ival(2)*dat.fs/1000);

if d == 3 || (d==2 && length(dat.chan)==1)
    if ival(1)<dat.ival(1) || ival(2)>dat.ival(end)
        warning('OpenBMI: Selected time interval is out of epoched interval')
        return
    end
    iv = [is:ie]-dat.ival(1)*dat.fs/1000+1;
    x = dat.x(iv,:,:);
    t = dat.t;
%     time = iv/dat.fs*1000; %% Revised hkkim
    time = dat.ival(iv);
    if a
        y_dec = dat.y_dec;
        y_logic = dat.y_logic;
        y_class = dat.y_class;
    end
elseif d == 2 && length(dat.chan)>1
    if ival(1)<0 || ival(2)/1000>size(dat.x,1)/dat.fs
        warning('OpenBMI: Selected time interval is out of time range')
        return
    end
    x = dat.x(is:ie,:);
    s = find((dat.t*1000/dat.fs)>=ival(1));
    e = find((dat.t*1000/dat.fs)<=ival(2));
    iv = s(1):e(end);
    t = dat.t(iv);
    if a
        y_dec = dat.y_dec;
        y_logic = dat.y_logic;
        y_class = dat.y_class;
    end
end
out = rmfield(dat,{'x','t'});
out.x = x;
out.t = t;
if isfield(out,'ival')
    out.ival = time;
end
if a
    out.y_dec = y_dec;
    out.y_logic = y_logic;
    out.y_class = y_class;
end

author: PatternRecognition

repo: OpenBMI


text:
function [Va_idx, gamma_idx, R_idx] = Td2idx(Td)
% Td2idx

Va_trim = Td(1);
gamma_trim = Td(2);
R_trim = Td(3);

if ~exist('P_Td', 'var')
    param_table_T
end
Va_idx = floor((Va_trim - P_Td.Va_min) / P_Td.Va_step)+1;
gamma_idx = floor((gamma_trim - P_Td.gamma_min) / P_Td.gamma_step)+1;
if ~isfinite(R_trim)
    R_idx = 1;
else
    R_idx = floor((R_trim - P_Td.R_min) / P_Td.R_step)+2;
end

if Va_idx > P_Td.Va_n_steps
    Va_idx = P_Td.Va_n_steps;
    str = ['Va_trim ' num2str(Va_trim) 'm/s too big.'];
    disp(str)
elseif Va_idx < 1
    Va_idx = 1;
    str = ['Va_trim ' num2str(Va_trim) 'm/s too small.'];
    disp(str)
end

if gamma_idx > P_Td.gamma_n_steps
    gamma_idx = P_Td.gamma_n_steps;
    str = ['gamma_trim ' num2str(rad2deg(gamma_trim)) ' too big.'];
    disp(str)
elseif gamma_idx < 1
    gamma_idx = 1;
    str = ['gamma_trim ' num2str(rad2deg(gamma_trim)) ' too small.'];
    disp(str)
end

if R_idx > P_Td.R_n_steps
    R_idx = P_Td.R_n_steps;
    str = ['R_trim ' num2str(R_trim) 'm too big.'];
    disp(str)
elseif R_idx < 1
    R_idx = 1;
    str = ['R_trim ' num2str(R_trim) 'm too small.'];
    disp(str)
end

end

author: lis-epfl

repo: swarmlab

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


text:
function probVal=bivarGaussRectangleCDF(rectMin,rectMax,mu,R)
%%BIVARGAUSSRECTANGLECDF Evaluate the an integral of the probability
%      density function of the bivariate normal distribution with a
%      specified mean and covariance matrix over a specified rectangular
%      region. In other words, it is the cumulative distribution function
%      (CDF) over a rectangular region.
%
%INPUTS: rectMin The 2X1 minimum bounds of the rectangular region.
%   rectMax The 2X1 maximum bounds of the rectangular region. Note that
%           rectMax>=rectMin.
%        mu The 2X1 mean of the distribution. If this is omitted or an
%           empty matrix is passed, then the default of [0;0] is used.
%         R The 2X2 covariance matrix of the distribution. If this is
%           omitted or an empty matrix is passed, then R=eye(2,2) is used.
%
%OUTPUTS: val The value of the integral.
%
%We can use the function bivarNormCDF to evaluate the probability of any
%region of the form Pr{x1<b(1), x2<b(2)). Consider 4 points of a rectangle
%dividing 2D space into 9 region.
%   a | b | c
% ----1---2---
%   d | e | f
% ----3---4---
%   g | h | i
%We would like to find the probability in region e. Well, if Pr{number} is
%the probability that everything is below and left of point specified by
%number (i.e. what one gets with the bivarNormCDF function), then one can
%write 4 equations in 4 unknowns:
% Pr{1}=Pr{d}+Pr{g}
% Pr{2}=Pr{d}+Pr{e}+Pr{g}+Pr{h}
% Pr{3}=Pr{g}
% Pr{4}=Pr{g}+Pr{h}
%where Pr{region name} is the probability of of that region. This is a
%system of 4 equation and 4 unknowns, and we solve it to find Pr{e}:
% Pr{e}=-Pr{1}+Pr{2}+Pr{3}-Pr{4}
%
%EXAMPLE 1:
%As a quick check, we know that the integral over the lower-left quadrant
%fof the bivarate normal distribution equals 0.25. That result is obtained
%here.
% mu=[0,;0];
% R=eye(2,2);
% rectMin=[-Inf;-Inf];
% rectMax=[0;0];
% P=bivarGaussRectangleCDF(rectMin,rectMax,mu,R)
%
%EXAMPLE 2:
%Here, we compare the result of this function with the lower-fidelity
%approximation from GaussianD.integralOverRegion. The values are in the
%general ballpark.
% mu=[1,;2];
% R=[3,-2;
%   -2, 4];
% rectMin=[-3;-2];
% rectMax=[2;1];
% P=bivarGaussRectangleCDF(rectMin,rectMax,mu,R)
% PApprox=GaussianD.integralOverRegion(mu,R,rectMin,rectMax)
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public releas

if(nargin<4||isempty(R))
    R=eye(2,2);
end

if(nargin<3||isempty(mu))
    mu=[0;0];
end

p1=[rectMin(1);rectMax(2)];
p2=[rectMax(1);rectMax(2)];
p3=[rectMin(1);rectMin(2)];
p4=[rectMax(1);rectMin(2)];

Pr1=bivarNormCDF(p1,mu,R);
Pr2=bivarNormCDF(p2,mu,R);
Pr3=bivarNormCDF(p3,mu,R);
Pr4=bivarNormCDF(p4,mu,R);

%The min and max is just to handle finite precision effects making the
%result slightly invalid.
probVal=min(1,max(0,-Pr1+Pr2+Pr3-Pr4));

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.


text:
function ierror = perm1_check ( n, p )

%*****************************************************************************80
%
%% PERM1_CHECK checks a permutation of (1,...,N).
%
%  Discussion:
%
%    The routine verifies that each of the integers from 1 to
%    to N occurs among the N entries of the permutation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, integer P(N), the array to check.
%
%    Output, integer IERROR:
%    0, P is a legal permutation of (1,...,N).
%    1, P is not a legal permutation of (1,...,N).
%
  ierror = 0;

  for value = 1 : n

    ierror = 1;

    for location = 1 : n
      if ( p(location) == value )
        ierror = 0;
        break;
      end
    end

    if ( ierror ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PERM1_CHECK - Warning!\n' );
      fprintf ( 1, '  Permutation is missing the value %d.\n', value );
      break;
    end

  end

  return
end


text:
function display_obj(obj,texture)
%
% function display_obj(obj,texture)
%
% displays an obj structure with texture
%
% INPUTS:  obj:     object data
%                   - obj.v:    vertices
%                   - obj.vt:   texture coordinates
%                   - obj.f.v:  face definition vertices
%                   - obj.f.vt: face definition texture
%
%       : texture -  texture image full path
%
% Author: Bernard Abayowa
% University of Dayton
% 6/16/08

texture = imread(texture);
texture_img = flipdim(texture,1);
[sy sx sz] = size(texture_img);
texture_img =  reshape(texture_img,sy*sx,sz);

% make image 3D if grayscale
if sz == 1
    texture_img = repmat(texture_img,1,3);
end

% select what texture correspond to each vertex according to face
% definition
[vertex_idx fv_idx] = unique(obj.f.v);
texture_idx = obj.f.vt(fv_idx);

x = abs(round(obj.vt(:,1)*(sx-1)))+1;
y = abs(round(obj.vt(:,2)*(sy-1)))+1;
xy = sub2ind([sy sx],y,x);
texture_pts = xy(texture_idx);
tval = double(texture_img(texture_pts,:))/255;

% display object
figure; patch('vertices',obj.v,'faces',obj.f.v,'FaceVertexCData', tval);
shading interp
colormap gray(256);
lighting phong;
camproj('perspective');
axis square;
axis off;
axis equal
axis tight;
cameramenu

author: Sable

repo: mcbench-benchmarks


text:
function d = fd04 ( p )

%*****************************************************************************80
%
%% FD04 is a signed distance function for problem 4.
%
%  Modified:
%
%    11 September 2005
%
%  Parameters:
%
%    Input, real P(N,3), one or more points.
%
%    Output, real D(N), the signed distance of each point to the boundary of the region.
%
  d = sqrt ( sum ( p .^ 2, 2 ) ) - 1.0;

  return
end


text:
function varargout = selection_gui(varargin)
% SELECTION_GUI M-file for selection_gui.fig
%      SELECTION_GUI, by itself, creates a new SELECTION_GUI or raises the existing
%      singleton*.
%
%      H = SELECTION_GUI returns the handle to a new SELECTION_GUI or the handle to
%      the existing singleton*.
%
%      SELECTION_GUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SELECTION_GUI.M with the given input arguments.
%
%      SELECTION_GUI('Property','Value',...) creates a new SELECTION_GUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before selection_gui_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to selection_gui_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help selection_gui

% Last Modified by GUIDE v2.5 04-Jan-2008 12:38:42

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @selection_gui_OpeningFcn, ...
                   'gui_OutputFcn',  @selection_gui_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before selection_gui is made visible.
function selection_gui_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to selection_gui (see VARARGIN)
handles.guifig = gcf;
movegui(handles.guifig,'center');

guidata(handles.guifig,handles);
% Choose default command line output for selection_gui
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes selection_gui wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = selection_gui_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
set(handles.guifig,'WindowStyle','Modal'); %make figure modal

uiwait; %wait till the figure is destroyed  or asked to resume

try %this statement is necessary if figure is destroyed , then output argument will be empty by default
    handles = guidata(handles.guifig);
    varargout{1} = handles.selection;
    closereq; % close the gui if OK is pressed
catch
    varargout{1} = [];
end

% --- Executes on button press in ok.
function ok_Callback(hObject, eventdata, handles)
% hObject    handle to ok (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'value');%if there is no selection of radio button and OK is pressed then the selection is Radio Button1 by default
    handles.selection = handles.value;
else
    handles.selection = 'Radio Button1';
end

guidata(hObject, handles);
guidata(handles.guifig, handles);

uiresume;

% --------------------------------------------------------------------
function uipanel1_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to uipanel1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

switch get(hObject,'Tag')   % Get Tag of selected object
    case 'radio1'
        handles.value = 'Radio Button1'; %if Radio Button1 is selected then update handles.value
    case 'radio2'
        handles.value = 'Radio Button2';%if Radio Button2 is selected then update handles.value
end
guidata(hObject, handles);

% --- Executes on button press in cancel.
function cancel_Callback(hObject, eventdata, handles)
% hObject    handle to cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
closereq;

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
delete(hObject);