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

% Edit the above text to modify the response to pushbutton_help DQ_Table_GUI

% Last Modified by GUIDE v2.5 26-Jan-2023 16:54:38

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

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

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

% Choose default command line output for DQ_Table_GUI
handles.output = [];

% check input ERPset DQ
ERP = varargin{1};
ALLERP = varargin{2};
current_ERP = varargin{3};
guiwin_num = varargin{4};

%GUI positions (for multiple DQ tables)

if guiwin_num ~= 1
    offset_GUI = 25 * randi([1 10],1); %offset by
    % %set(hObject, 'Units', 'normalized');
    figureposition = get(handles.figure1, 'Position');
    set(hObject, 'Position', ...
        [ figureposition(1)-offset_GUI figureposition(2)-offset_GUI  ...
        figureposition(3) figureposition(4)])
end

%check contents

try
    assert(exist('ERP','var')==1)
    assert(isempty(ERP)==0)
    assert(isfield(ERP,'dataquality'))
    assert(strcmp(ERP.dataquality(1).type,'empty')==0)
catch
    beep
    warning('Data Quality not present in current ERPset?')
    delete(hObject);
    return
end

% Update GUI with ERPSET info
n_erp = numel(ALLERP);

if n_erp == 0
    %for the case of showing aSME preavg
    erp_names = ERP.erpname;
    handles.text_ERPSET_title.String = 'Selected EEG Dataset:';
    handles.text_ERPSET_title.FontSize = 12;
    handles.text11.Visible = 0;
    handles.newerpwin.Visible = 0;
    handles.pushbutton_selERP.Visible = 0;
    handles.active_erpset.Style = 'text';
    if numel(erp_names) > 15
        handles.active_erpset.FontSize = 12.0;
        set(handles.active_erpset,'FontUnits', 'normalized');
    end
    %update handles
    guidata(hObject, handles);

else

    for i = 1:n_erp
        erp_names{i} = ALLERP(i).erpname;
    end

end

%ERPSET_title_str = ['ERPSET - ' erp_names(current_ERP)];
handles.active_erpset.String = erp_names;
handles.active_erpset.Value = current_ERP;

%AMS: Multiple ERPset options
handles.newerpwin.String = erp_names;
handles.newerpwin.Value = current_ERP;

n_dq = numel(ERP.dataquality);
for i=1:n_dq
    type_names{i} = ERP.dataquality(i).type;
end

handles.popupmenu_DQ_type.String = type_names;
%handles.popupmenu_DQ_type.Value = n_dq;
%default to aSME (if available) or end of list

where_aSME = strcmpi('aSME',type_names);
where_aSME_corr = strcmpi('aSME (Corrected)',type_names);

if any(where_aSME)
    indx_type = find(where_aSME);

elseif any(where_aSME_corr)
    indx_type = find(where_aSME_corr);
else
    indx_type = n_dq;
end
handles.popupmenu_DQ_type.Value = indx_type;

try
[n_elec, n_tw, n_bin] = size(ERP.dataquality(indx_type).data); %default to aSME (or last metric)
catch
[n_elec, n_tw, n_bin] = size(ERP.dataquality.data); %if only one measure
end

if n_elec == 0
    % if the dq data is of zero size, then the 3rd dim may be mis-sized.
    % So set to default of ERP.nbin instead.
    n_bin = ERP.nbin;
end

n_bin_names = numel(ERP.bindescr);
if n_bin_names == n_bin

    for i=1:n_bin
        bin_names{i} = ['BIN ' num2str(i) ' - ' ERP.bindescr{i}];
    end
else
    % if not every bin has a bin description, leave then off
    for i=1:n_bin
        bin_names{i} = ['BIN ' num2str(i)];
    end
end

handles.popupmenu_bin.String = bin_names;

% handles.popupmenu_DQ_type.Value saves the indx of the selected DQ Measure
selected_DQ_type = handles.popupmenu_DQ_type.Value;
selected_bin = handles.popupmenu_bin.Value;

table_data = ERP.dataquality(selected_DQ_type).data(:,:,selected_bin);

pSEM = 0;
if isempty(table_data) %check if potentially pointwise SEM
    if strcmp(ERP.dataquality(selected_DQ_type).type,'Point-wise SEM') && isempty(ERP.binerror) == 0
        table_data = ERP.binerror(:,:,selected_bin);
        pSEM = 1;

    elseif strcmp(ERP.dataquality(selected_DQ_type).type,'Point-wise SEM (Corrected)') && isempty(ERP.binerror) == 0
        table_data = ERP.binerror(:,:,selected_bin);
        pSEM = 1;

    end

end

handles.dq_table.Data = table_data;
handles.orig_data = table_data;

% electrode labels from ERPset, iff present and number matches
if isfield(ERP.chanlocs,'labels') && numel(ERP.chanlocs) == n_elec
    for i=1:n_elec
        elec_labels{i} = ERP.chanlocs(i).labels;
    end
else

    for i=1:ERP.nchan
        elec_labels{i} = i;
    end
end
handles.dq_table.RowName = elec_labels;
handles.orig_RowName = elec_labels;

% Time-window labels
if isfield(ERP.dataquality(selected_DQ_type),'time_window_labels') && isempty(ERP.dataquality(selected_DQ_type).time_window_labels) == 0  && handles.checkbox_text_labels.Value == 1
    tw_labels = ERP.dataquality(selected_DQ_type).time_window_labels;
elseif isempty(ERP.dataquality(selected_DQ_type).times)
    tw_labels = [];
else
    for i=1:n_tw
        tw_labels{i} = [num2str(ERP.dataquality(selected_DQ_type).times(i,1)) ' : ' num2str(ERP.dataquality(selected_DQ_type).times(i,2))];
    end
end
if n_tw == 0 %no ERP.dataquality.data?

    if pSEM == 1
        tw_labels = [];
    else
        tw_labels = 'No data here. Perhaps it was cleared?';
    end

end
handles.dq_table.ColumnName = tw_labels;

% Set font size of DQ Table
desired_fontsize = erpworkingmemory('fontsizeGUI');
handles.dq_table.FontSize = desired_fontsize;
handles.heatmap_on = 0;

%
% Prepare List of current Channels
%
nchan  = ERP.nchan; % Total number of channels
if ~isfield(ERP.chanlocs,'labels')
        for e=1:nchan
                ERP.chanlocs(e).labels = ['Ch' num2str(e)];
        end
end
listch = {''};
for ch =1:nchan
        listch{ch} = [num2str(ch) ' = ' ERP.chanlocs(ch).labels ];
end
handles.listch     = listch;

chanArray = 1:ERP.nchan; %default
handles.indxlistch = chanArray; %default
handles.orig_indxlistch = chanArray;
set(handles.chwindow, 'String', vect2colon(chanArray, 'Delimiter', 'off'));

% Set outliers text window off
set(handles.stdwindow,'Enable','Off');
set(handles.chwindow,'Enable','Off');
set(handles.chbutton,'Enable','Off');
set(handles.sdcorrection,'Enable','Off');
set(handles.sdcorrection,'Visible','Off');

%set outliers value in handles off initially
handles.outliers_on = 0;

% Update handles structure
handles.ERP = ERP;

if n_erp == 0
    %case of aSME preavg
    handles.ALLERP = ERP;
else
    handles.ALLERP = ALLERP;
end

guidata(hObject, handles);

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

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

% Get default command line output from handles structure
varargout{1} = [];

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

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

% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu_DQ_type contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu_DQ_type
% handles.popupmenu_DQ_type.Value saves the indx of the selected DQ Measure
selected_DQ_type = handles.popupmenu_DQ_type.Value;
selected_bin = handles.popupmenu_bin.Value;

% Check data exists, plot
if isempty(handles.ERP.dataquality(selected_DQ_type).data)

    % if pointwise SEM, use ERP.binerror
    if strcmp(handles.ERP.dataquality(selected_DQ_type).type,'Point-wise SEM') && isempty(handles.ERP.binerror) == 0
        table_data = handles.ERP.binerror(:,:,selected_bin);

    elseif strcmp(handles.ERP.dataquality(selected_DQ_type).type,'Point-wise SEM (Corrected)') && isempty(handles.ERP.binerror) == 0
        table_data = handles.ERP.binerror(:,:,selected_bin);
    else

        % Data empty here.
        table_data = nan(1);
        disp('DQ data not found for this DQ type and bin. Perhaps it has been cleared?');
    end

else
    % data is present and correct

    if strcmp(handles.ERP.dataquality(selected_DQ_type).type,'SD Across Trials') && isempty(handles.ERP.binerror) == 0
        table_data = handles.ERP.dataquality(selected_DQ_type).data.SD_bias(:,:,selected_bin); %SD (Divided by N-1)
        set(handles.sdcorrection,'Enable','On');
        set(handles.sdcorrection,'Visible','On');
    else
        table_data = handles.ERP.dataquality(selected_DQ_type).data(:,:,selected_bin);
        set(handles.sdcorrection,'Enable','Off');
        set(handles.sdcorrection,'Visible','Off');
    end

end

handles.dq_table.Data = table_data;
handles.orig_data = table_data;

% Time-window labels
clear tw_labels
ERP = handles.ERP;
n_dq = numel(ERP.dataquality);
if strcmp(handles.ERP.dataquality(selected_DQ_type).type,'SD Across Trials') && isempty(handles.ERP.binerror) == 0
    [n_elec, n_tw, n_bin] = size(ERP.dataquality(selected_DQ_type).data.SD_bias);
else
    [n_elec, n_tw, n_bin] = size(ERP.dataquality(selected_DQ_type).data);
end

if isfield(ERP.dataquality(selected_DQ_type),'time_window_labels') && isempty(ERP.dataquality(selected_DQ_type).time_window_labels) == 0 && handles.checkbox_text_labels.Value == 1
    tw_labels = ERP.dataquality(selected_DQ_type).time_window_labels;
elseif isempty(ERP.dataquality(selected_DQ_type).times)
    tw_labels = [];
elseif n_tw == 0  % some data problem, no tw here
    tw_labels{1} = 'No data here - measure cleared?';
else
    for i=1:n_tw
        tw_labels{i} = [num2str(ERP.dataquality(selected_DQ_type).times(i,1)) ' : ' num2str(ERP.dataquality(selected_DQ_type).times(i,2))];
    end
end
handles.dq_table.ColumnName = tw_labels;

if handles.heatmap_on
    redraw_heatmap(hObject, eventdata, handles);
end

if handles.outliers_on
   clear_heatmap(hObject, eventdata, handles);
   handles.heatmap_on = 0;

   redraw_outliers(hObject, eventdata, handles);
   set(handles.checkbox_outliers,'Value',1);
end

% Update handles structure
guidata(hObject, handles);

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

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

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

% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu_bin contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu_bin
% handles.popupmenu_DQ_type.Value saves the indx of the selected DQ Measure
selected_DQ_type = handles.popupmenu_DQ_type.Value;
selected_bin = handles.popupmenu_bin.Value;

% Check data exists, plot
if isempty(handles.ERP.dataquality(selected_DQ_type).data)

    % if pointwise SEM, use ERP.binerror
    if strcmp(handles.ERP.dataquality(selected_DQ_type).type,'Point-wise SEM') && isempty(handles.ERP.binerror) == 0
        table_data = handles.ERP.binerror(:,:,selected_bin);
    else

        % Data empty here.
        table_data = nan(1);
        disp('DQ data not found for this DQ type and bin. Perhaps it has been cleared?');
    end

else
    % data is present and correct
    table_data = handles.ERP.dataquality(selected_DQ_type).data(:,:,selected_bin);
end

handles.dq_table.Data = table_data;
handles.orig_data = table_data;

%if user switches bin and both heatmap & outliers on,
%use shortcircuit to clear outliers and redraw heatmap only
if handles.heatmap_on %|| handles.outliers_on
    clear_outliers(hObject, eventdata, handles);
    redraw_heatmap(hObject, eventdata, handles);
    %turn outliers off
    handles.outliers = 0;
    set(handles.checkbox_outliers,'Value',0);

end

%if only outliers on, then clear heatmap and redraw outliers
if handles.outliers_on
   clear_heatmap(hObject, eventdata, handles);
   handles.heatmap_on = 0;

   redraw_outliers(hObject, eventdata, handles);
   set(handles.checkbox_outliers,'Value',1);
end

% Update handles structure
guidata(hObject, handles);

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

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

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

% Hint: get(hObject,'Value') returns toggle state of checkbox_text_labels
selected_DQ_type = handles.popupmenu_DQ_type.Value;
selected_bin = handles.popupmenu_bin.Value;

% if strcmp(handles.ERP.dataquality(selected_DQ_type).type,'SD Across Trials') && isempty(handles.ERP.binerror) == 0
%     table_data = handles.ERP.dataquality(selected_DQ_type).data.SD_bias(:,:,selected_bin);
% else
%     table_data = handles.ERP.dataquality(selected_DQ_type).data(:,:,selected_bin);
%
% end
% handles.dq_table.Data = table_data;
% handles.orig_data = table_data;

% Time-window labels
clear tw_labels
ERP = handles.ERP;
n_dq = numel(ERP.dataquality);
%[n_elec, n_tw, n_bin] = size(ERP.dataquality(selected_DQ_type).data);
if strcmp(handles.ERP.dataquality(selected_DQ_type).type,'SD Across Trials') && isempty(handles.ERP.binerror) == 0
    [n_elec, n_tw, n_bin] = size(ERP.dataquality(selected_DQ_type).data.SD_bias);
else
    [n_elec, n_tw, n_bin] = size(ERP.dataquality(selected_DQ_type).data);
end
if isfield(ERP.dataquality(selected_DQ_type),'time_window_labels') && isempty(ERP.dataquality(selected_DQ_type).time_window_labels) == 0 && handles.checkbox_text_labels.Value == 1
    tw_labels = ERP.dataquality(selected_DQ_type).time_window_labels;
elseif isempty(ERP.dataquality(selected_DQ_type).times)
    tw_labels = [];
elseif n_tw == 0  % some data problem, no tw here
    tw_labels{1} = 'No data here - measure cleared?';
else
    for i=1:n_tw
        tw_labels{i} = [num2str(ERP.dataquality(selected_DQ_type).times(i,1)) ' : ' num2str(ERP.dataquality(selected_DQ_type).times(i,2))];
    end
end
handles.dq_table.ColumnName = tw_labels;

if handles.heatmap_on
    redraw_heatmap(hObject, eventdata, handles);
end

 % Update handles structure
 guidata(hObject, handles);

% --- Executes on button press in pushbutton_xls.
function pushbutton_xls_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_xls (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
empty_filename = [];
selected_DQ_type = handles.popupmenu_DQ_type.Value;
%sd_correction = handles.sdcorrection.Value;
save_data_quality(handles.ERP,empty_filename,'xlsx',selected_DQ_type)

% --- Executes on button press in pushbutton_mat.
function pushbutton_mat_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_mat (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
empty_filename = [];
selected_DQ_type = handles.popupmenu_DQ_type.Value;
save_data_quality(handles.ERP,empty_filename,'mat',selected_DQ_type)

% --- Executes on button press in pushbutton_help.
function pushbutton_help_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_help (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
web https://github.com/lucklab/erplab/wiki/ERPLAB-Data-Quality-Metrics -browser

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

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

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

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

% Hint: get(hObject,'Value') returns toggle state of checkbox_heatmap
heatmap_on = get(hObject,'Value');

if heatmap_on == 1

    %if heatmap_on, then outliers not possible so clear outliers
    clear_outliers(hObject, eventdata, handles);
    handles.outliers_on = 0;
    set(handles.checkbox_outliers, 'Value', 0);

    redraw_heatmap(hObject, eventdata, handles);
    handles.heatmap_on = 1;
else
    clear_heatmap(hObject, eventdata, handles);
    handles.heatmap_on = 0;

    %if outliers is on, then keep outliers active
    if handles.outliers_on == 1
        redraw_outliers(hOBject, eventdata,handles);
    end

end
% Update handles structure
guidata(hObject, handles);

function redraw_heatmap(hObject, eventdata, handles)

color_map = viridis;
data = handles.dq_table.Data;
try
    data_min = min(data(:));
catch
    disp('Cannot show text-labels while color heatmap is activated')
    return
end
range_val = max(data(:)) - data_min;
range_colormap = size(color_map,1); % as in, 256 shades?
val_increase_per_shade = range_val / range_colormap;

% Use this @ anonymous function to make HTML tag in box in loop below
colergen = @(color,text) ['<html><table border=0 width=400 bgcolor=',color,'><TR><TD>',text,'</TD></TR> </table>'];

for cell = 1:numel(data)
    data_here = data(cell);
    shades_up_here = round((data_here - data_min) / val_increase_per_shade);
    if shades_up_here < 1
        shades_up_here = 1;
    end
    if shades_up_here > range_colormap
        shades_up_here = range_colormap;
    end
    RGB_here = color_map(shades_up_here,:);
    hex_color_here = ['#' dec2hex(round(255*RGB_here(1)),2) dec2hex(round(255*RGB_here(2)),2) dec2hex(round(255*RGB_here(3)),2)];
    data2{cell} = colergen(hex_color_here,num2str(data_here));
end

data2 = reshape(data2,size(data));

handles.dq_table.Data = data2;

function clear_heatmap(hObject, eventdata, handles)
handles.dq_table.Data = handles.orig_data;

% --- Executes when entered data in editable cell(s) in dq_table.
function dq_table_CellEditCallback(hObject, eventdata, handles)
% hObject    handle to dq_table (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.CONTROL.TABLE)
%	Indices: row and column indices of the cell(s) edited
%	PreviousData: previous data for the cell(s) edited
%	EditData: string(s) entered by the user
%	NewData: EditData or its converted form set on the Data property. Empty if Data was not changed
%	Error: error string when failed to convert EditData to appropriate value for Data
% handles    structure with handles and user data (see GUIDATA)

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

% Hints: contents = cellstr(get(hObject,'String')) returns newerpwin contents as cell array
%        contents{get(hObject,'Value')} returns selected item from newerpwin

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

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

% --- Executes on button press in pushbutton_selERP.
function pushbutton_selERP_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_selERP (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
selected_erpset = handles.newerpwin.Value;
sel_ERP = handles.ALLERP(selected_erpset);
DQ_Table_GUI(sel_ERP, handles.ALLERP, selected_erpset,2);

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

% Hint: get(hObject,'Value') returns toggle state of checkbox_outliers

% Hint: get(hObject,'Value') returns toggle state of checkbox_heatmap
outliers_on = get(hObject,'Value'); % this is not the same as handles.outliers_on

if outliers_on == 1
    set(handles.stdwindow,'Enable','On');
    set(handles.chwindow,'Enable','On');
    set(handles.chbutton,'Enable','On');
else
    set(handles.stdwindow,'Enable','Off');
    set(handles.chwindow,'Enable','Off');
    set(handles.chbutton,'Enable','Off');
end

%updated the handles object for outliers to be in sync with hObject
if outliers_on == 1
    handles.outliers_on = 1;
else
    handles.outliers_on = 0;
end

if handles.outliers_on == 1

    %if heatmap is currently on, stop it and turn its states off
    %before applying outliers fxn

    if handles.heatmap_on == 1
        clear_heatmap(hObject, eventdata, handles);
        handles.heatmap_on = 0;
        set(handles.checkbox_heatmap, 'Value', 0);

        %handles.outliers_on = 1;
        redraw_outliers(hObject, eventdata, handles);

    else
        %handles.outliers_on = 1;
        clear_outliers(hObject, eventdata, handles);
        %guidata(hObject, handles);
        redraw_outliers(hObject, eventdata, handles);

    end

else %but if outliers is off, and heatmap is wanted, you can do it

%     if handles.heatmap_on ==1
%
%         clear_outliers(hObject, eventdata, handles);
%         redraw_heatmap(hObject, eventdata, handles);
%         handles.outliers_on = 0;
%
%     else
%         clear_outliers(hObject, eventdata, handles);
%         handles.outliers_on = 0;
%
%     end

    %also, redraw initial data/channels without outliers
    clear_outliers(hObject, eventdata, handles)

end

% %if heatmap is currently on, stop it and turn its states off
% %before applying outliers fxn
% if handles.heatmap_on == 1
%
%     clear_heatmap(hObject, eventdata, handles);
%     handles.heatmap_on = 0;
%     set(handles.checkbox_heatmap, 'Value', 0);
%
%     if outliers_on == 1
%         handles.outliers_on = 1;
%         redraw_outliers(hObject, eventdata, handles);
%     else
%
%         set(handles.stdwindow,'Enable','Off');
%         set(handles.chwindow,'Enable','Off');
%         set(handles.chbutton,'Enable','Off');
%         clear_outliers(hObject, eventdata, handles);
%         handles.outliers_on = 0;
%
%     end
%
% else
%
%      if outliers_on == 1
%         handles.outliers_on = 1;
%         redraw_outliers(hObject, eventdata, handles);
%      else
%         set(handles.stdwindow,'Enable','Off');
%         set(handles.chwindow,'Enable','Off');
%         set(handles.chbutton,'Enable','Off');
%         clear_outliers(hObject, eventdata, handles);
%         handles.outliers_on = 0;
%
%      end
%
% end
% Update handles structure
guidata(hObject, handles);

function redraw_outliers(hObject, eventdata, handles)
handles = guidata(hObject); %get the most updated handles struct
data = handles.dq_table.Data;

if iscell(data)
    %in case of new std from stdwindow
    data = handles.orig_data;
end

%either default(all) or from chbutton
chans_to_use = handles.indxlistch;
data = data(chans_to_use,:);
n_elec = length(chans_to_use);
ERPtouse = handles.ERP.chanlocs;

for i=1:n_elec
    curr_chanInd = chans_to_use(i);
    elec_labels{i} = ERPtouse(curr_chanInd).labels;
end

%save new labels to handles
handles.dq_table.RowName = elec_labels;

%compute outliers as X standard deviations from mean across channels

col_means = mean(data,1);
col_std = std(data,0,1);
chosen_std = str2double((handles.stdwindow.String));

% Use this @ anonymous function to make HTML tag in box in loop below
colergen = @(color,text) ['<html><table border=0 width=400 bgcolor=',color,'><TR><TD>',text,'</TD></TR> </table>'];

for cell = 1:numel(data)

    data_here = data(cell);
    [d_row,d_col] = ind2sub(size(data),cell);

    neg_thresh = col_means(d_col) - (chosen_std*col_std(d_col));
    pos_thresh = col_means(d_col) + (chosen_std*col_std(d_col));

    if data_here < neg_thresh || data_here > pos_thresh
        hex_color_here = ['#' dec2hex(255,2) dec2hex(0,2) dec2hex(0,2)];
        data2{cell} = colergen(hex_color_here,num2str(data_here));
    else
        data2{cell} = num2str(data_here);

    end

end

data2= reshape(data2,size(data));
handles.dq_table.Data = data2;

function clear_outliers(hObject, eventdata, handles)
%handles = guidata(hObject);
handles.dq_table.Data = handles.orig_data;
handles.dq_table.RowName = handles.orig_RowName;
handles.indxlistch = handles.orig_indxlistch;

% Re-Prepare List of current Channels
chanArray = handles.orig_indxlistch; %default
%handles.indxlistch = chanArray; %default
set(handles.chwindow, 'String', vect2colon(chanArray, 'Delimiter', 'off'));
 % Update handles structure
guidata(hObject, handles);

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

newstd = get(hObject,'String');
handles.stdwindow.String= newstd;
redraw_outliers(hObject, eventdata, handles);

guidata(hObject, handles);

% Hints: get(hObject,'String') returns contents of stdwindow as text
%        str2double(get(hObject,'String')) returns contents of stdwindow as a double

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

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

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

% Hints: get(hObject,'String') returns contents of chwindow as text
%        str2double(get(hObject,'String')) returns contents of chwindow as a double

ch = get(hObject,'string');
nch = length(handles.listch);

ch_val = eval(ch);

if ~isempty(ch)

%     testch = regexpi(ch,':') ;
%
%     if isempty(testch)
%
%         try
%             ch_val = eval(ch);
%             handles.listch(ch_val); %attempt to index
%
%         catch
%             msgboxText =  'Invalid channel input: please enter channels as integers, with a space between each number, or ch1:chN array syntax';
%             title = 'ERPLAB: channel GUI input';
%             errorfound(msgboxText, title);
%
%         end
%     else
%         ch_val = eval(ch);
%
%     end
%

    tf = checkchannels(ch_val, nch, 1);

    if tf
        return
    end

%     if length(ch_val) > length(handles.orig_indxlistch)
%         msgboxText =  'Exceeded Number of Available Channels';
%         title = 'ERPLAB: channel GUI input';
%         errorfound(msgboxText, title);
%         return
%     end
%
    set(handles.chwindow, 'String', vect2colon(ch_val, 'Delimiter', 'off'));
    handles.indxlistch = ch_val;
    % Update handles structure
    guidata(hObject, handles);
    redraw_outliers(hObject, eventdata, handles);
    guidata(hObject, handles);

else
    msgboxText =  'Not valid channel';
    title = 'ERPLAB: channel GUI input';
    errorfound(msgboxText, title);
    return
end

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

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

% --- Executes on button press in chbutton.
function chbutton_Callback(hObject, eventdata, handles)
% hObject    handle to chbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
listch = handles.listch; %true channels list as labels
indxlistch = handles.indxlistch; %true array of channel indexs
indxlistch = indxlistch(indxlistch<=length(listch)); % true array of channel index if less than expected labels
titlename = 'Select Channel(s)';
if get(hObject, 'Value')
        if ~isempty(listch)
                ch = browsechanbinGUI(listch, indxlistch, titlename);
                if ~isempty(ch)
                        set(handles.chwindow, 'String', vect2colon(ch, 'Delimiter', 'off'));
                        handles.indxlistch = ch;
                        % Update handles structure
                        guidata(hObject, handles);
                        redraw_outliers(hObject, eventdata, handles);
                        guidata(hObject, handles);

                else
                        disp('User selected Cancel')
                        return
                end
        else
                msgboxText =  'No channel information was found';
                title = 'ERPLAB: ploterp GUI input';
                errorfound(msgboxText, title);
                return
        end
end

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

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

newactive = get(hObject,'Value');
ERPs = handles.ALLERP;
ERP = ERPs(newactive);

n_dq = numel(ERP.dataquality);
for i=1:n_dq
    type_names{i} = ERP.dataquality(i).type;
end

handles.popupmenu_DQ_type.String = type_names;
handles.popupmenu_DQ_type.Value = n_dq;

[n_elec, n_tw, n_bin] = size(ERP.dataquality(n_dq).data);

if n_elec == 0
    % if the dq data is of zero size, then the 3rd dim may be mis-sized.
    % So set to default of ERP.nbin instead.
    n_bin = ERP.nbin;
end

n_bin_names = numel(ERP.bindescr);
if n_bin_names == n_bin

    for i=1:n_bin
        bin_names{i} = ['BIN ' num2str(i) ' - ' ERP.bindescr{i}];
    end
else
    % if not every bin has a bin description, leave then off
    for i=1:n_bin
        bin_names{i} = ['BIN ' num2str(i)];
    end
end

handles.popupmenu_bin.String = bin_names;

% handles.popupmenu_DQ_type.Value saves the indx of the selected DQ Measure
selected_DQ_type = handles.popupmenu_DQ_type.Value;
selected_bin = handles.popupmenu_bin.Value;

table_data = ERP.dataquality(selected_DQ_type).data(:,:,selected_bin);
handles.dq_table.Data = table_data;
handles.orig_data = table_data;

% electrode labels from ERPset, iff present and number matches
if isfield(ERP.chanlocs,'labels') && numel(ERP.chanlocs) == n_elec
    for i=1:n_elec
        elec_labels{i} = ERP.chanlocs(i).labels;
    end
else

    for i=1:ERP.nchan
        elec_labels{i} = i;
    end
end
handles.dq_table.RowName = elec_labels;

% Time-window labels
if isfield(ERP.dataquality(selected_DQ_type),'time_window_labels') && isempty(ERP.dataquality(selected_DQ_type).time_window_labels) == 0  && handles.checkbox_text_labels.Value == 1
    tw_labels = ERP.dataquality(selected_DQ_type).time_window_labels;
elseif isempty(ERP.dataquality(selected_DQ_type).times)
    tw_labels = [];
else
    for i=1:n_tw
        tw_labels{i} = [num2str(ERP.dataquality(selected_DQ_type).times(i,1)) ' : ' num2str(ERP.dataquality(selected_DQ_type).times(i,2))];
    end
end
if n_tw == 0
    tw_labels = 'No data here. Perhaps it was cleared?';
end
handles.dq_table.ColumnName = tw_labels;

% Set font size of DQ Table
desired_fontsize = erpworkingmemory('fontsizeGUI');
handles.dq_table.FontSize = desired_fontsize;
handles.heatmap_on = 0;

%
% Prepare List of current Channels
%
nchan  = ERP.nchan; % Total number of channels
if ~isfield(ERP.chanlocs,'labels')
        for e=1:nchan
                ERP.chanlocs(e).labels = ['Ch' num2str(e)];
        end
end
listch = {''};
for ch =1:nchan
        listch{ch} = [num2str(ch) ' = ' ERP.chanlocs(ch).labels ];
end
handles.listch     = listch;

chanArray = 1:ERP.nchan; %default
handles.indxlistch = chanArray; %default
set(handles.chwindow, 'String', vect2colon(chanArray, 'Delimiter', 'off'));

% Set outliers text window off
set(handles.stdwindow,'Enable','Off');
set(handles.chwindow,'Enable','Off');
set(handles.chbutton,'Enable','Off');

% Update handles structure
handles.ERP = ERP;
handles.ALLERP = ERPs;

%clear all perviously set options
set(handles.checkbox_outliers,'Value',0);
set(handles.checkbox_heatmap,'Value',0);
set(handles.checkbox_text_labels,'Value',0);
set(handles.stdwindow,'Enable','Off');
set(handles.chwindow,'Enable','Off');
set(handles.chbutton,'Enable','Off');

%check to see if outliers is on prior to switching active erpset?
%checkbox_outliers_Callback(hObject, eventdata, handles)

% outliers_on = get(handles.checkbox_outliers,'Value');
%
% if outliers_on == 1
%     handles.outliers_on = 1;
%     set(handles.stdwindow,'Enable','On');
%     set(handles.chwindow,'Enable','On');
%     set(handles.chbutton,'Enable','On');
%     redraw_outliers(hObject, eventdata, handles);
%
% else
%     handles.heatmap_on = 0;
%     set(handles.stdwindow,'Enable','Off');
%     set(handles.chwindow,'Enable','Off');
%     set(handles.chbutton,'Enable','Off');
%     clear_outliers(hObject, eventdata, handles);
%
%
% end

guidata(hObject, handles);

% Hints: contents = cellstr(get(hObject,'String')) returns active_erpset contents as cell array
%        contents{get(hObject,'Value')} returns selected item from active_erpset

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

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

%--------------------------------------------------------------------------
function tf = checkchannels(chx, nchan, showmsg)

if nargin<3
        showmsg = 1;
end
tf = 0; % no problem by default

if ~mod(chx, 1) == 0
    if showmsg
        msgboxText =  'Invalid channel indexing.';
        title = 'ERPLAB: basicfilterGUI() error:';
        errorfound(msgboxText, title);
    end
    tf = 1; %

end

if isempty(chx)
        if showmsg
                msgboxText =  'Invalid channel indexing.';
                title = 'ERPLAB: basicfilterGUI() error:';
                errorfound(msgboxText, title);
        end
        tf = 1; %
        return
end
if ~isempty(find(chx>nchan))
        if showmsg
                msgboxText =  ['You only have %g channels,\n'...
                        'so you cannot specify indices greater than this.'];
                title = 'ERPLAB: basicfilterGUI() error:';
                errorfound(sprintf(msgboxText, nchan), title);
        end
        tf = 1; %
        return
end
if ~isempty(find(chx<1))
        if showmsg
                msgboxText =  'You cannot use zero or a negative number as a channel indexing';
                title = 'ERPLAB: basicfilterGUI() error:';
                errorfound(msgboxText, title);
        end
        tf = 1; %
        return
end
if length(chx)>length(unique_bc2(chx))
        if showmsg
                msgboxText =  ['Repeated channels are not allowed.\n'...
                        'Therefore, ERPLAB will get rid of them.'];
                title = 'ERPLAB: basicfilterGUI() error:';
                errorfound(sprintf(msgboxText), title, [1 1 0], [0 0 0], 0)
        end
        tf = 0; %
        return
end

% --- Executes on button press in sdcorrection.
function sdcorrection_Callback(hObject, eventdata, handles)
% hObject    handle to sdcorrection (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
selected_DQ_type = handles.popupmenu_DQ_type.Value;
selected_bin = handles.popupmenu_bin.Value;

if handles.sdcorrection.Value == 1
    table_data = handles.ERP.dataquality(selected_DQ_type).data.SD_unbias(:,:,selected_bin); %SD (Gurland & Tripathi, 1971)
else
    table_data = handles.ERP.dataquality(selected_DQ_type).data.SD_bias(:,:,selected_bin); %SD (N-1)
end
handles.dq_table.Data = table_data;
handles.orig_data = table_data;

%update handles struct
guidata(hObject, handles);

% Hint: get(hObject,'Value') returns toggle state of sdcorrection

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function [D]=patchEdgeLengths(F,V)

% function [D]=patchEdgeLengths(F,V)
% -----------------------------------------------------------------------
% Computes the edge lengths (D) for the patch data specified by the faces
% (F) and vertices (V) arrays. If size(F,2)>2 it is assumed that F indeed
% represents faces. If however size(F,2)==2 it is instead assumed that F is
% an array representing edges. As such it skips the computation of the
% edges array. The edges array used is non-unique by default. See the
% |patchEdges| function for more details if the lengths of a unique set of
% edges is desired.
%
%
% See also: |patchEdges|
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2014/03/17
% 2019/04/23 Expanded for cell array input
%------------------------------------------------------------------------

%%

if isa(F,'cell')
    D=F; %Initialize as F
    for q=1:1:numel(F) %Loop over cell entries and call function for each cell entry
        D{q}=patchEdgeLengths(F{q},V); %parse cell entry
    end
else
    %Derive edge array
    if size(F,2)>2 %The input is assumed to represent faces hence an edge array is derived
        E=patchEdges(F,0);
    else %It is assumed that the input array represents an edges array
        E=F;
    end

    %Derive edge vertex arrays
    V_E1=V(E(:,1),:);
    V_E2=V(E(:,2),:);

    %Derive difference vectors
    VD=(V_E1-V_E2);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [Xgp Ygp]=projectToGroundPlane(Xi, Yi, sceneInfo)
%
% (C) Anton Andriyenko, 2012
%
% The code may be used free of charge for non-commercial and
% educational purposes, the only requirement is that this text is
% preserved within the derivative work. For any other purpose you
% must contact the authors for permission. This code may not be
% redistributed without written permission from the authors.

[F, N]=size(Xi);
Xgp=zeros(size(Xi));
Ygp=zeros(size(Xi));

for t=1:F
    extar=find(Xi(t,:));
    for id=extar
        [Xgp(t,id), Ygp(t,id), zw]=imageToWorld(Xi(t,id), Yi(t,id), sceneInfo.camPar);
    end
end

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function detector = acfTrain2( varargin )
% Train aggregate channel features object detector.
%
% Train aggregate channel features (ACF) object detector as described in:
%  P. Doll?, R. Appel, S. Belongie and P. Perona
%   "Fast Feature Pyramids for Object Detection", PAMI 2014.
% The ACF detector is fast (30 fps on a single core) and achieves top
% accuracy on rigid object detection. Please see acfReadme.m for details.
%
% Takes a set of parameters opts (described in detail below) and trains a
% detector from start to finish including performing multiple rounds of
% bootstrapping if need be. The return is a struct 'detector' for use with
% acfDetect.m which fully defines a sliding window detector. Training is
% fast (on the INRIA pedestrian dataset training takes ~10 minutes on a
% single core or ~3m using four cores). Taking advantage of parallel
% training requires launching matlabpool (see help for matlabpool). The
% trained detector may be altered in certain ways via acfModify(). Calling
% opts=acfTrain() returns all default options.
%
% (1) Specifying features and model: The channel features are defined by
% 'pPyramid'. See chnsCompute.m and chnsPyramid.m for more details. The
% channels may be convolved by a set 'filters' to remove local correlations
% (see our NIPS14 paper on LDCF), improving accuracy but slowing detection.
% If 'filters'=[wFilter,nFilter] these are automatically computed. The
% model dimensions ('modelDs') define the window height and width. The
% padded dimensions ('modelDsPad') define the extended region around object
% candidates that are used for classification. For example, for 100 pixel
% tall pedestrians, typically a 128 pixel tall region is used to make a
% decision. 'pNms' controls non-maximal suppression (see bbNms.m), 'stride'
% controls the window stride, and 'cascThr' and 'cascCal' are the threshold
% and calibration used for the constant soft cascades. Typically, set
% 'cascThr' to -1 and adjust 'cascCal' until the desired recall is reached
% (setting 'cascCal' shifts the final scores output by the detector by the
% given amount). Training alternates between sampling (bootstrapping) and
% training an AdaBoost classifier (clf). 'nWeak' determines the number of
% training stages and number of trees after each stage, e.g. nWeak=[32 128
% 512 2048] defines four stages with the final clf having 2048 trees.
% 'pBoost' specifies parameters for AdaBoost, and 'pBoost.pTree' are the
% decision tree parameters, see adaBoostTrain.m for details. Finally,
% 'seed' is the random seed used and makes results reproducible and 'name'
% defines the location for storing the detector and log file.
%
% (2) Specifying training data location and amount: The training data can
% take on a number of different forms. The positives can be specified using
% either a dir of pre-cropped windows ('posWinDir') or dirs of full images
% ('posImgDir') and ground truth labels ('posGtDir'). The negatives can by
% specified using a dir of pre-cropped windows ('negWinDir'), a dir of full
% images without any positives and from which negatives can be sampled
% ('negImgDir'), and finally if neither 'negWinDir' or 'negImgDir' are
% given negatives are sampled from the images in 'posImgDir' (avoiding the
% positives). For the pre-cropped windows all images must have size at
% least modelDsPad and have the object (of size exactly modelDs) centered.
% 'imreadf' can be used to specify a custom function for loading an image,
% and 'imreadp' are custom additional parameters to imreadf. When sampling
% from full images, 'pLoad' determines how the ground truth is loaded and
% converted to a set of positive bbs (see bbGt>bbLoad). 'nPos' controls the
% total number of positives to sample for training (if nPos=inf the number
% of positives is limited by the training set). 'nNeg' controls the total
% number of negatives to sample and 'nPerNeg' limits the number of
% negatives to sample per image. 'nAccNeg' controls the maximum number of
% negatives that can accumulate over multiple stages of bootstrapping.
% Define 'pJitter' to jitter the positives (see jitterImage.m) and thus
% artificially increase the number of positive training windows. Finally if
% 'winsSave' is true cropped windows are saved to disk as a mat file.
%
% USAGE
%  detector = acfTrain( opts )
%  opts = acfTrain()
%
% INPUTS
%  opts       - parameters (struct or name/value pairs)
%   (1) features and model:
%   .pPyramid   - [{}] params for creating pyramid (see chnsPyramid)
%   .filters    - [] [wxwxnChnsxnFilter] filters or [wFilter,nFilter]
%   .modelDs    - [] model height+width without padding (eg [100 41])
%   .modelDsPad - [] model height+width with padding (eg [128 64])
%   .pNms       - [..] params for non-maximal suppression (see bbNms.m)
%   .stride     - [4] spatial stride between detection windows
%   .cascThr    - [-1] constant cascade threshold (affects speed/accuracy)
%   .cascCal    - [.005] cascade calibration (affects speed/accuracy)
%   .nWeak      - [128] vector defining number weak clfs per stage
%   .pBoost     - [..] parameters for boosting (see adaBoostTrain.m)
%   .seed       - [0] seed for random stream (for reproducibility)
%   .name       - [''] name to prepend to clf and log filenames
%   (2) training data location and amount:
%   .posGtDir   - [''] dir containing ground truth
%   .posImgDir  - [''] dir containing full positive images
%   .negImgDir  - [''] dir containing full negative images
%   .posWinDir  - [''] dir containing cropped positive windows
%   .negWinDir  - [''] dir containing cropped negative windows
%   .imreadf    - [@imread] optional custom function for reading images
%   .imreadp    - [{}] optional custom parameters for imreadf
%   .pLoad      - [..] params for bbGt>bbLoad (see bbGt)
%   .nPos       - [inf] max number of pos windows to sample
%   .nNeg       - [5000] max number of neg windows to sample
%   .nPerNeg    - [25]  max number of neg windows to sample per image
%   .nAccNeg    - [10000] max number of neg windows to accumulate
%   .pJitter    - [{}] params for jittering pos windows (see jitterImage)
%   .winsSave   - [0] if true save cropped windows at each stage to disk
%
% OUTPUTS
%  detector   - trained object detector (modify only via acfModify)
%   .opts       - input parameters used for model training
%   .clf        - learned boosted tree classifier (see adaBoostTrain)
%   .info       - info about channels (see chnsCompute.m)
%
% EXAMPLE
%
% See also acfReadme, acfDetect, acfDemoInria, acfModify, acfTest,
% chnsCompute, chnsPyramid, adaBoostTrain, bbGt, bbNms, jitterImage
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.40
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% initialize opts struct
opts = initializeOpts( varargin{:} );
if(nargin==0), detector=opts; return; end

% load or initialize detector and begin logging
nm=[opts.name 'Detector.mat']; t=exist(nm,'file');
if(t), if(nargout), t=load(nm); detector=t.detector; end; return; end
t=fileparts(nm); if(~isempty(t) && ~exist(t,'dir')), mkdir(t); end
detector = struct( 'opts',opts, 'clf',[], 'info',[] );
startTrain=clock; nm=[opts.name 'Log.txt'];
if(exist(nm,'file')), diary(nm); diary('off'); delete(nm); end; diary(nm);
RandStream.setGlobalStream(RandStream('mrg32k3a','Seed',opts.seed));

% iterate bootstraping and training
for stage = 0:numel(opts.nWeak)-1
  diary('on'); fprintf([repmat('-',[1 75]) '\n']);
  fprintf('Training stage %i\n',stage); startStage=clock;

  % sample positives and compute info about channels
  if( stage==0 )
    [Is1,IsOrig1] = sampleWins( detector, stage, 1 );
    t=ndims(Is1); if(t==3), t=Is1(:,:,1); else t=Is1(:,:,:,1); end
    t=chnsCompute2(t,opts.pPyramid.pChns); detector.info=t.info;
  end

  % compute local decorrelation filters
  if( stage==0 && length(opts.filters)==2 )
    fs = opts.filters; opts.filters = [];
    X1 = chnsCompute1( IsOrig1, opts );
    fs = chnsCorrelation( X1, fs(1), fs(2) );
    opts.filters = fs; detector.opts.filters = fs;
  end

  % compute lambdas
  if( stage==0 && isempty(opts.pPyramid.lambdas) )
    fprintf('Computing lambdas... '); start=clock;
    ds=size(IsOrig1); ds(1:end-1)=1;
    IsOrig1=IsOrig1(:,:,1:3,:);
    IsOrig1=mat2cell2(IsOrig1,ds);
    ls=chnsScaling(opts.pPyramid.pChns,IsOrig1,0);
    ls=round(ls*10^5)/10^5; detector.opts.pPyramid.lambdas=ls;
    fprintf('done (time=%.0fs).\n',etime(clock,start));
  end

  % compute features for positives
  if( stage==0 )
    X1 = chnsCompute1( Is1, opts );
    X1 = reshape(X1,[],size(X1,4))';
    clear Is1 IsOrig1 ls fs ds t;
  end

  % sample negatives and compute features
  Is0 = sampleWins( detector, stage, 0 );
  X0 = chnsCompute1( Is0, opts ); clear Is0;
  X0 = reshape(X0,[],size(X0,4))';

  % accumulate negatives from previous stages
  if( stage>0 )
    n0=size(X0p,1); n1=max(opts.nNeg,opts.nAccNeg)-size(X0,1);
    if(n0>n1 && n1>0), X0p=X0p(randSample(n0,n1),:); end
    if(n0>0 && n1>0), X0=[X0p; X0]; end %#ok<AGROW>
  end; X0p=X0;

  % train boosted clf
  detector.opts.pBoost.nWeak = opts.nWeak(stage+1);
  detector.clf = adaBoostTrain(X0,X1,detector.opts.pBoost);
  detector.clf.hs = detector.clf.hs + opts.cascCal;

  % update log
  fprintf('Done training stage %i (time=%.0fs).\n',...
    stage,etime(clock,startStage)); diary('off');
end

% save detector
save([opts.name 'Detector.mat'],'detector');

% finalize logging
diary('on'); fprintf([repmat('-',[1 75]) '\n']);
fprintf('Done training (time=%.0fs).\n',...
  etime(clock,startTrain)); diary('off');

end

function opts = initializeOpts( varargin )
% Initialize opts struct.
dfs= { 'pPyramid',{}, 'filters',[], ...
  'modelDs',[100 41], 'modelDsPad',[128 64], ...
  'pNms',struct(), 'stride',4, 'cascThr',-1, 'cascCal',.005, ...
  'nWeak',128, 'pBoost', {}, 'seed',0, 'name','', 'posGtDir','', ...
  'posImgDir','', 'negImgDir','', 'posWinDir','', 'negWinDir','', ...
  'imreadf',@imread, 'imreadp',{}, 'pLoad',{}, 'nPos',inf, 'nNeg',5000, ...
  'nPerNeg',25, 'nAccNeg',10000, 'pJitter',{}, 'winsSave',0, 'trainSet', [] };
opts = getPrmDflt(varargin,dfs,1);
% fill in remaining parameters
% p=chnsPyramid([],opts.pPyramid); p=p.pPyramid;
p=chnsPyramid2([],opts.pPyramid); p=p.pPyramid;
p.minDs=opts.modelDs; shrink=p.pChns.shrink;
opts.modelDsPad=ceil(opts.modelDsPad/shrink)*shrink;
p.pad=ceil((opts.modelDsPad-opts.modelDs)/shrink/2)*shrink;
p=chnsPyramid([],p); p=p.pPyramid; p.complete=1;
p.pChns.complete=1; opts.pPyramid=p;
% initialize pNms, pBoost, pBoost.pTree, and pLoad
dfs={ 'type','maxg', 'overlap',.65, 'ovrDnm','min' };
opts.pNms=getPrmDflt(opts.pNms,dfs,-1);
dfs={ 'pTree',{}, 'nWeak',0, 'discrete',1, 'verbose',16 };
opts.pBoost=getPrmDflt(opts.pBoost,dfs,1);
dfs={'nBins',256,'maxDepth',2,'minWeight',.01,'fracFtrs',1,'nThreads',16};
opts.pBoost.pTree=getPrmDflt(opts.pBoost.pTree,dfs,1);
opts.pLoad=getPrmDflt(opts.pLoad,{'squarify',{0,1}},-1);
opts.pLoad.squarify{2}=opts.modelDs(2)/opts.modelDs(1);
end

function [Is,IsOrig] = sampleWins( detector, stage, positive )
% Load or sample windows for training detector.
opts=detector.opts; start=clock;
if( positive ), n=opts.nPos; else n=opts.nNeg; end
if( positive ), crDir=opts.posWinDir; else crDir=opts.negWinDir; end
if( exist(crDir,'dir') && stage==0 )
  % if window directory is specified simply load windows
  fs=bbGt2('getFiles',{crDir}); nImg=length(fs); assert(nImg>0);
  if(nImg>n), fs=fs(:,randSample(nImg,n)); end; n=nImg;
  for i=1:n, fs{i}=[{opts.imreadf},fs(i),opts.imreadp]; end
  Is=cell(1,n); parfor i=1:n, Is{i}=feval(fs{i}{:}); end
else
  % sample windows from full images using sampleWins1()
  hasGt=positive||isempty(opts.negImgDir); fs={opts.negImgDir};
  if(hasGt), fs={opts.posImgDir,opts.posImgDir,opts.posGtDir}; end
  % fs=bbGt('getFiles',fs);
  fs=bbGt2('getSubsetFiles',fs, opts.trainSet);
%   fp=fopen(opts.trainset, 'r'); info=textscan(fp, '%s%s%s', 'Delimiter', '/');  fclose(fp);
%   rgbFile=cellfun(@(x,y,z) fullfile(opts.posImgDir, x, y, 'visible', [z '.jpg']), info{1}, info{2}, info{3}, 'uniformoutput', false );
%   lwirFile=cellfun(@(x,y,z) fullfile(opts.posImgDir, x, y, 'lwir', [z '.jpg']), info{1}, info{2}, info{3}, 'uniformoutput', false );
%   gtFile=cellfun(@(x,y,z) fullfile(opts.posGtDir, x, y, [z '.txt']), info{1}, info{2}, info{3}, 'uniformoutput', false );
%   fs = cat(2, rgbFile, lwirFile, gtFile)';

  nImg=size(fs,2); assert(nImg>0);
  if(~isinf(n)), fs=fs(:,randperm(nImg)); end; Is=cell(nImg*1000,1);
  diary('off'); tid=ticStatus('Sampling windows',1,30); k=0; i=0; batch=64;
  while( i<nImg && k<n )
    batch=min(batch,nImg-i); Is1=cell(1,batch);
    parfor j=1:batch, ij=i+j;
%     for j=1:batch, ij=i+j;
      I = feval(opts.imreadf,fs{1,ij}, fs{2,ij},opts.imreadp{:}); %#ok<PFBNS>
      gt=[]; if(hasGt), [~,gt]=bbGt2('bbLoad',fs{3,ij},opts.pLoad); end
      Is1{j} = sampleWins1( I, gt, detector, stage, positive );
    end
    Is1=[Is1{:}]; k1=length(Is1); Is(k+1:k+k1)=Is1; k=k+k1;
    if(k>n), Is=Is(randSample(k,n)); k=n; end
    i=i+batch; tocStatus(tid,max(i/nImg,k/n));
  end
  Is=Is(1:k); diary('on');
  fprintf('Sampled %i windows from %i images.\n',k,i);
end
% optionally jitter positive windows
if(length(Is)<2), Is={}; return; end
nd=ndims(Is{1})+1; Is=cat(nd,Is{:}); IsOrig=Is;
if( positive && isstruct(opts.pJitter) )
  opts.pJitter.hasChn=(nd==4); Is=jitterImage(Is,opts.pJitter);
  ds=size(Is); ds(nd)=ds(nd)*ds(nd+1); Is=reshape(Is,ds(1:nd));
end
% make sure dims are divisible by shrink and not smaller than modelDsPad
ds=size(Is); cr=rem(ds(1:2),opts.pPyramid.pChns.shrink); s=floor(cr/2)+1;
e=ceil(cr/2); Is=Is(s(1):end-e(1),s(2):end-e(2),:,:); ds=size(Is);
if(any(ds(1:2)<opts.modelDsPad)), error('Windows too small.'); end
% optionally save windows to disk and update log
nm=[opts.name 'Is' int2str(positive) 'Stage' int2str(stage)];
if( opts.winsSave ), save(nm,'Is','-v7.3'); end
fprintf('Done sampling windows (time=%.0fs).\n',etime(clock,start));
diary('off'); diary('on');
end

function Is = sampleWins1( I, gt, detector, stage, positive )
% Sample windows from I given its ground truth gt.
opts=detector.opts; shrink=opts.pPyramid.pChns.shrink;
modelDs=opts.modelDs; modelDsPad=opts.modelDsPad;
if( positive ), bbs=gt; bbs=bbs(bbs(:,5)==0,:); else
  if( stage==0 )
    % generate candidate bounding boxes in a grid
    [h,w,~]=size(I); h1=modelDs(1); w1=modelDs(2);
    n=opts.nPerNeg; ny=sqrt(n*h/w); nx=n/ny; ny=ceil(ny); nx=ceil(nx);
    [xs,ys]=meshgrid(linspace(1,w-w1,nx),linspace(1,h-h1,ny));
    bbs=[xs(:) ys(:)]; bbs(:,3)=w1; bbs(:,4)=h1; bbs=bbs(1:n,:);
  else
    % run detector to generate candidate bounding boxes
    bbs=acfDetect2(I,detector); [~,ord]=sort(bbs(:,5),'descend');
    bbs=bbs(ord(1:min(end,opts.nPerNeg)),1:4);
  end
  if( ~isempty(gt) )
    % discard any candidate negative bb that matches the gt
    n=size(bbs,1); keep=false(1,n);
    for i=1:n, keep(i)=all(bbGt2('compOas',bbs(i,:),gt,gt(:,5))<.1); end
    bbs=bbs(keep,:);
  end
end
% grow bbs to a large padded size and finally crop windows
modelDsBig=max(8*shrink,modelDsPad)+max(2,ceil(64/shrink))*shrink;
r=modelDs(2)/modelDs(1); assert(all(abs(bbs(:,3)./bbs(:,4)-r)<1e-5));
r=modelDsBig./modelDs; bbs=bbApply('resize',bbs,r(1),r(2));
Is=bbApply('crop',I,bbs,'replicate',modelDsBig([2 1]));
end

function chns = chnsCompute1( Is, opts )
% Compute single scale channels of dimensions modelDsPad.
if(isempty(Is)), chns=[]; return; end
fprintf('Extracting features... '); start=clock; fs=opts.filters;
pChns=opts.pPyramid.pChns; smooth=opts.pPyramid.smooth;
dsTar=opts.modelDsPad/pChns.shrink; ds=size(Is); ds(1:end-1)=1;
Is=squeeze(mat2cell2(Is,ds)); n=length(Is); chns=cell(1,n);
% parfor i=1:n
for i=1:n
  C=chnsCompute2(Is{i},pChns); C=convTri(cat(3,C.data{:}),smooth);
  if(~isempty(fs)), C=repmat(C,[1 1 size(fs,4)]);
    for j=1:size(C,3), C(:,:,j)=conv2(C(:,:,j),fs(:,:,j),'same'); end; end
  if(~isempty(fs)), C=imResample(C,.5); shr=2; else shr=1; end
  ds=size(C); cr=ds(1:2)-dsTar/shr; s=floor(cr/2)+1; e=ceil(cr/2);
  C=C(s(1):end-e(1),s(2):end-e(2),:); chns{i}=C;
end; chns=cat(4,chns{:});
fprintf('done (time=%.0fs).\n',etime(clock,start));
end

function filters = chnsCorrelation( chns, wFilter, nFilter )
% Compute filters capturing local correlations for each channel.
fprintf('Computing correlations... '); start=clock;
[~,~,m,n]=size(chns); w=wFilter; wp=w*2-1;
filters=zeros(w,w,m,nFilter,'single');
for i=1:m
  % compute local auto-scorrelation using Wiener-Khinchin theorem
  mus=squeeze(mean(mean(chns(:,:,i,:)))); sig=cell(1,n);
  parfor j=1:n
    T=fftshift(ifft2(abs(fft2(chns(:,:,i,j)-mean(mus))).^2));
    sig{j}=T(floor(end/2)+1-w+(1:wp),floor(end/2)+1-w+(1:wp));
  end
  sig=double(mean(cat(4,sig{mus>1/50}),4));
  sig=reshape(full(convmtx2(sig,w,w)),wp+w-1,wp+w-1,[]);
  sig=reshape(sig(w:wp,w:wp,:),w^2,w^2); sig=(sig+sig')/2;
  % compute filters for each channel from sig (sorted by eigenvalue)
  [fs,D]=eig(sig); fs=reshape(fs,w,w,[]);
  [~,ord]=sort(diag(D),'descend');
  fs=flipdim(flipdim(fs,1),2); %#ok<DFLIPDIM>
  filters(:,:,i,:)=fs(:,:,ord(1:nFilter));
end
fprintf('done (time=%.0fs).\n',etime(clock,start));
end

author: SoonminHwang

repo: rgbt-ped-detection


text:
 function y = mtimes_block(ob, x, iblock, nblock)
%function y = mtimes_block(ob, x, iblock, nblock)
%	y = G(i'th block) * x	or y = G'(i'th block) * x
%	in either case the project data will be "small"
%	iblock is 1,...,nblock

% support 'exists' option for seeing if this routine is available
if nargin == 2 & ischar(x) & streq(x, 'exists')
	y = 1;
	return
end

if nargin ~= 4, ir_usage, end

if nblock ~= 1 | iblock ~= 1
	error 'no block for mtimes_block'
end
y = ob * x;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MATLAB Code for                                                  %
%                                                                   %
%  Multi-Objective Particle Swarm Optimization (MOPSO)              %
%  Version 1.0 - Feb. 2011                                          %
%                                                                   %
%  According to:                                                    %
%  Carlos A. Coello Coello et al.,                                  %
%  "Handling Multiple Objectives with Particle Swarm Optimization," %
%  IEEE Transactions on Evolutionary Computation, Vol. 8, No. 3,    %
%  pp. 256-279, June 2004.                                          %
%                                                                   %
%  Developed Using MATLAB R2009b (Version 7.9)                      %
%                                                                   %
%  Programmed By: S. Mostapha Kalami Heris                          %
%                                                                   %
%         e-Mail: sm.kalami@gmail.com                               %
%                 kalami@ee.kntu.ac.ir                              %
%                                                                   %
%       Homepage: http://www.kalami.ir                              %
%                                                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function pop=DetermineDomination(pop)

    npop=numel(pop);

    for i=1:npop
        pop(i).Dominated=false;
        for j=1:i-1
            if ~pop(j).Dominated
                if Dominates(pop(i),pop(j))
                    pop(j).Dominated=true;
                elseif Dominates(pop(j),pop(i))
                    pop(i).Dominated=true;
                    break;
                end
            end
        end
    end

end

author: sfvsfv

repo: Mathematical-modeling


text:
function engine = belprop_mrf2_inf_engine(mrf2, varargin)
% BELPROP_MRF2_INF_ENGINE Belief propagation for MRFs with discrete pairwise potentials
% engine = belprop_mrf2_inf_engine(mrf2, ...)
%
% This is like belprop_inf_engine, except it is designed for mrf2, so is much faster.
%
% [ ... ] = belprop_mrf2_inf_engine(..., 'param1',val1, 'param2',val2, ...)
% allows you to specify optional parameters as name/value pairs.
% Parameters modifying behavior of enter_evidence are below [default value in brackets]
%
% max_iter - max. num. iterations [ 5*nnodes]
% momentum - weight assigned to old message in convex combination
%            (useful for damping oscillations) [0]
% tol      - tolerance used to assess convergence [1e-3]
% verbose - 1 means print error at every iteration [0]
%
% Parameters can be changed later using set_params

% The advantages of pairwise potentials are
% (1) we can compute messages using vector-matrix multiplication
% (2) we can easily specify the parameters: one potential per edge
% In contrast, potentials on larger cliques are more complicated to deal with.

nnodes = length(mrf2.adj_mat);

[engine.max_iter, engine.momentum, engine.tol, engine.verbose] = ...
    process_options(varargin, 'max_iter', [], 'momentum', 0, 'tol', 1e-3, ...
		   'verbose', 0);

if isempty(engine.max_iter) % no user supplied value, so compute default
  engine.max_iter = 5*nnodes;
  %if acyclic(mrf2.adj_mat, 0) --- can be very slow!
  %  engine.max_iter = nnodes;
  %else
  %  engine.max_iter = 5*nnodes;
  %end
end

engine.bel = cell(1, nnodes); % store results of enter_evidence here
engine.mrf2 = mrf2;

engine = class(engine, 'belprop_mrf2_inf_engine');

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef ZXH_CF9 < PROBLEM
% <multi/many> <real> <large/none> <constrained>
% Constrained benchmark MOP proposed by Zhou, Xiang, and He

%------------------------------- Reference --------------------------------
% Y. Zhou, Y. Xiang, and X. He, Constrained multiobjective optimization:
% Test problem construction and performance evaluations, IEEE Transactions
% on Evolutionary Computation, 2021, 25(1): 172-186.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        k;  % Number of constrained variables
    end
	methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = obj.M+10;  end
            obj.lower    = zeros(1,obj.D) + 1e-10;
            obj.upper    = ones(1,obj.D)  - 1e-10;
            obj.encoding = ones(1,obj.D);
            if obj.M <= 3
                obj.k = obj.M - 1;
            elseif obj.M > 3 && obj.M <= 8
                obj.k = floor(obj.M/2);
            else
                obj.k = 3;
            end
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            PopDec = varargin{1};
            OptX   = 0.2;
            [N,D]  = size(PopDec);
            M      = obj.M;
            % Step 1: Compute cumsum
            Sx = cumsum(PopDec(:,1:M).^2,2,'reverse');
            % Step 2: Compute theta
            THETA = 2/pi*atan(sqrt(Sx(:,2:end))./PopDec(:,1:M-1));
            % Step 3: Calculate Ackley function
            h = 20 - 20 * exp(-0.2 * sqrt(sum((PopDec(:,M+1:end)-OptX).^2,2)/(D-M))) + exp (1) - exp(sum(cos(2 * pi .*(PopDec(:,M+1:end)-OptX)),2)/(D-M));
            % Step 4: Compute T_
            T = (1 - Sx(:,1)).^2 + h;
            % Step 5: Objectives (mixed)
            A      = 2; % number of segments
            G      = 1-[ones(N,1) cumprod(sin(pi/2*THETA),2)] .* [cos(pi/2*THETA) ones(N,1)];
            G(:,1) = THETA(:,1) - cos(2*pi*A*THETA(:,1)+pi/2)/2/A/pi;
            PopObj = G .* repmat((1+T),1,M);
            % Step 6: Constraints
            PopCon(:,1) = Sx(:,1) + h - 1;
            for i = 1 : obj.k
                PopCon(:,i+1) = min(THETA(:,i)-1/4,3/4-THETA(:,i));
            end
            Population = SOLUTION(varargin{1},PopObj,PopCon,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,obj.M);
            R = 1 - R./repmat(sqrt(sum(R.^2,2)),1,obj.M);
            V = asin(sqrt(sum((1-R(:,2:end)).^2,2)));
            R(:,1) = (2/pi).*(V+1/8*sin(8*V));
            T = zeros(size(R));
            for i = obj.M-1 : -1 : 2
                T(:,i) = atan((1-R(:,i+1))./(1-R(:,i))./cos(T(:,i+1)));
            end
            T(:,1) = asin((1-R(:,2))./cos(T(:,2)));
            THETA  = T(:,1:obj.k)*2/pi;
            Valid  = all(THETA<=1/4|THETA>=3/4,2);
            R      = R(Valid,:);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,pi/2,30)';
                R = {1-sin(a)*cos(a'),1-sin(a)*sin(a'),1-cos(a)*ones(size(a'))};
                V = asin(sqrt((1-R{2}).^2+(1-R{3}).^2));
                R{1} = (2/pi).*(V+1/8*sin(8*V));
                T2 = atan((1-R{3})./(1-R{2}));
                T1 = asin((1-R{2})./cos(T2));
                THETA = cat(3,T1,T2)*2/pi;
                Valid = all(THETA<=1/4|THETA>=3/4,3);
                R{1}(~Valid) = nan;
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%  INTERNAL FUNCTION: separate terms for conditional forecasting
%
%  ::
%
%    [Ty,Te,Tsig,C]=SEPARATE_TERMS(T,sstate,ny,nx,nshocks,k,h)
%
%  Args:
%
%     - **T** [cell array]: solution impact
%     - **sstate** [cell array]: steady state
%     - **state_cols** [vector]: location of state variables excluding shocks
%     - **k** [numeric]: number of forward steps (anticipation)
%     - **nshocks** [numeric]: number of shocks
%
%  Returns:
%     :
%
%     - **Ty** [cell array]: includes square matrices for the impact of
%       endogenous variables
%     - **Te** [cell array]: includes matrices for the impact of shocks
%     - **Tsig** [cell array]: includes vectors summing the impact of
%       uncertainty and the trend
%     - **C** [cell array]: includes vectors summing the impact of
%       uncertainty, the trend and the steady state
%
%

author: jmaih


text:
% feldkamp_example.m
% example of how to use feldkamp.m for cone-beam CT reconstruction
% Copyright 2004-8-28, Nicole Caparanis, Patty Laskowsky, Taka Masuda,
% and Jeff Fessler, The University of Michigan

if ~isvar('proj'), disp 'proj'
	down = 8;
	nh = 256/down;
	nv = 240/down;
	na = 224/down;
	ds = 1024/nh;
	dt = ds;
	dis_src_det = 949.075;
	dis_iso_det = 408.075;
	dis_src_iso = dis_src_det - dis_iso_det;
%	dis_foc_src = inf; % flat detector panel
	dis_foc_src = 0; % arc detector
	offset_det_h = 0.25; % quarter detector
	offset_det_v = 0.0;
	horiz = ([-(nh-1)/2:(nh-1)/2]' - offset_det_h) * ds;
	verti = ([-(nv-1)/2:(nv-1)/2]' - offset_det_v) * dt;
	printf('rmax=%g', dis_src_iso*sin(atan(max(abs(horiz)) / dis_src_det)))

	ell = [ ...
		[20 0 0  150 150 200 0 0.01]; % 30cm diam "cylinder"
		[80 0 0 50 50 30  0 0.01]; % bone-like inserts
		[0 0 75  40 40 40  0 0.01];
		[0 70 0  30 30 30  0 0.01];
	];

	proj = ellipsoid_proj(ell, horiz, verti, na, ...
		dis_src_iso, dis_iso_det, dis_foc_src);

	nx = 256/down;
	ny = 240/down;
	nz = 200/down;
	dx = 2*down; dy = dx; dz = dx; % 2mm voxels * down-sampling
	xtrue = ellipsoids(nx, ny, nz, ell, dx, dy, dz);

	% cone-beam system geometry, generalized from fan-beam geometry.
	% see ASPIRE users guide under tech. reports on web page for details.
	args = arg_pair('system', NaN, 'nx', nx, 'ny', ny, 'nz', nz, ...
		'nv', nv, 'nh', nh, 'na', na, 'support', 'all', ...
		'orbit', 360, 'orbit_start', 0, ...
		'pixel_size', dx, 'ray_spacing', ds, 'strip_width', 0, ...
		'dis_src_det', dis_src_det, ...
		'dis_iso_det', dis_iso_det, ...
		'dis_foc_src', dis_foc_src, ...
		'offset_source', 0, ...
		'offset_det_h', offset_det_h, ...
		'offset_det_v', offset_det_v);

	pl=330;
	im(pl+1, xtrue, 'x true'), cbar
	im(pl+4, proj, 'true projections'), cbar
	drawnow
prompt
end

% noisy data and estimated line integrals
if ~isvar('li_hat'), disp 'li_hat'
	% noisy data, if blank scan value has been specified.
	if isvar('bi') && isvar('ri')
		yb = bi .* exp(-proj) + ri;
		yi = poisson(yb);
		li_hat = -log((yi-ri) ./ bi);
		li_hat(yi-ri <= 0) = 0; % fix: need something better here...
	else
		li_hat = proj; % noiseless
	end
end

% FDK cone-beam reconstruction
if ~isvar('xfdk'), disp 'fdk'
	mask = true([nx ny nz]);
	xfdk = feldkamp_old(li_hat, 'ramp', mask, args);
end

if 1
	% show results (off-center slices worse than central slice)
	im(pl+2, xfdk, 'FDK recon'), cbar
	im(pl+3, xfdk - xtrue, 'FDK error'), cbar

	subplot(pl+5)
	ix = 1:nx; iy = ceil(ny/2); iz = ceil(nz/2);
	plot(ix, xtrue(ix,iy,iz), '-', ix, xfdk(ix,iy,iz), '--')
	axis([1 nx -0.005 0.025]), legend('true', 'FDK recon', 2)
	title 'middle slice', xlabel 'ix'

	subplot(pl+6)
	iz=1:nz; ix = 1+floor(nx/2); iy = 1+floor(ny/2);
	plot(iz, squeeze(xtrue(ix,iy,iz)), '-', iz, squeeze(xfdk(ix,iy,iz)), '--')
	axis([1 nz -0.005 0.025]), legend('true', 'FDK recon', 2), xlabel 'iz'
	title(sprintf('profile at (ix,iy)=(%g,%g)', ix,iy))

prompt
end


text:
function filename = dtiWriteBrainVoyagerFibers(fiberGroups, filename, coordsType)
% Anyone use this for BV?  Or should we rename for Tal only?
%
% filename = dtiWriteBrainVoyagerFibers(fiberGroups, [filename], [coordsType])
%
% fiberGroups(grpNum).name = 'some text string'
% fiberGroups(grpNum).visible = 0 or 1 (boolean)
% fiberGroups(grpNum).thickness = fiber thickness (diameter? in mm?)
% fiberGroups(grpNum).colorRgb = [R G B] (0-255)
% fiberGroups(grpNum).fibers is a cell array of Nx3 lists of points (in
% XYZ Talairach space if coordsType=TAL)
%
% Writes the fibers (assumed to be in Talairach coords) to the specified
% file. Prompts the user for the file if none is given.
%
% coordsType should be 'TAL', 'SYS' or 'BVI'. Default is TAL.
%
% HISTORY:
%   2003.07.02 RFD (bob@white.stanford.edu) wrote it.
%
% Bob (c) Stanford VISTASOFT team.  Probably should be removed.

if(~exist('filename','var') | isempty(filename))
    [f, p] = uiputfile('fibers.fib', 'Save BrainVoyager fibers file...');
    if(isnumeric(f))
        disp('User cancelled.');
        filename = '';
        return;
    end
    filename = fullfile(p, f);
end
if(~exist('coordsType','var') | isempty(coordsType))
    coordsType = 'TAL';
end
numGroups = length(fiberGroups);
fid = fopen(filename, 'wt');
fprintf(fid, 'FileVersion: 1.0\n');
fprintf(fid, 'CoordsType:  %s\n', coordsType);
fprintf(fid, 'NrOfGroups:  %d\n\n', numGroups);
for(gNum=1:numGroups)
    fprintf(fid, 'Name:        %s\n', fiberGroups(gNum).name);
    fprintf(fid, 'Visible:     %d\n', fiberGroups(gNum).visible);
    fprintf(fid, 'Animate:     0\n');
    fprintf(fid, 'Thickness:   %6.2f\n', fiberGroups(gNum).thickness);
    fprintf(fid, 'Color:       %03d %03d %03d\n',  fiberGroups(gNum).colorRgb);
    numFibers = length(fiberGroups(gNum).fibers);
    fprintf(fid, 'NrOfFibers:  %d\n', numFibers);
    for(fNum=1:numFibers)
        numPoints = size(fiberGroups(gNum).fibers{fNum},1);
        fprintf(fid, 'NrOfPoints:  %d\n', numPoints);
        for(p=1:numPoints)
            fprintf(fid, '%7.3f %7.3f %7.3f\n', fiberGroups(gNum).fibers{fNum}(1,p), ...
                                                fiberGroups(gNum).fibers{fNum}(2,p), ...
                                                fiberGroups(gNum).fibers{fNum}(3,p));
        end
        fprintf(fid, '\n');
    end
    fprintf(fid, '\n\n');
end
fclose(fid);
return;

% Code from BrainVoyager:
%
% 	float pos_x, pos_y, pos_z, thickness;
% 	int n_groups, n_fibers, n_points, cr, cg, cb, FileVersion;
% 	bool SystemCoords;
% 	QString strTemp, strTemp2, CoordsType;
% 	// file header
% 	ar << "FileVersion: " << FibersFileVersion << "\n";
% 	if(FibersSystemCoords == 2)
% 		ar << "CoordsType:  BVI\n";
% 	else if(FibersSystemCoords == 1)
% 		ar << "CoordsType:  SYS\n";
% 	else
% 		ar << "CoordsType:  TAL\n";
% 	ar << "NrOfGroups:  " << NrOfFiberGroups << "\n\n";
% 	int g, i, p;
% 	for(g=0; g<NrOfFiberGroups; g++){
% 		ar << "Name:        " << QString(FiberGroups[g].Label) << "\n";
% 		ar << "Visible:     " << (int)FiberGroups[g].visible << "\n";
% 		ar << "Animate:     " << (int)FiberGroups[g].animate << "\n";
% 		ar << "Thickness:   " << FiberGroups[g].thickness << "\n";    // %6.2f
% 		ar << "Color:       " << FiberGroups[g].r << " " << FiberGroups[g].g << " " << FiberGroups[g].b << "\n"; // %3i  %3i  %3i
%
% 		n_fibers = FiberGroups[g].n_fibers;
%
% 		ar << "NrOfFibers:  " << n_fibers << "\n";
%
% 		for(i=0; i<n_fibers; i++)
% 		{
% 			n_points = FiberGroups[g].fibers[i].n_points;
%
% 			ar << "NrOfPoints:  " << n_points << "\n";
%
% 			for(p=0; p<n_points; p++)
% 			{
% 				ar << FiberGroups[g].fibers[i].pos_x[p] << " " << FiberGroups[g].fibers[i].pos_y[p] << " " << FiberGroups[g].fibers[i].pos_z[p] << "\n";
% 				// fprintf(fp_fbr, "%7.3f %7.3f %7.3f\n", Fibers[i].pos_x[p], Fibers[i].pos_y[p], Fibers[i].pos_z[p]);
% 			}
% 			ar << "\n";
% 		}
% 		ar << "\n\n";
% 	}
%
% 	f.close();

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function Retorno_Puntero()
    global Puntero
        vect_Punt = getfield(Puntero, 'translation');
        Pos_actual = vect_Punt(1,2);
            despl_retor = Pos_actual:0.01:2.5;
        for r=1:length(despl_retor)
            Puntero.translation = [ 0 despl_retor(r) vect_Punt(1,3)];
            pause(0.01);
        end
end

author: Sable

repo: mcbench-benchmarks


text:
function t = typeof(a)
%TYPEOF       Type of a
%
%   t = typeof(a)
%
%For details, see intval\typeof and intval\typeadj.
%

% written  05/22/09     S.M. Rump  taylor added
%

% taylor\@taylor\typeof:  a  must be taylor
  if isa(a.t,'intval')
    t = 'taylorintval';
  else
    t = 'taylor';
  end

repo: OpenMAS


text:
function ind = find(S1G,x,epsilon)
% find close points
%
% Syntax
%   ind = find(S1G,x,epsilon) % find all points in distance epsilon
%   ind = find(S1G,x)         % find closest point
%
% Input
%  S1G     - @S1Grid
%  x       - double
%  epsilon - double
%
% Output
%  ind - int32

if S1G(1).periodic, p = S1G.max - S1G.min; else p = 0; end

if nargin == 2
  ind = S1Grid_find(S1G.points(:),S1G.min,p,x);
else
  ind = S1Grid_find_region(S1G.points(:),S1G.min,p,x,epsilon);
end

author: mtex-toolbox

repo: mtex

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


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Program to reconstruct Phantom-Head Model using Algebraic    %%%%%%%%%
%%%% Reconstruction Method.                                               %
%%%%     This code is implemented By :                                   %%                                                    %%
%%%%     AUTHOR: SAYEDALI A SHAIKH,                                     %%%
%%%%     M.Tech.(CSE)                                                  %%%%
%%%%     BVBCET,HUBLI-580031, E-mail Id:sayedalishaikh@gmail.com      %%%%%
%%%%     Date: 28/08/2009                                       %%%%%%%%%%%
% % %   Guided by: Mr Shrinivas D Desai                         %%%%%%%%%%%
%%%     Associate Prof, Dept of ISE, BVBCET Hubli - 580031      %%%%%%%%%%%
%%%%    sd_desai@bvb.edu    9845275066                          %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Call Function Inputs To Enter Following Parameters
[size1,rotation,incr,s]=inputart();
%Orginal Image
I=phantom(size1);
[r1,c1]=size(I);
figure,imshow(I);
title('Original Phantom-Head Model Image ');

%Guessed Image
G=zeros(size(I));
[r2,c2]=size(G);
G2=zeros(size(G));

% Call Function To Pad The Original Image
[r3,c3,padIMG]=padO(I);
%figure,imshow(padIMG);
T1=padIMG;

% Call Function To Pad The Guessed Image
[r4,c4,padGIMG]=padG(G);
%figure,imshow(padGIMG);
T2=padGIMG;

% Calculate The Correction Factor
%Call Function To Choose The Denomenator Value w.r.t Angle and Increment

[z]=chooseart(rotation,incr);
z1=0;
z2=incr;
T3=T2;
THETA=0:incr:rotation;
s1=length(THETA);
for a=1:s1,
    org1=imrotate(T1,z1,'bilinear','crop');
% Call Function To Calculate The Row Sum and Column Sum for Original Image
    [rsumO,csumO]=calc_sum(org1,r3,z);
% Call Function To Calculate The Row Sum and Column Sum for Guessed Image
    [rsumG,csumG,padGIMG]=corr_factor(T2,r4,rsumO,csumO,r3,c3,s,z);
    %figure,imshow(padGIMG);
    G2=padGIMG;
    %figure,imshow(G2);
    T3=T3+G2;
    figure,imshow(T3);
    T3=imrotate(T3,z2,'bilinear','crop');
    z1=z1+incr;
end


text:
 function [Af, Aa, dif] = test_adjoint(A, varargin)
%function [Af, Aa, dif] = test_adjoint(A, varargin)
%|
%| test that the adjoint of a 'matrix object' A really is its transpose
%|
%| in
%|	A		Fatrix or fatrix2 typically
%|
%| option
%|	'tolre'	dbl	tolerance (relative) for equality (default: 0)
%|	'big'	0|1	for big objects, use random data
%|	'tol'	dbl	tolerance for adjoint for 'big' option (mpd)
%|	'nrep'	int	multiple realizations for 'big' option
%|	'complex' 0|1	if 1, test with complex data (default: 0)
%|	'warn'	0|1	if 1, just print warning. 0 (default), fail if >tol
%|	'chat'	0|1	default 0
%|
%| Copyright 2005-8-2, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if streq(A, 'test'), test_adjoint_test, return, end

arg.big = 0;
arg.tolre = 0;
arg.tol = 1e-5;
arg.nrep = 1;
arg.chat = 0;
arg.warn = false;
arg.complex = 0;
arg = vararg_pair(arg, varargin);

test_adjoint_big(A, arg.nrep, arg.tol, arg.complex, arg.warn, arg.chat);
if arg.big, return, end % stop here if 'big'

[nd, np] = size(A);

try
	% must be very careful here because ob(:,:) calls full(ob) by default
	% which in turns calls fatrix2_subsref_colon() which is "too" clever
	% because it uses the smaller dimension, so transpose has no effect.
	% 2013-07-16 added 'col' option to 'full' to force using the column
	% dimension for the loop (for both ob and ob') thereby truly testing.

	if arg.chat, printm 'A(:,:)', end
	if isa(A, 'Fatrix')
		Af = A(:,:); % not for fatrix2! calls full()
	else
		if isa(A, 'fatrix2')
			Af = full(A, 'col'); % trick
		else
			Af = full(A);
		end
	end
	if arg.chat, printm 'done', end

	if arg.chat, printm 'A''(:,:)', end
	Aa = A';
	if isa(Aa, 'Fatrix')
		Aa = Aa(:,:); % not for fatrix2! calls full()
	else
		if isa(A, 'fatrix2')
			Aa = full(Aa, 'col'); % trick
		else
			Af = full(A);
		end
	end
	if arg.chat, printm 'done', end

catch
	warn 'A(:,:) or its adjoint failed!?'

	Af = zeros(nd,np);
	Aa = A';
	% forward projection
	for jj=1:np
		x = zeros(np,1);
		x(jj) = 1;
		Af(:,jj) = A * x(:);
	end

	% back projection
	for ii=1:nd
		y = zeros(nd, 1);
		y(ii) = 1;
		Aa(:,ii) = A' * y(:);
	end
end

dif = Af - Aa';
relerr = max(abs(dif(:))) / max(abs(Af(:)));
if relerr > arg.tolre
	printm('adjoint of %s is imperfect (relerr=%g):', inputname(1), relerr)
	printm('adjoint real minmax: %g %g', minmax(real(dif)).')
	printm('adjoint imag minmax: %g %g', minmax(imag(dif)).')
	printm('adjoint real max diff = %g%%', max_percent_diff(real(Af), real(Aa')))
	printm('adjoint imag max diff = %g%%', max_percent_diff(imag(Af), imag(Aa')))
	printm('adjoint class=%s range: %g %g', class(Aa), minmax(Aa(:)).')
else
	tmp = class(A);
	if streq(tmp, 'Fatrix') || streq(tmp, 'fatrix2')
		tmp = sprintf('%s:%s', tmp, A.caller);
	end
	if arg.tolre
		printm('adjoint of %s matches within %g < %g, %s', ...
			inputname(1), relerr, arg.tolre, tmp)
	else
		printm('adjoint of %s appears numerically exact, %s', ...
			inputname(1), tmp)
	end
end

end % test_adjoint

% test_adjoint_big()
% test for big operators using random vectors
function test_adjoint_big(A, nrep, tol, do_complex, do_warn, chat)
rng(0)
for ii=1:nrep
	x = rand(size(A,2),1) - 0.5;
	y = rand(size(A,1),1) - 0.5;
	if do_complex
		x = x + 1i * rand(size(A,2),1) - 0.5;
		y = y + 1i * rand(size(A,1),1) - 0.5;
	end
	Ax = A * x;
	if ~isreal(Ax) && ~do_complex
		fail 'must test complex systems with ''complex'' option'
	end
%	v1 = y' * (A * x);
	v1 = ir_dot_double(conj(y), A * x);
%	v2 = (x' * (A' * y))';
	v2 = ir_dot_double(conj(x), A' * y)';

	mpd = max_percent_diff(v1, v2);
	if mpd/100 > tol
		pr v1
		pr v2
		pr '[mpd/100 tol]'
		if do_warn
			warn 'adjoint mismatch'
		else
			fail 'adjoint mismatch'
		end

	elseif chat
		pr v1
		pr v2
	end
end

end % test_adjoint_big

function test_adjoint_test
	B = randn(4,3);
	test_adjoint(B);
	test_adjoint(B, 'big', true);
	B = B + 2i;
	test_adjoint(B, 'complex', true);
end


text:
function [ r, center ] = circle_dia2imp_2d ( p1, p2 )

%*****************************************************************************80
%
%% CIRCLE_DIA2IMP_2D converts a diameter to an implicit circle in 2D.
%
%  Discussion:
%
%    The diameter form of a circle is:
%
%      P1 and P2 are the endpoints of a diameter.
%
%    The implicit form of a circle in 2D is:
%
%      ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 = R^2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), two points that are the
%    endpoints of a diameter of the circle.
%
%    Output, real R, the radius of the circle.
%
%    Output, real CENTER(2,1), the center of the circle.
%
  r = 0.5 * sqrt ( sum ( ( p2(1:2,1) - p1(1:2,1) ).^2 ) );

  center(1:2,1) = 0.5 * ( p1(1:2,1) + p2(1:2,1) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [W,B] = updateWeight(obj,VOV,N)
% Update the weight vectors

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Tomoaki Takagi

    fmin = min(obj,[],1);
    fmax = max(obj,[],1);
    obj = unique([normalize(obj,'range');VOV],'rows');

    %% Select the representative objective vectors
    LpNormD = pdist2(obj,obj,'minkowski',0.5);
    Choose  = false(1,size(obj,1));
    % Select the extreme objective vectors
    [~,extreme]     = min(pdist2(obj,eye(size(obj,2)),'cosine'),[],1);
    Choose(extreme) = true;
    % Greedy inclusion distance-based subset selection
    while sum(Choose) < N
        Remain   = find(~Choose);
        [~, rho] = max(min(LpNormD(Remain,Choose),[],2));
        Choose(Remain(rho)) = true;
    end
    obj = obj(Choose,:);

    %% Update the weight vectors
    W = obj.*repmat(fmax-fmin,N,1);
    B = pdist2(W,W);
    [~,B] = sort(B,2);
    B = B(:,1:ceil(N/10));
    W = W./vecnorm(W,1,2);
end


text:
function [uOutput] = import_california_fm(nFunction, sFilename)
  % import NCEDC NCSN-FPFIT catalog from the UCB, CA, and the Hauksson relocated catalog for southern California.
% [uOutput] = import_california_fm(nFunction, sFilename)
%
% See:
% http://quake.geo.berkeley.edu
% http://www.data.scec.org/ftp/catalogs/hauksson/Socal_focal/
%
% updated: 17.05.2005, Danijel Schorlemmer
% schorlemmer@sed.ethz.ch

% Counter for double events
nCount = 0;

% Filter function switchyard
if nFunction == FilterOp.getDescription
  uOutput = 'California focal mechanism catalogs (NCEDC NCSN-FPFIT/SCSN Hauksson) raw format/multiple solutions eliminated';
elseif nFunction == FilterOp.importCatalog
  % Read formatted data
  mData = textread(sFilename, '%s', 'delimiter', '\n', 'whitespace', '');
  % Create empty catalog
  uOutput = zeros(length(mData), 24);
  % Loop thru all lines of catalog and convert them
  for i = 1:length(mData)
    if rem(i,100) == 0
      % Notify the user every 100 processed events
      disp([num2str(i) ' of ' num2str(length(mData)) ' events processed.']);
    end
    try
      % Replace whitespace in first column with '0'
      l = find(mData{i} == ' ');
      mData{i}(l) = '0';
      % Read values into catalog
      uOutput(i,1) = -str2num(mData{i}(29:32)) - str2num(mData{i}(34:38))/60 ;  % Longitude
      uOutput(i,2) =  str2num(mData{i}(20:22)) + str2num(mData{i}(24:28))/60 ;  % Latitude
      uOutput(i,3) = str2num(mData{i}(1:4));      % Year
      uOutput(i,4) = str2num(mData{i}(5:6));      % Month
      uOutput(i,5) = str2num(mData{i}(7:8));      % Day
      uOutput(i,6) = str2num(mData{i}(48:52));    % Magnitude
      uOutput(i,7) = str2num(mData{i}(39:45));    % Depth
      uOutput(i,8) = str2num(mData{i}(10:11));    % Hour
      uOutput(i,9) = str2num(mData{i}(12:13));    % Minute
      uOutput(i,10) = str2num(mData{i}(15:19));   % Second
      uOutput(i,11) = nan;                        % Reserved for cross-section values
      uOutput(i,13) = str2num(mData{i}(84:86));   % Dip direction (Plane 1)
      uOutput(i,12) = uOutput(i,13) - 90;         % Strike (Plane 1)
      uOutput(i,12) = mod(uOutput(i,12) + 360, 360);
      uOutput(i,14) = str2num(mData{i}(88:89));   % Dip (Plane 1)
      uOutput(i,15) = str2num(mData{i}(90:93));   % Rake (Plane 1)
      % Second nodal plane
      % uOutput(i,16) = Strike (Plane 2)
      % uOutput(i,17) = Dip direction (Plane 2)
      % uOutput(i,18) = Dip (Plane 2)
      % uOutput(i,19) = Rake (Plane 2)
      [uOutput(i,16), uOutput(i,18), uOutput(i,19), uOutput(i,17)] ...
        = ComputeSecondPlane(uOutput(i,12), uOutput(i,14), uOutput(i,15));
      uOutput(i,20) = str2num(mData{i}(72:74));   % Horizontal location error [km]
      uOutput(i,21) = str2num(mData{i}(77:79));   % Vertical location error [km]
      uOutput(i,22) = str2num(mData{i}(96:99));   % Solution misfit value. 0=perfect fit, 1=perfect misfit (never exceeds 0.5 in reality).
      uOutput(i,23) = str2num(mData{i}(132:141)); % Event ID # found in Hypoinverse archive file
      uOutput(i,24) = str2num(mData{i}(111:114)); % (STDR) Station distribution ratio (0-1). Lower numbers indicate data lie near nodal planes.
      uOutput(i,25) = str2num(mData{i}(101:103)); % Number of first motion observations used in solution
      uOutput(i,26) = str2num(mData{i}(122:123)); % Maximum half-width of 90% confidence range of strike
      uOutput(i,27) = str2num(mData{i}(125:126)); % Maximum half-width of 90% confidence range of dip
      uOutput(i,28) = str2num(mData{i}(128:129)); % Maximum half-width of 90% confidence range of rake
      % In case of multiple solutions, take only the best solution
      try
        if uOutput(i,23) == uOutput(i-1,23)       % Same event?
            nCount=nCount+1;
          if uOutput(i,24) > uOutput(i-1,24)      % Higher STDR?
            uOutput(i-1,1) = nan;                 % Invalidate the former event
          elseif uOutput(i,24) == uOutput(i-1,24) % Same STDR?
            if uOutput(i,22) < uOutput(i-1,22)    % Better solution misfit?
              uOutput(i-1,1) = nan;               % Invalidate the former event
            else                                  % Worse or equal solution misfit?
              uOutput(i,:) = uOutput(i-1,:);      % Copy the former event into current one
              uOutput(i-1,1) = nan;               % Invalidate the former event
            end
          else                                    % Lower STDR?
            uOutput(i,:) = uOutput(i-1,:);        % Copy the former event into current one
            uOutput(i-1,1) = nan;                 % Invalidate the former event
          end
        end
      catch
      end
      %Create decimal year
      uOutput(i,3) = decyear([uOutput(i,3) uOutput(i,4) uOutput(i,5) uOutput(i,8) uOutput(i,9) uOutput(i,10)]);
    catch
      msg.dbfprintf('Import: Problem in line %d of %s. Line ignored.\n',i, sFilename);
      uOutput(i,:) = nan;
    end
  end
  % Delete all invalidated events from catalog
  vSel = isnan(uOutput(:,1));
  uOutput(vSel,:) = [];
end

% --- Helper function

function [fStrike2, fDip2, fRake2, fDipDir2] = ComputeSecondPlane(fStrike1, fDip1, fRake1)

% Do not modify strike (Dip and rake need to be modified in order to
% minimize computational errors)
fDip1 = fDip1 + 0.000001;
fRake1 = mod((fRake1 + 0.000001) + 360, 360);
[fStrike2, fDip2, fRake2, fDipDir2, nError] = focal_pl2pl(fStrike1, fDip1, fRake1);

author: CelsoReyes


text:
function opt_vect = foptions()
% FOPTIONS Sets default parameters for optimisation routines
% For compatibility with MATLAB's foptions()
%
% Copyright (c) Dharmesh Maniyar, Ian T. Nabney (2004)

opt_vect      = zeros(1, 18);
opt_vect(2:3) = 1e-4;
opt_vect(4)   = 1e-6;
opt_vect(16)  = 1e-8;
opt_vect(17)  = 0.1;


text:
function [wordnet, counts] = findsenses(D, name)
%
% Looks at the wordnet field in D and return all the different senses that
% are returned when using 'name' for query.
%
% Example:
%  wordnet = findsenses(D, 'screen-device');
%    shutter. blind, screen. protective covering, protective cover, protection. covering. artifact, artefact. whole, unit. object, physical object. physical entity. entity.
%    windshield, windscreen. screen. protective covering, protective cover, protection. covering. artifact, artefact. whole, unit. object, physical object. physical entity. entity.
%
%  wordnet = findsenses(D, 'screen+device');
%    computer screen, computer display. screen, CRT screen. display, video display. electronic device. device. instrumentality, instrumentation. artifact, artefact. whole, unit. object, physical object. physical entity. entity.
%    screen, CRT screen. display, video display. electronic device. device. instrumentality, instrumentation. artifact, artefact. whole, unit. object, physical object. physical entity. entity.

if nargin == 2
    D = LMquery(D, 'object.synset', name);
end

Nannotations = length(D);

wordnet = [];
for i = 1:Nannotations
    if isfield(D(i).annotation, 'object')
        if isfield(D(i).annotation.object, 'synset')
            N = length(D(i).annotation.object);
            wordnet = [wordnet {D(i).annotation.object(:).synset}];
        end
    end
end

% Extract unique senses and counts
[wordnet, i, j] = unique(wordnet);
Nsenses = length(wordnet);
[counts, x] = hist(j, 1:Nsenses);

author: CSAILVision

repo: LabelMeToolbox