text:
function varargout = Read_FrictionmodeGUI(varargin)
% Last Modified by GUIDE v2.5 11-Feb-2019 18:03:10

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

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

% --- Executes just before Read_FrictionmodeGUI is made visible.
function Read_FrictionmodeGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% Choose default command line output for Read_FrictionmodeGUI
handles.output = hObject;
a=get(handles.global_slider,'Value');                                           %Get slider value of "global slider"
set(handles.global_text,'String',num2str(a));                                   %Set text above "global slider" to slider value
set(findall(handles.random_panel1, '-property', 'enable'),'enable', 'off');     %Disable Random_panel1
set(findall(handles.random_panel2, '-property', 'enable'),'enable', 'off');     %Disable Random_panel2
set(findall(handles.raceline_panel, '-property', 'enable'),'enable', 'off');    %Disable Raceline_panel
set(findall(handles.uipanel10, '-property', 'enable'),'enable', 'off');         %Disable Raceline percentage mode panel
handles.another_map_1=0;                                                        %Set default value of "another_map" to "0"
% Update handles structure
guidata(hObject, handles);

% --- Outputs from this function are returned to the command line.
function varargout = Read_FrictionmodeGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

% --- Executes on button press in global_button.
function global_button_Callback(hObject, eventdata, handles)
set(handles.raceline_pathname, 'String', '');                                   %Set raceline pathname to none
set(handles.random_text1,'String','');                                          %Set text above "random slider1" to none
set(handles.random_text2,'String','');                                          %Set text above "random slider2" to none
set(handles.slider1_text2,'String','');                                         %Set text above "raceline slider1" to none
set(handles.slider2_text2,'String','');                                         %Set text above "raceline slider2" to none
cla(handles.axes1);                                                             %Clear axes1
set(findall(handles.global_panel, '-property', 'enable'),'enable', 'on');       %Enable Global panel
set(handles.checkbox, 'Enable','on');                                           %Enable Checkbox of Raceline Percentage Mode
set(findall(handles.random_panel1, '-property', 'enable'),'enable', 'off');     %Disable Random_panel1
set(findall(handles.random_panel2, '-property', 'enable'),'enable', 'off');     %Disable Random_panel2
set(findall(handles.raceline_panel, '-property', 'enable'),'enable', 'off');    %Disable Raceline_panel
a=get(handles.global_slider,'Value');                                           %Get slider value of "global slider"
set(handles.global_text,'String',num2str(a));                                   %Set text above "global slider" to slider value

% --- Executes on button press in random_button.
function random_button_Callback(hObject, eventdata, handles)
set(handles.raceline_pathname, 'String', '');                                   %Set raceline pathname to none
set(handles.global_text,'String','');                                           %Set text above "global slider" to none
set(handles.slider1_text2,'String','');                                         %Set text above "raceline slider1" to none
set(handles.slider2_text2,'String','');                                         %Set text above "raceline slider2" to none
cla(handles.axes1);                                                             %Clear axes1
set(findall(handles.global_panel, '-property', 'enable'),'enable', 'off');      %Disable Global panel
set(handles.checkbox, 'Enable','on');                                           %Enable Checkbox of Raceline Percentage Mode
set(findall(handles.random_panel1, '-property', 'enable'),'enable', 'on');      %Enable Random_panel1
set(findall(handles.random_panel2, '-property', 'enable'),'enable', 'on');      %Enable Random_panel2
set(findall(handles.raceline_panel, '-property', 'enable'),'enable', 'off');    %Disable Raceline_panel
a=get(handles.random_slider1,'Value');                                          %Get slider value of "random slider1"
b=get(handles.random_slider2,'Value');                                          %Get slider value of "random slider2"
set(handles.random_text1,'String',num2str(a));                                  %Set text above "random slider1" to slider value
set(handles.random_text2,'String',num2str(b));                                  %Set text above "random slider2" to slider value

% --- Executes on button press in raceline_button.
function raceline_button_Callback(hObject, eventdata, handles)
set(findall(handles.global_panel, '-property', 'enable'),'enable', 'off');      %Disable Global panel
set(handles.checkbox, 'Value', 0);                                              %Set Checkbox Value to "0"
checkbox_Callback(handles.checkbox, eventdata, handles);                        %Run Checkbox_Callback
set(handles.checkbox, 'Enable','off');                                          %Disable Checkbox of Raceline Percentage Mode
set(handles.raceline_pathname_percent,'String','');                             %Set raceline_pathname_percent to none
set(findall(handles.random_panel1, '-property', 'enable'),'enable', 'off');     %Disable Random_panel1
set(findall(handles.random_panel2, '-property', 'enable'),'enable', 'off');     %Disable Random_panel2
set(findall(handles.raceline_panel, '-property', 'enable'),'enable', 'on');     %Enable Raceline_panel
set(handles.global_text,'String','');                                           %Set text above "global slider" to none
set(handles.random_text1,'String','');                                          %Set text above "random slider1" to none
set(handles.random_text2,'String','');                                          %Set text above "random slider2" to none
a=get(handles.raceline_slider1,'Value');                                        %Get slider value of "raceline slider1"
b=get(handles.raceline_slider2,'Value');                                        %Get slider value of "raceline slider2"
set(handles.slider1_text2,'String',num2str(a));                                 %Set text above "raceline slider1" to slider value
set(handles.slider2_text2,'String',num2str(b));                                 %Set text above "raceline slider2" to slider value
x=-6:0.1:6;                                                                     %Generate array
y=a*exp(-b*0.1*x.^2);                                                           %Generate function value
plot(handles.axes1,x,y)                                                         %Plot function
set(handles.axes1, 'Xlim', [-6,6]);                                             %Set x-limits to [-6,6]

% --- Executes on button press in checkbox.
function checkbox_Callback(hObject, eventdata, handles)
if get(handles.checkbox, 'Value')==1                                            %If checkbox is activated: Enable all features!
    set(findall(handles.uipanel10, '-property', 'enable'),'enable', 'on');      %Enable Raceline percentage mode panel
    set(handles.browse3_button, 'Enable', 'On');                                %Enable Raceline percentage mode browser button
    set(handles.percent_value, 'Enable', 'On');                                 %Enable Raceline percentage mode Value Text
elseif get(handles.checkbox, 'Value')==0                                        %If checkbox is not activated: Disable all features!
    set(findall(handles.uipanel10, '-property', 'enable'),'enable', 'off');     %Disable Raceline percentage mode panel
    set(handles.browse3_button, 'Enable', 'Off');                               %Disable Raceline percentage mode browser button
    set(handles.percent_value, 'Enable', 'Off');                                %Disable Raceline percentage mode Value Text
end

% --- Executes on slider movement.
function raceline_slider1_Callback(hObject, eventdata, handles)
a=get(handles.raceline_slider1,'Value');                                        %Get slider value of "raceline slider1"
b=get(handles.raceline_slider2,'Value');                                        %Get slider value of "raceline slider2"
a=round(a*100)/100;                                                             %Round value of "raceline slider1" to 1/100
set(handles.raceline_slider1,'Value', a);                                       %Set slider value of "raceline slider1" to rounded value
set(handles.slider1_text2,'String',num2str(a));                                 %Set text above "raceline slider1" to slider value
x=-6:0.1:6;                                                                     %Generate array
y=a*exp(-b*0.1*x.^2);                                                           %Generate function value
plot(handles.axes1,x,y)                                                         %Plot function
set(handles.axes1, 'Xlim', [-6,6]);                                             %Set x-limits to [-6,6]

% --- Executes during object creation, after setting all properties.
function raceline_slider1_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function raceline_slider2_Callback(hObject, eventdata, handles)
a=get(handles.raceline_slider1,'Value');                                        %Get slider value of "raceline slider1"
b=get(handles.raceline_slider2,'Value');                                        %Get slider value of "raceline slider2"
b=round(b*100)/100;                                                             %Round value of "raceline slider2" to 1/100
set(handles.raceline_slider2,'Value', b);                                       %Set slider value of "raceline slider2" to rounded value
set(handles.slider2_text2,'String',num2str(b));                                 %Set text above "raceline slider2" to slider value
x=-6:0.1:6;                                                                     %Generate array
y=a*exp(-b*0.1*x.^2);                                                           %Generate function value
plot(handles.axes1,x,y)                                                         %Plot function
set(handles.axes1, 'Xlim', [-6,6]);                                             %Set x-limits to [-6,6]

% --- Executes during object creation, after setting all properties.
function raceline_slider2_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function random_slider1_Callback(hObject, eventdata, handles)
a=get(handles.random_slider1,'Value');                                          %Get slider value of "random slider1"
a=round(a*50)/50;                                                               %Round value of "random slider1" to 1/50
set(handles.random_slider1,'Value', a);                                         %Set slider value of "random slider1" to rounded value
set(handles.random_text1,'String',num2str(a));                                  %Set text above "random slider1" to slider value

% --- Executes during object creation, after setting all properties.
function random_slider1_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function random_slider2_Callback(hObject, eventdata, handles)
a=get(handles.random_slider2,'Value');                                          %Get slider value of "random slider2"
a=round(a*50)/50;                                                               %Round value of "random slider2" to 1/50
set(handles.random_slider2,'Value', a);                                         %Set slider value of "random slider2" to rounded value
set(handles.random_text2,'String',num2str(a));                                  %Set text above "random slider2" to slider value

% --- Executes during object creation, after setting all properties.
function random_slider2_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function global_slider_Callback(hObject, eventdata, handles)
a=get(handles.global_slider,'Value');                                           %Get slider value of "global slider"
a=round(a*50)/50;                                                               %Round value of "global slider" to 1/50
set(handles.global_slider,'Value', a);                                          %Set slider value of "global slider" to rounded value
set(handles.global_text,'String',num2str(a));                                   %Set text above "global slider" to slider value

% --- Executes during object creation, after setting all properties.
function global_slider_CreateFcn(hObject, eventdata, handles) %#ok<*DEFNU,*INUSD>
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on button press in browse2_button.
function browse2_button_Callback(hObject, eventdata, handles)
[filename2, pathname2] = uigetfile('*.csv', 'Pick a .csv file');                    %Browse for file
handles.pathname_raceline = strcat(pathname2, filename2);                           %Connect path name and file name to one string
data_type=lower(filename2(end-2:end));                                              %Extract data type and convert to string
if data_type =='csv'                                                 %#ok<*BDSCA>
    set(handles.raceline_pathname, 'String', handles.pathname_raceline);            %Set GUI.full_path to string fullpathname
else
    set(handles.raceline_pathname, 'String', 'ERROR');                              %Display "ERROR"
end
guidata(hObject, handles);

% --- Executes on button press in browse3_button.
function browse3_button_Callback(hObject, eventdata, handles)
[filename3, pathname3] = uigetfile('*.csv', 'Pick a .csv file');                                %Browse for file
handles.raceline_pathname_percent_name = strcat(pathname3, filename3);                          %Connect path name and file name to one string
data_type=lower(filename3(end-2:end));                                                          %Extract data type and convert to string
if data_type =='csv'                                                 %#ok<*BDSCA>
    set(handles.raceline_pathname_percent, 'String', handles.raceline_pathname_percent_name);   %Set GUI.full_path to string fullpathname
else
    set(handles.raceline_pathname_percent, 'String', 'ERROR');                                  %Display "ERROR"
end
guidata(hObject, handles);

function percent_value_Callback(hObject, eventdata, handles)

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

% --- Executes on button press in create_2_map_button.
function create_2_map_button_Callback(hObject, eventdata, handles)
handles.another_map_1=1;                                                            %Set value of "another_map" to "1"
create_map_button_Callback(handles.create_map_button, eventdata, handles);          %Run create_map_button_Callback

% --- Executes on button press in create_map_button.
function create_map_button_Callback(hObject, eventdata, handles) %#ok<*INUSL>
if handles.another_map_1==1                                                         %if "another_map" =1
    assignin('base', 'another_map', 1);                                             %Handover "another_map"==1 to base workspace
else
    assignin('base', 'another_map', 0);                                             %Handover "another_map"==0 to base workspace
end
a=0;                                                                                %Initialize a consistency check variable
%Check Raceline_percentage_panel on consistency
if get(handles.checkbox, 'Value')==1                                                %If Checkbox for Raceline Percentage Mode is activated
    d=str2double(get(handles.percent_value, 'String'));                                 %Get Raceline Percentage Value
    if d>-99 && d<99                                                                    %If Raceline Percentage Value is consistent
        path_percent = get(handles.raceline_pathname_percent, 'String');                    %Get Raceline Percentage Pathname
        if strcmp(path_percent, "") || strcmp(path_percent,'ERROR')                         %If Pathname is not consistent
            set(handles.raceline_pathname_percent, 'String', 'ERROR');                          %Set Raceline Percentage Pathname to "ERROR"
        else                                                                                %Else (if pathname is consistent)
            assignin('base', 'raceline_percentage_mode', 1);                                    %Handover a positive ("1") raceline percentage mode
            assignin('base', 'raceline_percent_value', d);                                      %Handover the Raceline Percentage value
            assignin('base', 'path_raceline', path_percent);                                    %Handover Raceline Path
            a=1;                                                                                %Set Consistency Check Variable to "1"
        end
    else                                                                                %Else (if raceline percentage value is not consistent)
        set(handles.percent_value, 'String', 'ERROR')                                       %Set Raceline Percentage value to "ERROR"
    end
else                                                                                %Else (if checkbox for Raceline Percentage Mode is not activated)
    assignin('base', 'raceline_percentage_mode', 0);                                    %Handover a negative ("0") raceline percentage mode
    a=1;                                                                                %Set Consistency Check Variable to "1"
end
%End of Consistency Check: a=0: not consistent, a=1: consistent
% If everything is consistent, GUI continues
if a==1
    if get(handles.global_button,'Value')                                   %If Global_Button is activated
        global_value=get(handles.global_slider,'Value');                        %Get slider value of "global slider"
        assignin('base', 'global_value', global_value);                         %Handover "global value" to base workspace
        assignin('base', 'mode', 'global');                                     %Handover "mode" to base workspace
        close(Read_FrictionmodeGUI);                                                      %Close Read_FrictionmodeGUI
        Friction_Map_Creation();                                                %Run Friction_Map_Creation
    elseif get(handles.random_button,'Value')                               %If Random_Button is activated
        random_value_1=get(handles.random_slider1,'Value');                     %Get slider value of "random slider1"
        random_value_2=get(handles.random_slider2,'Value');                     %Get slider value of "random slider2"
        if random_value_1<random_value_2                                        %if slider values are consistent
            assignin('base', 'random_value_1', random_value_1);                     %Handover "random value1" to base workspace
            assignin('base', 'random_value_2', random_value_2);                     %Handover "random value2" to base workspace
            assignin('base', 'mode', 'random');                                     %Handover "mode" to base workspace
            close(Read_FrictionmodeGUI);                                                      %Close Read_FrictionmodeGUI
            Friction_Map_Creation();                                                %Run Friction_Map_Creation
        else                                                                    %else (if slider values are not consistent)
            set(handles.random_text1, 'String', '');                                %Set text above "random slider1" to none
            set(handles.random_text2, 'String', '');                                %Set text above "random slider2" to none
        end
    elseif get(handles.raceline_button,'Value')                             %If Raceline_Button is activated
        path = get(handles.raceline_pathname, 'String');                        %Get Raceline Pathname
        if strcmp(path, "")                                                     %If Raceline Pathname is none
            set(handles.raceline_pathname, 'String', 'ERROR');                      %Set Raceline Pathname to "ERROR"
        elseif strcmp(path,'ERROR')                                             %If Raceline Pathname is "ERROR"
            set(handles.raceline_pathname, 'String', 'ERROR');                      %Set Raceline Pathname to "ERROR"
        else                                                                    %else (if Raceline Pathname is consistent)
            a=get(handles.raceline_slider1,'Value');                                %Get slider value of "raceline slider1"
            b=get(handles.raceline_slider2,'Value');                                %Get slider value of "raceline slider2"
            assignin('base', 'raceline_value_1', a);                                %Handover "raceline value1" to base workspace
            assignin('base', 'raceline_value_2', b);                                %Handover "raceline value2" to base workspace
            assignin('base', 'mode', 'raceline');                                   %Handover "mode" to base workspace
            assignin('base','path_raceline', handles.pathname_raceline);            %Handover Path Raceline
            close(Read_FrictionmodeGUI);                                                      %Close Read_FrictionmodeGUI
            Friction_Map_Creation();                                                %Run Friction_Map_Creation
        end
    end
end

author: TUMFTM


text:
function varargout = process_simulate_ar_spectra( varargin )
% PROCESS_SIMULATE_AR_SPECTRA: Simulate source signals with an auto-regressive model
% Using this ARfit toolbox: https://climate-dynamics.org/software/#arfit
% The interactions in the AR model are defined as spectral peaks (frequency and magnitude)
% Coefficients are computed with the function DesignTwoPoles()
%
% USAGE:   OutputFiles = process_simulate_ar_spectra('Run', sProcess, sInputA)
%               signal = process_simulate_ar_spectra('Compute', b, A, C, nsamples)
%               sFiles = process_simulate_ar_spectra('Test')

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Raymundo Cassani, 2021-2022
%          Francois Tadel, 2022

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription()
    % Description the process
    sProcess.Comment     = 'Simulate AR signals';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Simulate';
    sProcess.Index       = 902;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity2#Simulated_data';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import'};
    sProcess.OutputTypes = {'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;

    % === SUBJECT NAME
    sProcess.options.subjectname.Comment = 'Subject name:';
    sProcess.options.subjectname.Type    = 'subjectname';
    sProcess.options.subjectname.Value   = 'Test';
    % === CONDITION NAME
    sProcess.options.condition.Comment = 'Condition name:';
    sProcess.options.condition.Type    = 'text';
    sProcess.options.condition.Value   = 'Simulation';
    % === NUMBER OF SAMPLES
    sProcess.options.samples.Comment = 'Number of time samples:';
    sProcess.options.samples.Type    = 'value';
    sProcess.options.samples.Value   = {12000, ' (Ntime)', 0};
    % === SAMPLING FREQUENCY
    sProcess.options.srate.Comment = 'Signal sampling frequency:';
    sProcess.options.srate.Type    = 'value';
    sProcess.options.srate.Value   = {1200, 'Hz', 2};
    % === PEAKS: FREQUENCY AND AMPLITUDE
    % === COEFFICIENT MATRIX
    sProcess.options.interactions.Comment = [
                                  'Interactions specification:
' ...
                                  'From, To  /  PeakFrequencies [Hz]  /  PeakMagnitudes [0-1]
' ...
                                  'Example:  1, 2  /  10, 20, 40  /  0.5, 0.2, 1.0
'];
    sProcess.options.interactions.Type    = 'textarea';
    sProcess.options.interactions.Value   = ['1, 1 / 10, 25 / 0.0, 1.0', 10 ...
                                             '2, 2 / 10, 25 / 1.0, 0.0', 10 ...
                                             '3, 3 / 10, 25 / 1.0, 0.2', 10 ...
                                             '1, 3 / 10, 25 / 0.0, 0.6'];
    % === DISPLAY SPECTRAL METRICS
    sProcess.options.display.Comment = {'process_simulate_ar_spectra(''DisplayMetrics'');', '
', 'View spectral metrics'};
    sProcess.options.display.Type    = 'button';
    sProcess.options.display.Value   = [];
    % === DISPLAY COEFFICIENT MATRIX
    sProcess.options.coeff.Comment = {'process_simulate_ar_spectra(''DisplayCoefficients'');', '
', 'Get coefficients matrix'};
    sProcess.options.coeff.Type    = 'button';
    sProcess.options.coeff.Value   = [];
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess)
    Comment = sProcess.Comment;
end

%% ===== GET COEFFICIENTS =====
function [A, b, C, errMsg] = GetCoefficients(sProcess)
    % Initialize returned values
    A = [];
    b = [];
    C = [];
    errMsg = '';
    % Parse interactions
    [interactions, errMsg] = process_tf_bands('ParseBands', sProcess.options.interactions.Value);
    if ~isempty(errMsg)
        return;
    end
    nInteractions = size(interactions, 1);
    ToFroms = zeros(nInteractions, 2);
    % Number of signals from interactions
    for i = 1: nInteractions
        FromTo = eval(['[' interactions{i,1} ']']);
        if length(FromTo) ~= 2
            errMsg = 'Direction must be a pair of values: From, To';
            return;
        end
        ToFroms(i, :) = fliplr(FromTo);
    end
    nSignals = max(ToFroms(:));
    % Parse definition of peaks: frequency and magnitude
    peaks = repmat(struct('freqs', [], 'mags', []), nSignals); % peaks(To, From)
    for i = 1: nInteractions
        peakFreqs = eval(['[' interactions{i,2} ']']);
        peakMags  = eval(['[' interactions{i,3} ']']);
        if length(peakFreqs) ~= length(peakMags)
            errMsg = 'Number of frequencies must match number of magnitudes';
            return;
        end
        peaks(ToFroms(i,1), ToFroms(i,2)).freqs = peakFreqs;
        peaks(ToFroms(i,1), ToFroms(i,2)).mags  = peakMags;
    end
    % Find the maximum order requested (2 * maximun_requested_peaks)
    nOrder = max(max(arrayfun(@(x) length(x.freqs),peaks))) * 2;
    % Signal sampling frequency [Hz]
    sfreq = sProcess.options.srate.Value{1};

    % Compute coefficients
    A = zeros(nSignals, nSignals, nOrder); % To, From, Order
    for iFrom = 1 : nSignals
        for iTo = 1 : nSignals
            frqs = peaks(iTo, iFrom).freqs;
            mags = peaks(iTo, iFrom).mags;
            if ~isempty(frqs)
                [~, a] = DesignTwoPoles(sfreq, frqs, mags, nOrder);
                A(iTo, iFrom, :) = -a(2:end);
            end
        end
    end
    A = reshape(A, nSignals, nSignals * nOrder);
    b = zeros(1, nSignals);
    C = eye(nSignals, nSignals);
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    OutputFiles = {};
    % === GET OPTIONS ===
    % Get subject name
    SubjectName = file_standardize(sProcess.options.subjectname.Value);
    if isempty(SubjectName)
        bst_report('Error', sProcess, sInputs, 'Subject name is empty.');
        return
    end
    % Get condition name
    Condition = file_standardize(sProcess.options.condition.Value);
    % Get signal options
    nsamples = sProcess.options.samples.Value{1};
    srate   = sProcess.options.srate.Value{1};

    % Get coefficients
    [A, b, C, errMsg] = GetCoefficients(sProcess);
    if ~isempty(errMsg)
        bst_report('Error', sProcess, [], errMsg);
        return;
    end
    % Generate the signal
    try
        Data = Compute(b, A, C, nsamples);
    catch
        e = lasterr();
        bst_report('Error', sProcess, [], e);
        return;
    end

    % ===== GENERATE FILE STRUCTURE =====
    % Create empty matrix file structure
    FileMat = db_template('matrixmat');
    FileMat.Value       = Data;
    FileMat.Time        = (0:nsamples-1) ./ srate;
    FileMat.Comment     = sprintf('Simulated AR (%dx%d)', size(Data,1), nsamples);
    FileMat.Description = cell(size(Data,1),1);
    for i = 1:size(Data,1)
        FileMat.Description{i} = ['s', num2str(i)];
    end
    % Add history entry
    FileMat = bst_history('add', FileMat, 'process', 'Simulate AR signals:');
    FileMat = bst_history('add', FileMat, 'process', ['   ' strrep(sProcess.options.interactions.Value, char(10), ' ')]);

    % === OUTPUT CONDITION ===
    % Get subject
    sSubject = bst_get('Subject', SubjectName);
    % Create subject if it does not exist yet
    if isempty(sSubject)
        sSubject = db_add_subject(SubjectName);
    end
    % Default condition name
    if isempty(Condition)
        Condition = 'Simulation';
    end
    % Get condition asked by user
    [sStudy, iStudy] = bst_get('StudyWithCondition', bst_fullfile(SubjectName, Condition));
    % Condition does not exist: create it
    if isempty(sStudy)
        iStudy = db_add_condition(SubjectName, Condition, 1);
        sStudy = bst_get('Study', iStudy);
    end

    % === SAVE FILE ===
    % Output filename
    OutputFiles{1} = bst_process('GetNewFilename', bst_fileparts(sStudy.FileName), 'matrix_simar');
    % Save file
    bst_save(OutputFiles{1}, FileMat, 'v6');
    % Register in database
    db_add_data(iStudy, OutputFiles{1}, FileMat);
end

%% ===== GENERATE SIGNALS =====
% Generate synthetic AR signals
function yOut = Compute(b,A,C,samples)
    yOut = process_simulate_ar('arsim', b',A,C,samples)';
end

%% ===== DISPLAY SPECTRAL METRIC =====
function DisplayMetrics()
    % Get current process structure
    sProcess = panel_process_select('GetCurrentProcess');
    % Get options
    sfreq = sProcess.options.srate.Value{1}; % Signal sampling frequency [Hz]
    % Get coefficients
    [A,~,~,errMsg] = GetCoefficients(sProcess);
    if isempty(A)
        bst_error(['Error: Cannot compute coefficients.' 10 10 errMsg], 'Display metrics', 0);
        return;
    end
    A = reshape(A, size(A,1), size(A,1), []);
    % Display spectral metrics
    hFig = process_simulate_ar('HDisplayMetrics', A, sfreq);
end

%% ===== DISPLAY COEFFICIENTS =====
function DisplayCoefficients()
    % Get current process structure
    sProcess = panel_process_select('GetCurrentProcess');
    % Get coefficients
    [A,~,~,errMsg] = GetCoefficients(sProcess);
    if ~isempty(errMsg)
        bst_error(['Error: Cannot compute coefficients.' 10 10 errMsg], 'Display coefficients', 0);
        return;
    end

    % Titles
    strTitleA = [...
        '
Using the ARSIM function from ARFIT Toolbox (T.Schneider):
' ...
        '    X = b + A1.X1 + ... + Ap.Xp + noise

' ...
        'Coefficient matrix A =  [A1 ... Ap]
' ...
        '     - p is the order of the autorregresive model
' ...
        '     - Ai: [Nsignals_to x Nsignals_from]
'];
    strTitleB = 'Intercept b:   [1 x Nsignals]';
    strTitleC = 'Noise covariance matrix C:   [Nsignals x Nsignals]
';
    % Reshape A to [nSignals, nSignals, nOrder]
    nSignals = size(A, 1);
    A = reshape(A, nSignals, nSignals, []);
    nOrder = size(A,3);
    textA = [];
    textAtot = 'A = [';
    for iOrder = 1 : nOrder
        tmpA = A(:,:,iOrder);
        tmpTextA = evalc('disp(tmpA)');
        tmptextA = ['A', num2str(iOrder), ' = [' 10, tmpTextA(1:end-1), '];', 10];
        textA = [textA, tmptextA];
        textAtot  = [textAtot, 'A', num2str(iOrder), ' '];
    end
    textA = [textA, textAtot(1:end-1), '];'];

    % Text B
    textB = 'b = [';
    for iSignal = 1 : nSignals
        textB = [textB, '0 '];
    end
    textB = [textB(1:end-1), '];'];

    % Text C
    textC = ['C = eye(', num2str(nSignals), ', ', num2str(nSignals), ');'];

    % Using uifigure
    if exist('uitextarea', 'file')
        % Get existing specification figure
        hFig = findall(0, 'Type', 'Figure', 'Tag', 'CoeffMatrix');
        % If the figure doesn't exist yet: create it
        if isempty(hFig)
            hFig = uifigure(...
                'Name',   'AR Coefficients', ...
                'Tag',    'CoeffMatrix', ...
                'Units',  'Pixels');
        % Figure already exists: re-use it
        else
            clf(hFig);
            figure(hFig);
        end
        figpos = get(hFig, 'Position');
        wLabel = figpos(3) - 60;

        % Display A
        labelA = uilabel(hFig, ...
            'Position', [30, 260, wLabel, 220], ...
            'Text', strTitleA, ...
            'Interpreter', 'html', ...
            'FontSize', 12);
        textareaA = uitextarea(hFig, ...
            'Position',[30, 165, wLabel, 150], ...
            'Value', textA, ...
            'WordWrap', 'off');
        % Display B
        labelB = uilabel(hFig, 'Position', [30, 135, wLabel, 20], ...
                 'Text', strTitleB, ...
                 'Interpreter', 'html', 'FontSize', 12);
        textareaB = uitextarea(hFig, 'Position',[30, 105, wLabel, 30], ...
                                     'Value', textB, ...
                                     'WordWrap', 'off');
        % Display C
        labelC = uilabel(hFig, 'Position', [30, 75, wLabel, 20], ...
                 'Text', strTitleC, ...
                 'Interpreter', 'html', 'FontSize', 12);
        textareaC = uitextarea(hFig, 'Position',[30, 45, wLabel, 30], ...
                                     'Value', textC, ...
                                     'WordWrap', 'off');
    % Otherwise, just display as text
    else
        strFigure = [...
            str_striptag(strrep(strTitleA, '
', char(10))), ...
            '-------------------------------------------------' 10 ...
            textA 10, ...
            '-------------------------------------------------' 10 10, ...
            str_striptag(strrep(strTitleB, '
', char(10))), ...
            10 '-------------------------------------------------' 10 ...
            textB 10, ...
            '-------------------------------------------------' 10 10, ...
            str_striptag(strrep(strTitleC, '
', char(10))), ...
            '-------------------------------------------------' 10 ...
            textC 10 ...
            '-------------------------------------------------' 10];
        view_text(strFigure, 'AR Coefficients');
    end
end

%% ===== DESIGN TWO POLES FILTER =====
function [b, a, n_order, rads, bws, rel_gains] = DesignTwoPoles(fs, freqs, gains, n_order)
% Computes a filter with two-poles peak
% A peak is defined by its frequency and its gain
% FREQS are peak frequencies in Hz
% GAINS are relative gains in [0,1] range,
%   GAINS are relative with 1 = to gain for FS/4 peak with BW = 1Hz
%
% References
% [1] Ifeachor, EC (1993), Digital signal processing: A practical approach.
%     Chapter: Design of IIR digital filters ISBN: 020154413X
% [2] Smith, JO (2007), Introduction to digital filters: With audio applications
%     Appendix: Elementary Audio Digital Filters ISBN: 0974560715

    % Compute gain for frequency peak at Fs/4 @ BW = 1 Hz;
    bw_ref = 1;
    radius_ref = exp(-pi * bw_ref / fs);
    theta_ref  = (fs/4) / fs * 2 * pi;
    gain_ref   = 1 / ((1 - radius_ref) * sqrt(1 - (2 * radius_ref * cos(2 * theta_ref))  + (radius_ref .^ 2)));

    % Number of peaks and order of filter
    n_peaks = length(freqs);
    if isempty(n_order)
        n_order = 2 * n_peaks;
    end

    % Find radius for each peak, and compute filter
    for i_peak = 1 : n_peaks
        % Estimate radius for given peak
        peak_theta = (freqs(i_peak) / fs) * 2 * pi;
        % Range to explore
        radius_range = 0 : 0.001 : 1;
        gain_range = (1 / gain_ref) ./ ((1 - radius_range) .* sqrt(1 - (2 .* radius_range .* cos(2 * peak_theta))  + (radius_range .^ 2)));
        [~, ix] = min(abs(gain_range - gains(i_peak)));
        rads(i_peak) = radius_range(ix);
        rel_gains(i_peak) = gain_range(ix);
        bws(i_peak) = -log(rads(i_peak)) * fs / pi;
        % Two poles
        poles_peak = rads(i_peak) .* exp(1i * peak_theta);
        poles_peak = transpose([poles_peak, conj(poles_peak)]);
        % Compute one SOS
        [bs{i_peak},as{i_peak}] = zp2tf([], poles_peak, 1);
    end

    % Convolve coefficients
    a = as{1};
    b = 1;
    for i_peak = 2 : n_peaks
        a = conv(a, as{i_peak});
    end

    % Complete order if required
    if length(a) ~= n_order +1
        n_miss = n_order - (length(a) - 1);
        a = conv([1, zeros(1, n_miss)], a);
    end
end

%% ===== TEST FUNCTION =====
function sFiles = Test() %#ok<DEFNU>
    % Three-variable AR model
    sFiles = bst_process('CallProcess', 'process_simulate_ar_spectra', [], [], ...
        'subjectname', 'Test', ...
        'condition',   'Simulation', ...
        'samples',     12000, ...
        'srate',       100, ...
        'interactions' , ['1, 1 / 10, 25 / 0.0, 1.0', 10 ...
                          '2, 2 / 10, 25 / 1.0, 0.0', 10 ...
                          '3, 3 / 10, 25 / 1.0, 0.2', 10 ...
                          '1, 3 / 10, 25 / 0.0, 0.6']);
end

author: brainstorm-tools


text:
function result = gaussint(P,a,b,ngp)
% Integration using Gauss Quadratures
% by Manuel Diaz 2013.03.20

% compute gauss locations and weights
[w,psi] = gauss1d(ngp);

% compute the change the integration range to -1 to 1
x = (b+a)/2 + (b-a)/2*psi;

% perform quadrature
result  = (b-a)/2*sum(w.*P(x));

repo: Aero-matlab


text:
function [MPa] = dynpcm22MPa(dynpcm2)
% Convert pressure from dynes per square centimeter to megapascals
% Chad Greene 2012
MPa = dynpcm2*1.00000e-7;

author: Sable

repo: mcbench-benchmarks


text:
function msm_to_mm_coordinate_integer_symmetric ( output_filename, a )

%*****************************************************************************80
%
%% MSM_TO_MM_COORDINATE_INTEGER_SYMMETRIC writes a "matrix coordinate integer symmetric" Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
  [ nrow, ncol ] = size ( a );
  nnzeros = nnz ( a );

  fid = fopen ( output_filename, 'wt+' );

  if ( fid < 0 );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM_COORDINATE_INTEGER_SYMMETRIC - Fatal error!\n' );
    fprintf ( 1, '  Cannot open the output file.\n' );
    error ( 'MSM_TO_MM_COORDINATE_INTEGER_SYMMETRIC - Fatal error!');
  end;

  fprintf ( fid, '%%%%MatrixMarket matrix coordinate integer symmetric\n');
  fprintf ( fid, '%%%%  Created by MSM_TO_MM_COORDINATE_INTEGER_SYMMETRIC.M\n' );
  fprintf ( fid, '  %d  %d  %d\n', nrow, ncol, nnzeros );

  for j = 1 : ncol
    [ rows, temp, vals ] = find ( a(:,j) );
    sz = size ( rows );
    for k2 = 1 : sz
      if ( j <= rows(k2) )
        fprintf ( fid, '  %d  %d  %d\n', rows(k2), j, vals(k2) );
      end
    end

  end

  fclose ( fid );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function c=addlogs(a,b)
%ADDLOGS Add numbers represented by their logarithms.
%
%  Description
%    C=ADDLOGS(A,B) computes C=log(exp(A)+exp(B)) in such a fashion
%    that it works even when A and B have large magnitude.

% Copyright (c) 2003 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if a>b
  c = a + log(1+exp(b-a));
else
  c = b + log(1+exp(a-b));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

%*****************************************************************************80
%
%% RULE01 returns the rule of degree 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer N, the number of nodes.
%
%    Output, real X(2,N), the coordinates of the nodes.
%
%    Output, real W(N), the weights.
%
  xs = [ ...
       0.00000000000000000E+00 ];
  ys = [ ...
       0.00000000000000000E+00 ];
  ws = [ ...
       0.28284271247461904E+01 ];

  x(1,1:n) = xs(1:n);
  x(2,1:n) = ys(1:n);
  w(1:n) = ws(1:n);

  return
end


text:
function grid = GetSearchGrid(model, material, fix, fixedvals)
% Returns a list of parameter combinations for the specified model
% for the GridSearch function to use in search for the combination
% that best fits a set of measurements.
%
% grid=GetSearchGrid(model, material, fix, fixedvals)
% returns a list of parameter combinations for the model that
% are realistic for a particular material type.
%
% model is a string specifying the model.
%
% material is a string specifying the type of material, which
% determines the range of each parameter in the grid.
% Options are:
% invivo
% invivopreterm
% invivowhitematter
% fixedwhitematter
%
% fix is a list of binary numbers specifying which
% model parameters have fixed values.  By default the
% list is all zeros so no parameters are fixed.
%
% fixedvals is an array the same size as fix that specifies
% the fixed values of any fixed parameters.  Entries in
% fixedvals in locations where fix has value zero are not
% used.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%         Gary Hui Zhang     (gary.zhang@ucl.ac.uk)
%

if(nargin<3)
    fix=zeros(6);
    fixedvals = zeros(6);
end

if(strcmp(material, 'invivo'))
    fs = [0.0 0.25 0.5 0.75 1.0];
    dpars = [13.0 17.0 21.0]*1E-10;
    disos = [10.0 20.0 30.0 50.0]*1E-10;
    Rs=[0.5 1 2 4]*1E-6;
    irfracs=[0 0.1 0.2 0.3];
    fisos=[0.0 0.25 0.5 0.75 1.0];
    kappas = [0.5 1 2 4 8];
    fic = [0.3 0.5 0.7];
elseif(strcmp(material, 'exvivo'))
    fs = [0.0 0.25 0.5 0.75 1.0];
    dpars = [3.0 4.5 6.0 7.5]*1E-10;
    disos = [5.0 10.0 15.0]*1E-10;
    Rs=[1 2 4 8]*1E-6;
    irfracs=[0 0.1 0.2 0.3];
    fisos=[0.0 0.25 0.5 0.75 1.0];
    kappas = [0.5 1 2 4 8];
    fic = [0.3 0.5 0.7];
elseif(strcmp(material, 'invivopreterm'))
    fs = [0.0 0.1 0.2 0.3];
    dpars = [13.0 17.0 21.0]*1E-10;
    disos = [10.0 20.0 30.0 50.0]*1E-10;
    Rs=[1 2 4 8]*1E-6;
    irfracs=[0 0.1 0.2 0.3];
    fisos=[0.0 0.25 0.5 0.75 1.0];
    kappas = [0.5 1 2 4 8];
    fic = [0.3 0.5 0.7];
elseif(strcmp(material, 'invivowhitematter'))
    fs = [0.5 0.7 0.9];
    dpars = [10.0 13.0 15.0 17.0 19.0 21.0 23.0 25.0]*1E-10;
    disos = [10.0 20.0 30.0 50.0]*1E-10;
    Rs=[1 2 4 8]*1E-6;
    irfracs=[0 0.1 0.2 0.3];
    fisos=[0 0.2 0.4];
    kappas = [4 8 16 32 64 128];
    fic = [0.3 0.5 0.7];
elseif(strcmp(material, 'postmortemwhitematter'))
    fs = [0.5 0.7 0.9];
    dpars = [2.0 3.0 4.0 5.0 6.0]*1E-10;
    disos = [5.0 10.0 15.0]*1E-10;
    Rs=[1 2 4 8]*1E-6;
    irfracs=[0 0.1 0.2 0.3];
    fisos=[0 0.2 0.4];
    kappas = [4 8 16 32 64 128];
    fic = [0.3 0.5 0.7];
else
    error(['Unknown material: ', tissue]);
end

% Adjust for fixed parameters.  The first two parameters are the same for
% all models.
if(fix(1))
    fs = [fixedvals(1)];
end
if(fix(2))
    dpars = [fixedvals(2)];
end

if(strcmp(model, 'CylSingleRadTortGPD'))
    if(fix(3))
        Rs = [fixedvals(3)];
    end

    numCombs = length(Rs)*length(fs)*length(dpars);
    grid = zeros(3, numCombs);
    ind = 1;
    for i=1:length(Rs)
        for j=1:length(fs)
            for k=1:length(dpars)
                pars = [fs(j) dpars(k) Rs(i)];
                grid(:,ind) = pars;
                ind = ind + 1;
            end
        end
    end
elseif(strcmp(model, 'CylSingleRadGPD'))
    if(fix(4))
        Rs = [fixedvals(4)];
    end

    numCombs = length(Rs)*length(fs)*length(dpars);
    grid = zeros(4, numCombs);
    ind = 1;
    for i=1:length(Rs)
        for j=1:length(fs)
            for k=1:length(dpars)
                % Set dperp using the standard tortuosity model for
                % randomly placed cylinders unless fixed.
                dperp = dpars(k)*(1-fs(j));
                if(fix(3))
                    dperp = fixedvals(3);
                end
                pars = [fs(j) dpars(k) dperp Rs(i)];
                grid(:,ind) = pars;
                ind = ind + 1;
            end
        end
    end
elseif(strcmp(model, 'CylSingleRadTortIsoGPD'))
    if(fix(3))
        Rs = [fixedvals(3)];
    end
    if(fix(4))
        fisos = [fixedvals(4)];
    end

    numCombs = length(Rs)*length(fs)*length(dpars)*length(fisos);
    grid = zeros(4, numCombs);
    ind = 1;
    for i=1:length(Rs)
        for j=1:length(fs)
            for k=1:length(dpars)
                for l=1:length(fisos)
                    pars = [fs(j) dpars(k) Rs(i) fisos(l)];
                    grid(:,ind) = pars;
                    ind = ind + 1;
                end
            end
        end
    end
elseif(strcmp(model, 'CylSingleRadIsoV_GPD') || strcmp(model, 'CylSingleRadIsoV_GPD_B0'))
    if(fix(4))
        Rs = [fixedvals(4)];
    end
    if(fix(5))
        fisos = [fixedvals(5)];
    end
    if(fix(6))
        disos = [fixedvals(6)];
    end

    numCombs = length(Rs)*length(fs)*length(dpars)*length(fisos)*length(disos);
    grid = zeros(6, numCombs);
    ind = 1;
    for i=1:length(Rs)
        for j=1:length(fs)
            for k=1:length(dpars)
                for l=1:length(fisos)
                    for m=1:length(disos)
                        % Set dperp using the standard tortuosity model for
                        % randomly placed cylinders unless fixed.
                        dperp = dpars(k)*(1-fs(j));
                        if(fix(3))
                            dperp = fixedvals(3);
                        end
                        pars = [fs(j) dpars(k) dperp Rs(i) fisos(l) disos(m)];
                        grid(:,ind) = pars;
                        ind = ind + 1;
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'CylSingleRadTortIsoV_GPD') || strcmp(model, 'CylSingleRadTortIsoV_GPD_B0'))
    if(fix(3))
        Rs = [fixedvals(3)];
    end
    if(fix(4))
        fisos = [fixedvals(4)];
    end
    if(fix(5))
        disos = [fixedvals(5)];
    end

    numCombs = length(Rs)*length(fs)*length(dpars)*length(fisos)*length(disos);
    grid = zeros(5, numCombs);
    ind = 1;
    for i=1:length(Rs)
        for j=1:length(fs)
            for k=1:length(dpars)
                for l=1:length(fisos)
                    for m=1:length(disos)
                        pars = [fs(j) dpars(k) Rs(i) fisos(l) disos(m)];
                        grid(:,ind) = pars;
                        ind = ind + 1;
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'CylSingleRadIsoDotGPD'))
    if(fix(4))
        Rs = [fixedvals(4)];
    end
    if(fix(5))
        irfracs = [fixedvals(5)];
    end
    numCombs = length(Rs)*length(fs)*length(dpars)*length(irfracs);
    grid = zeros(5, numCombs);
    ind = 1;
    for i=1:length(Rs)
        for j=1:length(fs)
            for k=1:length(dpars)
                for l=1:length(irfracs)
                    dperp = dpars(k)*(1-fs(j));
                    pars = [fs(j) dpars(k) dperp Rs(i) irfracs(l)];
                    grid(:,ind) = pars;
                    ind = ind + 1;
                end
            end
        end
    end
elseif(strcmp(model, 'CylSingleRadIsoResTortIsoV_GPD') || strcmp(model, 'CylSingleRadIsoResTortIsoV_GPD_B0')...
        || strcmp(model, 'CylSingleRadIsoStickTortIsoV_GPD') || strcmp(model, 'CylSingleRadIsoStickTortIsoV_GPD_B0')...
        || strcmp(model, 'CylSingleRadIsoSphereTortIsoV_GPD') || strcmp(model, 'CylSingleRadIsoSphereTortIsoV_GPD_B0')...
        || strcmp(model, 'CylSingleRadIsoDotTortIsoV_GPD') || strcmp(model, 'CylSingleRadIsoDotTortIsoV_GPD_B0'))
    if(fix(3))
        Rs = [fixedvals(3)];
    end
    if(fix(4))
        irfracs = [fixedvals(4)];
    end
    if(fix(5))
        fisos = [fixedvals(5)];
    end
    if(fix(6))
        disos = [fixedvals(6)];
    end

    numCombs = length(Rs)*length(fs)*length(dpars)*length(fisos)*length(disos);
    grid = zeros(6, numCombs);
    ind = 1;
    for i=1:length(Rs)
        for j=1:length(fs)
            for k=1:length(dpars)
                for l=1:length(fisos)
                    for m=1:length(disos)
                        for n=1:length(irfracs)
                            pars = [fs(j) dpars(k) Rs(i) irfracs(n) fisos(l) disos(m)];
                            grid(:,ind) = pars;
                            ind = ind + 1;
                        end
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'Stick'))
    numCombs = length(fs)*length(dpars);
    grid = zeros(3, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
            % Set dperp using the standard tortuosity model for
            % randomly placed cylinders unless fixed.
            dperp = dpars(j)*(1-fs(i));
            if(fix(3))
                dperp = fixedvals(3);
            end
            pars = [fs(i) dpars(j) dperp];
            grid(:,ind) = pars;
            ind = ind + 1;
        end
    end
elseif(strcmp(model, 'StickIsoV_B0'))
    if(fix(4))
        fisos = [fixedvals(4)];
    end
    if(fix(5))
        disos = [fixedvals(5)];
    end

    numCombs = length(fs)*length(dpars)*length(fisos)*length(disos);
    grid = zeros(5, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
            for k=1:length(fisos)
                for l=1:length(disos)
                    % Set dperp using the standard tortuosity model for
                    % randomly placed cylinders unless fixed
                    dperp = dpars(j)*(1-fs(i));
                    if(fix(3))
                        dperp = fixedvals(3);
                    end
                    pars = [fs(i) dpars(j) dperp fisos(k) disos(l)];
                    grid(:,ind) = pars;
                    ind = ind + 1;
                end
            end
        end
    end
elseif(strcmp(model, 'StickTortIsoV_B0'))
    if(fix(3))
        fisos = [fixedvals(3)];
    end
    if(fix(4))
        disos = [fixedvals(4)];
    end

    numCombs = length(fs)*length(dpars)*length(fisos)*length(disos);
    grid = zeros(4, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
            for k=1:length(fisos)
                for l=1:length(disos)
                    pars = [fs(i) dpars(j) fisos(k) disos(l)];
                    grid(:,ind) = pars;
                    ind = ind + 1;
                end
            end
        end
    end
elseif(strcmp(model, 'WatsonStick') || strcmp(model, 'WatsonSHStick'))
    numCombs = length(fs)*length(dpars)*length(kappas);
    grid = zeros(4, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
		  for k=1:length(kappas)
                % Set dperp using the standard tortuosity model for
                % randomly placed cylinders unless fixed.
                dperp = dpars(j)*(1-fs(i));
                if(fix(3))
                    dperp = fixedvals(3);
                end
                pars = [fs(i) dpars(j) dperp kappas(k)];
                grid(:,ind) = pars;
                ind = ind + 1;
            end
        end
    end
elseif(strcmp(model, 'WatsonSHStickIsoV_B0'))
    if(fix(4))
        kappas = [fixedvals(4)];
    end
    if(fix(5))
        fisos = [fixedvals(5)];
    end
    if(fix(6))
        disos = [fixedvals(6)];
    end

    numCombs = length(fs)*length(dpars)*length(fisos)*length(disos)*length(kappas);
    grid = zeros(6, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
            for k=1:length(fisos)
                for l=1:length(disos)
                    for m=1:length(kappas)
                      dperp = dpars(j)*(1-fs(i));
                      if(fix(3))
                          dperp = fixedvals(3);
                      end
                      pars = [fs(i) dpars(j) dperp kappas(m) fisos(k) disos(l)];
                      grid(:,ind) = pars;
                      ind = ind + 1;
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'WatsonSHStickIsoVIsoDot_B0'))
    if(fix(5))
        fisos = [fixedvals(5)];
    end
    if(fix(6))
        disos = [fixedvals(6)];
    end
    if(fix(7))
        irfracs = [fixedvals(7)];
    end

    numCombs = length(fs)*length(dpars)*length(fisos)*length(disos)*length(kappas)*length(irfracs);
    grid = zeros(7, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
            for k=1:length(fisos)
                for l=1:length(disos)
                    for m=1:length(kappas)
                        for n=1:length(irfracs)
                          dperp = dpars(j)*(1-fs(i));
                          if(fix(3))
                             dperp = fixedvals(3);
                          end
                          pars = [fs(i) dpars(j) dperp kappas(m) fisos(k) disos(l) irfracs(n)];
                          grid(:,ind) = pars;
                          ind = ind + 1;
                        end
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'WatsonStickTort') || strcmp(model, 'WatsonSHStickTort'))
    numCombs = length(fs)*length(dpars)*length(kappas);
    grid = zeros(3, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
		  for k=1:length(kappas)
                pars = [fs(i) dpars(j) kappas(k)];
                grid(:,ind) = pars;
                ind = ind + 1;
            end
        end
    end
elseif(strcmp(model, 'WatsonSHStickTortIsoV_B0'))
    if(fix(4))
        fisos = [fixedvals(4)];
    end
    if(fix(5))
        disos = [fixedvals(5)];
    end

    numCombs = length(fs)*length(dpars)*length(fisos)*length(disos)*length(kappas);
    grid = zeros(5, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
            for k=1:length(fisos)
                for l=1:length(disos)
                    for m=1:length(kappas)
                      pars = [fs(i) dpars(j) kappas(m) fisos(k) disos(l)];
                      grid(:,ind) = pars;
                      ind = ind + 1;
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'WatsonSHStickTortIsoVIsoDot_B0'))
    if(fix(4))
        fisos = [fixedvals(4)];
    end
    if(fix(5))
        disos = [fixedvals(5)];
    end
    if(fix(6))
        irfracs = [fixedvals(6)];
    end

    numCombs = length(fs)*length(dpars)*length(fisos)*length(disos)*length(kappas)*length(irfracs);
    grid = zeros(6, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
            for k=1:length(fisos)
                for l=1:length(disos)
                    for m=1:length(kappas)
                        for n=1:length(irfracs)
                          pars = [fs(i) dpars(j) kappas(m) fisos(k) disos(l) irfracs(n)];
                          grid(:,ind) = pars;
                          ind = ind + 1;
                        end
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'WatsonSHCylSingleRadTortIsoV_GPD') || strcmp(model, 'WatsonSHCylSingleRadTortIsoV_GPD_B0'))
    if(fix(3))
        Rs = [fixedvals(3)];
    end
    if(fix(5))
        fisos = [fixedvals(5)];
    end
    if(fix(6))
        disos = [fixedvals(6)];
    end

    numCombs = length(Rs)*length(fs)*length(dpars)*length(fisos)*length(disos)*length(kappas);
    grid = zeros(6, numCombs);
    ind = 1;
    for i=1:length(Rs)
        for j=1:length(fs)
            for k=1:length(dpars)
                for l=1:length(fisos)
                    for m=1:length(disos)
                      for n=1:length(kappas)
                        pars = [fs(j) dpars(k) Rs(i) kappas(n) fisos(l) disos(m)];
                        grid(:,ind) = pars;
                        ind = ind + 1;
                      end
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'BinghamCylSingleRadTortIsoV_GPD_B0'))
    if(fix(3))
        Rs = [fixedvals(3)];
    end
    if(fix(7))
        fisos = [fixedvals(7)];
    end
    if(fix(8))
        disos = [fixedvals(8)];
    end

    numCombs = length(Rs)*length(fs)*length(dpars)*length(fisos)*length(disos)*length(kappas);
    grid = zeros(8, numCombs);
    ind = 1;
    for i=1:length(Rs)
        for j=1:length(fs)
            for k=1:length(dpars)
                for l=1:length(fisos)
                    for m=1:length(disos)
                        for n=1:length(kappas)
                            pars = [fs(j) dpars(k) Rs(i) kappas(n) 0 0 fisos(l) disos(m)];
                            grid(:,ind) = pars;
                            ind = ind + 1;
                        end
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'BinghamStickTortIsoV_B0'))
    if(fix(6))
        fisos = [fixedvals(6)];
    end
    if(fix(7))
        disos = [fixedvals(7)];
    end

    numCombs = length(fs)*length(dpars)*length(fisos)*length(disos)*length(kappas);
    grid = zeros(7, numCombs);
    ind = 1;
    for i=1:length(fs)
        for j=1:length(dpars)
            for k=1:length(fisos)
                for l=1:length(disos)
                    for m=1:length(kappas)
                        pars = [fs(i) dpars(j) kappas(m) 0 0 fisos(k) disos(l)];
                        grid(:,ind) = pars;
                        ind = ind + 1;
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'ExCrossingCylSingleRadGPD'))
    R1s = Rs;
    R2s = Rs;
    if(fix(4))
        R1s = [fixedvals(4)];
    end

    if(fix(5))
        R2s = [fixedvals(5)];
    end

    numCombs = length(R1s)*length(R2s)*length(fs)*length(dpars)*length(fic);
    grid = zeros(6, numCombs);
    ind = 1;
    for i=1:length(R1s)
        for j=1:length(R2s)
            for k=1:length(fs)
                for l=1:length(dpars)
                    for m=1:length(fic)
                        % Set dperp using the standard tortuosity model for
                        % randomly placed cylinders unless fixed.
                        dperp = dpars(l)*(1-fs(k));
                        if(fix(3))
                            dperp = fixedvals(3);
                        end
                        pars = [fs(k) dpars(l) dperp R1s(i) R2s(j) fic(m)];
                        grid(:,ind) = pars;
                        ind = ind + 1;
                    end
                end
            end
        end
    end
elseif(strcmp(model, 'ExCrossingCylSingleRadIsoDotTortIsoV_GPD_B0'))
    R1s = Rs;
    R2s = Rs;
    if(fix(3))
        R1s = [fixedvals(3)];
    end

    if(fix(4))
        R2s = [fixedvals(4)];
    end

    if(fix(6))
        irfracs = [fixedvals(6)];
    end

    if(fix(7))
        fisos = [fixedvals(7)];
    end

    if(fix(8))
        disos = [fixedvals(8)];
    end

    numCombs = length(R1s)*length(R2s)*length(fs)*length(dpars)*length(fic)*length(irfracs)*length(fisos)*length(disos);
    grid = zeros(8, numCombs);
    ind = 1;
    for i=1:length(R1s)
        for j=1:length(R2s)
            for k=1:length(fs)
                for l=1:length(dpars)
                    for m=1:length(fic)
								 for n=1:length(irfracs)
										for p=1:length(fisos)
											 for q=1:length(disos)
						  				  		  pars = [fs(k) dpars(l) R1s(i) R2s(j) fic(m) irfracs(n) fisos(p) disos(q)];
						  				  		  grid(:,ind) = pars;
						  				  		  ind = ind + 1;
									 		 end
										end
								 end
                    end
                end
            end
        end
    end
else
    error(['Starting combinations not implemented for model: ', model]);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [inchiString] = getInchiStringFromHMDB(HMDBID)
% This function retrieves the inchiString from HMDB (online) for a given
% HMDB ID.
%
% INPUT
% HMDBID    Human metabolome database (HMDB) ID
%
% OUTPUT
% inchiString   Retrieved inchiString
%
% Ines Thiele, 09/2021
%

% get inchi from HMDB
try
    url=strcat('https://hmdb.ca/metabolites/',HMDBID);
    syst = urlread(url);
    [tok] = split(syst,'InChI=1S');
    tok2= split(tok{2},'INCHI_KEY');
    tok3 = regexprep(tok2{1},'\n\n\> \<','');
    inchiString = ['InChI=1S' tok3];
catch
    inchiString =[];
end

author: opencobra

repo: cobratoolbox


text:
function f = local ( m, x )

%*****************************************************************************80
%
%% LOCAL computes the local function.
%
%  Discussion:
%
%    This function has a local minimizer:
%
%      X* = ( 0.2858054412..., 0.2793263206...), F(X*) = 5.9225...
%
%    and a global minimizer:
%
%      X* = ( -21.02667179..., -36.75997872...), F(X*) = 0.
%
%    Suggested starting point for local minimizer:
%
%      X = ( 1, 1 ), F(X) = 3.33 * 10^6.
%
%    Suggested starting point for global minimizer:
%
%      X = ( -15, -35), F(X) = 1.49 * 10^8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Himmelblau,
%    Applied Nonlinear Programming,
%    McGraw Hill, 1972,
%    ISBN13: 978-0070289215,
%    LC: T57.8.H55.
%
%  Parameters:
%
%    Input, integer M, the number of variables.
%
%    Input, real X(M), the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  f = ( x(1)^2 + 12.0 * x(2) - 1.0 )^2 ...
    + ( 49.0 * x(1)^2 + 49.0 * x(2)^2 + 84.0 * x(1) + 2324.0 * x(2) - 681.0 )^2;

  return
end


text:
function r=ksrmv(x,y,hx,z)
% KSRMV   Multivariate kernel smoothing regression
%
% r=ksrmv(x,y) returns the Gaussian kernel regression in structure r such that
%   r.f(r.x) = y(x) + e
% The bandwidth and number of samples are also stored in r.h and r.n
% respectively.
%
% r=ksrmv(x,y,h) performs the regression using the specified bandwidth, h.
%
% r=ksrmv(x,y,h,z) calculates the regression at location z (default z=x).
%
% Algorithm
% The kernel regression is a non-parametric approach to estimate the
% conditional expectation of a random variable:
%
% E(Y|X) = f(X)
%
% where f is a non-parametric function. Based on the kernel density
% estimation, this code implements the Nadaraya-Watson kernel regression
% using the Gaussian kernel as follows:
%
% f(x) = sum(kerf((x-X)/h).*Y)/sum(kerf((x-X)/h))
%
% See also gkde, ksdensity, ksr

% Example 1: smoothing a noised logo
%{
L = 40*membrane(1,25)+randn(51);
[x,y]=meshgrid(0:50);
r=ksrmv([x(:) y(:)],L(:));
Lr=L;
Lr(:)=r.f;
subplot(121), surf(x,y,L)
subplot(122), surf(x,y,Lr)
%}
% Example 2: smoothing noised peaks with 20% missing data
%{
L = 10*peaks(50)+randn(50);
I = ceil(rand(500,1)*2500);
L(I) = NaN;
[x,y]=meshgrid(1:50);
r=ksrmv([x(:) y(:)],L(:));
Lr=L;
Lr(:)=r.f;
subplot(121), surf(x,y,L)
subplot(122), surf(x,y,Lr)
%}

% By Yi Cao at Cranfield University on 20 March 2008.
%

% Check input and output
error(nargchk(2,4,nargin));
error(nargoutchk(0,1,nargout));
y=y(:);
if size(x,1)~=size(y,1)
    error('x and y have different rows.');
end
d=size(x,2);

% Default parameters
if nargin<4
    z=x;
elseif size(z,2)~=d
    error('z must have the same number of columns as x.')
end
r.x=z;
N=size(z,1);

% clean missing or invalid data points
inv=(y~=y);
x(inv,:)=[];
y(inv)=[];
r.n=numel(y);

if nargin<3
    % optimal bandwidth suggested by Bowman and Azzalini (1997) p.31
    hy=median(abs(y-median(y)))/0.6745*(4/(d+2)/r.n)^(1/(d+4));
    hx=median(abs(x-repmat(median(x),r.n,1)))/0.6745*(4/(d+2)/r.n)^(1/(d+4));
    hx=sqrt(hy*hx);
elseif size(hx,2)~=d
    error('h must be a scalar.')
end
r.h=hx;
%
% % Vectorization, Not suitable for large data set
% % Firstly, diagnal matrix for bandwidth
% H=diag(1./hx);
%
% % Then scale X by the bendwidth
% X=r.x*H;
%
% % Square of X
% X2=sum(X.*X,2)/2;
%
% % Scale Xi
% C=H*x';
%
% % Square of Xi
% C2=sum(C.*C)/2;
%
% % D = ((X - Xi) / H)^T ((X - Xi) / H)
% D=C2(ones(N,1),:)+X2(:,ones(1,r.n))-X*C;
%
% % only elements whos values < 12 need to evaluated since exp(-12) = 6.1e-6
% idx=D<12;
%
% % prepare local variable
% Z=zeros(N,r.n);
%
% % Z = exp(-D)
% Z(idx)=exp(-D(idx));
%
% % f(x) = sum Z*Yi / sum Z
% r.f=sum(Z.*y(:,ones(1,N))',2)./sum(Z,2);

% Improved efficient code

% Scaling first
H=diag(1./hx);
x=x*H;
x1=r.x*H;

% Gaussian kernel function
kerf=@(z)exp(-sum(z.*z,2)/2);

% allocate memory
r.f=zeros(N,1);

% Loop through each regression point
for k=1:N
    % scaled deference from regression point
    xx=abs(x-x1(k+zeros(r.n,1),:));
    % select neighbours using exp(-5^2/2)<5e-6
    idx=all(xx<5,2);
    % kernel function
    z=kerf(xx(idx,:));
    % regression
    r.f(k)=sum(z.*y(idx))/sum(z);
end


text:
function [IM] = InterpMatrix2D(rout, sout)

% function [IM] = InterpMatrix2D(rout, sout)
% Purpose: Compute local elemental interpolation matrix

Globals2D;

% compute Vandermonde at (rout,sout)
Vout = Vandermonde2D(N, rout, sout);

% build interpolation matrix
IM = Vout*invV;
return


text:
% Applies a Hanning window of specified size in sec, to the log-likelihood (LLt) resulting in llt_smooth and then
% converts it to probability of models v_smooth.
% Inputs:
%   sRate          - Sampling rate of data (EEG.srate)
%   LLt            - Log-likelihood of data
%   smoothlength   - Length of smoothing window in sec.

function [v_smooth llt_smooth] = smooth_amica_prob(sRate,LLt, smoothlength)

if nargin<2
    smoothlength = 2; % default smoothing length (in sec.);
end

% Hanning smoothing window generation
smoothwnd = hann(round(smoothlength*sRate)); smoothwnd = smoothwnd/sum(smoothwnd);

if size(LLt,3) == 1

    llt_smooth = filtfilt_fast(fastif(smoothlength == 0, 1,smoothwnd),1,LLt')';
    v_smooth = LLt2v(llt_smooth);

else
    llt = reshape(LLt,size(LLt,1),size(LLt,2)*size(LLt,3));
    llt_smooth = filtfilt_fast(fastif(smoothlength == 0, 1,smoothwnd),1,llt')';
    v_smooth = LLt2v(llt_smooth);
    llt_smooth = reshape(llt_smooth,size(LLt,1),size(LLt,2),size(LLt,3));
    v_smooth = reshape(v_smooth,size(LLt,1),size(LLt,2),size(LLt,3));
end

author: PatternRecognition

repo: OpenBMI


text:
function fixedtrx2ann(trx,oldannname,newannname)

% open the annotation file for reading and writing
fin = fopen(oldannname,'rb');
if exist(newannname,'file'),
  b = questdlg(sprintf('File %s exists. Overwrite?',newannname));
  if ~strcmpi(b,'yes'),
    return;
  end
end
fout = fopen(newannname,'wb');

% find the end of the header
while true,
  [s,param,value] = read_ann_line(fin);
  if ~ischar(s) && s == -1,
    error('End of file encountered before end of header found');
  end
  if ~ischar(value) && length(value) > 1,
    fprintf(fout,'%s\n',s);
    fwrite(fout,value,'uint8');
  else
    fprintf(fout,'%s\n',s);
  end
  if strcmp(s,'end header'),
    break;
  end
end

if ftell(fin) ~= ftell(fout),
  error('input and output file do not match');
end

fclose(fin);

nframes = max([trx.endframe]);
nflies = length(trx);

for i = 1:nframes,
  if mod(i,100) == 0,
    fprintf('Writing frame %d\n',i);
  end
  for j = 1:nflies,
    if trx(j).firstframe > i || trx(j).endframe < i,
      continue;
    end
    k = i - trx(j).firstframe + 1;
    fprintf(fout,'%f\t%f\t%f\t%f\t%f\t%d\t',...
      trx(j).x(k)-1,trx(j).y(k)-1,trx(j).b(k),trx(j).a(k),trx(j).theta(k),j-1);
  end
  fprintf(fout,'\n');
end

fclose(fout);

if 0,
  fclose('all');
  fin = fopen(oldannname,'rb');
  fout = fopen(newannname,'rb');
  while true,
[sin,paramin,valuein] = read_ann_line(fin);
[sout,paramout,valueout] = read_ann_line(fout);
if ~strcmp(sin,sout),
fprintf('sin does not match sout\n');
keyboard;
else
fprintf([sin,'\n']);
end
if strcmp(sin,'end header'),
   if ~strcmp(sout,'end header'),
     fprintf('end of headerin but not headerout\n');
     keyboard;
   end
   break;
elseif strcmp(sout,'end header'),
  fprintf('end of headerout but not headerin\n');
  keyboard;
  break;
end
if isempty(paramin) && isempty(paramout),
  fprintf('both empty\n');
  continue;
end
if ~strcmp(paramin,paramout),
fprintf('paramin does not match paramout\n');
keyboard;
else
fprintf([paramin,'\n']);
end
if ischar(valuein),
if ~ischar(valueout),
fprintf('ischar(valuein) but not valueout\n');
keyboard;
end
if ~strcmp(valuein,valueout),
fprintf('valuein does not match valueout\n');
end
elseif length(valuein) == 1,
if ischar(valueout),
fprintf('ischar(valueout) but not valuein\n');
keyboard;
end
if length(valueout) ~= 1,
fprintf('length(valuein) == 1, but not valueout\n');
keyboard;
end
if valuein ~= valueout,
fprintf('valuein and valueout do not match\n');
keyboard;
else
fprintf([num2str(value),'\n']);
end
else
if ischar(valueout),
fprintf('ischar(valueout) but not valuein\n');
keyboard;
end

if length(valueout) == 1,
fprintf('length(valueout) == 1 but not valuein\n');
keyboard;
end
if length(valuein) ~= length(valueout),
fprintf('lengths do not match\n');
keyboard;
end
if ~all(valuein == valueout),
fprintf('not all values match\n');
keyboard;
else
fprintf('value vectors of length %d match\n',length(valuein));
end
end
end
end

author: kristinbranson

repo: JAABA


text:
%% Quasi symmetry
%
%%
% By quasi crystals we mean materials with a non crystallographic symmetry
% group, e.g. including a pentagonal or decagonal symmetry axis. Such
% symmetry groups do not lead to perfect translational invariant lattices
% and hence are not crystals.
%
% Setting up a non standard symmetry group can be done by specifying all
% its symmetry elements. To illustrate this definition lets consider the
% symmetry group corresponding to the dodecaeder which is generated by a 5
% fold and a 3 fold symmetry axis.

rot2 = rotation.byAxisAngle(vector3d.Z,180*degree);
rot5 = rotation.byAxisAngle(vector3d.byPolar(31.7171*degree,0*degree),72*degree);
%rot3 = rotation.byAxisAngle(vector3d.byPolar(20.9054*degree,90*degree),120*degree);

cs = crystalSymmetry.byElements([rot5,rot2])

plot(cs,'symbolSize',0.5,'projection','eangle','grid','on')

%%
% Many of the MTEX methods do work also with quasi symmetries. E.g. we can
% identify the fundamental sector by

hold on
plot(cs.fundamentalSector,'color','red')
hold off

%%
% the fundamental region is orientation space, which is exactly a
% dodecaeder,

oR = cs.fundamentalRegion;
plot(oR)
axis off

%%
% or define an approbiate color key

ipfKey = ipfHSVKey(cs);

plot(ipfKey,'complete','upper','resolution',0.5*degree)
hold on
plot(cs,'SymbolSize',0.6,'linewidth',1)
hold off

%%
% We may also choose a different setup where the 5 fold axis is aligned
% parallel to z

rot5 = rotation.byAxisAngle(zvector,72*degree);
rot3 = rotation.byAxisAngle(vector3d('polar',37.377*degree,0) ,120*degree);

cs = crystalSymmetry.byElements([rot5,rot3])

plot(cs,'symbolSize',0.5,'projection','eangle','grid','on')

%%
% or include the inversion

cs = crystalSymmetry.byElements([rot5,rot3,rotation.inversion])

ipfKey = ipfHSVKey(cs);

plot(ipfKey,'complete','upper','resolution',0.5*degree)

hold on
plot(cs,'symbolSize',0.6,'linewidth',2)
hold off

author: mtex-toolbox

repo: mtex

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


text:
% Make new eyeCatch library

% find ori_map on eyeCatch webiste, filename:
% 'eyeChannelWeightNormalizedpart1.mat' and 'eyeChannelWeightNormalizedpart1.mat'

function new_map = interpo_lib(ori_map, downsize)

[sampleNumber, oldsize] = size(ori_map);
if downsize > sqrt(oldsize)
    error('Downsize map should be smaller than original map.');
end

oldsize = sqrt(oldsize);

% reshape old_map
ori_map = reshape(ori_map, [sampleNumber, oldsize, oldsize]);
new_map = zeros(sampleNumber, downsize^2);

%% resize
for i = 1:sampleNumber
    new_map(i,:) = reshape(imresize(squeeze(ori_map(i,:,:)),[32 32],'nearest'), [1, downsize^2]);
end

% save to mat
filename = 'libEyeCatch';

% for oversize library
% k = 1;
% filename = ['new_lib',int2str(k)];
% while exist([filename,'.mat'],'file')
%     k = k+1;
%     filename = ['new_lib.mat', int2str(k)];
% end

save(filename, 'new_map');

end

repo: REST


text:
function [G]=gsp_2dgrid(N,M)
%GSP_2dgrid  Initialize a 2 dimentional grid graph
%   Usage:  G=gsp_path(N);
%
%   Input parameters:
%         N     : Number of vertices along the first dimention (default 16)
%         M     : Number of vertices along the second dimention (default N)
%   Output parameters:
%         G     : Graph structure.
%
%   The 2d grid graph correspond the  graph used for the DCT2????
%
%   Example:::
%
%          G = gsp_2dgrid(16);
%          param.show_edges = 1;
%          gsp_plot_graph(G,param);
%
%   See also: gsp_ring, gsp_path
%
%   References: strang1999discrete

% Author: Nathanael Perraudin
% Date: 15 March 2014
% Testing: test_graph

if nargin <1
   N = 16;
end

if nargin < 2
   M = N;
end

G.N=N*M;

% Create weighted adjancency matrix
K = 2*(N-1);
J = 2*(M-1);
i_inds = zeros(K*M+J*N,1);
j_inds = zeros(K*M+J*N,1);
for ii = 1:M
    i_inds((ii-1)*K+(1:K)) = (ii-1)*N+[1:N-1,2:N]';
    j_inds((ii-1)*K+(1:K)) = (ii-1)*N+[2:N,1:N-1]';
end

for ii = 1:M-1
    i_inds(K*M+(ii-1)*2*N+(1:2*N)) = [((ii-1)*N+(1:N)),(ii*N+(1:N))]';
    j_inds(K*M+(ii-1)*2*N+(1:2*N)) = [(ii*N+(1:N)),((ii-1)*N+(1:N))]';
end

G.W = sparse(i_inds,j_inds,ones(K*M+J*N,1),N*M,N*M);

% Create coordinates
G.coords = [repmat((0:(N-1))'/N,M,1),...
        reshape(repmat((0:(M-1))/M,N,1),M*N,1) ];
G.plotting.limits = [-1/N, 1+1/N, -1/M, 1+1/M];

G.type = '2d-grid';
G.plotting.vertex_size = 30;

G = gsp_graph_default_parameters(G);

end

repo: gspbox


text:
function libdir = opti_GetLibPath()
% Return architecture dependent general library path

switch(computer)
    case 'PCWIN'
        libdir = 'lib\win32\ ';
    case 'PCWIN64'
        libdir = 'lib\win64\ ';
    otherwise
        error('This function is only setup for Windows PCs');
end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function obj = add_matrix(obj, mat3d, mat, imgno)
% Add 3d matrix image vol to slice overlay
% FORMAT obj = add_matrix(obj, mat3d, mat, imgno)
%
% Inputs
% obj          - object
% mat3d        - 3D matrix to add as img
% mat          - optional 4x4 voxel->world translation
% imgno        - optional img no to add to (defaults to last in object)
%
% Ouputs
% obj          - modified object
%__________________________________________________________________________

% Matthew Brett
% $Id: add_matrix.m 6623 2015-12-03 18:38:08Z guillaume $

if nargin < 2
    return
end
if nargin < 3
    mat = [];
end
if nargin < 4
    imgno = [];
end
if isempty(imgno)
    imgno = length(obj.img);
end
if ~isempty(mat3d)
    obj.img(imgno).vol = pr_matrix2vol(mat3d, mat);
end

author: spm


text:
function bs = bct_toolbox_undirected_graph_metrics(bs, thresh, varargin)
% bs = bct_toolbox_undirected_graph_metrics(bs)
%
% thresh = proportional link density threshold; 0 to 1 value. (.1 is a common value)
%
% Method for brainpathway_multisubject that extracts graph metrics using
% the Sporns lab (et al.) BCT toolbox.
% Uses bct_toolbox_undirected_graph_metrics
%
% Stores results in bs.graph_properties.regions
% - Nodes x subjects matrix for each graph property
% - Can do group stats/t-tests/etc on (matrices)'
%
% Examples:
% ------------------------------------------------------------------
% Load brainpathway_data_gsr_censoring_pipeline.mat % from OLP/Yoni
% bs = bct_toolbox_undirected_graph_metrics(bs);
% mean_degree = nanmean(bs.graph_properties.regions.degree')'; % mean degree for each node
%
% % Plot montage of mean degree on brain
% fmri_dat_obj = brainpathway2fmri_data(bs, bs.graph_properties.regions.degree);
% montage(mean(fmri_dat_obj));

r = double(bs.connectivity.regions.r);

n = size(r, 3); % subjects
n_nodes = size(r, 1); % nodes

bs.graph_properties.regions = table(); % clear out -- overwrite

fprintf('Completed  subject        ');

warning off
for i = 1:n

    [graph_prop, graph_prop_glob] = bct_toolbox_undirected_graph_metrics(r(:, :, i), thresh, varargin{:});

    % nodal properties
    vnames = graph_prop.Properties.VariableNames;
    for j = 1:length(vnames)

        bs.graph_properties.regions.(vnames{j})(i, 1:n_nodes) = graph_prop.(vnames{j});

    end

    % global properties
    vnames = graph_prop_glob.Properties.VariableNames;
    for j = 1:length(vnames)

        bs.graph_properties.regions.(vnames{j})(i, 1) = graph_prop_glob.(vnames{j});

    end

    % print to update status
    fprintf('\b\b\b\b\b\b%5d\n',i)

end % subject loop
warning on

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function sgmga_index_test ( dim_num, importance, level_weight, ...
  level_max_min, level_max_max, rule, growth, np, p, tol )

%****************************************************************************80
%
%% SGMGA_INDEX_TEST tests SGMGA_INDEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real IMPORTANCE(DIM_NUM), the anisotropic importance of
%    each dimension.
%
%    Input, real LEVEL_WEIGHT(DIM_NUM), the anisotropic weight
%    for each dimension.
%
%    Input, integer LEVEL_MAX_MIN, LEVEL_MAX_MAX, the minimum and
%    maximum values of LEVEL_MAX.
%
%    Input, integer RULE(DIM_NUM), the rule in each dimension.
%     1, "CC",  Clenshaw Curtis, Closed Fully Nested.
%     2, "F2",  Fejer Type 2, Open Fully Nested.
%     3, "GP",  Gauss Patterson, Open Fully Nested.
%     4, "GL",  Gauss Legendre, Open Weakly Nested.
%     5, "GH",  Gauss Hermite, Open Weakly Nested.
%     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested.
%     7, "LG",  Gauss Laguerre, Open Non Nested.
%     8, "GLG", Generalized Gauss Laguerre, Open Non Nested.
%     9, "GJ",  Gauss Jacobi, Open Non Nested.
%    10, "HGK", Hermite Genz-Keister, Open Fully Nested.
%    11, "UO",  User supplied Open, presumably Non Nested.
%    12, "UC",  User supplied Closed, presumably Non Nested.
%
%    Input, integer GROWTH(DIM_NUM), the growth in each dimension.
%    0, "DF", default growth associated with this quadrature rule;
%    1, "SL", slow linear, L+1;
%    2  "SO", slow linear odd, O=1+2((L+1)/2)
%    3, "ML", moderate linear, 2L+1;
%    4, "SE", slow exponential;
%    5, "ME", moderate exponential;
%    6, "FE", full exponential.
%
%    Input, integer NP(DIM_NUM), the number of parameters used by each rule.
%
%    Input, real P(*), the parameters needed by each rule.
%
%    Input, real TOL, a tolerance for point equality.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_INDEX_TEST\n' );
  fprintf ( 1, '  SGMGA_INDEX returns index and order vectors that\n' );
  fprintf ( 1, '  identify each point in a multidimensional sparse grid with mixed factors.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Each sparse grid is of spatial dimension DIM_NUM,\n' );
  fprintf ( 1, '  and is made up of product grids of levels up to LEVEL_MAX.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  IMPORTANCE:    ');
  for dim = 1 : dim_num
    fprintf ( 1, '  %14f', importance(dim) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL_WEIGHT:  ');
  for dim = 1 : dim_num
    fprintf ( 1, '  %14f', level_weight(dim) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ' Dimension      Rule       Growth     Parameters\n' );
  fprintf ( 1, '\n' );

  p_index = 1;

  for dim = 1 : dim_num

    if ( rule(dim) == 1 )
      fprintf ( 1, '  %8d  %8d  %11d\n', dim, rule(dim), growth(dim) );
    elseif ( rule(dim) ==  2 )
      fprintf ( 1, '  %8d  %8d  %11d\n', dim, rule(dim), growth(dim) );
    elseif ( rule(dim) ==  3 )
      fprintf ( 1, '  %8d  %8d  %11d\n', dim, rule(dim), growth(dim) );
    elseif ( rule(dim) ==  4 )
      fprintf ( 1, '  %8d  %8d  %11d\n', dim, rule(dim), growth(dim) );
    elseif ( rule(dim) ==  5 )
      fprintf ( 1, '  %8d  %8d  %11d\n', dim, rule(dim), growth(dim) );
    elseif ( rule(dim) ==  6 )
      fprintf ( 1, '  %8d  %8d  %11d  %14f\n', dim, rule(dim), growth(dim), p(p_index) );
    elseif ( rule(dim) ==  7 )
      fprintf ( 1, '  %8d  %8d  %11d\n', dim, rule(dim), growth(dim) );
    elseif ( rule(dim) ==  8 )
      fprintf ( 1, '  %8d  %8d  %11d  %14f\n', dim, rule(dim), growth(dim), p(p_index) );
    elseif ( rule(dim) ==  9 )
      fprintf ( 1, '  %8d  %8d  %11d  %14f  %14f\n', ...
        dim, rule(dim), growth(dim), p(p_index), p(p_index+1) );
    elseif ( rule(dim) == 10 )
      fprintf ( 1, '  %8d  %8d  %11d\n', dim, rule(dim), growth(dim) );
    elseif ( rule(dim) ==  11 )
      fprintf ( 1, '  %8d  %8d  %11d\n', dim, rule(dim), growth(dim) );
    elseif ( rule(dim) ==  12 )
      fprintf ( 1, '  %8d  %8d  %11d\n', dim, rule(dim), growth(dim) );
    end

    p_index = p_index + np(dim);

  end

  for level_max = level_max_min : level_max_max

    point_total_num = sgmga_size_total ( dim_num, level_weight, level_max, ...
      rule, growth );

    point_num = sgmga_size ( dim_num, level_weight, level_max, rule, growth, ...
      np, p, tol );

    sparse_unique_index = sgmga_unique_index ( dim_num, level_weight, ...
      level_max, rule, growth, np, p, tol, point_num, point_total_num );

    [ sparse_order, sparse_index ] = sgmga_index ( dim_num, level_weight, ...
      level_max, rule, growth, point_num, point_total_num, sparse_unique_index );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  For LEVEL_MAX = %d\n', level_max );
    fprintf ( 1, '\n' );
    for point = 1 : point_num
      fprintf ( 1, '  %4d', point );
      for dim = 1 : dim_num
        fprintf ( 1, '  %3d / %3d', sparse_index(dim,point), sparse_order(dim,point) );
      end
      fprintf ( 1, '\n' );
    end

  end

  return
end


text:
function ns = count_squares_fisheye_distorted(I,x1,y1,x2,y2,win, fg, cg, kg);

[ny,nx] = size(I);

if ((x1-win <= 0) || (x1+win >= nx) || (y1-win <= 0) || (y1+win >= ny) || ...
        (x2-win <= 0) || (x2+win >= nx) || (y2-win <= 0) || (y2+win >= ny))
    ns = -1;
    return;
end;

if ((x1 - x2)^2+(y1-y2)^2) <  win,
    ns = -1;
    return;
end;

nX = round(sqrt((x1-x2)^2 + (y1-y2)^2));
alpha_x = (0:nX)/nX;

pt1n = normalize_pixel_fisheye([x1;y1]-1,fg,cg,kg,0);
pt2n = normalize_pixel_fisheye([x2;y2]-1,fg,cg,kg,0);

ptsn = repmat(pt1n,[1 nX+1]) + (pt2n - pt1n)*alpha_x;

pts = apply_fisheye_distortion(ptsn,kg);
pts(1,:) = fg(1)*pts(1,:) + cg(1);
pts(2,:) = fg(2)*pts(2,:) + cg(2);

% Check that the curve is within the image:
good_line = (min(pts(1,:))-win > 0) && (max(pts(1,:))+win < (nx-1)) && ...
    (min(pts(2,:))-win > 0) && (max(pts(2,:))+win <(ny-1));

if ~good_line,
    ns = -1;
    return;
end;

% Deviate the trajectory orthogonally:
lambda = [y1 - y2 ; x2 - x1];
lambda = lambda / sqrt(sum(lambda.^2));

Np = size(pts,2);
xs_mat = ones(2*win + 1,1)*pts(1,:);
ys_mat = ones(2*win + 1,1)*pts(2,:);
win_mat = (-win:win)'*ones(1,Np);
xs_mat2 = round(xs_mat - win_mat * lambda(1));
ys_mat2 = round(ys_mat - win_mat * lambda(2));
ind_mat = (xs_mat2) * ny + ys_mat2 + 1;
ima_patch = zeros(2*win + 1,Np);
ima_patch(:) = I(ind_mat(:));

filtk = [ones(win,Np);zeros(1,Np);-ones(win,Np)];
out_f = sum(filtk.*ima_patch);
out_f_f = conv2(out_f,[1/4 1/2 1/4],'same');
out_f_f = out_f_f(win+1:end-win);
ns = length(find(((out_f_f(2:end)>=0)&(out_f_f(1:end-1)<0)) | ((out_f_f(2:end)<=0)&(out_f_f(1:end-1)>0))))+1;

repo: ekfmonoslam


text:
function value = r8_scinvchi_sample ( df, s )

%*****************************************************************************80
%
%% R8_SCINVCHI_SAMPLE: sample a scaled inverse chi-squared distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real DF, the degrees of freedom.
%    0.0 < DF.
%
%    Input, real S, the scale factor.
%    0.0 < S.
%
%    Input, real VALUE, a sample value.
%
  a = 0.5 * df * s;
  b = 0.5 * df;

  value = r8_gamma_sample ( a, b );

  if ( value ~= 0.0 )
    value = 1.0 / value;
  end

  return
end


text:
function [tout,yout] = pendulumDE(tVals,omega0,startAngle)

%   Copyright 2008-2009 The MathWorks, Inc.

% Equations of motion
eqnsOfMotion = @(t,x) [x(2); -omega0^2*sin(x(1))];

[tout,yout] = ode45(eqnsOfMotion,tVals,[startAngle; 0]);


text:
function [image_roidb, bbox_means, bbox_stds] = fast_rcnn_prepare_image_roidb(conf, imdbs, roidbs, bbox_means, bbox_stds)
% [image_roidb, bbox_means, bbox_stds] = fast_rcnn_prepare_image_roidb(conf, imdbs, roidbs, cache_img, bbox_means, bbox_stds)
%   Gather useful information from imdb and roidb
%   pre-calculate mean (bbox_means) and std (bbox_stds) of the regression
%   term for normalization
% --------------------------------------------------------
% Fast R-CNN
% Reimplementation based on Python Fast R-CNN (https://github.com/rbgirshick/fast-rcnn)
% Copyright (c) 2015, Shaoqing Ren
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

    if ~exist('bbox_means', 'var')
        bbox_means = [];
        bbox_stds = [];
    end

    if ~iscell(imdbs)
        imdbs = {imdbs};
        roidbs = {roidbs};
    end

    imdbs = imdbs(:);
    roidbs = roidbs(:);

    image_roidb = ...
        cellfun(@(x, y) ... // @(imdbs, roidbs)
                arrayfun(@(z) ... //@([1:length(x.image_ids)])
                        struct('image_path', x.image_at(z), 'image_id', x.image_ids{z}, 'im_size', x.sizes(z, :), 'imdb_name', x.name, ...
                        'overlap', y.rois(z).overlap, 'boxes', y.rois(z).boxes, 'class', y.rois(z).class, 'image', [], 'bbox_targets', []), ...
                [1:length(x.image_ids)]', 'UniformOutput', true),...
        imdbs, roidbs, 'UniformOutput', false);

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

    % enhance roidb to contain bounding-box regression targets
    [image_roidb, bbox_means, bbox_stds] = append_bbox_regression_targets(conf, image_roidb, bbox_means, bbox_stds);
end

function [image_roidb, means, stds] = append_bbox_regression_targets(conf, image_roidb, means, stds)
    % means and stds -- (k+1) * 4, include background class

    num_images = length(image_roidb);
    % Infer number of classes from the number of columns in gt_overlaps
    num_classes = size(image_roidb(1).overlap, 2);
    valid_imgs = true(num_images, 1);
    for i = 1:num_images
       rois = image_roidb(i).boxes;
       [image_roidb(i).bbox_targets, valid_imgs(i)] = ...
           compute_targets(conf, rois, image_roidb(i).overlap);
    end
    if ~all(valid_imgs)
        image_roidb = image_roidb(valid_imgs);
        num_images = length(image_roidb);
        fprintf('Warning: fast_rcnn_prepare_image_roidb: filter out %d images, which contains zero valid samples\n', sum(~valid_imgs));
    end

    if ~(exist('means', 'var') && ~isempty(means) && exist('stds', 'var') && ~isempty(stds))
        % Compute values needed for means and stds
        % var(x) = E(x^2) - E(x)^2
        class_counts = zeros(num_classes, 1) + eps;
        sums = zeros(num_classes, 4);
        squared_sums = zeros(num_classes, 4);
        for i = 1:num_images
           targets = image_roidb(i).bbox_targets;
           for cls = 1:num_classes
              cls_inds = find(targets(:, 1) == cls);
              if ~isempty(cls_inds)
                 class_counts(cls) = class_counts(cls) + length(cls_inds);
                 sums(cls, :) = sums(cls, :) + sum(targets(cls_inds, 2:end), 1);
                 squared_sums(cls, :) = squared_sums(cls, :) + sum(targets(cls_inds, 2:end).^2, 1);
              end
           end
        end

        means = bsxfun(@rdivide, sums, class_counts);
        stds = (bsxfun(@minus, bsxfun(@rdivide, squared_sums, class_counts), means.^2)).^0.5;

        % add background class
        means = [0, 0, 0, 0; means];
        stds = [0, 0, 0, 0; stds];
    end

    % Normalize targets
    for i = 1:num_images
        targets = image_roidb(i).bbox_targets;
        for cls = 1:num_classes
            cls_inds = find(targets(:, 1) == cls);
            if ~isempty(cls_inds)
                image_roidb(i).bbox_targets(cls_inds, 2:end) = ...
                    bsxfun(@minus, image_roidb(i).bbox_targets(cls_inds, 2:end), means(cls+1, :));
                image_roidb(i).bbox_targets(cls_inds, 2:end) = ...
                    bsxfun(@rdivide, image_roidb(i).bbox_targets(cls_inds, 2:end), stds(cls+1, :));
            end
        end
    end
end

function [bbox_targets, is_valid] = compute_targets(conf, rois, overlap)

    overlap = full(overlap);

    [max_overlaps, max_labels] = max(overlap, [], 2);

    % ensure ROIs are floats
    rois = single(rois);

    bbox_targets = zeros(size(rois, 1), 5, 'single');

    % Indices of ground-truth ROIs
    gt_inds = find(max_overlaps == 1);

    if ~isempty(gt_inds)
        % Indices of examples for which we try to make predictions
        ex_inds = find(max_overlaps >= conf.bbox_thresh);

        % Get IoU overlap between each ex ROI and gt ROI
        ex_gt_overlaps = boxoverlap(rois(ex_inds, :), rois(gt_inds, :));

        assert(all(abs(max(ex_gt_overlaps, [], 2) - max_overlaps(ex_inds)) < 10^-4));

        % Find which gt ROI each ex ROI has max overlap with:
        % this will be the ex ROI's gt target
        [~, gt_assignment] = max(ex_gt_overlaps, [], 2);
        gt_rois = rois(gt_inds(gt_assignment), :);
        ex_rois = rois(ex_inds, :);

        [regression_label] = fast_rcnn_bbox_transform(ex_rois, gt_rois);

        bbox_targets(ex_inds, :) = [max_labels(ex_inds), regression_label];
    end

    % Select foreground ROIs as those with >= fg_thresh overlap
    is_fg = max_overlaps >= conf.fg_thresh;
    % Select background ROIs as those within [bg_thresh_lo, bg_thresh_hi)
    is_bg = max_overlaps < conf.bg_thresh_hi & max_overlaps >= conf.bg_thresh_lo;

    % check if there is any fg or bg sample. If no, filter out this image
    is_valid = true;
    if ~any(is_fg | is_bg)
        is_valid = false;
    end
end

author: ShaoqingRen


text:
function [ panoSegment ] = gbPanoSegment( img, sigma, k, minSz )
%GBPANOSEGMENT Summary of this function goes here
%   Detailed explanation goes here
% global config

[height, width, ~] = size(img);
img_smooth = smooth(img, sigma);

%% pixel segmentation, finish
% edges = zeros(width*height*4,3);
% [gridX, gridY] = meshgrid(1:width,1:height);
% num = 0;
%
% x = gridX(:);
% y = gridY(:);
%
% % xyz = uv2xyzN(coords2uv([x y], 2048, 1024));
% % angleNorm = 2*pi/width;
%
% vector = [1 0; 0 1; 1 1; 1 -1];
% inda = sub2ind([height width], y, x);
%
% for vid = 1:size(vector,1)
%     xv = x+vector(vid,1);
%     yv = y+vector(vid,2);
%     valid = xv>=1 & xv<=width & yv>=1 & yv<=height;
%     indav = inda(valid);
%     indbv = sub2ind([height width], yv(valid), xv(valid));
%     diff = (img_smooth(indav)-img_smooth(indbv)).^2 ...
%          + (img_smooth(indav+width*height)-img_smooth(indbv+width*height)).^2 ...
%          + (img_smooth(indav+2*width*height)-img_smooth(indbv+2*width*height)).^2;
% %     edges(num+1:num+sum(valid),:) = [x(valid)-1 + (y(valid)-1)*width ...
% %                                      xv(valid)-1 + (yv(valid)-1)*width ...
% %                                      sqrt(diff)];
% %     angleWeight = sqrt(1-dot(xyz(indav,:), xyz(indbv,:), 2).^2)/angleNorm;
%     angleWeight = ones(length(indav),1);
%     edges(num+1:num+sum(valid),:) = [indav indbv sqrt(diff).*angleWeight];
%     num = num + sum(valid);
% end
%
% % connect left to right
% XY(1).xyLHS = [ones(height,1) (1:height)'];
% XY(1).xyRHS = [width*ones(height,1) (1:height)'];
% XY(2).xyLHS = [ones(height-1,1) (1:height-1)'];
% XY(2).xyRHS = [width*ones(height-1,1) (2:height)'];
% XY(3).xyLHS = [ones(height-1,1) (2:height)'];
% XY(3).xyRHS = [width*ones(height-1,1) (1:height-1)'];
% %
% % combos = combntns(1:width,2);
% % XY(4).xyLHS = [combos(:,1) ones(size(combos,1),1)];
% % XY(4).xyRHS = [combos(:,2) ones(size(combos,1),1)];
% % XY(5).xyLHS = [combos(:,1) height*ones(size(combos,1),1)];
% % XY(5).xyRHS = [combos(:,2) height*ones(size(combos,1),1)];
% %
% for sid = 1:length(XY)
%     xyLHS = XY(sid).xyLHS;
%     xyRHS = XY(sid).xyRHS;
%     indav = sub2ind([height width], xyLHS(:,2), xyLHS(:,1));
%     indbv = sub2ind([height width], xyRHS(:,2), xyRHS(:,1));
%     diff = (img_smooth(indav)-img_smooth(indbv)).^2 ...
%              + (img_smooth(indav+width*height)-img_smooth(indbv+width*height)).^2 ...
%              + (img_smooth(indav+2*width*height)-img_smooth(indbv+2*width*height)).^2;
%     addNum = size(xyLHS,1);
% %     edges(num+1:num+addNum,:) = [xyLHS(:,1)-1+(xyLHS(:,2)-1)*width ...
% %                                  xyRHS(:,1)-1+(xyRHS(:,2)-1)*width ...
% %                                  sqrt(diff)];
% %     angleWeight = sqrt(1-dot(xyz(indav,:), xyz(indbv,:), 2).^2)/angleNorm;
%     angleWeight = ones(length(indav),1);
%     edges(num+1:num+addNum,:) = [indav indbv sqrt(diff).*angleWeight];
%     num = num + addNum;
% end
%
% edges = edges';
% segment = segmentGraphMex_edge(width*height, num, edges, k, minSz);
% L = unique(segment);
% temp = zeros(height, width);
% for i = 1:length(L)
%     temp(segment==L(i)) = i;
% end
% panoSegment = temp;
%
% % segment = segmentGraphMex(width, height, num, edges, k, minSz);
% % L = unique(segment);
% %
% % temp = zeros(size(segment));
% % for i = 1:length(L)
% %     temp(segment==L(i)) = i;
% % end

%% uniformly sample vectors on sphere and segment, test later
% global coor;
% global tri;
% load('./rectangleDetector/segmentation/uniformvector_lvl8.mat');
[coor, tri] = getUniformVector(8);
% load('./region_based_hypothesis/SketchTokens-master/models/forest/modelSmall.mat');
% st = stDetect( img, model );
% E = stToEdges( st, 1 );
[ E ] = getSketchTokenEdgemap( img );

% h = [0 -1 0; -1 4 -1; 0 -1 0];
% E = filter2(h,rgb2gray(img));
[EE, Ix, Iy] = dt2(double(E), 0.1, 0, 0.1, 0 );
% EE = zeros(1024, 2048);
% EE(747:846,1641:1697) = 1;

% [coor,tri] = icosahedron2sphere(level);
xySubs = uv2coords(xyz2uvN(coor), width, height);
xyinds = sub2ind([height width], xySubs(:,2), xySubs(:,1));
offset = width*height;

edges = [tri(:,1) tri(:,2); tri(:,2) tri(:,3); tri(:,3) tri(:,1)];
invert = edges(:,2)<edges(:,1);
edges(invert,:) = edges(invert,[2 1]);

uniEdges = unique(edges, 'rows');
weight = (img_smooth(xyinds(uniEdges(:,1)))-img_smooth(xyinds(uniEdges(:,2)))).^2 ...
       + (img_smooth(xyinds(uniEdges(:,1))+offset)-img_smooth(xyinds(uniEdges(:,2))+offset)).^2 ...
       + (img_smooth(xyinds(uniEdges(:,1))+2*offset)-img_smooth(xyinds(uniEdges(:,2))+2*offset)).^2;
gdweight = (EE(xyinds(uniEdges(:,1)))+EE(xyinds(uniEdges(:,2))))/2;
panoEdge = [uniEdges(:,1)'; uniEdges(:,2)'; sqrt(weight)'+10*double(gdweight)'];

maxID = size(coor,1);
num = size(uniEdges,1);

edgeLabel = segmentGraphMex_edge(maxID, num, panoEdge, k, minSz);

L = unique(edgeLabel);
temp = zeros(size(edgeLabel));
% center = zeros(length(L),2);
[gridX, gridY] = meshgrid(1:width, 1:height);
for i = 1:length(L)
%     fprintf('%d\n',i);
    temp(edgeLabel==L(i)) = i;
%     center(i,:) = [mean(gridX(edgeLabel==L(i))) mean(gridY(edgeLabel==L(i)))];
end

pixelvector = uv2xyzN(coords2uv([gridX(:) gridY(:)], width, height));

k = 1;
% [nnidx, dists] = annsearch( coor', pixelvector', k);
[nnidx, dists] = knnsearch( coor, pixelvector);

panoSegment = reshape(temp(nnidx), height, width);
% [ colors, pc ] = pcsel( center, 6 );
% C = [1 0 0;0 1 0;0 0 1;1 1 0;1 0 1;0 1 1];
% figure; imshow(label2rgb(panoSegment, C(colors,:)));

% panoSegment = zeros(height, width);
% for i = 1:length(L)
%     panoSegment(nnidx
% end

end

author: zouchuhang

repo: LayoutNet