text:
function [C, sigma] = dataset3Params(X, y, Xval, yval)
%DATASET3PARAMS returns your choice of C and sigma for Part 3 of the exercise
%where you select the optimal (C, sigma) learning parameters to use for SVM
%with RBF kernel
%   [C, sigma] = DATASET3PARAMS(X, y, Xval, yval) returns your choice of C and
%   sigma. You should complete this function to return the optimal C and
%   sigma based on a cross-validation set.
%

% You need to return the following variables correctly.
C = 1;
sigma = 0.3;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the optimal C and sigma
%               learning parameters found using the cross validation set.
%               You can use svmPredict to predict the labels on the cross
%               validation set. For example,
%                   predictions = svmPredict(model, Xval);
%               will return the predictions on the cross validation set.
%
%  Note: You can compute the prediction error using
%        mean(double(predictions ~= yval))
%

params = [0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30];

min_error = 0;

for c=params
    for s=params
        model = svmTrain(X, y, c, @(x1, x2) gaussianKernel(x1, x2, s)); % 模仿ex6的调用方式
        predictions = svmPredict(model, Xval); % 查看svmPredict
        error = mean(double(predictions ~= yval));

        if min_error == 0 || error < min_error
            min_error = error;
            C = c;
            sigma = s;
        end
    end
end

% =========================================================================

end

author: xjwhhh

repo: AndrewNgMachineLearning


text:
% FactorMarginalization Sums given variables out of a factor.
%   B = FactorMarginalization(A,V) computes the factor with the variables
%   in V summed out. The factor data structure has the following fields:
%       .var    Vector of variables in the factor, e.g. [1 2 3]
%       .card   Vector of cardinalities corresponding to .var, e.g. [2 2 2]
%       .val    Value table of size prod(.card)
%
%   The resultant factor should have at least one variable remaining or this
%   function will throw an error.
%
%   See also FactorProduct.m, IndexToAssignment.m, and AssignmentToIndex.m
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function B = FactorMarginalization(A, V)

% Check for empty factor or variable list
if (isempty(A.var) || isempty(V)), B = A; return; end;

% Construct the output factor over A.var \ V (the variables in A.var that are not in V)
% and mapping between variables in A and B
[B.var, mapB] = setdiff(A.var, V);

% Check for empty resultant factor
if isempty(B.var)
    %error('Error: Resultant factor has empty scope');
    B.var = [];
    B.card = [];
    B.val = [];
    return;
end;

% Initialize B.card and B.val
B.card = A.card(mapB);
B.val = zeros(1,prod(B.card));

% Compute some helper indices
% These will be very useful for calculating C.val
% so make sure you understand what these lines are doing
assignments = IndexToAssignment(1:length(A.val), A.card);
indxB = AssignmentToIndex(assignments(:, mapB), B.card);

for i = 1:length(A.val),
    B.val(indxB(i)) = B.val(indxB(i)) + A.val(i);
end;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
%% Simulating Iron Ore Prices
% This example demonstrates calibrating an geomtric brownian motion &
% mean reverting stochastic model from historical data of iron ore prices.
% The model is then used to simulate the prices into the future using
% the Stochastic Differential Equation Simulation engine in Econometrics
% Toolbox.
%
%% Import Historical Data
% Data is saved as a excel file based on historical iron ore prices from
% May 10, 2012 to June 22, 2012
% http://www.indexmundi.com/commodities/?commodity=iron-ore&months=240&currency=aud

% Database Version
importirondb
dates = data.Month;
Date = datenum(dates,'yyyy-mm-dd HH:MM:SS.FFF');
S = data.Price;

% Excel
% [data,txt] = xlsread('ironoremonthly.xlsx');
% dates = txt(3:end,1);
% Date = datenum(dates,'dd/mm/yyyy');
% S = data(:,1);
SR = S(2:end)./S(1:end-1) - 1;
%% The Models
% Return Model (gbm)
% The returns based model used for simulation Iron Orce Prices is geometric
% brownian motion.  It can be described as:
%
% $$ dX_t = \mu(t)X_t dt + D(t,X_t)V(t)dW_t $$
%
% $$ \mu = \textrm{ Instantaneous Rate of Return} $$
%
% $$ \sigma = \textrm { Volatility Rate} $$
%
% Mean Reversion Model (hwv)
% The mean reversion model used for simulating the Iron Ore prices is an Ornstein-Uhlenbeck
% brownian motion with mean reverting drift. This model is fit to the prices
% iron ore prices. The discrete-time equation of this model can be written as,
%
% $$ \Delta x_t = \alpha(\mu - x_t)\Delta t + \sigma dz_t, \textrm{ where } dz_t \sim N(0, \sqrt\Delta_t) $$
%
% $$ \alpha = \textrm{ Mean Reversion Rate}$$
%
% $$ \mu = \textrm{ Mean level} $$
%
% $$ \sigma = \textrm { Volatility } $$

%% Calibrate Parameters

modelsel = 'gbm'; %hwv - mean reverting, %gbm - based on returns
startdate = '2010-08-01';
i = find(datenum(startdate,'yyyy-mm-dd') == Date);
S2 = S(i:end);
Date2 = Date(i:end);
SR2 = SR(i+1:end);

switch modelsel
    case 'gbm'
        dt = 1;
        sigma = std(SR2);
        mu = mean(SR2);
        model = gbm(mu,sigma,'StartState', S(end));
    case 'hwv'
        % hwv - mean reversion model
        % The reversion rate and mean level can be calculated from the coefficients
        % of a linear fit between the prices and their first difference scaled
        % by the time interval parameter. All quantities are specified on an annual
        % scale.
        regressors = [ones(length(S2) - 1, 1) S2(1:end-1)];
        [coefficients, intervals, residuals] = regress(diff(S2), regressors);
        dt    = 1;
        speed = -coefficients(2)/dt;
        level = -coefficients(1)/coefficients(2);
        sigma =  std(residuals)/sqrt(dt);
        model = hwv(speed,level,sigma,  'StartState', S(end));
end

%% Monte-Carlo Simulation
% The model defined above can be simulated with the simulate method of the
% SDE object to generate multiple log price paths. These are exponentiated
% to compute the simulated natural gas prices. The plot below shows 1000
% paths simulated daily, 8 years into the future.

NTrials = 1000;
years = 8;
ann = 12;
NSteps  = years*ann;

tic
Xsim = simulate(model, NSteps, 'NTrials', NTrials, 'DeltaTime', dt);
Xsim = squeeze(Xsim); % Remove redundant dimension
toc

Ssim = Xsim;
SYear = Ssim(13:12:end,:);

% Visualize first 20 prices of 1000 paths
plot(Date, S, Date(end)+(0:30:(years*ann+1)*30-1), Ssim(:,1:100));
datetick; xlabel('Date'); ylabel('Iron Ore Price $US/tonne');
title('Historical & Simulated Iron Ore Prices')
%% Save Model
% The calibrated model is saved in a MAT-file for later use.

save ironoremodel model dt

%% Visual Analysis of Simulated Price Paths
% Instead of plotting a number of paths at once, we can plot longer single paths
% against the observed historical data to visually validate the simulated
% paths. This can serve as a final sanity check.

path = 15;
figure
plot(Date, S, 'b', Date(end)+(0:30:NSteps*30), Ssim(:,path), 'r');
title(['Historical & Simulated Prices, Path ' int2str(path)]);
datetick('x','keeplimits');
xlabel('Date'); ylabel('Iron Ore Price $US/tonne');
%% Histogram results of price paths at the end of each year

figure

for i = 1:years
    subplot(years/2,2,i)
    SsimM = mean(Ssim(i*12+1,:));
    hist(Ssim(i*12,:),20)
    hold on
    hm = plot(SsimM,0,'ro');
    title(['Year: ',num2str(i),' of Iron Ore Prices'])
    xlabel(['Iron Ore Price: Mean =',num2str(round(100*(SsimM))/100)])
    zlabel('Simulation Count')
end

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DeepSqueak 2.7.0                                                        %
% Copyright (c) 2021 Kevin Coffey , Ruby Marx, & Robert Ciszek         %
%                                                                         %
% Licensed under the BSD 3-Clause Licence:                                %
% https://opensource.org/licenses/BSD-3-Clause                            %
%                                                                         %
% Unless required by applicable law or agreed to in writing, software     %
% distributed under the License is distributed on an "AS IS" BASIS,       %
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = DeepSqueak(varargin)

set(groot,'defaultFigureVisible','on');
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @DeepSqueak_OpeningFcn, ...
    'gui_OutputFcn',  @DeepSqueak_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

% --- Executes just before DeepSqueak is made visible.
function DeepSqueak_OpeningFcn(hObject, eventdata, handles, varargin)
% Very Important Logo (Mouse from hjw)
disp '                                                                                .---.'
disp '                                                                                /  .  \  '
disp '                                                                                ) _     _               |\_/|   |'
disp '                                                                                ( (^)-~-(^)              |   |   |'
disp '    ._______________________________________________________________________________________________,-.\_( 0 0 )__,-.__________|___|__,|'
disp '   /  .-.                                                                                ''M''   \   /   ''M''                  |'
disp '  |  /   \                                                                                >o<                         |'
disp '  | |\_.  |                                                                                |  '
disp '  |\|  | /|        ,---,                                    .--.--.                                                         ,-.        |  '
disp '  | `---'' |      .''  .'' `\                      ,-.----.   /  /    ''.   ,----.                                          ,--/ /|        |'
disp '  |       |    ,---.''     \                     \    /  \ |  :  /`. /  /   /  \-.         ,--,                        ,--. :/ |        |'
disp '  |       |    |   |  .`\  |                    |   :    |;  |  |--`  |   :    :|       ,''_ /|                        :  : '' /         |'
disp '  |       |    :   : |  ''  |   ,---.     ,---.  |   | .\ :|  :  ;_    |   | .\  .  .--. |  | :    ,---.     ,--.--.   |  ''  /          |'
disp '  |       |    |   '' ''  ;  :  /     \   /     \ .   : |: | \  \    `. .   ; |:  |,''_ /| :  . |   /     \   /       \  ''  |  :          |'
disp '  |       |    ''   | ;  .  | /    /  | /    /  ||   |  \ :  `----.   \''   .  \  ||  '' | |  . .  /    /  | .--.  .-. | |  |   \         | '
disp '  |       |    |   | :  |  ''.    '' / |.    '' / ||   : .  |  __ \  \  | \   `.   ||  | '' |  | | .    '' / |  \__\/: . . ''  : |. \        |'
disp '  |       |    ''   : | /  ; ''   ;   /|''   ;   /|:     |`-'' /  /`--''  /  `--''""| |:  | : ;  ; | ''   ;   /|  ," .--.; | |  | '' \ \       |'
disp '  |       |    |   | ''` ,/  ''   |  / |''   |  / |:   : :   ''--''.     /     |   | |''  :  `--''   \''   |  / | /  /  ,.  | ''  : |--''        |'
disp '  |       |    ;   :  .''    |   :    ||   :    ||   | :     `--''---''      |   | ::  ,      .-./|   :    |;  :   .''   \;  |,''           |'
disp '  |       |    |   ,.''       \   \  /  \   \  / `---''.|                   `---''.| `--`----''     \   \  / |  ,     .-./''--''             |'
disp '  |       |    ''---''          `----''    `----''    `---`                     `---`                `----''   `--`---''                     |'
disp '  |       |                                                                                |'
disp '  \       |____________________________________________________________________________________________________________________________/'
disp '   \     /'
disp '    `---'''
disp '  '
disp '  '
disp '  '

% Set Handles
hFig = hObject;
handles.hFig=hFig;

% % Fullscreen
% warning ('off','all');
% pause(0.00001);
% frame_h = get(handle(gcf),'JavaFrame');
% set(frame_h,'Maximized',1);

% Create a class to hold the data
squeakfolder = fileparts(mfilename('fullpath'));

% Add to MATLAB path and check for toolboxes
if ~isdeployed
    % Add DeepSqueak to the path
    addpath(squeakfolder);
    addpath(genpath(fullfile(squeakfolder, 'Functions')));
    savepath

    %% Display error message if running on matlab before 2017b or toolboxes not found
    if verLessThan('matlab','9.9')
        errordlg(['Warning, DeepSqueak V3 requires MATLAB 20201 or later. It looks like you are use MATLAB ' version('-release')],'upgrade your matlab')
    end

    try
        verLessThan('nnet','1');
    catch
        warning('Deep Learning Toolbox not found')
    end

    try
        verLessThan('curvefit','1');
    catch
        warning('Curve Fitting Toolbox not found')
    end

    try
        verLessThan('vision','1');
    catch
        warning('Computer Vision System Toolbox not found')
    end

    try
        verLessThan('images','1');
    catch
        warning('Image Processing Toolbox not found')
    end

    try
        verLessThan('parallel','1');
    catch
        warning('Parallel Computing Toolbox not found')
    end
end

handles.data = squeakData(squeakfolder);

set ( hFig, 'Color', [.1 .1 .1] );
handles.output = hObject;
cd(handles.data.squeakfolder);

%% Display version
try % Read the current changelog and find the version (## number)
    fid = fopen(fullfile(handles.data.squeakfolder,'CHANGELOG.md'));
    changelog = fscanf(fid,'%c');
    fclose(fid);
    tokens = regexp(changelog, '## ([.\d])+', 'tokens');
    handles.DSVersion =  tokens{1}{:};
catch
    handles.DSVersion = '? -- can''t read CHANGELOG.md. Make sure you have the latest version!';
end
fprintf(1,'%s %s\n\n', 'DeepSqueak version', handles.DSVersion);
try % Check if a new version is avaliable by comparing changelog to whats online
    WebChangelog = webread('https://raw.githubusercontent.com/DrCoffey/DeepSqueak/master/CHANGELOG.md');
    [changes, tokens] = regexp(WebChangelog, '## ([.\d])+', 'start', 'tokens');
    WebVersion = tokens{1}{:};
    if ~strcmp(WebVersion, handles.DSVersion)
        fprintf(1,'%s\n%s\n\n%s\n',...
            'A new version of DeepSqueak is avaliable.',...
            'Download it here!',...
            WebChangelog(1:changes(2)-1))
    end
catch
    fprintf(1,'Can''t check for a updates online right now\n');
end

% set(handles.spectogramWindow,'Visible', 'off');
% set(handles.epochSpect,'Visible', 'off');
% set(handles.topRightButton, 'Visible', 'off');
% set(handles.topLeftButton, 'Visible', 'off');

if ~(exist(fullfile(handles.data.squeakfolder,'Background.png'), 'file')==2)
    disp('Background image not found')
    background = zeros(280);
    fonts = listTrueTypeFonts;
    background = insertText(background,[10 8],'DeepSqueak','Font',char(datasample(fonts,1)),'FontSize',30);
    background = insertText(background,[10 80],'DeepSqueak','Font',char(datasample(fonts,1)),'FontSize',30);
    background = insertText(background,[10 150],'DeepSqueak','Font',char(datasample(fonts,1)),'FontSize',30);
    background = insertText(background,[10 220],'DeepSqueak','Font',char(datasample(fonts,1)),'FontSize',30);
    handles.background = background;
else
    handles.background=imread('Background.png');
end
if ~(exist(fullfile(handles.data.squeakfolder,'DeepSqueak.fig'), 'file')==2)
    errordlg('"DeepSqueak.fig" not found');
end

% Cool Background Image
imshow(handles.background, 'Parent', handles.focusWindow);
set(handles.focusWindow,'Color',[0.1 0.1 0.1],'YColor',[1 1 1],'XColor',[1 1 1]);
set(handles.focusWindow,'XTick',[]);
set(handles.focusWindow,'YTick',[]);
update_folders(hObject, eventdata, handles);
handles = guidata(hObject);  % Get newest version of handles

% Set the sliders to the saved values
% set(handles.TonalitySlider, 'Value', handles.data.settings.EntropyThreshold);

% Set the page and focus window dropdown boxes to the values defined in
% squeakData, and set the current value to the one closest to the save value.
handles.epochWindowSizePopup.String = compose('%gs', handles.data.pageSizes);
[~, handles.epochWindowSizePopup.Value] =  min(abs(handles.data.pageSizes - handles.data.settings.pageSize));
handles.focusWindowSizePopup.String = compose('%gs', handles.data.focusSizes);
[~, handles.focusWindowSizePopup.Value] =  min(abs(handles.data.focusSizes - handles.data.settings.focus_window_size));

guidata(hObject, handles);

set(handles.contourWindow,'Color',[0.1 0.1 0.1],'YColor',[1 1 1],'XColor',[1 1 1],'Box','off','Clim',[0,1]);
set(handles.contourWindow,'XTickLabel',[]);
set(handles.contourWindow,'XTick',[]);
set(handles.contourWindow,'YTick',[]);

set(handles.waveformWindow,'Color',[0.1 0.1 0.1],'YColor',[1 1 1],'XColor',[1 1 1],'Box','off','Clim',[0,1]);
set(handles.waveformWindow,'XTickLabel',[]);
set(handles.waveformWindow,'XTick',[]);
set(handles.waveformWindow,'YTick',[]);

C = spatialPattern([1000,10000],-2);
imagesc(C(1:900,1:10000),'Parent', handles.spectogramWindow);
colormap(handles.spectogramWindow,inferno);
set(handles.spectogramWindow,'Color',[0.1 0.1 0.1],'YColor',[1 1 1],'XColor',[1 1 1]);
set(handles.spectogramWindow,'XTickLabel',[]);
set(handles.spectogramWindow,'XTick',[]);
set(handles.spectogramWindow,'YTick',[]);

% imagesc(C(900:1000,1:10000),'Parent', handles.detectionAxes);
% colormap(handles.detectionAxes,inferno);
set(handles.detectionAxes,'Color',[64/255 10/255 103/255],'YColor',[1 1 1],'XColor',[1 1 1]);
set(handles.detectionAxes,'XTickLabel',[]);
set(handles.detectionAxes,'XTick',[]);
set(handles.detectionAxes,'YTick',[]);
set(handles.spectogramWindow,'Parent',handles.hFig);

% Set the list of colormaps
handles.popupmenuColorMap.String = {
    'inferno'
    'magma'
    'plasma'
    'viridis'
    'cubehelix'
    'gray'
    'jet'
    'turbo'
    'hot'
    'parula'
    'hsv'
    'cool'
    'spring'
    'summer'
    'autumn'
    'winter'
    'bone'
    'copper'
    'pink'};

function varargout = DeepSqueak_OutputFcn(hObject, eventdata, handles)
shg;
varargout{1} = handles.output;

% --- Executes on button press in PlayCall.
function PlayCall_Callback(hObject, eventdata, handles)
% Play the sound within the boxs
audio = handles.data.AudioSamples(...
    handles.data.calls.Box(handles.data.currentcall, 1),...
    handles.data.calls.Box(handles.data.currentcall, 1) + handles.data.calls.Box(handles.data.currentcall, 3));
playbackRate = handles.data.audiodata.SampleRate * handles.data.settings.playback_rate; % set playback rate
audio = resample(audio, 192000, playbackRate);
audio = audio - mean(audio);
% Use a window funtion to smooth the beginning and end to remove clicks
w = hamming(2000);
audio(1:1000) = audio(1:1000) .* w(1:1000);
audio(end-999:end) = audio(end-999:end) .* w(1001:2000);
% Bandpass Filter
% audio = bandpass(audio,[handles.data.calls.Box(handles.data.currentcall, 2), handles.data.calls.Box(handles.data.currentcall, 2) + handles.data.calls.Box(handles.data.currentcall, 4)] * 1000,handles.data.audiodata.SampleRate);
soundsc(audio,192000);

% --- Executes on button press in NextCall.
function NextCall_Callback(hObject, eventdata, handles)
if handles.data.currentcall < height(handles.data.calls) % If not the last call
    handles.data.currentcall=handles.data.currentcall+1;
    handles.data.focusCenter = handles.data.calls.Box(handles.data.currentcall,1) + handles.data.calls.Box(handles.data.currentcall,3)/2;
end
handles.data.current_call_valid = true;
update_fig(hObject, eventdata, handles);

% --- Executes on button press in PreviousCall.
function PreviousCall_Callback(hObject, eventdata, handles)
if handles.data.currentcall > 1 % If not the first call
    handles.data.currentcall=handles.data.currentcall-1;
    handles.data.focusCenter = handles.data.calls.Box(handles.data.currentcall,1) + handles.data.calls.Box(handles.data.currentcall,3)/2;
end
handles.data.current_call_valid = true;
update_fig(hObject, eventdata, handles);

% --- Executes on selection change in Networks Folder Pop up.
function neuralnetworkspopup_Callback(hObject, eventdata, handles)
guidata(hObject, handles);

% --- Executes during Networks Folder Pop up.
function neuralnetworkspopup_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in Audio Folder Pop up.
function AudioFilespopup_Callback(hObject, eventdata, handles)
guidata(hObject, handles);

% --- Executes during Audio Folder Pop up.
function AudioFilespopup_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in Detection Folder Pop up.
function popupmenuDetectionFiles_Callback(hObject, eventdata, handles)
guidata(hObject, handles);

% --- Executes during Detection Folder Pop up.
function popupmenuDetectionFiles_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in AcceptCall.
function AcceptCall_Callback(hObject, eventdata, handles)
handles.data.calls.Accept(handles.data.currentcall) = true;
handles.update_position_axes = 1;
NextCall_Callback(hObject, eventdata, handles)

% --- Executes on button press in RejectCall.
function RejectCall_Callback(hObject, eventdata, handles)
handles.data.calls.Accept(handles.data.currentcall) = false;
handles.update_position_axes = 1;
NextCall_Callback(hObject, eventdata, handles)

% --- Executes during MAIN AXES CREATION
function focusWindow_CreateFcn(hObject, eventdata, handles)

function slide_focus(focus_offset, hObject, eventdata, handles)
% Move the focus window one unit over
new_position = handles.data.focusCenter + focus_offset;
new_position = min(new_position, handles.data.audiodata.Duration - handles.data.settings.focus_window_size ./ 2);
new_position = max(new_position, handles.data.settings.focus_window_size ./ 2);
handles.data.focusCenter = new_position;

if new_position > handles.data.windowposition + handles.data.settings.pageSize
    forwardButton_Callback(hObject, eventdata, handles);
elseif new_position < handles.data.windowposition
    backwardButton_Callback(hObject, eventdata, handles);
else
    calls_within_window = find(...
        handles.data.calls.Box(:,1) > new_position - handles.data.settings.focus_window_size/2 & ...
        handles.data.calls.Box(:,1) < new_position + handles.data.settings.focus_window_size/2 | ...
        sum(handles.data.calls.Box(:,[1,3]),2) > new_position - handles.data.settings.focus_window_size/2 &...
        sum(handles.data.calls.Box(:,[1,3]),2) < new_position + handles.data.settings.focus_window_size/2);
    if ~isempty(calls_within_window)
        handles.data.currentcall = calls_within_window(1);

        handles.data.current_call_valid = true;
    end

    guidata(hObject,handles);
    update_fig(hObject, eventdata, handles);
end

function score_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function score_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function status_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function status_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on figure1 or any of its controls.
function figure1_WindowKeyPressFcn(hObject, eventdata, handles)
%disp(eventdata);
switch eventdata.Character
    case 'p'
        PlayCall_Callback(hObject, eventdata, handles)
    case {'e', 29} % char(29) is right arrow key
        NextCall_Callback(hObject, eventdata, handles)
    case {'q', 28} % char(28) is left arrow key
        PreviousCall_Callback(hObject, eventdata, handles)
    case 'a'
        AcceptCall_Callback(hObject, eventdata, handles)
    case 'r'
        RejectCall_Callback(hObject, eventdata, handles)
    case 'd'
        rectangle_Callback(hObject, eventdata, handles)
    case 127 % Delete key
        handles.data.calls(handles.data.currentcall,:) = [];
        SortCalls(hObject, [], handles, 'time', 0, handles.data.currentcall - 1);
    case 30 % char(30) is up arrow key
        slide_focus(+ handles.data.settings.focus_window_size, hObject, eventdata, handles)
    case 31 % char(31) is down arrow key
        slide_focus(- handles.data.settings.focus_window_size, hObject, eventdata, handles)
    case 32 % 'space'
        forwardButton_Callback(hObject, eventdata, handles);
    case handles.data.labelShortcuts
        %% Update the call labels
        % Index of the shortcut
        idx = contains(handles.data.labelShortcuts, eventdata.Character);
        handles.data.calls.Type(handles.data.currentcall) = categorical(handles.data.settings.labels(idx));
        update_fig(hObject, eventdata, handles);
end
% drawnow

function figure1_KeyPressFcn(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function popupmenuColorMap_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes during object creation, after setting all properties.
function freq_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in rectangle.
function rectangle_Callback(hObject, eventdata, handles)
current_box = drawrectangle( 'Parent',handles.focusWindow,...
                            'FaceAlpha',0,...
                            'LineWidth',1 );
% Don't do anything if the new box is empty
if current_box.Position(3) == 0 || current_box.Position(4) == 4
    delete(current_box)
    return
end
new_box = table();
new_box.Box = current_box.Position;
new_box.Score = 1;
new_box.Type = categorical({'USV'});
new_box.Accept = true;
handles.data.calls = [handles.data.calls; new_box];

%Now delete the roi and render the figure. The roi will be rendered along
%with the existing boxes.
handles.data.current_call_valid = true;
SortCalls(hObject, eventdata, handles, 'time', 0, -1);
delete(current_box)

% --------------------------------------------------------------------
function select_audio_Callback(hObject, eventdata, handles)
% Find audio in folder
path=uigetdir(handles.data.settings.audiofolder,'Select Audio File Folder');
if isnumeric(path);return;end
handles.data.settings.audiofolder = path;
handles.data.saveSettings();
update_folders(hObject, eventdata, handles);

% --------------------------------------------------------------------
function load_networks_Callback(hObject, eventdata, handles)
% Find networks
path=uigetdir(handles.data.settings.networkfolder,'Select Network Folder');
if isnumeric(path);return;end
handles.data.settings.networkfolder = path;
handles.data.saveSettings();
update_folders(hObject, eventdata, handles);

function load_detectionFolder_Callback(hObject, eventdata, handles)
% Find audio in folder
path=uigetdir(handles.data.settings.detectionfolder,'Select Detection File Folder');
if isnumeric(path);return;end
handles.data.settings.detectionfolder = path;
handles.data.saveSettings();
update_folders(hObject, eventdata, handles);

% --------------------------------------------------------------------
function folders_Callback(hObject, eventdata, handles)

% --------------------------------------------------------------------
function export_raven_Callback(hObject, eventdata, handles)
% Export current file as a txt file for viewing in Raven
% http://www.birds.cornell.edu/brp/raven/RavenOverview.html
raventable = [{'Selection'} {'View'} {'Channel'} {'Begin Time (s)'} {'End Time (s)'} {'Low Freq (Hz)'} {'High Freq (Hz)'} {'Delta Time (s)'} {'Delta Freq (Hz)'} {'Avg Power Density (dB FS)'} {'Annotation'} {'Begin Path'} {'File Offset'}];
View = 'Spectrogram 1';
Channel = 1;
for i = 1:height(handles.data.calls)
    if handles.data.calls.Accept(i)
        Selection = i;
        BeginTime = handles.data.calls.Box(i, 1);
        EndTime = sum(handles.data.calls.Box(i ,[1, 3]));
        LowFreq = handles.data.calls.Box(i, 2) * 1000;
        HighFreq = sum(handles.data.calls.Box(i, [2, 4])) * 1000;
        DeltaTime = EndTime - BeginTime;
        DeltaFreq = HighFreq - LowFreq;
        AvgPwr = 1;
        Annotation = handles.data.calls.Accept(i);
        BeginPath = handles.data.audiodata.Filename;
        FileOffset = handles.data.calls.Box(i, 1);
        raventable = [raventable; {Selection} {View} {Channel} {BeginTime} {EndTime} {LowFreq} {HighFreq} {DeltaTime} {DeltaFreq} {AvgPwr} {Annotation} {BeginPath} {FileOffset}];
    end
end
a  = cell2table(raventable);
handles.current_file_id = get(handles.popupmenuDetectionFiles,'Value');
current_detection_file = handles.detectionfiles(handles.current_file_id).name;
ravenname=[strtok(current_detection_file,'.') '_Raven.txt'];
[FileName,PathName] = uiputfile(ravenname,'Save Raven Truth Table (.txt)');
writetable(a,[PathName FileName],'delimiter','\t','WriteVariableNames',false);
guidata(hObject, handles);

% --------------------------------------------------------------------
function export_Callback(hObject, eventdata, handles)

function training_Callback(hObject, eventdata, handles)

% --------------------------------------------------------------------
function customlabels_Callback(hObject, eventdata, handles)
% Define call categories
prompt = {
    'Label 1  --- Key 1'
    'Label 2  --- Key 2'
    'Label 3  --- Key 3'
    'Label 4  --- Key 4'
    'Label 5  --- Key 5'
    'Label 6  --- Key 6'
    'Label 7  --- Key 7'
    'Label 8  --- Key 8'
    'Label 9  --- Key 9'
    'Label 10  --- Key 0'
    'Label 11  --- Key -'
    'Label 12  --- Key ='
    'Label 13  --- Key !'
    'Label 14  --- Key @'
    'Label 15  --- Key #'
    'Label 16  --- Key $'
    'Label 17  --- Key %'
    'Label 18  --- Key ^'
    'Label 19  --- Key &'
    'Label 20  --- Key *'
    'Label 21  --- Key ('
    'Label 22  --- Key )'
    'Label 23  --- Key _'
    'Label 24  --- Key +'
    };

dlg_title = 'Set Custom Label Names';
num_lines=[1,60]; options.Resize='off'; options.WindowStyle='modal'; options.Interpreter='none';
old_labels = handles.data.settings.labels;
new_labels = inputdlgcol(prompt,dlg_title,num_lines,old_labels,options,3);
if ~isempty(new_labels)
    handles.data.settings.labels = new_labels;
    handles.data.saveSettings();
    update_folders(hObject, eventdata, handles);
end
guidata(hObject, handles);

% --------------------------------------------------------------------
function Change_Playback_Rate_Callback(hObject, eventdata, handles)
prompt = {'Playback Rate: (default = 0.0.5)'};
dlg_title = 'Change Playback Rate';
num_lines=1; options.Resize='off'; options.WindowStyle='modal'; options.Interpreter='tex';
defaultans = {num2str(handles.data.settings.playback_rate)};
rate = inputdlg(prompt,dlg_title,num_lines,defaultans);
if isempty(rate); return; end

[newrate,~,errmsg] = sscanf(rate{1},'%f',1);
disp(errmsg);
if ~isempty(newrate)
    handles.data.settings.playback_rate = newrate;
    handles.data.saveSettings();
    update_folders(hObject, eventdata, handles);
end
guidata(hObject, handles);

% --------------------------------------------------------------------
function Help_Callback(hObject, eventdata, handles)

% --------------------------------------------------------------------
function Untitled_3_Callback(hObject, eventdata, handles)

% --------------------------------------------------------------------
function CallClassification_Callback(hObject, eventdata, handles)

% --------------------------------------------------------------------
function ChangeContourThreshold_Callback(hObject, eventdata, handles)
% Change the contour threshold
prompt = {'Entropy Threshold: (range = 0-1, default = 0.215)', 'Amplitude Percentile Threshold: (range = 0-1, default = 0.825)'};
dlg_title = 'New Contour Threshold:';
num_lines=[1 50]; options.Resize='off'; options.WindowStyle='modal'; options.Interpreter='tex';
defaultans = {num2str(handles.data.settings.EntropyThreshold),num2str(handles.data.settings.AmplitudeThreshold)};
threshold = inputdlg(prompt,dlg_title,num_lines,defaultans);
if isempty(threshold); return; end

[EntropyThreshold,~,errmsg] = sscanf(threshold{1},'%f',1);
disp(errmsg);
[AmplitudeThreshold,~,errmsg] = sscanf(threshold{2},'%f',1);
disp(errmsg);

if ~isempty(EntropyThreshold) && ~isempty(AmplitudeThreshold)

    if AmplitudeThreshold < .001 || AmplitudeThreshold > .999
        disp('Warning! Amplitude Percentile Threshold Must be (0 > 1), Reverting to Default (.825)');
        AmplitudeThreshold = handles.data.defaultSettings.AmplitudeThreshold;
    end
    if EntropyThreshold < .001 || EntropyThreshold > .999
        disp('Warning! Entropy Threshold Must be (0 > 1), Reverting to Default (.215)');
        EntropyThreshold = handles.data.defaultSettings.EntropyThreshold;
    end

    handles.data.settings.EntropyThreshold = EntropyThreshold;
    handles.data.settings.AmplitudeThreshold = AmplitudeThreshold;
    handles.data.saveSettings();
    update_folders(hObject, eventdata, handles);
    try
        update_fig(hObject, eventdata, handles);
    catch
        disp('Could not update figure. Is a call loaded?')
    end
end
guidata(hObject, handles);

% --------------------------------------------------------------------
function ViewManual_Callback(hObject, eventdata, handles)
web('https://github.com/DrCoffey/DeepSqueak/wiki','-browser');

% --------------------------------------------------------------------
function AboutDeepSqueak_Callback(hObject, eventdata, handles)
title = 'About DeepSqueak';

d = dialog('Position',[300 350 500  600],'Name',title,'WindowStyle','Normal','Visible', 'off','Color', [0,0,0]);
movegui(d,'center');
ha = axes(d,'Units','Normalized','Position',[0,0,1,1]);

A = zeros(128);
A = insertText(A,[64,20],'Coffey & Marx, 2019','TextColor','white','BoxColor','Black','AnchorPoint','Center','FontSize',11);
A = insertText(A,[64,64],'DeepSqueak','TextColor','white','BoxColor','Black','AnchorPoint','Center');
A = insertText(A,[64,80],['Version ' handles.DSVersion],'TextColor','white','BoxColor','Black','AnchorPoint','Center','FontSize',11);

A = A(:,:,1);

P = [64, 64];
D = 3;
T = [1,0,-1,0;0,1,0,-1];	% 4 directions
k = 0;

handle_image = imshow(A,[0,1],'parent',ha);

btn = uicontrol('Parent',d,...
    'Units','Normalized',...
    'Position',[.42 .01 .16 .06],...
    'String','Okay',...
    'Callback','delete(gcf)');
set(d,'Visible','on')

while isvalid(handle_image)
    k = k+1;
    a = A(P(1),P(2));
    A(P(1),P(2)) = ~a;
    if ( a )
        D = mod(D+1,4);
    else
        D = mod(D-1,4);
    end
    P = P+T(:,D+1);
    handle_image.CData = A;
    pause(.01)
end

% --- Executes on slider movement.
function TonalitySlider_Callback(hObject, eventdata, handles)
handles.data.settings.EntropyThreshold=(get(hObject,'Value'));
handles.data.saveSettings();
update_fig(hObject, eventdata, handles);

% --- Executes during object creation, after setting all properties.
function TonalitySlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to TonalitySlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --------------------------------------------------------------------
function submit_a_bug_Callback(hObject, eventdata, handles)
% hObject    handle to submit_a_bug (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
web('https://github.com/UEFepilepsyAIVI/DeepSqueak/issues','-browser');

% --- Executes on slider movement.
function optimization_slider_Callback(hObject, eventdata, handles)
hObject.Value = round(hObject.Value);

% --- Executes on button press in backwardButton.
function backwardButton_Callback(hObject, eventdata, handles)
% handles.data.windowposition = max(0, handles.data.windowposition - handles.data.settings.pageSize);
handles.data.focusCenter = max(0, handles.data.windowposition - handles.data.settings.focus_window_size ./ 2);
% get_closest_call_to_focus(hObject, eventdata, handles);

jumps = floor(handles.data.focusCenter / handles.data.settings.pageSize);
handles.data.windowposition = jumps*handles.data.settings.pageSize;

calls_within_window = find(handles.data.calls.Box(:,1) < handles.data.windowposition + handles.data.settings.pageSize, 1, 'last');
if ~isempty(calls_within_window)
    handles.data.currentcall = calls_within_window;
    handles.data.current_call_valid = true;
end

update_fig(hObject, eventdata, handles);

% --- Executes on button press in forwardButton.
function forwardButton_Callback(hObject, eventdata, handles)
% handles.data.windowposition = min(handles.data.audiodata.Duration - handles.data.settings.pageSize, handles.data.windowposition + handles.data.settings.pageSize);
handles.data.focusCenter = handles.data.windowposition + handles.data.settings.pageSize + handles.data.settings.focus_window_size ./ 2;
handles.data.focusCenter = min(handles.data.focusCenter, handles.data.audiodata.Duration - handles.data.settings.focus_window_size ./ 2);

jumps = floor(handles.data.focusCenter / handles.data.settings.pageSize);
handles.data.windowposition = jumps*handles.data.settings.pageSize;

% handles.data.focusCenter = max(0, handles.data.windowposition - handles.data.settings.focus_window_size ./ 2);
% get_closest_call_to_focus(hObject, eventdata, handles);

calls_within_window = find(handles.data.calls.Box(:,1) > handles.data.windowposition, 1);
if ~isempty(calls_within_window)
    handles.data.currentcall = calls_within_window;
    handles.data.current_call_valid = true;
end

update_fig(hObject, eventdata, handles);
% guidata(hObject, handles);

function get_closest_call_to_focus(hObject, eventdata, handles)
calls_within_window = find(...
    handles.data.calls.Box(:,1) > handles.data.windowposition &...
    sum(handles.data.calls.Box(:,[1,3]),2) < handles.data.windowposition + handles.data.settings.pageSize);
if ~isempty(calls_within_window)
    callMidpoints = handles.data.calls.Box(calls_within_window,1) + handles.data.calls.Box(calls_within_window,3)/2;
    [~, closestCall] = min(abs(callMidpoints - handles.data.focusCenter));
    handles.data.currentcall = calls_within_window(closestCall);
    handles.data.current_call_valid = true;
end
update_fig(hObject, eventdata, handles);

% --- Executes during object creation, after setting all properties.
function epochWindowSizePopup_CreateFcn(hObject, eventdata, handles)
% hObject    handle to epochWindowSizePopup (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
padding = cellstr(get(hObject,'String'));
seconds = regexp(padding{get(hObject,'Value')},'([\d*.])*','match');
seconds = str2num(seconds{1});
handles.data.settings.pageSize = seconds;
guidata(hObject, handles);

% --- Executes on selection change in spectogramScalePopup.
function spectogramScalePopup_Callback(hObject, eventdata, handles)
% hObject    handle to spectogramScalePopup (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns spectogramScalePopup contents as cell array
%        contents{get(hObject,'Value')} returns selected item from spectogramScalePopup
update_fig(hObject, [], handles);

% --- Executes during object creation, after setting all properties.
function spectogramScalePopup_CreateFcn(hObject, eventdata, handles)
% hObject    handle to spectogramScalePopup (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function SpectogramMax_Callback(hObject, eventdata, handles)
% hObject    handle to SpectogramMax (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    set(hObject, 'String', num2str(get_spectogram_max(hObject, handles)));
    guidata(hObject,handles);
    update_fig(hObject, [], handles);

% --- Executes during object creation, after setting all properties.
function SpectogramMax_CreateFcn(hObject, eventdata, handles)
% hObject    handle to SpectogramMax (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes during object creation, after setting all properties.
function focusWindowSizePopup_CreateFcn(hObject, eventdata, handles)
% hObject    handle to focusWindowSizePopup (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

padding = cellstr(get(hObject,'String'));
seconds = regexp(padding{get(hObject,'Value')},'([\d*.])*','match');
seconds = str2double(seconds{1});
handles.data.settings.focus_window_size = seconds;
guidata(hObject, handles);

% --- Executes on button press in topLeftButton.
function topLeftButton_Callback(hObject, eventdata, handles)
% hObject    handle to topLeftButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    slide_focus(handles.data.settings.focus_window_size*(-1), hObject, eventdata, handles)

% --- Executes on button press in topRightButton.
function topRightButton_Callback(hObject, eventdata, handles)
% hObject    handle to topRightButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    slide_focus(handles.data.settings.focus_window_size*(1), hObject, eventdata, handles)

% --- Executes on button press in loadAudioFile.
function loadAudioFile_Callback(hObject, eventdata, handles)
h = waitbar(0,'Loading Audio Please wait...');
update_folders(hObject, eventdata, handles);
handles = guidata(hObject);

if nargin == 3 % if "Load Calls" button pressed
    if isempty(handles.audiofiles)
        close(h);
        errordlg(['No valid audio files in current audio folder. Select a folder containing audio with '...
            '"File -> Select Audio Folder", then choose the desired file in the "Audio Files" dropdown box.'])
        return
    end
    handles.current_file_id = get(handles.AudioFilespopup,'Value');
    handles.current_audio_file = handles.audiofiles(handles.current_file_id).name;
end

handles.data.audiodata = audioinfo(fullfile(handles.data.settings.audiofolder,handles.current_audio_file));

Calls = table(zeros(0,4),[],[],[], 'VariableNames', {'Box', 'Score', 'Type', 'Accept'});
% Calls.Box = [0 0 1 1];
% Calls.Score = 0;
% Calls.Type = categorical({'NA'});
% Calls.Power = 1;
% Calls.Accept = false;
handles.data.calls = Calls;
% Position of the focus window
handles.data.focusCenter = handles.data.settings.focus_window_size ./ 2;
initialize_display(hObject, eventdata, handles);
close(h);

% --- Executes during object creation, after setting all properties.
function detectionAxes_CreateFcn(hObject, eventdata, handles)
% hObject    handle to detectionAxes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% --------------------------------------------------------------------
function contributorToolsMenu_Callback(hObject, eventdata, handles)
web('https://gitter.im/DeepSqueak_Community/General')

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

% --- Executes on button press in invert_cmap.
function invert_cmap_Callback(hObject, eventdata, handles)
colormap(handles.spectogramWindow, flipud(colormap(handles.spectogramWindow)))
colormap(handles.focusWindow, flipud(colormap(handles.focusWindow)))

% --- Executes during object creation, after setting all properties.
function waveformWindow_CreateFcn(hObject, eventdata, handles)
% hObject    handle to waveformWindow (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate waveformWindow

% --------------------------------------------------------------------
function Keyboard_Shortcuts_Callback(hObject, eventdata, handles)
% Display a list of keyboard shortcuts
Keyboard_Shortcuts = [
    "Save file", "ctrl + s"
    "Next call", "e, right arrow"
    "Previous call", "q, left arrow"
    "Accept call", "a"
    "Reject call", "r"
    "Delete call", "delete, right click"
    "Redraw box", "d"
    "Play call audio", "p"
    "Set call label", "See ""Add Custom Labels"", double click"
    "Slide foucs forward", "up arrow"
    "Slide foucs back", "down arrow"
    "Next page", "space"
    ];
CreateStruct.Interpreter = 'tex';
CreateStruct.WindowStyle = 'modal';
msgbox(['\fontname{Courier}\fontsize{12}' sprintf('%-20s |   %s\n', Keyboard_Shortcuts')], 'Keyboard Shortcuts', CreateStruct)

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function r8lib_test127 ( )

%*****************************************************************************80
%
%% R8LIB_TEST127 tests R8VEC_MAX_INDEX and R8VEC_MIN_INDEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8LIB_TEST127\n' );
  fprintf ( 1, '  For an R8VEC:\n' );
  fprintf ( 1, '  R8VEC_MAX_INDEX: index of maximum entry;\n' );
  fprintf ( 1, '  R8VEC_MIN_INDEX: index of minimum entry;\n' );

  b = -n;
  c = n;

  seed = 123456789;

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

  r8vec_print ( n, a, '  Input vector:' );

  fprintf ( 1, '\n' );

  ival = r8vec_max_index ( n, a );
  fprintf ( 1, '  Maximum index: %d\n', ival );

  ival = r8vec_min_index ( n, a );
  fprintf ( 1, '  Minimum index: %d\n', ival );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function recon_im = sense_recon(smap, reduced_fft, reduced_dim, np, regularizer, mask, figs_on)
% function recon_im = sense_recon(smap, reduced_fft, reduced_dim, np,
%       regularizer, figs_on)
% in:
%   smap          [N M nc]                complex sensitivity maps
%   reduced_fft   [N/np M] or [N M/np]    undersampled fft
%   reduced_dim   1 or 2                  dimension of undersampling
%   np            degree of undersampling
%   regularizer   string, 'none' or 'tikhonov'
%   mask          [N M] boolean
%   figs_on       boolean that shows intermediate images for debugging
% out:
%   recon_im      [N M]                   SENSE reconstructed image
% This implementation assumes that the fft will be reduced in such a way
% that the DC term is still included (i.e. even DFT coefficients retained)
% 2012-06-07 Mai Le, University of Michigan

assert(isequal(size(mask), size(squeeze(smap(:,:,1)))), 'mask size does not match image size');

% dimensions of original image
dims = [size(smap,1) size(smap,2)];

% number of coils
nc = size(smap,3);

% create aliased intermediate images
aliased_im = zeros(size(reduced_fft,1), size(reduced_fft,2), nc);
for ii = 1:nc
     aliased_im(:,:,ii) = ifft2(reduced_fft(:,:,ii));
end

% plot aliased images
if (figs_on)
    figure;
    for ii = 1:nc
        subplot(2,2,ii); imshow(abs(aliased_im(:,:,ii)),[]); colorbar;
    end
end

% calculate regularization parameter
% I found empirically that somewhere between 0.5% and 1% of the max SVD
% value works well for beta
C = Cdiff1(np);
[U, SIG, V] = svd(C'*C);
% beta = 0.01*max(SIG(:));
beta = 0.005*max(SIG(:));

% reconstruct image by finding (regularized) LS solution for each set of pixels
recon_im = zeros(dims);
indeces = (0:np-1)*dims(reduced_dim)/np;
if (reduced_dim == 1)
    for ii = 1:dims(1)/np
        for jj = 1:dims(2)
            if any(mask(ii+indeces,jj))
%                 S = constructS(reduced_dim, smap, ii+indeces, jj, np);
                ivect = ii+indeces;
                S = squeeze([smap(ivect(mask(ii+indeces,jj)),jj,:)]).';
                if ~all(mask(ii+indeces,jj))
                    S = S.';
                end
                a = squeeze([aliased_im(ii,jj,1:nc)]);
                v = recon_pixels(S,a,regularizer, beta);
                [recon_im] = place_pixels(recon_im, v, reduced_dim, dims, ii, jj, np, mask(ii+indeces,jj));
            end
        end
    end
else
    for ii = 1:dims(1)
        for jj = 1:dims(2)/np
            if any(mask(ii,jj+indeces))
%                 S = constructS(reduced_dim, smap, ii, jj+indeces, np, mask(ii,jj+indeces));
                jvect = jj+indeces;
                S = squeeze([smap(ii,jvect(mask(ii,jj+indeces)),:)]).';
                if ~all(mask(ii,jj+indeces))
                    S = S.';
                end
                a = squeeze([aliased_im(ii,jj,1:nc)]);
                v = recon_pixels(S,a,regularizer, beta);
                [recon_im] = place_pixels(recon_im, v, reduced_dim, dims, ii, jj, np, mask(ii,jj+indeces));
            end
        end
    end
end
% figure; subplot(121); imshow(test_fill,[]); colorbar;
% subplot(122); imshow(real(recon_im),[]); colorbar;
% keyboard;

% selects values from sensitivity map relevant for specific pixels
% depending on reduced_dim, either ivect or jvect is a 1x2 vector
% function S = constructS(reduced_dim, smap, ivect, jvect, np, mask_vector)
%
%
% S = squeeze([smap(ivect,jvect(mask_vector),:)]).';

% places each of the reconstructed pixels in the reconstructed image
function [recon_im] = place_pixels(recon_im, v, reduced_dim, dims, ii, jj, np, mask)
indeces = (0:np-1)*dims(reduced_dim)/np;
if (reduced_dim == 1)
    mk = 1;
    vk = 1;
    while (mk <= np)
        if (mask(mk))
            recon_im(ii+indeces(mk),jj) = v(vk);
            vk = vk + 1;
        end
        mk = mk + 1;
    end
else
    vk = 1;
    mk = 1;
    while (mk <= np)
        if (mask(mk))
            recon_im(ii,jj+indeces(mk)) = v(vk);
            vk = vk + 1;
        end
        mk = mk + 1;
    end
end

% reconstruct pixels with or without Tikhonov regularization
function pixels = recon_pixels(S,a,regularizer, beta)
switch regularizer
    case 'none'
        pixels = S\a;
    case 'tikhonov'
        M = min(size(S));
        pixels = inv(S'*S+(beta^2)*eye(M))*S'*a;
    otherwise
        pixels = S\a;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


% Get the size of the long and short cyclic prefixes based on the sample rate
%
% @param sample_rate Sample rate (in Hz)
% @return long_cp_len Long cyclic prefix length (in samples)
% @return short_cp_len Short cyclic prefix length (in samples)
function [long_cp_len, short_cp_len] = get_cyclic_prefix_lengths(sample_rate)
    long_cp_len = round(1/192000 * sample_rate);
    short_cp_len = round(0.0000046875 * sample_rate);
end


text:
function fx1 = p16_fx1 ( x )

%*****************************************************************************80
%
%% P16_FX1 evaluates the derivative of the function for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX1, the first derivative of the function at X.
%
  e = 0.8;

  fx1 = ( pi / 180.0 ) - e * pi * cos ( pi * x / 180.0  ) / 180.0;

  return
end


text:
function planC = updatePlanIVH(planD)
% updatePlanIVH
%
% Updates the CERR planC with IVH (intensity volume histogram) field.
%
% Usuage
% planC = updatePlanIVH(planC);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

planC = initializeCERR;
indexS = planD{end};
fields = fieldnames(indexS);
indexSNew = planC{end};

for i = 1:length(fields)
    field = fields{i};
    if ~strcmpi(field,'indexS')
        try
            planC{indexSNew.(field)} = planD{indexS.(field)};
        catch
            warning(['Unknown field "' field '"in CERR plan. Rejecting this field']);
            continue
        end
    end
end
planC{end}=indexSNew;
clear planD

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function SR = superresolve_irwsr(slidingWindow, magFactor)

    model = SRModel;
    model.magFactor = magFactor;
    model.psfWidth = 0.4;
    model.motionParams = slidingWindow.flowToReference;
    model.SR = imresize(slidingWindow.referenceFrame, magFactor);

    model.imagePrior = SRPrior('function', @btvPriorWeighted, 'gradient', @btvPriorWeighted_grad, 'weight', [], 'parameters', {size(model.SR), 1, 0.5, []});
    reweightedOptimParams = getReweightedOptimizationParams;
    SR = superresolve(slidingWindow.frames, model, 'reweightedSR', reweightedOptimParams);

author: thomas-koehler

repo: SupER


text:
function traj = ScrewTrajectory(Xstart, Xend, Tf, N, method)
% *** CHAPTER 9: TRAJECTORY GENERATION ***
% Takes Xstart: The initial end-effector configuration,
%       Xend: The final end-effector configuration,
%       Tf: Total time of the motion in seconds from rest to rest,
%       N: The number of points N > 1 (Start and stop) in the discrete
%          representation of the trajectory,
%       method: The time-scaling method, where 3 indicates cubic
%               (third-order polynomial) time scaling and 5 indicates
%               quintic (fifth-order polynomial) time scaling.
% Returns traj: The discretized trajectory as a list of N matrices in SE(3)
%               separated in time by Tf/(N-1). The first in the list is
%               Xstart and the Nth is Xend .
% This function calculates a trajectory corresponding to the screw motion
% about a space screw axis.
% Example Input:
%
% clear; clc;
% Xstart = [[1 ,0, 0, 1]; [0, 1, 0, 0]; [0, 0, 1, 1]; [0, 0, 0, 1]];
% Xend = [[0, 0, 1, 0.1]; [1, 0, 0, 0]; [0, 1, 0, 4.1]; [0, 0, 0, 1]];
% Tf = 5;
% N = 4;
% method = 3;
% traj = ScrewTrajectory(Xstart, Xend, Tf, N, method)
%
% Output:
% traj =
%    1.0000         0         0    1.0000
%         0    1.0000         0         0
%         0         0    1.0000    1.0000
%         0         0         0    1.0000
%
%    0.9041   -0.2504    0.3463    0.4410
%    0.3463    0.9041   -0.2504    0.5287
%   -0.2504    0.3463    0.9041    1.6007
%         0         0         0    1.0000
%
%    0.3463   -0.2504    0.9041   -0.1171
%    0.9041    0.3463   -0.2504    0.4727
%   -0.2504    0.9041    0.3463    3.2740
%         0         0         0    1.0000
%
%   -0.0000    0.0000    1.0000    0.1000
%    1.0000   -0.0000    0.0000   -0.0000
%    0.0000    1.0000   -0.0000    4.1000
%         0         0         0    1.0000

timegap = Tf / (N - 1);
traj = cell(1, N);
for i = 1: N
    if method == 3
        s = CubicTimeScaling(Tf, timegap * (i - 1));
    else
        s = QuinticTimeScaling(Tf, timegap * (i - 1));
    end
    traj{i} = Xstart * MatrixExp6(MatrixLog6(TransInv(Xstart) * Xend) * s);
end
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function s = collapse(t,dims,fun)
%COLLAPSE Collapse sparse tensor along specified dimensions.
%
%  S = COLLAPSE(T,DIMS) sums the entries of T along all dimensions
%  specified in DIMS. If DIMS is negative, then T is summed across
%  all dimensions *not* specified by -DIMS.
%
%  S = COLLAPSE(T) is shorthand for S = COLLAPSE(T,1:ndims(T)).
%
%  S = COLLAPSE(T,DIMS,FUN) accumulates the entries of T using the
%  accumulation function @FUN.
%
%  Examples
%  subs = [1 1 1; 1 1 3; 2 2 4; 4 4 4]
%  vals = [10.5; 1.5; 2.5; 3.5]
%  X = sptensor(subs,vals,[4 4 4]);
%  Y = collapse(X,[2 3]) %<-- sum of entries in each mode-1 slice
%  Y = collapse(ones(X),[1 2]) %<-- nnz in each mode-3 slide
%  Y = collapse(ones(X),[1 2],@max) %<-- 1 if mode-3 has any entry
%  Y = collapse(ones(X),-3,@max); %<-- equivalent
%
%  See also SPTENSOR, SPTENSOR/SCALE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ~exist('fun', 'var')
    fun = @sum;
end

if ~exist('dims', 'var')
    dims = 1:ndims(t);
end

dims = tt_dimscheck(dims,ndims(t));
remdims = setdiff(1:ndims(t),dims);

% Check for the case where we accumulate over *all* dimensions
if isempty(remdims)
    s = fun(t.vals);
    return;
end

% Calculate the size of the result
newsiz = size(t,remdims);

% Check for the case where the result is just a dense vector
if numel(remdims) == 1
    if ~isempty(t.subs)
        s = accumarray(t.subs(:,remdims), t.vals, [newsiz 1], fun);
    else
        s = zeros(newsiz,1);
    end
    return;
end

% Create the result
if ~isempty(t.subs)
  s = sptensor(t.subs(:,remdims), t.vals, newsiz, fun);
else
  s = sptensor([],[],newsiz);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function metsBenchmark = evaluateTendMark( allMets, world )

% Aggregate scores from all sequences
MT = 0; PT = 0; ML = 0; FRA = 0;
falsepositives = 0; missed = 0; idswitches = 0;
Fgt = 0; distsum = 0; Ngt = 0; sumg = 0;
Nc = 0;
numGT = 0; numPRED = 0; IDTP = 0; IDFP = 0; IDFN = 0;

for ind = 1:length(allMets)
    if isempty(allMets(ind).m)
        fprintf('\n\nResults missing for sequence #%d\n', ind)
        continue;
    end
    numGT = numGT + allMets(ind).IDmeasures.numGT;
    numPRED = numPRED + allMets(ind).IDmeasures.numPRED;
    IDTP = IDTP + allMets(ind).IDmeasures.IDTP;
    IDFN = IDFN + allMets(ind).IDmeasures.IDFN;
    IDFP = IDFP + allMets(ind).IDmeasures.IDFP;

    MT = MT + allMets(ind).additionalInfo.MT;
    PT = PT + allMets(ind).additionalInfo.PT;
    ML = ML + allMets(ind).additionalInfo.ML;
    FRA = FRA + allMets(ind).additionalInfo.FRA;
    Fgt = Fgt + allMets(ind).additionalInfo.Fgt;
    Ngt = Ngt + allMets(ind).additionalInfo.Ngt;
    Nc = Nc + sum(allMets(ind).additionalInfo.c);
    sumg = sumg + sum(allMets(ind).additionalInfo.g);
    falsepositives = falsepositives + sum(allMets(ind).additionalInfo.fp);
    missed = missed + sum(allMets(ind).additionalInfo.m);
    idswitches = idswitches + sum(allMets(ind).additionalInfo.mme);
    dists = allMets(ind).additionalInfo.d;
    td = allMets(ind).additionalInfo.td;
    distsum = distsum + sum(sum(dists));
end

IDPrecision = IDTP / (IDTP + IDFP);
IDRecall = IDTP / (IDTP + IDFN);
IDF1 = 2*IDTP/(numGT + numPRED);
if numPRED==0, IDPrecision = 0; end
IDP = IDPrecision * 100;
IDR = IDRecall * 100;
IDF1 = IDF1 * 100;

FAR = falsepositives / Fgt;
MOTP = (1-distsum/Nc) * 100;
if world, MOTP = MOTP / td; end
if isnan(MOTP), MOTP = 0; end
MOTAL=(1-(missed+falsepositives+log10(idswitches+1))/sumg)*100;
MOTA=(1-(missed+falsepositives+idswitches)/sumg)*100;
recall=Nc/sumg*100;
precision=Nc/(falsepositives+Nc)*100;

metsBenchmark = [IDF1, IDP, IDR, recall, precision, FAR, Ngt, MT, PT, ML, falsepositives, missed, idswitches, FRA, MOTA, MOTP, MOTAL];

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function H= plot_channel(epo, clab, varargin)
%PLOT_CHANNEL - Plot the classwise averages of one channel of 1D (time OR
% frequency) or 2D data (time x frequency)
%
%Synopsis:
% H= plot_channel(EPO, CLAB, <OPT>)
%
%Input:
% EPO:  Struct of epoched signals, see makeEpochs
% CLAB: Name (or index) of the channel to be plotted.
% OPT:  Struct or property/value list of optional
% properties, used to generate a plot of 1D or 2D data:
%  ----------------------------------
%   (1) 1D channel data plot:
%  .Butterfly    - Butterfly plot. Recommended when CLAB =  '*'. Only
%                  first class is considered for plotting. If ColorOrder is
%                  default, the current colormap is used.
%  .PlotStat     - plot additional statistic: 'std' for standard deviation,
%                  'sem' for standard error of the mean, 'perc' for percentiles
%                  or 'none' for nothing.
%  .Legend       - show Class legend (1, default), or not (0).
%  .LegendPos    - position of the legend, see help of function 'legend'.
%  .XUnit        - unit of x axis, default 'ms'
%  .YUnit        - unit of y axis, default epo.yUnit if this field
%                  exists, 'a.u.' otherwise
%  .YDir         - 'normal' (negative down) or 'reverse' (negative up)
%  .RefCol       -  Color of patch indicating the baseline interval
%  .RefVSize     - (value, default 0.05) controls the Height of the patch
%                  marking the Reference interval.
%  .ColorOrder   - specifies the Colors for drawing the curves of the
%                  different Classes. If not given the ColorOrder
%                  of the current axis is taken. As special gimmick
%                  you can use 'rainbow' as ColorOrder.
%  .LineWidthOrder - (numerical vector) analog to ColorOrder.
%  .LineStyleOrder - (cell array of strings) analog to ColorOrder.
%  .LineSpecOrder - (cell array of cell array) analog to ColorOrder,
%                  giving full control on the appearance of the curves.
%                  (If LineSpecOrder is defined, LineWidthOrder and
%                  LineStyleOrder are ignored, but ColorOrder is not.)
%  .XGrid, ...   - many axis properties can be used in the usual
%                  way
%  .YLim         - Define the y limits. If empty (default) an automatic
%                  selection according to 'YLimPolicy' is performed.
%  .YLimPolicy   - policy of how to select the YLim (only if 'YLim' is
%                  empty resp unspecified): 'auto' uses the
%                  usual mechanism of Matlab; 'tightest' selects YLim as the
%                  exact data range; 'tight' (default) takes the data range,
%                  adds a little border and selects 'nice' limit values.
%  .Title        - title of the plot to be displayed above the axis.
%                  If OPT.title equals 1, the channel label is used.
%  .Title*       - with * in {'Color', 'FontWeight', 'FontSize'}
%                  selects the appearance of the title.
%  .XZeroLine    - draw an axis along the x-axis at y=0
%  .YZeroLine    - draw an axis along the y-axis at x=0
%  .ZeroLine*    - with * in {'Color','Style'} selects the
%                  drawing style of the axes at x=0/y=0
%  .AxisTitle    - (string) title to be displayed *within* the axis.
%  .AxisTitle*   - with * in {'Color', 'HorizontalAlignment',
%                  'VerticalAlignment', 'FontWeight', 'FontSize'}
%                  selects the appearance of the subplot titles.
%  ----------------------------------
%   (2) 2D channel data plot:
%  .XUnit  - unit of x axis, default 'ms'
%  .YUnit  - unit of y axis, default epo.unit if this field
%                     exists, 'Hz' otherwise
%  .YDir     - 'normal' (low FreqLimuencies at bottom) or 'reverse'
%  .FreqLim     - A vector giving lowest and highest Frequency. If not specified,
%              [1 size(epo.x,2)] is taken.
%  .PlotRef  - if 1 plot Reference interval (default 0)
%  .RefYPos     - y position of Reference line
%  .RefWhisker - length of whiskers (vertical lines)
%  .Ref*     - with * in {'LineStyle', 'LineWidth', 'Color'}
%              selects the appearance of the Reference interval.
%  .Colormap - specifies the colormap for depicting amplitude. Give either
%              a string or a x by 3 Color matrix (default 'jet')
%  .CLim   - Define the Color (=amplitude) limits. If empty (default),
%              limits correspond to the data limits.
%  .CLimPolicy - if 'sym', Color limits are symmetric (so that 0
%              corresponds to the middle of the colormap) (default
%              'normal')
%  .XGrid, ... -  many axis properties can be used in the usual
%                 way
%  .GridOverPatches - if 1 plot grid (default 0)
%  .Title   - title of the plot to be displayed above the axis.
%             If OPT.title equals 1, the channel label is used.
%  .Title*  - with * in {'Color', 'FontWeight', 'FontSize'}
%             selects the appearance of the title.
%  .YTitle  - if set, the title is displayed within the axis, with its
%             Y position corresponding to Ytitle (default [])
%  .ZeroLine  - draw an axis along the y-axis at x=0
%  .ZeroLine*  - with * in {'Color','Style'} selects the
%                drawing style of the axes at x=0/y=0

%Output:
% H - Handle to several graphical objects.
%
% See also plotutil_channel1D, plotutil_channel2D

if nargin==0,
  H= opt_catProps(plotutil_channel1D, plotutil_channel2D);
  return
end

if nargin<2,
  clab= 1;
end

if util_getDataDimension(epo)==1
  if ~isempty(varargin)
    opt1D= plotutil_channel1D;
    opt= opt_structToProplist(opt_substruct(opt_proplistToStruct(varargin{:}),opt1D(:,1)));
  else
    opt=varargin;
  end
  H= plotutil_channel1D(epo, clab, opt{:});
else
  if ~isempty(varargin)
    opt2D= plotutil_channel2D;
    opt=opt_structToProplist( opt_substruct(opt_proplistToStruct(varargin{:}),opt2D(:,1)));
  else
    opt=varargin;
  end
  H= plotutil_channel2D(epo, clab, opt{:});
end

author: PatternRecognition

repo: OpenBMI


text:
%% isGlobalSurfDirOutward
% Below is a demonstration of the features of the |isGlobalSurfDirOutward| function

%%
clear; close all; clc;

%% Syntax
% |[L]=isGlobalSurfDirOutward(F,V);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function OutputToCSV(whatData,writeTimeSeriesData,writeMasterFeatures)
% OutputToCSV     Outputs data to csv for external analysis
%
%---INPUTS:
% whatData, which hctsa .mat file to use (default: HCTSA.mat)
% writeTimeSeriesData, (logical)[false] whether to also output time-series data to file
% writeMasterFeatures, (logical)[false] whether to also output master operation info

% ------------------------------------------------------------------------------
% 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 work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Check inputs and set defaults
%-------------------------------------------------------------------------------
if nargin < 1 || isempty(whatData)
    whatData = 'raw';
end
if nargin < 2 || isempty(writeTimeSeriesData)
    writeTimeSeriesData = false;
end
if nargin < 3
    writeMasterFeatures = false;
end
theDelimiter = ',';
%-------------------------------------------------------------------------------

% Load in the hctsa data:
[TS_DataMat,TimeSeries,Operations,whatData] = TS_LoadData(whatData);

% Get the master operations
if writeMasterFeatures
    MasterOperations = TS_GetFromData(whatData,'MasterOperations');
end

% Get the quality info:
TS_Quality = TS_GetFromData(whatData,'TS_Quality');

% Put NaNs in TS_DataMat where good quality is lacking (ignoring the code):
TS_DataMat(TS_Quality~=0) = NaN;

%-------------------------------------------------------------------------------
% Output data matrix to file:
%-------------------------------------------------------------------------------
fileName = 'hctsa_datamatrix.csv';
dlmwrite(fileName,TS_DataMat,'delimiter',theDelimiter);
fprintf(1,'Wrote feature matrix to %s.\n',fileName);

%-------------------------------------------------------------------------------
% Output time-series data to file:
%-------------------------------------------------------------------------------
if writeTimeSeriesData
    numTimeSeries = height(TimeSeries);
    fprintf(1,'Writing time-series data for %u time series...\n',numTimeSeries);
    fileName = 'hctsa_timeseries-data.csv';
    fid = fopen(fileName,'w');
    for i = 1:numTimeSeries
        x = TimeSeries.Data{i};
        L = length(x);
        for t = 1:L-1
            fprintf(fid,'%.6g,',x(t)); % up to last value with commas
        end
        fprintf(fid,'%.6g\n',x(L)); % final value with new line
    end
    fclose(fid);
    fprintf(1,'Wrote comma-delimited time-series data (6-decimal precision) to %s.\n',fileName);
end

%-------------------------------------------------------------------------------
% Output time-series info to .csv file:
%-------------------------------------------------------------------------------
fileName = 'hctsa_timeseries-info.csv';
if strcmp(theDelimiter,',')
    doQuoteStrings = true;
    % warning(['The comma-delimited Keywords variable may cause trouble for ',...
        % 'the comma-delimited %s file: quoting strings should help'],fileName);
else
    doQuoteStrings = false;
end

% We don't want the data to be written out in this info file:
v = version;
if str2double(v(1:3)) < 9.4
    % 2017 or older MATLAB version
    TimeSeries(:,'Data') = [];
else
    TimeSeries = removevars(TimeSeries,{'Data'});
end
writetable(TimeSeries,fileName,'FileType','text','WriteVariableNames',true,...
                'Delimiter',theDelimiter,'QuoteStrings',doQuoteStrings)
fprintf(1,'Wrote the TimeSeries table to %s.\n',fileName);

%-------------------------------------------------------------------------------
% Output feature info to file:
%-------------------------------------------------------------------------------
fileName = 'hctsa_features.csv';
writetable(Operations,fileName,'FileType','text','WriteVariableNames',true,...
                            'Delimiter',theDelimiter)
fprintf(1,'Wrote the Operations table to %s.\n',fileName);

%-------------------------------------------------------------------------------
% Output master features info to file:
%-------------------------------------------------------------------------------
if writeMasterFeatures
    fileName = 'hctsa_masterfeatures.csv';
    writetable(MasterOperations,fileName,'FileType','text','WriteVariableNames',true,...
                                'Delimiter',theDelimiter)
    fprintf(1,'Wrote the MasterOperations table to %s.\n',fileName);
end

%-------------------------------------------------------------------------------
% Output time-series group info to file:
% (it's accessible from the time-series table)
%-------------------------------------------------------------------------------
% if ismember('Group',TimeSeries.Properties.VariableNames)
%     fileName = 'hctsa_grouplabel-info.csv';
%     fid = fopen(fileName,'w');
%     % Header:
%     fprintf(fid,'%s%s%s\n','Name',theDelimiter,'Group');
%     for i = 1:height(TimeSeries)
%         fprintf(fid,'%s%s%u\n',TimeSeries.Name{i},theDelimiter,TimeSeries.Group(i));
%     end
%     fclose(fid);
%     fprintf(1,'Wrote time-series info to %s.\n',fileName);
% end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [out] = rainfall_1(In,T,p1,varargin)
%rainfall_1

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Rainfall based on temperature threshold
% Constraints:  -
% @(Inputs):    p1   - temperature threshold above which rainfall occurs [oC]
%               T    - current temperature [oC]
%               In   - incoming precipitation flux [mm/d]
%               varargin(1) - smoothing variable r (default 0.01)

if size(varargin,2) == 0
    out = In.*(1-smoothThreshold_temperature_logistic(T,p1));
elseif size(varargin,2) == 1
    out = In.*(1-smoothThreshold_temperature_logistic(T,p1,varargin(1)));
end

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
% [template t2 valid] = PPG_SQI(wave,anntime,annot,template_ahead,windowlen,Fs)
%
% PPG_SQI.m - PPG SQI based on beat template correlation.
% (as an advice, the algorithm get 30 beats at each call and run in loop)
% by Qiao Li 30 Mar 2011
%
% PPG sampling frequency is Fs
%
% input:
%     wave:       PPG data;
%     anntime:    PPG annotation time (samples), read from ple annot file,
%                 But the ann-time is the OFFSET based on wave(1)
%     annot:      Annotation of beats, read from from ple annot file
%                 directly
%     template:   Last PPG beat template
%     windowlen:  length of window to calculate template(default: 30s)
%     Fs       :  sampling frequency (defatult: 125 to work with pervious code)
% output:
%     annot:      ppg sqi annotation
%                     annot.typeMnemonic: E - excellent beat;
%                                         A - acceptable beat;
%                                         Q - unacceptable beat
%                     annot.subtype: SQI based on Direct compare
%                     annot.chan:    SQI based on Linear resampling
%                     annot.num:     SQI based on Dynamic time warping
%                     annot.auxInfo: SQI based on Clipping detection
%     template:   Current PPG beat template
%     valid:      1 or greater for valid template,
%                 0 for invalid template
%
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%
% 12-01-2107 Modified by Giulia Da Poian: replace fixed samplig frequency
% (125) with a variable Fs

function [annot template valid] = PPG_SQI(wave,anntime,annot,template,windowlen,Fs)

    if nargin < 6
        Fs =125;
    end
    if nargin < 5
        windowlen=30*Fs;
    end
    if nargin < 4
        template=[];
    end
    if nargin < 3
        sprintf('Error: must provide wave, anntime and annot');
        annot=[];
        template=[];
        valid=0;
        return;
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 19/09/2011 ADD baseline wander filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    wave = PPGmedianfilter(wave, Fs,Fs);
    % get PPG template
    [t t2 v]=template_pleth(wave(1:min(windowlen,length(wave))), anntime(find(anntime<min(windowlen,length(wave)))), 0,Fs);

    if v<1 && length(template)<1 % Current template invalid && no previous template available
        for j=1:length(annot)
            annot(j).typeMnemonic='Q'; % Unacceptable
        end
        t=[];
    else
        % Using previous template as the template
        if v<1
            t=template;
        end
        % Using t2 if available
        if v>1
            t=t2;
        end

            % Calculate the PLA of template for dynamic time warping
            d1=t;
            d1=(d1-min(d1))/(max(d1)-min(d1)).*100;
            [y1 pla1]=PLA(d1,1,1);

% Main Loop
            for j=1:length(annot)-1

%             SQI1: Direct compare

%             Calculate correlation coefficients based on the template
%             length
                beatbegin=anntime(j);
                beatend=anntime(j+1);
% 07/11/2011 ADD max beat length <= 3s detection
		if beatend-beatbegin>3*Fs
			beatend=beatbegin+3*Fs;
		end
		templatelength=length(t);
		complength=min(templatelength,beatend-beatbegin-1);
                if beatbegin+complength-1 > length(wave) || beatend > length(wave) || beatbegin < 1
                    continue;
                end
                currentb=j;
                cc=corrcoef(t(1:complength),wave(beatbegin:beatbegin+complength-1));
                c1(j)=cc(1,2);
                if (c1(j)<0)
                    c1(j)=0;
                end
                annot(currentb).subtype=int8(c1(j)*100);

%             SQI2: Linear resampling

%             Calculate correlation coefficients based on the linear
%             resampling (interp1)

                y=interp1(1:beatend-beatbegin, wave(beatbegin:beatend-1),1:(beatend-beatbegin-1)/(templatelength-1):(beatend-beatbegin),'spline');
                y(isnan(y))=0;
                cc=corrcoef(t,y);
                c2(j)=cc(1,2);
                if (c2(j)<0)
                    c2(j)=0;
                end
                annot(currentb).chan=int8(c2(j)*100);

%             SQI3: Dynamic Time Warping

%             Calculate correlation coefficients based on the dynamic time
%             warping

                d2=wave(beatbegin:beatend-1);

                % if beat too long, set SQI=0;
                if (length(d2)>length(d1)*10)
                    c3(j)=0;
                else
                    d2=(d2-min(d2))/(max(d2)-min(d2)).*100;
                   [y2 pla2]=PLA(d2,1,1);

                   [w ta tb] = simmx_dtw(y1,pla1,y2,pla2);
                   [p,q,Dm] = dp_dwt2(w,ta,tb);
                   [ym1 ym2 yout1]=draw_dtw(y1,pla1,p,y2,pla2,q);
                    cc=corrcoef(y1,ym2);
                    c3(j)=cc(1,2);
                    if (c3(j)<0)
                        c3(j)=0;
                    end
                end
                annot(currentb).num=int8(c3(j)*100);

%             SQI4: Clipping detection

                d2=wave(beatbegin:beatend-1);
                y=diff(d2);
                clipthreshold=0.5;
                c4(j)=int8(length(find(abs(y)>clipthreshold))/length(y)*100);

%             SQI: Combined

                sqibuf=[annot(currentb).subtype annot(currentb).chan annot(currentb).num c4(j)];
                if min(sqibuf)>=90
                    annot(currentb).typeMnemonic='E'; % Excellent
                else
                    if length(find(sqibuf>=90))>=3 || (median(sqibuf(1:3))>=80 && sqibuf(1)>=50 && sqibuf(4)>=70) || min(sqibuf)>=70
                        annot(currentb).typeMnemonic='A'; % Acceptable
                    else
                        annot(currentb).typeMnemonic='Q'; % Unacceptable
                    end
                end
                annot(currentb).auxInfo=num2str(int8(c4(j)));

            end

    end
    template=t;
    valid=v;

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [mrk, why, channel] = find_artifacts(dat, channels, criteria, amplitude, maxmin, low)
%[mrk, why] = find_artifacts(dat, channels, criteria)
% or
%[mrk, why] = find_artifacts(dat, channels, gradient, amplitude, maxmin, low)
%
% Find artifacts and gives position back. If why is asked for, the
% reason are given back. If chan is asked for chan are given back
% dat are a struct with .x as datas and .fs as used frequence and
% .clab as known channels.
% channels are the channels where the criteria must be
% right. (cell array) (empty = allChannels)
% criteria can be a struct with fields
%  gradient (value for maximal allowed voltage step)
%  amplitude (struct with field min and max (minimal and maximal
%  allowed amplitude), or a twodimenisional vector)
%  maxmin (value with the maximal allowed absolute difference in a
%  segment)
%  low struct with the two fields activity (the lowest allowed
%  activity) and length (the interval length (not set the whole
%  segment)) or twodimensional vector with this values or one value
%  for the lowest allowed activity
%  In this case amplitude, maxmin, low are only used if this
%  values are not determined in criteria.
% criteria can be a value, then equal to criteria.gradient
%  amplitude is then a struct or a twodimensional vector (equal to
%  criteria.amplitude)
%  maxmin a value (equal to criteria.maxmin)
%  low a struct or a twodimensional vector (equal to criteria.low)
%  or a value

misc_checkType(dat,'STRUCT(x fs)');
misc_checkTypeIfExists('channels','CELL{CHAR}|CHAR');

if exist('channels') && ~isempty(channels)
  dat = proc_selectChannels(dat, channels);
end
freq = dat.fs;
chan = dat.clab;
dat = dat.x;

n = nargout;

if ~exist('criteria')
  mrk = [];
  return
end

if ~isstruct(criteria) && ~isempty(criteria)
  criteria.gradient = criteria;
else
  criteria.dummy = [];
end

if ~isfield(criteria,'amplitude') && exist('amplitude') && ~isempty(amplitude)
  criteria.amplitude = amplitude;
end

if ~isfield(criteria,'maxmin') && exist('maxmin') && ~isempty(maxmin)
  criteria.maxmin = maxmin;
end

if ~isfield(criteria,'low') && exist('low') && ~isempty(low)
  criteria.low = low;
end

mrkgr = [];
if isfield(criteria,'gradient')
  gr = dat(2:end,:,:)-dat(1:end-1,:,:);
  gr = abs(gr)>criteria.gradient;
  gro = sum(gr,1);
  gr = sum(gro,2);
  mrkgr = find(gr>0);
  if n > 1
     whygr = cell(size(dat,3),1);
     channelgr = cell(size(dat,3),1);

     for i =mrkgr'
       whygr{i,1} = 'High gradient in channels';
       ch = find(gro(1,:,i)>0);
       for j=ch
         whygr{i,1} = [whygr{i,1}, ' ', chan{j}];
       end
       whygr{i,1} = [whygr{i,1}, '\n'];
       channelgr{i} = {chan{ch}};
     end
  end
end

mrkamp = [];
if isfield(criteria,'amplitude')
  if ~isstruct(criteria.amplitude)
     amp = criteria.amplitude;
     criteria.amplitude.min = amp(1);
     criteria.amplitude.max = amp(2);
  end
  mrkmi = [];
  if isfield(criteria.amplitude,'min')
    mi = dat<criteria.amplitude.min;
    mii = sum(mi,1);
    mi = sum(mii,2);
    mrkmi = find(mi>0);
    if n>1
      whymi = cell(size(dat,3),1);
      channelmi = cell(size(dat,3),1);
      for i =mrkmi'
        whymi{i,1} = 'Too minimal amplitude in channels';
        ch = find(mii(1,:,i)>0);
        for j=ch
          whymi{i,1} = [whymi{i,1}, ' ', chan{j}];
        end
        whymi{i,1} = [whymi{i,1}, '\n'];
	channelmi{i} = {chan{ch}};
      end
    end

  end
  mrkma = [];
  if isfield(criteria.amplitude,'max')
    ma = dat>criteria.amplitude.max;
    maa = sum(ma,1);
    ma = sum(maa,2);
    mrkma = find(ma>0);
    if n>1
      whyma = cell(size(dat,3),1);
      channelma = cell(size(dat,3),1);
      for i =mrkma'
        whyma{i,1} = 'Too maximal amplitude in channels';
        ch = find(maa(1,:,i)>0);
        for j=ch
          whyma{i,1} = [whyma{i,1}, ' ', chan{j}];
        end
        whyma{i,1} = [whyma{i,1}, '\n'];
	channelma{i} = {chan{ch}};
      end
    end
   end
  mrkamp = union(mrkma,mrkmi,'legacy');
end

mrkmaxmin = [];
if isfield(criteria,'maxmin')
  maxmin = max(dat,[],1)-min(dat,[],1);
  maxmini = maxmin>criteria.maxmin;
  maxmin = sum(maxmini,2);
  mrkmaxmin = find(maxmin>0);
  if n > 1
     whymax = cell(size(dat,3),1);
     channelmax = cell(size(dat,3),1);
     for i =mrkmaxmin'
       whymax{i,1} = 'Too maximal difference in channels';
       ch = find(maxmini(1,:,i)>0);
       for j=ch
         whymax{i,1} = [whymax{i,1}, ' ', chan{j}];
       end
       whymax{i,1} = [whymax{i,1}, '\n'];
       channelmax{i} = {chan{ch}};
     end
  end
end

mrklow = [];
if isfield(criteria,'low')
  if ~isstruct(criteria.low)
     lo = criteria.low;
     criteria.low.activity = lo(1);
     if length(lo)>1
     criteria.low.length = lo(2);
     end
  end
  if ~isfield(criteria.low,'length')
      num = size(dat,1);
  else
    num = criteria.low.length *freq/1000;
  end
  if num>=size(dat,1)
    steps = 1;
    num = size(dat,1);
  else
    steps = floor(size(dat,1) - num);
  end
  data= dat;
  for i=1:steps
    lowac = max(data(i:i+num-1,:,:),[],1)-min(data(i:i+num-1,:,:), ...
                                              [],1);
    lowa{i} = lowac < criteria.low.activity;
    lowac = sum(lowa{i});
    mrklow = union(mrklow,find(lowac>0),'legacy');
  end
  if n>1
    whylow = cell(size(dat,3),1);
    channellow = cell(size(dat,3),1);
    for i =mrklow'
      whylow{i,1} = 'Too low activity in channels';
      ch = [];
      for k =1:steps
        ch = union(ch,find(lowa{k}(1,:,i)>0),'legacy');
      end
      for j=ch
        whylow{i,1} = [whylow{i,1}, ' ', chan{j}];
      end
      whylow{i,1} = [whylow{i,1}, '\n'];
      channellow{i} = {chan{ch}};
    end
  end
end

mrk = union(union(mrkgr,mrklow,'legacy'),union(mrkamp,mrkmaxmin,'legacy'),'legacy');
if n>1
  why  = cell(size(dat,3),1);
  for i=1:size(dat,3)
    if exist('whygr')
      why{i,1} = [why{i,1}, whygr{i,1}];
    end
    if exist('whymi')
      why{i,1} = [why{i,1}, whymi{i,1}];
    end
    if exist('whyma')
      why{i,1} = [why{i,1}, whyma{i,1}];
    end
    if exist('whymax')
      why{i,1} = [why{i,1}, whymax{i,1}];
    end
    if exist('whylow')
      why{i,1} = [why{i,1}, whylow{i,1}];
    end
    if ~isempty(why{i,1})
      why{i,1} = sprintf(why{i,1});
    end
  end
end

if n>2
  channel= cell(size(dat,3),1);
  for i=1:size(dat,3)
    channel{i} = {};
    if exist('channelgr') && ~isempty(channelgr{i})
      channel{i} = {channel{i}{:}, channelgr{i}{:}};
    end
    if exist('channelmi') && ~isempty(channelmi{i})
      channel{i} = {channel{i}{:}, channelmi{i}{:}};
    end
    if exist('channelma') && ~isempty(channelma{i})
      channel{i} = {channel{i}{:}, channelma{i}{:}};
    end
    if exist('channelmax') && ~isempty(channelmax{i})
      channel{i} = {channel{i}{:}, channelmax{i}{:}};
    end
    if exist('channellow') && ~isempty(channellow{i})
      channel{i} = {channel{i}{:}, channellow{i}{:}};
    end
  end
end

author: bbci


text:
function [T, roffset, V, P, W] = rendertransform(ax)
% Render Transform
%   [T, offset] = rendertransform(AX) returns the 4-by-4 render transform of axes AX
%   and the render offset.
%   The render transform maps data coordinates to pixels locations in the
%   axes parent coordinate system (with flipped y direction).
%   The render offset is the [x y z] values to subtract from the data inputs
%   before transforming by the render transform.
%
%   [..., V, P, W] = rendertransform(AX) returns additionally the view,
%   projection and viewport transforms.

%   Copyright 2014 The MathWorks, Inc.

newgraphics = ishandle(ax) && isprop(handle(ax), 'XTickLabelRotation');
cpos = get(ax,'CameraPosition');
ctar = get(ax,'CameraTarget');
cup = get(ax,'CameraUpVector');
cva = get(ax,'CameraViewAngle');
cvaAuto = strcmp(get(ax,'CameraViewAngleMode'), 'auto');
oldunits = get(ax,'Units');
set(ax,'Units','pixels')
axpos = get(ax,'Position');
set(ax,'Units',oldunits)
if newgraphics
    bottomleft = floor(axpos(1:2)) - 1;
    topright = ceil(axpos(1:2) + axpos(3:4));
else
    bottomleft = round(axpos(1:2)) - 1;
    topright = round(axpos(1:2) - 1 + axpos(3:4));
end
% round position to integer pixels
viewport = [bottomleft (topright - bottomleft)];
pxscale = max(viewport(3), 1);
pyscale = max(viewport(4), 1);

pbar = get(ax,'PlotBoxAspectRatio');
warp = strcmp(get(ax,'PlotBoxAspectRatioMode'),'auto') && ...
    strcmp(get(ax,'DataAspectRatioMode'),'auto') && ...
    strcmp(get(ax,'CameraViewAngleMode'),'auto');
if warp && newgraphics
    pbar = [1 1 1];
end
perspective = strcmp(get(ax,'Projection'), 'perspective');
xreverse = strcmp(get(ax,'XDir'),'reverse');
yreverse = strcmp(get(ax,'YDir'),'reverse');
zreverse = strcmp(get(ax,'ZDir'),'reverse');
dsx = xlim;
dsy = ylim;
dsz = zlim;
dsrange = [diff(dsx) diff(dsy) diff(dsz)];
roffset = [dsx(1) dsy(1) dsz(1)];
ascale = dsrange;
dsoffset = roffset;
if any(isinf([dsx dsy dsz]))
    error('infinite limits not supported')
end
if any(isinf(dsrange))
    error('limits too large')
end
if strcmp(get(ax,'XScale'),'log') || ...
        strcmp(get(ax,'YScale'),'log') || ...
        strcmp(get(ax,'ZScale'),'log')
    error('log scale not supported')
end

% normalize the camera properties
if xreverse
    roffset(1) = dsx(2);
    cpos(1) = 2*mean(dsx) - cpos(1);
    ctar(1) = 2*mean(dsx) - ctar(1);
    ascale(1) = -ascale(1);
    cup(1) = -cup(1);
end
if yreverse
    roffset(2) = dsy(2);
    cpos(2) = 2*mean(dsy) - cpos(2);
    ctar(2) = 2*mean(dsy) - ctar(2);
    ascale(2) = -ascale(2);
    cup(2) = -cup(2);
end
if zreverse
    roffset(3) = dsz(2);
    cpos(3) = 2*mean(dsz) - cpos(3);
    ctar(3) = 2*mean(dsz) - ctar(3);
    ascale(3) = -ascale(3);
    cup(3) = -cup(3);
end
ncpos = (cpos - dsoffset).*pbar./dsrange;
nctar = (ctar - dsoffset).*pbar./dsrange;
ncup  = cup.*pbar./dsrange;
ncdir = nctar - ncpos;

% make View
dot1 = ncdir*ncdir.';
dot2 = ncdir*ncup.';
u = dot1*ncup - dot2*ncdir;
lat = cross(ncdir, u);
ncdir = ncdir/norm(ncdir);
u = u/norm(u);
lat = lat/norm(lat);
p = ncpos.';
offset = [-lat*p ; -u*p ; -ncdir*p ; 1];
V = [lat.*pbar ; u.*pbar ; ncdir.*pbar ; 0 0 0];
V = [V offset];

% make Projection
P = eye(4);
if perspective
    P(4,4) = 0;
    P(4,3) = 1;
end

% compute size
T = P*V;
cube = [ 0 0 0 0 1 1 1 1   % x
         0 0 1 1 0 0 1 1   % y
         0 1 0 1 0 1 0 1   % z
         1 1 1 1 1 1 1 1]; % w
cube2 = T*cube;
cube3 = cube2;
cube3(1,:) = cube3(1,:) ./ cube2(4,:);
cube3(2,:) = cube3(2,:) ./ cube2(4,:);
width = 2*max(abs(max(cube3(1,:))), abs(min(cube3(1,:))));
height = 2*max(abs(max(cube3(2,:))), abs(min(cube3(2,:))));

% view angle
xConstraint = false;
if cvaAuto
    if abs(pxscale*width) < abs(pyscale*height)
        xConstraint = true;
    end
end

% scale projection
fov = cva*pi/360;
if perspective
    scale = 1/(2*tan(fov));
else
    v = ncpos - nctar;
    len = norm(v);
    scale = 1/(2*tan(fov)*len);
end
P(1,1) = P(1,1)*scale;
P(2,2) = P(2,2)*scale;

% viewport transform
T = P*V;
cube2 = T*cube;
cube3 = cube2;
cube3(1,:) = cube3(1,:) ./ cube2(4,:);
cube3(2,:) = cube3(2,:) ./ cube2(4,:);
width = max(cube3(1,:)) - min(cube3(1,:));
height = max(cube3(2,:)) - min(cube3(2,:));
pxoffset = bottomleft(1);
% flip y in the parent reference frame
parent = get(ax,'Parent');
oldunits = get(parent,'Units');
set(parent,'Units','pixels')
parentpos = get(parent,'Position');
set(parent,'Units',oldunits)
pyoffset = parentpos(4) - bottomleft(2);
ov = [pxoffset+pxscale/2 pyoffset-pyscale/2 0];
if ~warp
    if cvaAuto
        if xConstraint
            scale = pxscale;
        else
            scale = pyscale;
        end
    else
        scale = min(pxscale, pyscale);
    end
    sv = [scale -scale];
else
    xscale = pxscale/width;
    yscale = pyscale/height;
    sv = [xscale -yscale];
end
SC = diag([sv 1 1]);
W = makehgtform('translate',ov) * SC;

NORM = diag([1./ascale 1]);

T = W * P * V * NORM;
roffset = roffset.';

author: brainstorm-tools


text:
function rx = rxOpenTSeriesFig(rx);
% rx = rxOpenTSeriesFig(rx);
%
% Open a figure to navigate across
% the currently-loaded tSeries.
%
% ras 03/05.
if ieNotDefined('rx')
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');
end

if ~isfield(rx,'tSeries')
    msg = 'You must load a tSeries first! See File | Load | mrVista tSeries.';
    myWarnDlg(msg);
    return
end

javaFigs=feature('javafigures');
feature('javafigures', 0);

mrGlobals;
loadSession;

% make a tSeries navigation figure
rx.ui.tSeriesFig = figure('Color','w',...
                          'Units','Normalized',...
                          'Position',[.35 .1 .3 .12],...
                          'Name','Navigate tSeries Frames');

nFrames = size(rx.tSeries,4);

labelTxt = sprintf('Session %s, Data Type %s \n Scans: %s',...
                    mrSESSION.sessionCode,dataTYPES(dt).name,num2str(scans));
rx.ui.tSeriesLabel = uicontrol('Style','text','Units','Normalized',...
                               'Position',[.1 .7 .8 .3],'String',labelTxt,...
                               'BackgroundColor',get(gcf,'Color'),...
                               'FontWeight','bold','FontAngle','italic',...
                               'FontSize',10);

rx.ui.tSeriesVol = rxMakeSlider('Current (target) frame:',[1 nFrames],[.2 .4 .6 .2],1,1);

rx.ui.tSeriesRef = rxMakeSlider('Reference (base) frame:',[1 nFrames],[.2 .1 .6 .2],1,1);

%%%%%%%%%%%%%%%%%
% set callbacks %
%%%%%%%%%%%%%%%%%
% Slider callback string (vol)
cb = 'slider = get(gcbo,''UserData'');';
cb = sprintf('%s \n val = rxSetSlider(slider);',cb);
cb = sprintf('%s \n rxSet([],''targetframe'',val);',cb);
set(rx.ui.tSeriesVol.sliderHandle,'Callback',cb);

% Edit field callback string (vol):
cb = ['val = str2num(get(gcbo,''String'')); ' cb];
set(rx.ui.tSeriesVol.editHandle,'Callback',cb);

% Slider callback string (ref)
cb = 'slider = get(gcbo,''UserData'');';
cb = sprintf('%s \n val = rxSetSlider(slider);',cb);
cb = sprintf('%s \n rxSet([],''baseframe'',val);',cb);
set(rx.ui.tSeriesRef.sliderHandle,'Callback',cb);

% Edit field callback string (ref):
cb = ['val = str2num(get(gcbo,''String'')); ' cb];
set(rx.ui.tSeriesRef.editHandle,'Callback',cb);

feature('javafigures',javaFigs);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function spline_test001 ( )

%*****************************************************************************80
%
%% TEST001 tests PARABOLA_VAL2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  ndim = 1;
  ndata = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST001\n' );
  fprintf ( 1, '  PARABOLA_VAL2 evaluates parabolas through\n' );
  fprintf ( 1, '    3 points in a table\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Our data tables will actually be parabolas:\n' );
  fprintf ( 1, '    A: 2*x**2 + 3 * x + 1.\n' );
  fprintf ( 1, '    B: 4*x**2 - 2 * x + 5.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ndata

    xval = 2.0 * i;
    xdata(i) = xval;
    ydata(1,i) = 2.0 * xval * xval + 3.0 * xval + 1.0;
    zdata(i) = 4.0 * xval * xval - 2.0 * xval + 5.0;

    fprintf ( 1, '%6d  %10f  %10f  %10f\n', i, xdata(i), ydata(1,i), zdata(i) );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Interpolated data:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEFT, X, Y1, Y2\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 5

    xval = 2 * i - 1;
    left = i;

    if ( ndata - 2 < left )
      left = ndata - 2;
    end

    if ( left < 1 )
      left = 1;
    end

    yval = parabola_val2 ( ndim, ndata, xdata, ydata, left, xval );

    zval = parabola_val2 ( ndim, ndata, xdata, zdata, left, xval );

    fprintf ( 1, '%6d  %10f  %10f  %10f\n', left, xval, yval(1), zval(1) );

  end

  return
end


text:
function [cm] = pm2cm(pm)
% Convert length from picometers to centimeters.
% Chad A. Greene 2012
cm = pm*1e-10;

author: Sable


text:
function s = y_to_s_republican ( y, s )

%*****************************************************************************80
%
%% Y_TO_S_REPUBLICAN writes a Republican year into a string.
%
%  Format:
%
%    YearNumber ER
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, string S, a representation of the year.
%
  s = sprintf ( '%d ER', y );

  return
end