text:
function varargout = panel_cluster(varargin)
% PANEL_CLUSTER: Create a panel to add/remove/edit clusters of sensors attached to a given 3DViz figure.
%
% USAGE:  bstPanelNew = panel_cluster('CreatePanel')
%                       panel_cluster('UpdatePanel')
%                       panel_cluster('CurrentFigureChanged_Callback')
%
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2023

eval(macro_method);
end

%% ===== CREATE PANEL =====
function bstPanelNew = CreatePanel() %#ok<DEFNU>
    panelName = 'Cluster';
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    import org.brainstorm.icon.*;
    % Create tools panel
    jPanelNew = gui_component('Panel');
    TB_DIM = java_scaled('dimension',25,25);
    % Font size for the lists
    fontSize = round(11 * bst_get('InterfaceScaling') / 100);

    % ===== TOOLBAR =====
    jMenuBar = gui_component('MenuBar', jPanelNew, BorderLayout.NORTH);
        jToolbar = gui_component('Toolbar', jMenuBar);
        jToolbar.setPreferredSize(TB_DIM);
        jToolbar.setOpaque(0);
        % First buttons
        gui_component('ToolbarButton', jToolbar,[],[], IconLoader.ICON_NEW_SEL, 'Create new cluster: mouse selection.
Use sensors selected in any time series or topography figure.', @(h,ev)CreateNewCluster('Selection'));
        gui_component('ToolbarButton', jToolbar,[],[], IconLoader.ICON_NEW_IND, 'Create new cluster: list of indices/names.
Type the list of sensors indices, names or types you want in the new cluster.', @(h,ev)CreateNewCluster('Indices'));
        gui_component('ToolbarButton', jToolbar,[],[], IconLoader.ICON_TS_DISPLAY, 'Display cluster time series    [ENTER]', @(h,ev)view_clusters());

        % === MENU EDIT ===
        jMenu = gui_component('Menu', jMenuBar, [], 'Edit', IconLoader.ICON_MENU, [], [], 11);
        jMenu.setBorder(BorderFactory.createEmptyBorder(0,2,0,2));
            % Set color
            gui_component('MenuItem', jMenu, [], 'Set color', IconLoader.ICON_COLOR_SELECTION, [], @(h,ev)bst_call(@EditClusterColor));
            % Menu: Set cluster function
            jMenuTs = gui_component('Menu', jMenu, [], 'Set function', IconLoader.ICON_PROPERTIES, [], []);
                gui_component('MenuItem', jMenuTs, [], 'Mean',    [], [], @(h,ev)SetClusterFunction('Mean'));
                gui_component('MenuItem', jMenuTs, [], 'Mean+Std',    [], [], @(h,ev)SetClusterFunction('Mean+Std'));
                gui_component('MenuItem', jMenuTs, [], 'Mean+StdErr',    [], [], @(h,ev)SetClusterFunction('Mean+StdErr'));
                gui_component('MenuItem', jMenuTs, [], 'PCA',     [], [], @(h,ev)SetClusterFunction('PCA'));
                gui_component('MenuItem', jMenuTs, [], 'FastPCA', [], [], @(h,ev)SetClusterFunction('FastPCA'));
                gui_component('MenuItem', jMenuTs, [], 'Max',     [], [], @(h,ev)SetClusterFunction('Max'));
                gui_component('MenuItem', jMenuTs, [], 'Power',   [], [], @(h,ev)SetClusterFunction('Power'));
                gui_component('MenuItem', jMenuTs, [], 'All',     [], [], @(h,ev)SetClusterFunction('All'));
            jMenu.addSeparator();
            gui_component('MenuItem', jMenu, [], 'Rename   [Double-click]', IconLoader.ICON_EDIT,    [], @(h,ev)EditClusterLabel);
            gui_component('MenuItem', jMenu, [], 'Remove   [DEL]',          IconLoader.ICON_DELETE,  [], @(h,ev)RemoveClusters);
            gui_component('MenuItem', jMenu, [], 'Deselect all  [ESC]',     IconLoader.ICON_RELOAD,  [], @(h,ev)SetSelectedClusters(0));
            jMenu.addSeparator();
            gui_component('MenuItem', jMenu, [], 'Export to Matlab', IconLoader.ICON_MATLAB_EXPORT, [], @(h,ev)bst_call(@ExportClustersToMatlab));
            gui_component('MenuItem', jMenu, [], 'Import from Matlab', IconLoader.ICON_MATLAB_IMPORT, [], @(h,ev)bst_call(@ImportClustersFromMatlab));
            jMenu.addSeparator();
            gui_component('MenuItem', jMenu, [], 'Copy to other folders', IconLoader.ICON_COPY,  [], @(h,ev)CopyClusters('AllConditions'));
            gui_component('MenuItem', jMenu, [], 'Copy to other subjects', IconLoader.ICON_COPY,  [], @(h,ev)CopyClusters('AllSubjects'));

    % ===== PANEL MAIN =====
    jPanelMain = java_create('javax.swing.JPanel');
    jPanelMain.setLayout(BoxLayout(jPanelMain, BoxLayout.Y_AXIS));
    jPanelMain.setBorder(BorderFactory.createEmptyBorder(7,7,7,7));
        % ===== FIRST PART =====
        jPanelFirstPart = gui_component('Panel');
        jPanelFirstPart.setPreferredSize(java_scaled('dimension', 100,180));
        jPanelFirstPart.setMaximumSize(java_scaled('dimension', 500,180));
            % ===== Vertical Toolbar =====
            jToolbar2 = gui_component('Toolbar', jPanelFirstPart, BorderLayout.EAST);
            jToolbar2.setOrientation(jToolbar2.VERTICAL);
                gui_component('ToolbarButton', jToolbar2,[],[], IconLoader.ICON_FOLDER_OPEN, 'Load clusters file', @LoadClusters);
                gui_component('ToolbarButton', jToolbar2,[],[], IconLoader.ICON_SAVE, 'Save clusters file', @SaveClusters);

            % ===== Clusters list =====
             jPanelList = gui_component('Panel');
                jBorder = java_scaled('titledborder', 'Available clusters');
                jPanelList.setBorder(jBorder);

                jListClusters = java_create('org.brainstorm.list.BstClusterList');
                jListClusters.setCellRenderer(java_create('org.brainstorm.list.BstClusterListRenderer', 'I', fontSize));
                java_setcb(jListClusters, 'ValueChangedCallback', @ListValueChanged_Callback, ...
                                          'KeyTypedCallback',     @ListKeyTyped_Callback, ...
                                          'MouseClickedCallback', @ListClick_Callback);
                jPanelScrollList = java_create('javax.swing.JScrollPane');
                jPanelScrollList.getLayout.getViewport.setView(jListClusters);
                jPanelScrollList.setBorder([]);
                jPanelList.add(jPanelScrollList);
            jPanelFirstPart.add(jPanelList, BorderLayout.CENTER);
        jPanelMain.add(jPanelFirstPart);

        % ===== Clusters options panel =====
        jPanelOptions = gui_river([0,3], [0,5,10,3], 'Options');
            % Add extra space when not on a Mac
            if strncmp(computer,'MAC',3)
                strSpace = '';
            else
                strSpace = '   ';
            end
            % OPTIONS : Cluster size in number of sensors
            gui_component('Label', jPanelOptions, [], 'Number of sensors:');
            jLabelClusterSize = gui_component('Label', jPanelOptions, 'hfill', '');
            % OPTIONS : Overlay clusters/conditions
            gui_component('Label', jPanelOptions, 'br', ['Overlay:' strSpace]);
            jCheckOverlayClusters   = gui_component('CheckBox', jPanelOptions, [], 'Cluster', {Insets(0,0,0,0)});
            jCheckOverlayConditions = gui_component('CheckBox', jPanelOptions, [], 'Files',   {Insets(0,0,0,0)});
        panelPrefSize = jPanelOptions.getPreferredSize();
        jPanelOptions.setMaximumSize(Dimension(32000, panelPrefSize.getHeight()));
        jPanelMain.add(jPanelOptions);
        jPanelMain.add(Box.createVerticalGlue());
    jPanelNew.add(jPanelMain);

    % Create the BstPanel object that is returned by the function
    % => constructor BstPanel(jHandle, panelName, sControls)
    bstPanelNew = BstPanel(panelName, ...
                           jPanelNew, ...
                           struct('jPanelList',       jPanelList, ...
                                  'jPanelOptions',    jPanelOptions, ...
                                  'jMenuBar',         jMenuBar, ...
                                  'jToolbar',         jToolbar, ...
                                  'jToolbar2',        jToolbar2, ...
                                  'jLabelClusterSize',       jLabelClusterSize, ...
                                  'jCheckOverlayClusters',   jCheckOverlayClusters, ...
                                  'jCheckOverlayConditions', jCheckOverlayConditions, ...
                                  'jListClusters',           jListClusters));
end

%% =================================================================================
%  === CONTROLS CALLBACKS  =========================================================
%  =================================================================================

%% ===== LIST SELECTION CHANGED CALLBACK =====
function ListValueChanged_Callback(h, ev)
    if ~ev.getValueIsAdjusting()
        % Update number of sensors
        UpdateProperties();
        % Update mouse selection in all the Datasets
        UpdateChannelSelection();
    end
end

%% ===== LIST KEY TYPED CALLBACK =====
function ListKeyTyped_Callback(h, ev)
    switch(uint8(ev.getKeyChar()))
        % DELETE
        case {ev.VK_DELETE, ev.VK_BACK_SPACE}
            RemoveClusters();
        case ev.VK_ENTER
            tree_view_clusters();
        case ev.VK_ESCAPE
            SetSelectedClusters([]);
    end
end

%% ===== LIST CLICK CALLBACK =====
function ListClick_Callback(h, ev)
    % If DOUBLE CLICK
    if (ev.getClickCount() == 2)
        % Rename selection
        EditClusterLabel();
    end
end

%% =================================================================================
%  === EXTERNAL PANEL CALLBACKS  ===================================================
%  =================================================================================
%% ===== CURRENT FIGURE CHANGED =====
function CurrentFigureChanged_Callback(hFig)
    % Update list of clusters in the panel
    UpdatePanel();
    % Update sensor selection
    UpdateChannelSelection();
end

%% ===== UPDATE CALLBACK =====
function UpdatePanel()
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Cluster');
    if isempty(ctrl)
        return;
    end
    % Update clusters JList
    UpdateClustersList();
end

%% ===== UPDATE CLUSTERS LIST =====
function UpdateClustersList()
    import org.brainstorm.list.*;
    % Get "Cluster" panel controls
    ctrl = bst_get('PanelControls', 'Cluster');
    if isempty(ctrl)
        return;
    end
    % Get clusters
    sClusters = GetClusters();

    % Remove temporarily the list callback
    callbackBak = java_getcb(ctrl.jListClusters, 'ValueChangedCallback');
    java_setcb(ctrl.jListClusters, 'ValueChangedCallback', []);
    % Get selected clusters
    SelNames = {};
    selValues = ctrl.jListClusters.getSelectedValues();
    for i = 1:length(selValues)
        SelNames{end+1} = char(selValues(i));
    end
    if ~isempty(SelNames) && ~isempty(sClusters)
        iSelected = find(ismember({sClusters.Label}, SelNames));
    else
        iSelected = [];
    end
    % Create a new empty list
    listModel = java_create('javax.swing.DefaultListModel');
    % Get font with which the list is rendered
    fontSize = round(11 * bst_get('InterfaceScaling') / 100);
    jFont = java.awt.Font('Dialog', java.awt.Font.PLAIN, fontSize);
    tk = java.awt.Toolkit.getDefaultToolkit();
    % Add an item in list for each cluster
    Wmax = 0;
    for i = 1:length(sClusters)
        itemType  = sClusters(i).Function;
        itemText  = sClusters(i).Label;
        itemColor = sClusters(i).Color;
        listModel.addElement(BstListItem(itemType, [], itemText, i, itemColor(1), itemColor(2), itemColor(3)));
        % Get longest string
        W = tk.getFontMetrics(jFont).stringWidth(sClusters(i).Label);
        if (W > Wmax)
            Wmax = W;
        end
    end
    % Update list model
    ctrl.jListClusters.setModel(listModel);
    % Update cell rederer based on longest channel name
    ctrl.jListClusters.setCellRenderer(java_create('org.brainstorm.list.BstClusterListRenderer', 'II', fontSize, Wmax + 28));
    % Select previously selected clusters
    if ~isempty(iSelected)
        ctrl.jListClusters.setSelectedIndices(iSelected - 1);
    end
    % Reset cluster comments
    ctrl.jLabelClusterSize.setText('');

    % Restore callback
    drawnow;
    java_setcb(ctrl.jListClusters, 'ValueChangedCallback', callbackBak);
end

%% ===== UPDATE CLUSTER PROPERTIES DISPLAY =====
function UpdateProperties()
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Cluster');
    if isempty(ctrl)
        return;
    end
    % Get selected clusters
    sClusters = GetSelectedClusters();
    % Add all the selected clusters to compute
    nbSensors = 0;
    for i = 1:length(sClusters)
        nbSensors = nbSensors + length(sClusters(i).Sensors);
    end
    % Format results
    if (nbSensors == 0)
        strSize = '';
    else
        strSize = sprintf(' [ %d ]', nbSensors);
    end
    % Update panel
    ctrl.jLabelClusterSize.setText(strSize);
end

%% ===== UPDATE CHANNEL SELECTION =====
function UpdateChannelSelection()
    % Get selected clusters
    [sClusters, iClusters] = GetSelectedClusters();
    if isempty(sClusters)
        return;
    end
    % Join all the clusters contents
    SensorsNames = unique([sClusters.Sensors]);
    % Update sensors list
    bst_figures('SetSelectedRows', SensorsNames, 0);
end

%% ===== GET CLUSTERS =====
% USAGE:  [sClusters, iDSall, iFigall, hFigall] = panel_cluster('GetClusters')
function [sClusters, iDSall, iFigall, hFigall] = GetClusters()
    global GlobalData;
    % Get current figure
    [hFigall,iFigall,iDSall] = bst_figures('GetCurrentFigure');
    % Check if there are electrodes defined for this file
    if isempty(hFigall) || isempty(GlobalData.DataSet(iDSall).Clusters) || isempty(GlobalData.DataSet(iDSall).Clusters)
        sClusters = [];
        return;
    end
    % Return all the available electrodes
    sClusters = GlobalData.DataSet(iDSall).Clusters;
    ChannelFile = GlobalData.DataSet(iDSall).ChannelFile;
    % Get all the figures that share this channel file
    for iDS = 1:length(GlobalData.DataSet)
        % Skip if not the correct channel file
        if ~file_compare(GlobalData.DataSet(iDS).ChannelFile, ChannelFile)
            continue;
        end
        % Get all the figures
        for iFig = 1:length(GlobalData.DataSet(iDS).Figure)
            if ((iDS ~= iDSall(1)) || (iFig ~= iFigall(1))) && ismember(GlobalData.DataSet(iDS).Figure(iFig).Id.Type, {'DataTimeSeries', '3DViz', 'Topography'})
                iDSall(end+1) = iDS;
                iFigall(end+1) = iFig;
                hFigall(end+1) = GlobalData.DataSet(iDS).Figure(iFig).hFigure;
            end
        end
    end
end

%% ===== SET CLUSTER =====
% USAGE:  iClusters = SetClusters(iClusters=[], sClusters)
%         iClusters = SetClusters('Add', sClusters)
function iClusters = SetClusters(iClusters, sClusters)
    global GlobalData;
    % Parse input
    isAdd = ~isempty(iClusters) && ischar(iClusters) && strcmpi(iClusters, 'Add');
    % Get dataset
    [sClustersOld, iDSall] = GetClusters();
    % If there is no selected dataset
    if isempty(iDSall)
        return;
    end
    % Perform operations only once per dataset
    iDSall = unique(iDSall);
    for iDS = iDSall
        % Reset clusters list
        if isempty(sClusters)
            GlobalData.DataSet(iDS).Clusters(iClusters) = [];
        % Set clusters
        else
            % Add new clusters
            if isAdd
                iClusters = length(GlobalData.DataSet(iDS).Clusters) + (1:length(sClusters));
                % Add clusters
                for i = 1:length(sClusters)
                    % Default cluster name
                    if isempty(sClusters(i).Label)
                        sClusters(i).Label = sprintf('c%d', iClusters(i));
                    end
                    % Make new cluster names unique
                    if ~isempty(GlobalData.DataSet(iDS).Clusters)
                        sClusters(i).Label = file_unique(sClusters(i).Label, {GlobalData.DataSet(iDS).Clusters.Label, sClusters(1:i-1).Label});
                    end
                end
            end
            % Set clusters in global structure
            if isempty(GlobalData.DataSet(iDS).Clusters)
                GlobalData.DataSet(iDS).Clusters = sClusters;
            else
                GlobalData.DataSet(iDS).Clusters(iClusters) = sClusters;
            end
        end
        % Add color if not defined yet
        for i = 1:length(GlobalData.DataSet(iDS).Clusters)
            if isempty(GlobalData.DataSet(iDS).Clusters(i).Color)
                ColorTable = panel_scout('GetScoutsColorTable');
                iColor = mod(i-1, length(ColorTable)) + 1;
                GlobalData.DataSet(iDS).Clusters(i).Color = ColorTable(iColor,:);
            end
        end
    end
    % Mark channel file as modified (only in first dataset)
    GlobalData.DataSet(iDSall(1)).isChannelModified = 1;
end

%% ===== GET SELECTED CLUSTERS =====
% NB: Returned indices are indices in Clusters array
function [sSelClusters, iSelClusters] = GetSelectedClusters()
    sSelClusters = [];
    iSelClusters = [];
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Cluster');
    if isempty(ctrl)
        return;
    end
    % Get current clusters
    sClusters = GetClusters();
    if isempty(sClusters)
        return
    end
    % Get JList selected indices
    iSelClusters = uint16(ctrl.jListClusters.getSelectedIndices())' + 1;
    if isempty(iSelClusters)
        return
    end
    sSelClusters = sClusters(iSelClusters);
end

%% ===== SET SELECTED CLUSTERS =====
% WARNING: Input indices are references in the Clusters array, not in the JList
function SetSelectedClusters(iSelClusters, isUpdateMouseSel)
    if (nargin < 2) || isempty(isUpdateMouseSel)
        isUpdateMouseSel = 1;
    end
    % === CHECK FOR MODIFICATIONS ===
    iSelClusters = iSelClusters - 1;
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Cluster');
    if isempty(ctrl)
        return;
    end
    % Get previous selection
    iPrevItems = ctrl.jListClusters.getSelectedIndices();
    % If selection did not change: exit
    if isequal(iPrevItems, iSelClusters) || (isempty(iPrevItems) && isequal(iSelClusters, -1))
        return
    end

    % === UPDATE SELECTION ===
    % Temporality disables JList selection callback
    jListCallback_bak = java_getcb(ctrl.jListClusters, 'ValueChangedCallback');
    java_setcb(ctrl.jListClusters, 'ValueChangedCallback', []);
    % Select items in JList
    if isempty(iSelClusters)
        ctrl.jListClusters.setSelectedIndices(-1);
    else
        ctrl.jListClusters.setSelectedIndices(iSelClusters);
    end
    % Restore JList callback
    java_setcb(ctrl.jListClusters, 'ValueChangedCallback', jListCallback_bak);
    % Update panel "Clusters" fields
    UpdateProperties();
    if isUpdateMouseSel
        % Update selected channels for all figures
        UpdateChannelSelection();
    end
end

%% ===== GET CLUSTER DISPLAY TYPE =====
% ClustersOptions:
%    |- overlayClusters   : {0, 1}
%    |- overlayConditions : {0, 1}
function ClustersOptions = GetClusterOptions() %#ok<DEFNU>
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Cluster');
    if isempty(ctrl)
        ClustersOptions = [];
        return;
    end
    % Get current clusters
    sClusters = GetClusters();
    % If at least one of the cluster functions is "All", ignore the overlay checkboxes
    if any(strcmpi({sClusters.Function}, 'All'))
        ClustersOptions.overlayClusters   = 0;
        ClustersOptions.overlayConditions = 0;
    else
        % Overlay
        ClustersOptions.overlayClusters   = ctrl.jCheckOverlayClusters.isSelected();
        ClustersOptions.overlayConditions = ctrl.jCheckOverlayConditions.isSelected();
    end
end

%% ===== SET CLUSTER DISPLAY TYPE =====
function SetClusterOptions(overlayClusters, overlayConditions) %#ok<DEFNU>
    % Overlay
    if ~isempty(overlayClusters)
        ctrl.jCheckOverlayClusters.setSelected(overlayClusters);
    end
    if ~isempty(overlayConditions)
        ctrl.jCheckOverlayConditions.setSelected(overlayConditions);
    end
end

%% ===== SET CLUSTER FUNCTION =====
% USAGE:  SetClusterFunction(Function)
function SetClusterFunction(Function)
    % Get select clusters
    [sClusters, iClusters] = GetSelectedClusters();
    if isempty(iClusters)
        return
    end
    % Set function
    [sClusters.Function] = deal(Function);
    % Save clusters
    SetClusters(iClusters, sClusters);
    % Update JList
    UpdateClustersList();
    % Select edited clusters (selection was lost during update)
    SetSelectedClusters(iClusters);
end

%% ===== GET CHANNELS IN CLUSTER =====
function [iChannel, Modality] = GetChannelsInCluster(sCluster, Channel, ChannelFlag)
    % Get channels
    iChannel = zeros(1, length(sCluster.Sensors));
    for ic = 1:length(iChannel)
        tmp = find(strcmpi(sCluster.Sensors{ic}, {Channel.Name}));
        if isempty(tmp)
            iChannel = [];
            bst_error(['Sensor "' sCluster.Sensors{ic} '" does not exist.'], 'Clusters', 0);
            return;
        end
        iChannel(ic) = tmp;
    end
    % Check cluster homogeneity
    Modality = unique({Channel(iChannel).Type});
    if (length(Modality) > 1)
        iChannel = [];
        bst_error('Sensors in a cluster must be of the same type.', 'Clusters', 0);
        return;
    end
    Modality = Modality{1};
    % Remove channels that are "BAD"
    iBadChannels = find(ChannelFlag == -1);
    iChannel = setdiff(iChannel, iBadChannels);
    % Check if there are still some channels in the cluster
    if isempty(iChannel)
        iChannel = [];
        bst_error('All the channels in the cluster are BAD.', 'Clusters', 0);
        return;
    end
end

%% ===== CREATE NEW CLUSTER =====
% Usage:  [sCluster, iCluster] = CreateNewCluster(Sensors)      : Cell array of sensors names
%         [sCluster, iCluster] = CreateNewCluster('Selection')  : Get the current sensors selected by the user with his mouse
%         [sCluster, iCluster] = CreateNewCluster('Indices')    : Ask user to type the list of sensors indices
function [sCluster, iCluster] = CreateNewCluster(Sensors)
    global GlobalData;
    % Update current figure
    bst_figures('CheckCurrentFigure');

    % === GET SENSORS NAMES ===
    if ischar(Sensors)
        % Get selected figure
        [hFig,iFig,iDS] = bst_figures('GetCurrentFigure');
        if isempty(hFig)
            bst_error('No figure available.', 'Create new cluster', 0);
            return
        end

        % SELECTION: Get mouse selected sensors
        if strcmpi(Sensors, 'Selection')
            % Get the mouse selected sensors (generic interface)
            [SelChan, iSelChan] = bst_figures('GetSelectedChannels', iDS);
            % Selected sensors = intersection(mouse selection, figure sensors)
            iSensors = intersect(iSelChan, GlobalData.DataSet(iDS).Figure(iFig).SelectedChannels);
            % If no selection in current figure: error
            if isempty(iSensors)
                bst_error('Please select at least one sensor in the current figure.', 'Create new cluster', 0);
                return;
            end
        % INDICES: Ask user to give a list of sensors indices
        elseif strcmpi(Sensors, 'Indices')
            % Ask user to enter manually the indices of the sensors
            res = java_dialog('input', ['Enter the list of channels separated with commas.' 10 10 ...
                                        'Channels can be selected by:' 10 ...
                                        '     1) Types: "MEG, Misc"' 10 ...
                                        '     2) Names: "EEG021, EEG023"' 10 ...
                                        '     3) Indices: "5, 12, 13"' 10 10], 'Create cluster');
            if isempty(res)
                return;
            end
            % Get the channels indices
            iSensors = channel_find(GlobalData.DataSet(iDS).Channel, res);
            % If invalid sensors list
            if isempty(iSensors)
                bst_error('Invalid sensors selection.', 'Create new cluster', 0);
                return;
            end
            % Order sensors names and make unique
            iSensors = unique(iSensors);
        end
        % Check that the sensors are all of the same type
        uniqueType = unique({GlobalData.DataSet(iDS).Channel(iSensors).Type});
        if (length(uniqueType) > 1)
            java_dialog('warning', ['Warning: Sensors in a cluster must be of the same type.' 10 10 ...
                                    'Different types in this selection: ' uniqueType{1}, ', ' uniqueType{2}, 10 10], 'Create cluster');
            return
        end
        % Get the names of the mouse selected sensors
        Sensors = {GlobalData.DataSet(iDS).Channel(iSensors).Name};
    end

    % === NEW CLUSTER ===
    % New cluster structure
    sCluster  = db_template('Cluster');
    % Store current cluster sensors list
    sCluster.Sensors = Sensors;
    % Add clusters
    iCluster = SetClusters('Add', sCluster);

    % Update clusters list
    UpdateClustersList();
    % Select new cluster
    SetSelectedClusters(iCluster);
end

%% ===============================================================================
%  ====== CLUSTERS OPERATIONS ====================================================
%  ===============================================================================
%% ===== REMOVE CLUSTERS =====
% Usage : RemoveClusters()  : remove the clusters selected in the JList
function RemoveClusters()
    global GlobalData;
    % Get dataset
    [sClusters, iDSall, iFigall] = GetClusters();
    if isempty(iDSall)
        return;
    end
    % Get selected clusters
    [sClusters, iClusters] = GetSelectedClusters();
    % Check whether a cluster is selected
    if isempty(sClusters)
        java_dialog('warning', 'No cluster selected.', 'Remove cluster');
        return
    end
    % Ask for confirmation
    if (length(sClusters) == 1)
        strConfirm = ['Delete cluster "' sClusters(1).Label '"?'];
    else
        strConfirm = ['Delete ' num2str(length(sClusters)) ' clusters?'];
    end
    if ~java_dialog('confirm', strConfirm)
        return;
    end
    % Remove clusters definitions from global data structure
    for iDS = unique(iDSall)
        GlobalData.DataSet(iDS).Clusters(iClusters) = [];
        GlobalData.DataSet(iDS).isChannelModified = 1;
    end
    % Update Clusters list
    UpdateClustersList();
    % Reset list of selected sensors
    bst_figures('SetSelectedRows', [], 0);
end

%% ===== EDIT CLUSTER LABEL =====
% Usage : EditClusterLabel()  : Interactive edition of cluster name
function EditClusterLabel()
    % Get all clusters
    sClustersAll = GetClusters();
    % Get selected clusters
    [sCluster, iCluster] = GetSelectedClusters();
    % Warning message if no cluster selected
    if isempty(sCluster)
        java_dialog('warning', 'No cluster selected.', 'Rename selected cluster');
        return;
    % If more than one cluster selected: keep only the first one
    elseif (length(sCluster) > 1)
        iCluster = iCluster(1);
        sCluster = sCluster(1);
        SetSelectedClusters(iCluster);
    end
    % Ask user for a new Cluster Label
    newLabel = java_dialog('input', sprintf('Please enter a new label for cluster "%s":', sCluster.Label), ...
                           'Rename selected cluster', [], sCluster.Label);
    if isempty(newLabel) || strcmpi(newLabel, sCluster.Label)
        return
    elseif ~isempty(sClustersAll) && any(strcmpi({sClustersAll.Label}, newLabel))
        java_dialog('warning', 'Cluster name already exists.', 'Rename selected cluster');
        return;
    end

    % Update cluster definition
    sCluster.Label = newLabel;
    SetClusters(iCluster, sCluster);
    % Update JList
    UpdateClustersList();
    % Select edited clusters (selection was lost during update)
    SetSelectedClusters(iCluster);
end

%% ===== EDIT CLUSTER COLOR =====
function EditClusterColor(newColor)
    % Get selected clusters
    [sSelClusters, iSelClusters] = GetSelectedClusters();
    if isempty(iSelClusters)
        java_dialog('warning', 'No cluster selected.', 'Edit cluster color');
        return
    end
    % If color is not specified in argument : ask it to user
    if (nargin < 1)
        newColor = java_dialog('color');
        if (length(newColor) ~= 3) || all(sSelClusters(1).Color(:) == newColor(:))
            return
        end
    end
    % Update cluster color
    for i = 1:length(sSelClusters)
        sSelClusters(i).Color = newColor;
    end
    % Save clusters
    SetClusters(iSelClusters, sSelClusters);
    % Update clusters list
    UpdateClustersList();
end

%% ===== SAVE CLUSTERS =====
function SaveClusters(varargin)
    % Get selected clusters
    sClusters = GetSelectedClusters();
    if isempty(sClusters)
        return
    end
    % Default folder name
    ClusterDir = GetDefaultClusterDir();
    % Build a default file name
    ClusterFile = bst_fullfile(ClusterDir, ['cluster', sprintf('_%s', sClusters.Label), '.mat']);
    % Get filename where to store the filename
    ClusterFile = java_getfile( 'save', 'Save selected clusters', ClusterFile, ...
                                'single', 'files', ...
                                {{'_cluster'}, 'Brainstorm sensor clusters (*cluster*.mat)', 'BST'}, 1);
    if isempty(ClusterFile)
        return;
    end
    % Make sure that filename contains the 'cluster' tag
    if isempty(strfind(ClusterFile, '_cluster')) && isempty(strfind(ClusterFile, 'cluster_'))
        [filePath, fileBase, fileExt] = bst_fileparts(ClusterFile);
        ClusterFile = bst_fullfile(filePath, ['cluster_' fileBase fileExt]);
    end
    % Save file
    FileMat.Clusters = sClusters;
    bst_save(ClusterFile, FileMat, 'v7');
end

%% ===== GET DEFAULT CLUSTER DIR =====
function ClusterDir = GetDefaultClusterDir()
    global GlobalData;
    % Get subject of current figure
    [hFig,iFig,iDS] = bst_figures('GetCurrentFigure');
    if ~isempty(iDS) && ~isempty(GlobalData.DataSet(iDS).SubjectFile)
        SubjectFile = GlobalData.DataSet(iDS).SubjectFile;
    % Get current subject in the database
    else
        sSubject = bst_get('Subject');
        if ~isempty(sSubject)
            SubjectFile = sSubject.FileName;
        else
            SubjectFile = [];
        end
    end
    % Get subject anatomy folder
    if ~isempty(SubjectFile)
        ClusterDir = bst_fileparts(file_fullpath(SubjectFile));
    else
        ClusterDir = '';
    end
end

%% ===== LOAD CLUSTER =====
function LoadClusters(varargin)
    % Get default cluster folder
    ClusterDir = GetDefaultClusterDir();
    % Ask user which are the files to be loaded
    [ClusterFiles, FileFormat] = java_getfile(...
        'open', 'Import clusters', ClusterDir, ...
        'multiple', 'files', ...
        bst_get('FileFilters', 'clusterin'), 1);
    if isempty(ClusterFiles)
        return
    end

    % ==== CREATE AND DISPLAY ====
    iNewClusters = [];
    bst_progress('start', 'Load clusters', 'Load cluster file');
    % Load all files selected by user
    for iFile = 1:length(ClusterFiles)
        % Load clusters
        sClusters = in_clusters(ClusterFiles{iFile}, FileFormat);
        % Add to current clusters
        iNewClusters = [iNewClusters, SetClusters('Add', sClusters)];
    end
    % Update cluster list
    UpdateClustersList();
    % Select first cluster
    if isempty(iNewClusters)
        SetSelectedClusters(iNewClusters(1));
    end
    bst_progress('stop');
end

%% ===== COPY CLUSTERS =====
function CopyClusters(Target)
    global GlobalData;
    % Get loaded clusters
    [sClusters, iDSall] = GetSelectedClusters();
    if isempty(sClusters)
        bst_error('No clusters selected.', 'Copy clusters', 0);
        return
    end
    % Copy clusters
    db_set_clusters(GlobalData.DataSet(iDSall(1)).ChannelFile, Target, sClusters);
end

%% ===== EXPORT CLUSTERS TO MATLAB =====
function ExportClustersToMatlab()
    % Get selected clusters
    sClusters = GetSelectedClusters();
    % If nothing selected, take all clusters
    if isempty(sClusters)
        sClusters = GetClusters();
    end
    % If nothing: exit
    if isempty(sClusters)
        return;
    end
    % Export to the base workspace
    export_matlab(sClusters, []);
    % Display in the command window the selected clusters
    disp([10 'List of sensors for each cluster:']);
    for i = 1:length(sClusters)
        disp(['   ' sClusters(i).Label ': ' sprintf('%s ', sClusters(i).Sensors{:})]);
    end
    disp(' ');
end

%% ===== IMPORT CLUSTERS FROM MATLAB =====
function ImportClustersFromMatlab()
    % Export to the base workspace
    sClusters = in_matlab_var([], 'struct');
    if isempty(sClusters)
        return;
    end
    % Check structure
    sTemplate = db_template('cluster');
    if isempty(sClusters) || ~isequal(fieldnames(sClusters), fieldnames(sTemplate))
        bst_error('Invalid clusters structure.', 'Import from Matlab', 0);
        return;
    end
    % Save new cluster
    iNewCluster = SetClusters('Add', sClusters);
    % Update "Clusters Manager" panel
    UpdateClustersList();
    % Select last cluster in list (new cluster)
    SetSelectedClusters(iNewCluster);
end

author: brainstorm-tools


text:
disp('make some variables...')
x = 2.3
y = rand(4, 3, 2)
z = 'durp durp durp'

disp('import the matlab module')
pymex('import matlab')

disp('"pull" the variables from matlab to python')
pymex('x, y, z = matlab.pull("x", "y", "z")')

disp('print the variables as seen by python')
pymex('matlab.mex_print(x, y, z)')

disp('"push" the variables from python to matlab')
pymex('matlab.push("x2", x, "y2", y, "z2", z)')

disp('test that y is equal to y2')
y == y2

author: Sable

repo: mcbench-benchmarks


text:
% run Demo_Get_PCA_matrix.m first

ksize = 15;

load PCA_P.mat;

% show the PCA basis

for i = 1:size(P,1)

    kernel = reshape(P(i,:),ksize,ksize);

    subplot 121
    imagesc(kernel);
    title([int2str(i)])
    axis square;

    subplot 122
    surf(kernel);
    title([int2str(i)])
    view(45,55)
    xlim([1 15]);
    ylim([1 15]);
    axis square;
    pause(2)
end

author: cszn

repo: SRMD

save_path: github-repos/MATLAB/cszn-SRMD


text:
function [EEG, choppedFront, choppedBack] = chop(EEG, frontChop, backChop)
% Deletes all but specified seconds from outside first and last event
%
% This removes bad segments from the beginning and end of the data
% set, which often correspond to bad times for continuous EEG
%
% Parameters:
%    EEG = EEGLAB structure
%    frontChop =
choppedFront = 0;
choppedBack = 0;
if size(EEG.data, 3) > 1
    warning('chop:DataNotContinuous', 'Only works on continuous data');
    return;
elseif nargin < 3
    error('chop:NotEnoughArguments', 'Requires 3 arguments');
elseif isempty(EEG.event)
    warning('chop:NoEvents', 'EEG has no events, so no chop is done');
    return;
end
srate = EEG.srate;
firstEventTime = double(EEG.event(1).latency - 1)/srate;
lastEventTime = double(EEG.event(end).latency - 1)/srate;
lastDataTime = double(EEG.pnts-1)/srate;
choppedFront = max(0, firstEventTime - frontChop);
endTime = min(lastDataTime, lastEventTime + backChop);
choppedBack = lastDataTime - endTime;
EEG = pop_select(EEG, 'time', [choppedFront, endTime]);
boundaryEvents = strcmpi({EEG.event.type}, 'boundary');
if sum(strcmpi(EEG.event(1).type, 'boundary')) > 0
    fprintf('Removing %g boundary events\n', sum(boundaryEvents));
    EEG.event(boundaryEvents) = [];
end

author: VisLab

repo: EEG-Clean-Tools


text:
function Vf = fillMissingGrid(xGridV, yGridV, zGridV, Vf)
% function matrixOut = fillMissingGrid(xGridV, yGridV, zGridV, Vf)
%
% This function takes in a matrix that contains NaN values and replaces
% them with values that are linearly interpolated from the surrounding.
%
% APA, 10/03/2012

% Find NaN locations
indNaN = isnan(Vf(:,:,:,1));

% Find Not NaN locations
indNotNaN = ~isnan(Vf(:,:,:,1));

% Create x,y,z 3D grid
[xM,yM,zM] = meshgrid(xGridV, yGridV, zGridV);

% Get x,y,z of NaN points
xNaNv = xM(indNaN);
yNaNv = yM(indNaN);
zNaNv = zM(indNaN);

% Get x,y,z of Not-NaN points
xNotNaNv = xM(indNotNaN);
yNotNaNv = yM(indNotNaN);
zNotNaNv = zM(indNotNaN);

% Get values at Not-NaN points
xVf = Vf(:,:,:,1);
xVfv = xVf(indNotNaN);
yVf = Vf(:,:,:,2);
yVfv = yVf(indNotNaN);
zVf = Vf(:,:,:,3);
zVfv = zVf(indNotNaN);

Fx = TriScatteredInterp([xNotNaNv yNotNaNv zNotNaNv],double(xVfv));
Fy = TriScatteredInterp([xNotNaNv yNotNaNv zNotNaNv],double(yVfv));
Fz = TriScatteredInterp([xNotNaNv yNotNaNv zNotNaNv],double(zVfv));

xVfNaN = Fx([xNaNv yNaNv zNaNv]);
yVfNaN = Fy([xNaNv yNaNv zNaNv]);
zVfNaN = Fz([xNaNv yNaNv zNaNv]);

% xVfNaN = zeros(1,length(xNaNv));
% yVfNaN = xVfNaN;
% zVfNaN = xVfNaN;
% numVoxels = length(xNaNv);
%
% for i=1:numVoxels
%     disp([num2str(i),' / ', num2str(numVoxels)])
%     %     % Find distance of NaN from all the Not-NaN points
%     %     distV = (xNaNv(i) - xNotNaNv).^2 + (yNaNv(i) - yNotNaNv).^2 + (zNaNv(i) - zNotNaNv).^2;
%     %     % Sort by distance
%     %     [~,indV] = sort(distV);
%
%     %     indKrig = indV(1:20);
%     indKrig = [];
%     % Search in neighborhood
%     searchRadius = 0.2;
%     while sum(indKrig>0) < 5
%         indKrig = (xNaNv(i) - xNotNaNv).^2 + (yNaNv(i) - yNotNaNv).^2 + (zNaNv(i) - zNotNaNv).^2 <= searchRadius^2;
%         searchRadius = searchRadius + 0.1;
%     end
%     % Kriging
%     dmodelX = dacefit([xNotNaNv(indKrig) yNotNaNv(indKrig) zNotNaNv(indKrig)], double(xVfv(indKrig)), @regpoly0, @correxp, 10, 1e-1, 20);
%     xVfKrig = predictor([xNaNv(i) yNaNv(i) zNaNv(i)], dmodelX);
%     dmodelY = dacefit([xNotNaNv(indKrig) yNotNaNv(indKrig) zNotNaNv(indKrig)], double(yVfv(indKrig)), @regpoly0, @correxp, 10, 1e-1, 20);
%     yVfKrig = predictor([xNaNv(i) yNaNv(i) zNaNv(i)], dmodelY);
%     dmodelZ = dacefit([xNotNaNv(indKrig) yNotNaNv(indKrig) zNotNaNv(indKrig)], double(zVfv(indKrig)), @regpoly0, @correxp, 10, 1e-1, 20);
%     zVfKrig = predictor([xNaNv(i) yNaNv(i) zNaNv(i)], dmodelZ);
%
%     xNotNaNv = [xNotNaNv; xNaNv(i)];
%     yNotNaNv = [yNotNaNv; yNaNv(i)];
%     zNotNaNv = [zNotNaNv; zNaNv(i)];
%     xVfv = [xVfv; xVfKrig];
%     yVfv = [yVfv; yVfKrig];
%     zVfv = [zVfv; zVfKrig];
%
%     xVfNaN(i) = xVfKrig;
%     yVfNaN(i) = yVfKrig;
%     zVfNaN(i) = zVfKrig;
%
% end

xVf(indNaN) = xVfNaN;
yVf(indNaN) = yVfNaN;
zVf(indNaN) = zVfNaN;
Vf(:,:,:,1) = xVf;
Vf(:,:,:,2) = yVf;
Vf(:,:,:,3) = zVf;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [s, lambda_k] = tr_subsolver(problem, w, grad, tr_radius, sub_hess_indices, ...
                                                successful_flag, lambda_k, subproblem_solver,...
                                                exact_tol, krylov_tol)
    if strcmp(subproblem_solver, 'cauchy_point')
        %Hg = Hv_f(w, new_X, new_Y, grad,alpha);
        Hg = problem.hess_vec(w, grad, sub_hess_indices);
        gBg = grad'*Hg;
        tau = 1;
        if gBg > 0  % if model is convex quadratic the unconstrained minimizer may be inside the TR
            tau = min([(norm(grad))^3 / (tr_radius * gBg), 1]);
        end
        pc = - tau * tr_radius * (grad./norm(grad));
        s = pc;
        lamda_k = 0;
    elseif strcmp(subproblem_solver, 'dog_leg')
        %H = hessian_f(w, new_X, new_Y, alpha);
        H = problem.hess(w, sub_hess_indices);
        gBg = grad'*(H'*grad);
        if gBg <= 0
            error('dog_leg requires H to be positive definite in all steps!') ;
        end

        % Compute the Newton Point and return it if inside the TR
        L = chol(H,'lower');
        y = L\grad;
        pn = -L'\y;
%         cholesky_B = chol(H,'lower');
%         pn = - (cholesky_B\grad);
        if (norm(pn) < tr_radius)
            s = pn;
            lambda_k = 0;
            return;
        end

        % Compute the 'unconstrained Cauchy Point'
        pc = -((grad'*grad)/gBg) * grad;
        pc_norm = norm(pc);

        % if it is outside the TR, return the point where the path intersects the boundary
        if pc_norm >= tr_radius
            p_boundary = pc * (tr_radius / pc_norm);
            s = p_boundary;
            lambda_k = 0;
            return;
        end

        % else, give intersection of path from pc to pn with tr_radius.
        [t_lower, t_upper] = solve_quadratic_equation(pc, pn, tr_radius);
        p_boundary = pc + t_upper * (pn - pc);
        s = p_boundary;
        lambda_k = 0;

    elseif strcmp(subproblem_solver, 'cg')
        grad_norm = norm(grad);
        p_start = zeros(size(grad,1),size(grad,2));

        if grad_norm < min([sqrt(norm(grad)) * norm(grad), krylov_tol])
            s = p_start;
            lambda_k = 0;
            return;
        end

        % initialize
        z = p_start;
        r = grad;
        d = -r;
        k = 0;

        while true
            %Bd = Hv_f(w, new_X, new_Y, d, alpha);
            Bd = problem.hess_vec(w, d, sub_hess_indices);
            dBd = d'*Bd;
            % terminate when encountering a direction of negative curvature with lowest boundary point along current search direction
            if dBd <= 0
                [t_lower, t_upper] = solve_quadratic_equation(z, d, tr_radius);
                p_low = z + t_lower * d;
                p_up = z + t_upper * d;
                %m_p_low = loss_f(w + p_low, X, Y, alpha) + grad'*p_low + 0.5 * p_low'*(H'*p_low);
                %m_p_up = loss_f(w + p_up, X, Y, alpha) + grad'*p_up + 0.5 * p_up'*(H'*p_up);

                H = problem.hess(w, sub_hess_indices); % Added by HK (Need to be checked.)
                m_p_low = problem.cost(w+p_low) + grad'*p_low + 0.5 * p_low'*(H'*p_low);
                m_p_up = problem.cost(w+p_up) + grad'*p_up + 0.5 * p_up'*(H'*p_up);
                problem.cost(w);
                if m_p_low < m_p_up
                    s = p_low;
                    lambda_k = 0;
                    return;
                else
                    s = p_up;
                    lambda_k = 0;
                    return;
                end
            end

            alpha = (r'*r) / dBd;
            z_next = z + alpha * d;
            % terminate if z_next violates TR bound
            if norm(z_next) >= tr_radius
                % return intersect of current search direction w/ boud
                [t_lower, t_upper] = solve_quadratic_equation(z, d, tr_radius);
                s = z + t_upper * d;
                lambda_k = 0;
                return;
            end

            r_next = r + alpha * Bd;
            if norm(r_next) < min([sqrt(norm(grad)) * norm(grad),krylov_tol])
                s = z_next;
                lambda_k = 0;
                return;
            end

            beta_next = (r_next'*r_next) / (r'*r);
            d_next = -r_next + beta_next * d;
            % update iterates
            z = z_next;
            r = r_next;
            d = d_next;
            k = k + 1;
        end
    elseif strcmp(subproblem_solver, 'GLTR')
        g_norm = norm(grad);
        s = zeros(size(grad,1), size(grad,2));

        if g_norm == 0
            % escape along the direction of the leftmost eigenvector as far as tr_radius permits
            fprintf('zero gradient encountered');
            %H = hessian_f(w, new_X, new_Y, alpha);
            H = problem.hess(w, sub_hess_indices);
            [s, lambda_k] = exact_TR_suproblem_solver(grad, H, tr_radius, exact_tol, successful_flag, lambda_k);
        else
            % initialize
            g = grad;
            p = -g;
            gamma = g_norm;
            T = zeros(1, 1);
            alpha_k = [];
            beta_k = [];
            interior_flag = true;
            k = 0;

            while true
                %Hp = Hv_f(w, new_X, new_Y, p, alpha);
                Hp = problem.hess_vec(w, p, sub_hess_indices);
                pHp = p'*Hp;
                alpha = g'*g / pHp;

                alpha_k = [alpha_k; alpha];

                %Lanczos Step 1: Build up subspace
                % a) Create g_lanczos = gamma*e_1

                e_1 = zeros(k + 1, 1);
                e_1(1) = 1.0;
                g_lanczos = gamma .* e_1;

                % b) Create T for Lanczos Model
                T_new = zeros(k + 1, k + 1);
                if k == 0
                    T(k+1, k+1) = 1.0/alpha;
                    T_new(1:k+1, 1:k+1) = T;
                else
                    T_new(1:size(T,1), 1:size(T,2)) = T;
                    T_new(k+1, k+1) = 1.0 / alpha + beta_/ alpha_k(k);
                    T_new(k, k+1) = sqrt(beta_) / abs(alpha_k(k));
                    T_new(k+1, k) = sqrt(beta_) / abs(alpha_k(k));
                    T = T_new;
                end

                if (interior_flag == true && alpha < 0) || norm(s + alpha * p) >= tr_radius
                    interior_flag = false;
                end

                if interior_flag == true
                    s = s + alpha * p;
                else
                    % Lanczos Step 2: solve problem in subspace
                    [h, lambda_k] = exact_TR_suproblem_solver(g_lanczos, T, tr_radius, exact_tol, successful_flag, lambda_k);
                end
                g_next = g + alpha * Hp;

                % test for convergence
                e_k = zeros(k + 1,1);
                e_k(k+1) = 1.0;

                if interior_flag == true && norm(g_next) < min([sqrt(norm(grad)) * norm(grad), krylov_tol])
                    break;
                end
                if interior_flag == false && norm(g_next) * abs(h'*e_k) < min([sqrt(norm(grad)) * norm(grad), krylov_tol])
                    break;
                end

                if k == problem.d
                    fprintf('Krylov dimensionality reach full space! Breaking out..\n');
                    break;
                end
                beta_ = (g_next'*g_next) / (g'*g);
                beta_k = [beta_k; beta_];
                p = -g_next + beta_ * p;
                g = g_next;
                k = k + 1;
            end

            if interior_flag == false
                % Recover Q by building up the lanczos space, TBD: keep storable Qs in memory
                n = size(grad,1);
                Q1 = zeros(n, k + 1);

                g = grad;
                p = -g;

                for j = 0 : k
                    gn = norm(g);
                    if j == 0
                        sigma = 1;
                    else
                        sigma = -sign(alpha_k(j)) * sigma;
                    end
                    Q1(:, j+1) = sigma * g / gn;

                    if ~ (j == k)
                        %Hp = Hv_f(w, new_X, new_Y, p, alpha);
                        Hp = problem.hess_vec(w, p, sub_hess_indices);
                        g = g + alpha_k(j+1) * Hp;
                        p = -g + beta_k(j+1) * p;
                    end
                end

                % compute final step in R^n
                s = Q1*h;
            end
        end
    elseif strcmp(subproblem_solver, 'exact')
        %H = hessian_f(w, new_X, new_Y, alpha);
        Hw = problem.hess(w, sub_hess_indices);
        [s, lambda_k] = exact_TR_suproblem_solver(grad, H, tr_radius, exact_tol, successful_flag, lambda_k);
    else
    	error('solver unknown\n');
    end
end

function [s, lambda_j] = exact_TR_suproblem_solver(grad, H, tr_radius, exact_tol, successful_flag, lambda_k)

    s = zeros(size(grad,1),size(grad,2));
    % Step 0: initialize safeguards
    H_ii_min = min(diag(H));
    absH = abs(H);
    H_max_norm = sqrt((size(H,1))^2) * max(absH(:));
    H_fro_norm = norm(H, 'fro');
    list_l = [];
    list_u = [];
    for i = 1 : length(H)
        l = H(i, i) + sum(abs(H(i, :))) - abs(H(i, i));
        u = -H(i, i) + sum(abs(H(i, :))) - abs(H(i, i));
        list_l = [list_l l];
        list_u = [list_u u];
    end
    gerschgorin_l = max(list_l);
    gerschgorin_u = max(list_u);

    lambda_lower = max([0, -H_ii_min, norm(grad) / tr_radius - min([H_fro_norm, H_max_norm, gerschgorin_l])]);
    lambda_upper = max([0, norm(grad) / tr_radius + min([H_fro_norm, H_max_norm, gerschgorin_u])]);

    if successful_flag == false && (lambda_lower <= lambda_k) && (lambda_k <= lambda_upper) % reinitialize at previous lambda in case of unscuccesful iterations
        lambda_j = lambda_k;
    elseif lambda_lower == 0  % allow for fast convergence in case of inner solution
        lambda_j = lambda_lower;
    else
        lambda_j = lambda_lower + (lambda_upper-lambda_lower)*rand(1,1);
    end

    i = 0;
    % Root Finding
    while true
        i = i + 1;
        lambda_in_N = false;
        lambda_plus_in_N = false;
        B = H + lambda_j * eye(size(H,1), size(H,2));
        try
            % 1 Factorize B
            L = chol(B, 'lower');
            % 2 Solve LL^Ts=-g
            Li = inv(L);
            s = - (Li'*Li)*grad;
            sn = norm(s);
            % 2.1 Termination: Lambda in F, if q(s(lamda))<eps_opt q(s*) and sn<eps_tr tr_radius -> stop. By Conn: Lemma 7.3.5:
            phi_lambda = 1.0 / sn - 1.0 / tr_radius;
            %if (abs(sn - tr_radius) <= exact_tol * tr_radius):
            if (abs(phi_lambda)<=exact_tol) %
                break;
            end

            % 3 Solve Lw=s
                w = Li*s;
                wn = norm(w);

            % Step 1: Lambda in L
            if lambda_j > 0 && (phi_lambda) < 0
                % print ('lambda: ',lambda_j, ' in L')
                lambda_plus = lambda_j + ((sn - tr_radius) / tr_radius) * ((sn^2) / (wn^2));
                lambda_j = lambda_plus;

            % Step 2: Lambda in G    (sn<tr_radius)
            elseif (phi_lambda) > 0 && lambda_j > 0 && any(grad ~= 0) % TBD: remove grad
                % print ('lambda: ',lambda_j, ' in G')
                lambda_upper = lambda_j;
                lambda_plus = lambda_j + ((sn - tr_radius) / tr_radius) * ((sn^2) / (wn^2));

                % Step 2a: If factorization succeeds: lambda_plus in L
                if lambda_plus > 0
                    try
                        % 1 Factorize B
                        B_plus = H + lambda_plus * eye(size(H,1), size(H,2));
                        L = chol(B_plus, 'lower');
                        lambda_j = lambda_plus;
                        % print ('lambda+', lambda_plus, 'in L')
                    catch
                        lambda_plus_in_N = true;
                    end
                end

                % Step 2b/c: If not: Lambda_plus in N
                if lambda_plus <= 0 || lambda_plus_in_N == true
                    % 1. Check for interior convergence (H pd, phi(lambda)>=0, lambda_l=0)
                    try
                        U = chol(H, 'upper');
                        H_pd = true;
                    catch
                        H_pd = false;
                    end

                    if lambda_lower == 0 && H_pd == true && phi_lambda >= 0 %cannot happen in ARC!
                        lambda_j = 0;
                        % print ('inner solution found');
                        break;
                    % 2. Else, choose a lambda within the safeguard interval
                    else
                        % print ('lambda_plus', lambda_plus, 'in N')
                        lambda_lower = max([lambda_lower, lambda_plus]);  % reset lower safeguard
                        lambda_j = max([sqrt(lambda_lower * lambda_upper),lambda_lower + 0.01 * (lambda_upper - lambda_lower)]);
                        lambda_upper = single(lambda_upper);

                        if lambda_lower == lambda_upper
                            lambda_j = lambda_lower;
                            % Hard case
                            [ev, ew] = eig(H);
                            d = ev(:, 1);
                            dn = norm(d);
                            assert((ew == -lambda_j), 'Ackward: in hard case but lambda_j != -lambda_1');
                            [tao_lower, tao_upper] = mitternachtsformel(1, 2*(s'*d), (s'*s)-tr_radius^2);
                            s = s + tao_lower * d;
                            fprintf('hard case resolved inside');
                        end
                    end
                end

            elseif (phi_lambda) == 0
                break;
            else % TBD:  move into if lambda+ column #this only happens for Hg=0 -> s=(0,..,0)->phi=inf -> lambda_plus=nan -> hard case (e.g. at saddle)
                lambda_in_N = true;
            end
        % Step 3: Lambda in N
        catch
            lambda_in_N = true;
        end

        if lambda_in_N == true
            % print ('lambda: ',lambda_j, ' in N')
            try
                U = chol(H, 'upper');
                H_pd = true;
            catch
                H_pd = false;
            end

            % 1. Check for interior convergence (H pd, phi(lambda)>=0, lambda_l=0)
            if lambda_lower == 0 && H_pd == true && phi_lambda >= 0
                lambda_j = 0;
                %print ('inner solution found')
                break;
            % 2. Else, choose a lambda within the safeguard interval
            else
                lambda_lower = max([lambda_lower, lambda_j]);  % reset lower safeguard
                lambda_j = max([sqrt(lambda_lower * lambda_upper),lambda_lower + 0.01 * (lambda_upper - lambda_lower)]);  % eq 7.3.14
                lambda_upper = single(lambda_upper);
                % Check for Hard Case:
                if lambda_lower == lambda_upper
                    lambda_j = lambda_lower;
                    [ev, ew] = eig(H);
                    d = ev(:, 1);
                    dn = norm(d);
                    assert((ew == -lambda_j), 'Ackward: in hard case but lambda_j != -lambda_1');
                    [tao_lower, tao_upper] = mitternachtsformel(1, 2*(s'*d), (s'*s)-tr_radius^2);
                    s = s + tao_lower * d;

                    fprintf('hard case resolved outside');
                end
            end
        end
    end

    % compute final step
    B = H + lambda_j * eye(size(H,1), size(H,2));
    % 1 Factorize B
    L = chol(B, 'lower');
    % 2 Solve LL^Ts=-g
    Li = inv(L);
    s = - Li'*Li*grad;
    %print (i,' exact solver iterations')

end

% Auxiliary Functions
function [t_lower, t_upper] = mitternachtsformel(a, b, c)
    sqrt_discriminant = sqrt(b * b - 4 * a * c);
    t_lower = (-b - sqrt_discriminant) / (2 * a);
    t_upper = (-b + sqrt_discriminant) / (2 * a);
end

function [t_lower, t_upper] = solve_quadratic_equation(pc, pn, tr_radius)
    % solves ax^2+bx+c=0
    a = (pn - pc)'*(pn - pc);
    b = 2 * (pc'*(pn - pc));
    c = (pc'*pc) - tr_radius^2;
    sqrt_discriminant = sqrt(b * b - 4 * a * c);
    t_lower = (-b - sqrt_discriminant) / (2 * a);
    t_upper = (-b + sqrt_discriminant) / (2 * a);
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function output = ...
  ConvertToObservationProbatility(f0_structure, frame_shift)
% Convert wavelet filter bank output to ovserbation probability map
%
%   output = ConvertToObservationProbatility(f0_structure, frame_shift)
%
% Arguments
%
% f0_structure: structure with the following fields.
%     Note that only used fieds are listed below.
%   temporal_positions: vector: center location of the analysis window.
%     represented in terms of second
%   sampling_frequency: scalar: sampling frequency (Hz)
%   center_frequency_list: vector: center frequency of each wavelet filter
%   (Hz)
%   inst_freq_map: matrix: instantaneous frequency of each filter output
%   (Hz)
%   residual_map: matrix: relative level of random component (power ratio)
%   amplitude_map: matrix: amplitude of filter output (absolute value)
% frame_shift: optional variable, scalar, real value in second
%   Interval between adjacent frame center positions
%
% Return value
%
% output: structure with the following fields.
%   temporal_positions: vector: center location of the analysis window.
%   frame_shift: scalar: interval between frame locations
%   inst_freq_map: matrix: instantaneous frequency of each channel (Hz)
%   residual_map: relative resitual of each channel
%           (relative power, calibrated for Nuttall-11 window.)
%   observation_probability_map: matrix: probability time series of each
%           channel
%   amplitude_map: matrix: amplitude of filter output (absolute value)
%   channels_in_octave: scalar: number of wavelet channels in one octave
%   center_frequency_list: vector: list of center frequencies of wavelet
%           filters (Hz)

% Copyright 2016 Google Inc. All Rights Reserved
% Author: hidekik@google.com (Hideki Kawahara)

if nargin == 1, frame_shift = 0.005; end;

startTic = tic;
calibration_constant = 64; % This is for NuttallWindows(~, 11)
time_axis = f0_structure.temporal_positions;
fs = f0_structure.sampling_frequency;
temporal_positions = (0:frame_shift:time_axis(end))';
frame_index = ...
  min(f0_structure.sample_length,max(1, round(fs * temporal_positions)));
number_of_channels = length(f0_structure.center_frequency_list);
inst_freq_map = f0_structure.inst_freq_map(:, frame_index);
residual_map = f0_structure.residual_map(:, frame_index);
amplitude_map = f0_structure.amplitude_map(:, frame_index);
tmp_variance_map = ...
  ConvertPowerToDecibel(residual_map) + calibration_constant;
snr_map = 10.0 .^ (tmp_variance_map / 20);
best_weight = snr_map * 0;
for ii = 1 : length(temporal_positions)
  best_weight(:, ii) = GetBestMixingWeight(snr_map(:, ii));
end;
observation_probability_map = tmp_variance_map * 0;
center_frequency_list = f0_structure.center_frequency_list;
tmp_half_bw = sqrt(center_frequency_list(2) / center_frequency_list(1));
fixed_best_weight = best_weight .* snr_map;
fc_in_cent_U = 1200 * log2(center_frequency_list * tmp_half_bw);
fc_in_cent_L = 1200 * log2(center_frequency_list / tmp_half_bw);
for ii = 1:length(temporal_positions)
  for kk = 1 : number_of_channels
    fqi_cent = 1200 * log2(max(20, inst_freq_map(kk, ii)));
    tmpU = erf((fc_in_cent_U - fqi_cent) / snr_map(kk, ii) / sqrt(2));
    tmpL = erf((fc_in_cent_L - fqi_cent) / snr_map(kk, ii) / sqrt(2));
    tmpProb = (tmpU - tmpL) / sum(tmpU - tmpL);
    observation_probability_map(:, ii) = ...
      observation_probability_map(:, ii) ...
      + fixed_best_weight(kk, ii) * tmpProb(:);
  end;
  observation_probability_map(:, ii) = observation_probability_map(:, ii) ...
    / sum(observation_probability_map(:, ii));
end;
output = struct('temporal_positions', temporal_positions, ...
                'frame_shift', frame_shift', ...
                'inst_freq_map', inst_freq_map, ...
                'residual_map', residual_map, ...
                'observation_probability_map', observation_probability_map, ...
                'amplitude_map', amplitude_map, ...
                'channels_in_octave', f0_structure.channels_in_octave', ...
                'center_frequency_list', center_frequency_list);
end

author: google


text:
function square_symq_gnuplot ( n, x, header )

%*****************************************************************************80
%
%% SQUARE_SYMQ_GNUPLOT: GNUPLOT plot of the symmetric square quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    11 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, real N, the number of nodes.
%
%    Input, real X(2,N), the coordinates of the nodes.
%
%    Input, string HEADER, a string to be used to identify
%    the files created.
%

%
%  Create the vertex file.
%
  vertex_filename = strcat ( header, '_vertices.txt' );
  vertex_unit = fopen ( vertex_filename, 'wt' );
  fprintf ( vertex_unit, '%g  %g\n', -1.0, -1.0 );
  fprintf ( vertex_unit, '%g  %g\n', +1.0, -1.0 );
  fprintf ( vertex_unit, '%g  %g\n', +1.0, +1.0 );
  fprintf ( vertex_unit, '%g  %g\n', -1.0, +1.0 );
  fprintf ( vertex_unit, '%g  %g\n', -1.0, -1.0 );
  fclose ( vertex_unit );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created vertex file "%s"\n', vertex_filename );
%
%  Create node file.
%
  node_filename = strcat ( header, '_nodes.txt' );
  node_unit = fopen ( node_filename, 'wt' );
  for j = 1 : n
    fprintf ( node_unit, '%g  %g\n', x(1:2,j) );
  end
  fclose ( node_unit );
  fprintf ( 1, '  Created node file "%s"\n', node_filename );
%
%  Create graphics command file.
%
  command_filename = strcat ( header, '_commands.txt' );
  command_unit = fopen ( command_filename, 'wt' );
  fprintf ( command_unit, '# %s\n', command_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, '# Usage:\n' );
  fprintf ( command_unit, '#  gnuplot < %s\n', command_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set term png\n' );
  plot_filename = strcat ( header, '.png' );
  fprintf ( command_unit, 'set output "%s"\n', plot_filename );
  fprintf ( command_unit, 'set xlabel "<--- X --->"\n' );
  fprintf ( command_unit, 'set ylabel "<--- Y --->"\n' );
  fprintf ( command_unit, 'set title "%s"\n', header );
  fprintf ( command_unit, 'set grid\n' );
  fprintf ( command_unit, 'set key off\n' );
  fprintf ( command_unit, 'set size ratio -1\n' );
  fprintf ( command_unit, 'set style data lines\n' );
  fprintf ( command_unit, 'set timestamp\n' );
  fprintf ( command_unit, 'plot "%s" with lines lw 3, \\\n', vertex_filename );
  fprintf ( command_unit, '     "%s" with points pt 7 lt 0\n', node_filename );
  fclose ( command_unit );

  fprintf ( 1, '  Created command file "%s"\n', command_filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
if 0
  % functionality test
  repmat([1 2; 3 4],2,4,2,2)
  repmat([1 2; 3 4],[2 4 2 2])
  repmat([1 2; 3 4],2)
  repmat(j,3,2)
  repmat('hello',3,2)
  repmat({7},3,2)
  repmat(sparse(7),3,2)
end

x = rand(300,1);
x = rand(10,1);
% run it once to load the definition
repmat(1,1,1);
niter = 100000/prod(size(x));
n = 100;
fprintf('repmat(x,1,n)\n');
tic; for i = 1:niter xrepmat(x,1,n); end; t0=toc;
fprintf('old repmat: %g\n',t0);
tic; for i = 1:niter repmat(x,1,n); end; t=toc;
fprintf('new repmat: %g (%g times faster)\n',t,t0/t);

if 0
  % repmat is faster than ones
  tic; for i = 1:niter ones(300,1000); end; toc
  tic; for i = 1:niter xrepmat(ones(300,1),1,1000); end; toc
  tic; for i = 1:niter repmat(1,300,1000); end; toc
end

if 0
  % zeros is faster than repmat (as expected)
  tic; for i = 1:niter zeros(300,1000); end; toc
  tic; for i = 1:niter repmat(zeros(300,1),1,1000); end; toc
  tic; for i = 1:niter xrepmat(0,300,1000); end; toc
end

if 0
  fprintf('new repmat:');
  tic; for i = 1:niter repmat(x',1000,1); end; toc
  fprintf('old repmat:');
  tic; for i = 1:niter xrepmat(x',1000,1); end; toc
end

if 0
  fprintf('new repmat:');
  tic; for i = 1:niter repmat(x,1000,1); end; toc
  fprintf('old repmat:');
  tic; for i = 1:niter xrepmat(x,1000,1); end; toc
end

author: Cloud-CV

repo: object-proposals


text:
function spectra = loadallspectra(file)
% load all spectra files of a Dubna meassurement cycle
%
% Input
%  filename -
%
% Output
%  spectra - []
%

dir = file(1:max(strfind(file,'/')));
maxrot = findmax(dir);
fprintf('load detector: ');

for det=1:19
	fprintf([char(96+det),', ']);

	for rot=0:maxrot

		s = int2str(rot);
		if rot<10, s = ['0',s];end %#ok<AGROW>

		spec = loadspectra([file,'.',char(96+det),s]);
		if ~isempty(spec),  spectra(:,rot+1,det)=spec;end %#ok<AGROW>
	end
end
disp('.');

function m = findmax(directory)

l = dir(directory);
if isempty(l), error('directory not found or empty');end
for i=3:length(l), m(i-2) = max([0,str2num(l(i).name(end-1:end))]);end %#ok<ST2NM,AGROW>
m = max(m);

author: mtex-toolbox

repo: mtex

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


text:
function dist = polygon_point_dist ( n, v, p )

%*****************************************************************************80
%
%% POLYGON_POINT_DIST: distance ( polygon, point ).
%
%  Discussion:
%
%    Thanks to Stefano Zappacosta for pointing out that the arguments
%    passed to SEGMENT_POINT_DIST_2D needed to be transposed,
%    27 June 2006.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 June 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of vertices.
%
%    Input, real V(2,N), the vertices of the polygon.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, real DIST, the distance from the point to the polygon.
%

%
%  Find the distance to each of the line segments.
%
  dist = inf;

  for j = 1 : n

    jp1 = j + 1;
    if ( n < jp1 ) then
      jp1 = jp1 - n;
    end

    dist2 = segment_point_dist ( v(1:2,j), v(1:2,jp1), p );

    if ( dist2 < dist )
      dist = dist2;
    end

  end

  return
end