text:
function varargout = panel_dipoles(varargin)
% PANEL_DIPOLES: Create a panel to add/remove/edit scouts attached to a given 3DViz figure.
%
% USAGE:  bstPanelNew = panel_dipoles('CreatePanel')
%                       panel_dipoles('UpdatePanel')
%                       panel_dipoles('CurrentFigureChanged_Callback')
%                       panel_dipoles('SetGoodness', Goodness)

% @=============================================================================
% 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: Elizabeth Bock, 2010-2017
%          Francois Tadel, 2010-2020

eval(macro_method);
end

%% ===== CREATE PANEL =====
function bstPanelNew = CreatePanel() %#ok<DEFNU>
    panelName = 'Dipoles';
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    import org.brainstorm.list.*;

    % The the better JList height
    LABEL_WIDTH      = java_scaled('value', 40);
    SLIDER_WIDTH     = java_scaled('value', 5);
    DEFAULT_HEIGHT   = java_scaled('value', 22);
    LABEL_FONT = java.awt.Font('Arial', java.awt.Font.PLAIN, 9);
    % Font size for the lists
    fontSize = round(11 * bst_get('InterfaceScaling') / 100);

    % Main panel container
    jPanelNew = java_create('javax.swing.JPanel');
    jPanelNew.setLayout(BoxLayout(jPanelNew, BoxLayout.PAGE_AXIS));
    jPanelNew.setBorder(BorderFactory.createEmptyBorder(10,10,0,10));

    % ===== PANEL: LOADED DIPOLES =====
    jPanelList = gui_component('Panel');
    jPanelList.setPreferredSize(java_scaled('dimension', 100, 150));
    jPanelList.setMaximumSize(java_scaled('dimension', 500, 150));
        jBorder = java_scaled('titledborder', 'Loaded dipoles');
        jPanelList.setBorder(jBorder);
        % Groups
        jListDipoles = JList();
        jListDipoles.setCellRenderer(BstStringListRenderer(fontSize));
        java_setcb(jListDipoles, 'ValueChangedCallback', @DipolesListValueChanged_Callback, ...
                                 'KeyPressedCallback',   @DipolesListKeyTyped_Callback);
        jPanelScrollList = JScrollPane();
        jPanelScrollList.getLayout.getViewport.setView(jListDipoles);
        jPanelScrollList.setHorizontalScrollBarPolicy(jPanelScrollList.HORIZONTAL_SCROLLBAR_NEVER);
        jPanelScrollList.setBorder([]);
        % Subsets
        jListSubsets = JList();
        jListSubsets.setCellRenderer(BstStringListRenderer(fontSize));
        java_setcb(jListSubsets, 'ValueChangedCallback', @SubsetsListValueChanged_Callback, ...
                                 'MouseClickedCallback', @SubsetsListValueChanged_Callback, ...
                                 'KeyPressedCallback',   @DipolesListKeyTyped_Callback);

        jPanelScrollSubset = JScrollPane();
        jPanelScrollSubset.getLayout.getViewport.setView(jListSubsets);
        jPanelScrollSubset.setHorizontalScrollBarPolicy(jPanelScrollList.HORIZONTAL_SCROLLBAR_NEVER);
        jPanelScrollSubset.setBorder([]);

    jSplitDipoles = JSplitPane(JSplitPane.HORIZONTAL_SPLIT, jPanelScrollSubset, jPanelScrollList );
    jSplitDipoles.setResizeWeight(0.3);
    jSplitDipoles.setDividerSize(4);
    jSplitDipoles.setBorder([]);

    jPanelList.add(jSplitDipoles, BorderLayout.CENTER);
    jPanelList.setVisible(0);
    jPanelNew.add(jPanelList);

    % ===== PANEL: FILTER =====
    jPanelFilter = gui_river([0,1], [0,0,0,0], 'Filter dipoles');
        % Goodness title
        jTitleGoodness = gui_component('Label', jPanelFilter, 'br', 'Goodness:');
        % Goodness slider
        jSliderGoodness = JSlider(0, 100, 0);
        jSliderGoodness.setPreferredSize(Dimension(SLIDER_WIDTH, DEFAULT_HEIGHT));
        java_setcb(jSliderGoodness, 'MouseReleasedCallback', @SliderGoodness_Callback,  'KeyPressedCallback', @SliderGoodness_Callback);
        jPanelFilter.add('tab hfill', jSliderGoodness);
        % Goodness label
        jLabelGoodness = gui_component('Label', jPanelFilter, [], '     0%', {JLabel.RIGHT, Dimension(LABEL_WIDTH, DEFAULT_HEIGHT)});
        jLabelGoodness.setFont(LABEL_FONT);

        % Confidence Volume title
        jTitleConfVol = gui_component('Label', jPanelFilter, 'br', 'ConfVol:');
        % Confidence volume slider
        jSliderConfVol = JSlider(0, 100, 100);
        jSliderConfVol.setPreferredSize(Dimension(SLIDER_WIDTH, DEFAULT_HEIGHT));
        jSliderConfVol.setInverted(1);
        java_setcb(jSliderConfVol, 'MouseReleasedCallback', @SliderConfVol_Callback, ...
                                   'KeyPressedCallback',    @SliderConfVol_Callback);
        jPanelFilter.add('tab hfill', jSliderConfVol);
        % Confidence volume label
        jLabelConfVol = gui_component('Label', jPanelFilter, [], '         ', {JLabel.RIGHT, Dimension(LABEL_WIDTH, DEFAULT_HEIGHT)});
        jLabelConfVol.setFont(LABEL_FONT);

        % Chi-square title
        jTitleKhi2 = gui_component('Label', jPanelFilter, 'br', 'Chi-sqr:');
        % Chi-square slider
        jSliderKhi2 = JSlider(0, 100, 100);
        jSliderKhi2.setPreferredSize(Dimension(SLIDER_WIDTH, DEFAULT_HEIGHT));
        jSliderKhi2.setInverted(1);
        java_setcb(jSliderKhi2, 'MouseReleasedCallback', @SliderKhi2_Callback, ...
                                'KeyPressedCallback',    @SliderKhi2_Callback);
        jPanelFilter.add('tab hfill', jSliderKhi2);
        % Chi-square label
        jLabelKhi2 = gui_component('Label', jPanelFilter, [], '         ', {JLabel.RIGHT, Dimension(LABEL_WIDTH, DEFAULT_HEIGHT)});
        jLabelKhi2.setFont(LABEL_FONT);
    jPanelNew.add(jPanelFilter);

    % ===== PANEL: DISPLAY OPTIONS =====
    jPanelOptions = gui_river([0,1], [2,4,4,0], 'Display filters');
        % Options
        jToggleMaxGoodness = gui_component('Checkbox', jPanelOptions, 'br', 'Show max goodness of fit', [], 'Show only the dipole with the maximal goodness of fit', @(h,ev)FireUpdateDisplayOptions);
        jToggleAllTimes    = gui_component('Checkbox', jPanelOptions, 'br', 'Show all time', [], 'Show all time points', @(h,ev)FireUpdateDisplayOptions);
        jToggleSelTimes    = gui_component('Checkbox', jPanelOptions, 'br', 'Show only preferred times', [], 'Show only the dipoles at custom selected times (default time=0)', @(h,ev)FireUpdateDisplayOptions);
        jToggleAllTimes.setSelected(1);
        jButtonSetSel = gui_component('Button', jPanelOptions, '', 'Set', Insets(1,5,1,5), 'Set the current time as the preferred time for the loaded group(s)', @SetSelectedTime_Callback);
        jButtonSetSel.setFocusable(0);
    jPanelNew.add(jPanelOptions);

    % ===== PANEL: COLOR =====
    jPanelColor = gui_river([0,1], [2,4,4,0], 'Color dipoles according to:');
        colorButtonGroup = javax.swing.ButtonGroup();
        jToggleColorTime     = gui_component('Radio', jPanelColor, 'br', 'Time', colorButtonGroup, 'Set dipole colors according to the time instant', @CheckColor_Callback);
        jToggleColorGoodness = gui_component('Radio', jPanelColor, 'br', 'Goodness of fit', colorButtonGroup, 'Set dipole colors according to the goodness of fit', @CheckColor_Callback);
        jToggleColorGroup    = gui_component('Radio', jPanelColor, 'br', 'Group', colorButtonGroup, 'Set dipole colors according to the dipole group', @CheckColor_Callback);
        jToggleColorTime.setSelected(1);

        %dipole size
        jTitleDipSize = gui_component('Label', jPanelColor, 'br', 'Point display size:');
        jSliderDipSize = JSlider(5, 20, 8);
        jSliderDipSize.setPreferredSize(Dimension(SLIDER_WIDTH, DEFAULT_HEIGHT));
        java_setcb(jSliderDipSize, 'MouseReleasedCallback', @FireUpdateDisplayOptions, ...
                                   'KeyPressedCallback',    @FireUpdateDisplayOptions);
        jPanelColor.add('tab hfill', jSliderDipSize);
        %dipole tail size
        jTitleTailWidth = gui_component('Label', jPanelColor, 'br', 'Tail display width:');
        jSliderTailWidth = JSlider(4, 40, 8);
        jSliderTailWidth.setPreferredSize(Dimension(SLIDER_WIDTH, DEFAULT_HEIGHT));
        java_setcb(jSliderTailWidth, 'MouseReleasedCallback', @FireUpdateDisplayOptions, ...
                                   'KeyPressedCallback',    @FireUpdateDisplayOptions);
        jPanelColor.add('tab hfill', jSliderTailWidth);
    jPanelNew.add(jPanelColor);

    % Set max panel sizes
    drawnow;
    jPanelList.setMaximumSize(java.awt.Dimension(jPanelList.getMaximumSize().getWidth(), jPanelList.getPreferredSize().getHeight()));
    jPanelOptions.setMaximumSize(java.awt.Dimension(jPanelOptions.getMaximumSize().getWidth(), jPanelOptions.getPreferredSize().getHeight()));
    jPanelFilter.setMaximumSize(java.awt.Dimension(jPanelFilter.getMaximumSize().getWidth(), jPanelFilter.getPreferredSize().getHeight()));
    jPanelColor.setMaximumSize(java.awt.Dimension(jPanelColor.getMaximumSize().getWidth(), jPanelColor.getPreferredSize().getHeight()));
    % Add an extra glue at the end, so that panel stay small
    jPanelNew.add(Box.createVerticalGlue());

    % Create the BstPanel object that is returned by the function
    % => constructor BstPanel(jHandle, panelName, sControls)
    bstPanelNew = BstPanel(panelName, ...
                           jPanelNew, ...
                           struct('jPanelList',            jPanelList, ...
                                  'jPanelOptions',         jPanelOptions, ...
                                  'jPanelFilter',          jPanelFilter, ...
                                  'jListDipoles',          jListDipoles, ...
                                  'jListSubsets',          jListSubsets, ...
                                  'jToggleAllTimes',       jToggleAllTimes, ...
                                  'jToggleColorGroup',     jToggleColorGroup, ...
                                  'jToggleColorTime',      jToggleColorTime, ...
                                  'jToggleColorGoodness',  jToggleColorGoodness, ...
                                  'jToggleMaxGoodness',    jToggleMaxGoodness, ...
                                  'jTitleGoodness',        jTitleGoodness, ...
                                  'jLabelGoodness',        jLabelGoodness, ...
                                  'jSliderGoodness',       jSliderGoodness, ...
                                  'jTitleKhi2',            jTitleKhi2, ...
                                  'jLabelKhi2',            jLabelKhi2, ...
                                  'jSliderKhi2',           jSliderKhi2, ...
                                  'jTitleConfVol',         jTitleConfVol, ...
                                  'jLabelConfVol',         jLabelConfVol, ...
                                  'jSliderConfVol',        jSliderConfVol, ...
                                  'jToggleSelTimes',       jToggleSelTimes, ...
                                  'jButtonSetSel',         jButtonSetSel, ...
                                  'jTitleDipSize',         jTitleDipSize, ...
                                  'jSliderDipSize',        jSliderDipSize, ...
                                  'jTitleTailWidth',       jTitleTailWidth, ...
                                  'jSliderTailWidth',      jSliderTailWidth));

%% ===== INTERNAL CALLBACKS =====
    function SliderGoodness_Callback(varargin)
        % Update text field
        val = jSliderGoodness.getValue();
        jLabelGoodness.setText(sprintf('%d%%', val));
        % Update display options
        FireUpdateDisplayOptions();
    end

    function SliderKhi2_Callback(varargin)
        % Update display options
        FireUpdateDisplayOptions();
    end

    function SliderConfVol_Callback(varargin)
        % Update display options
        FireUpdateDisplayOptions();
    end
end

%% =================================================================================
%  === CONTROLS CALLBACKS  =========================================================
%  =================================================================================
%% ===== CHECKBOX: 3D COLOR =====
function CheckColor_Callback(varargin)
    % Get current figure
    hFig = bst_figures('GetCurrentFigure', '3D');
    if isempty(hFig)
        return
    end
    % Update display
    FireUpdateDisplayOptions();
    % Update colormap
    FireUpdateSurfaceColormap();
end

%% ===== LIST SELECTION CHANGED CALLBACK =====
function DipolesListValueChanged_Callback(h, ev)
    if ~ev.getValueIsAdjusting()
        FireUpdateDisplayOptions();
    end
end

%% ===== LIST KEY TYPED CALLBACK =====
function DipolesListKeyTyped_Callback(h, ev)
    import java.awt.event.KeyEvent;
    if ismember(ev.getKeyCode(), [KeyEvent.VK_LEFT, KeyEvent.VK_RIGHT])
        panel_time('TimeKeyCallback', ev);
    end
end

%% ===== SUBSET SELECTION CHANGED CALLBACK =====
function SubsetsListValueChanged_Callback(h, ev)
    % Get panel controls
    ctrl = bst_get('PanelControls','Dipoles');
    if isempty(ctrl)
        return;
    end
    % Get selected subset
    iSubset = double(ctrl.jListSubsets.getSelectedIndices())' + 1;

    % Remove JList callback
    bakCallback = java_getcb(ctrl.jListDipoles, 'ValueChangedCallback');
    java_setcb(ctrl.jListDipoles, 'ValueChangedCallback', []);

    hFig = bst_figures('GetCurrentFigure');
    DipolesInfo = GetDipolesForFigure(hFig);
    iNames = [];
    % Find all DipoleNames that contain the subset number
    for i=1:length(iSubset)
        strSubset = ['(' num2str(iSubset(i)) ')'];
        x=strfind([DipolesInfo.DipoleNames], strSubset);
        iNames = cat(2,iNames,find(~cellfun(@isempty,x))-1);
    end
    % Select the DipoleNames, or All, if no match...
    if isempty(iNames)
        iNames = 0:length(DipolesInfo.DipoleNames)-1;
    end
    ctrl.jListDipoles.setSelectedIndices(iNames);
    java_setcb(ctrl.jListDipoles, 'ValueChangedCallback', bakCallback);

    % Display/hide dipoles
    FireUpdateDisplayOptions()
end
%% ===== SET TIME SELECTION CALLBACK =====
function SetSelectedTime_Callback(h,ev)
    global GlobalData
    hFig = bst_figures('GetCurrentFigure');
    % Get the dipole info structure
    DipolesInfo = GetDipolesForFigure(hFig);
    % Get selected dipole(s)
    iSel = GetDisplayedDipoles(hFig);
    % Set the current time as the selected time for the displayed groups
    DipolesInfo.PreferredTimes(iSel) = GlobalData.UserTimeWindow.CurrentTime;

    % Save the times back to the dipole file
    DipolesFileInfo = getappdata(hFig, 'Dipoles');
    DipolesFile = file_fullpath(DipolesFileInfo.FileName);
    save(DipolesFile,'-struct','DipolesInfo','PreferredTimes','-append')

    % Update the loaded dipoles
    [iDS, iDipoles] = bst_memory('LoadDipolesFile', DipolesFile);
    GlobalData.DataSet(iDS).Dipoles(iDipoles).PreferredTimes = DipolesInfo.PreferredTimes;

    % Display/hide dipoles
    FireUpdateDisplayOptions()
end
%% =================================================================================
%  === EXTERNAL PANEL CALLBACKS  ===================================================
%  =================================================================================
%% ===== UPDATE CALLBACK =====
function UpdatePanel() %#ok<DEFNU>
    hFig = bst_figures('GetCurrentFigure');
    UpdateDipolesList(hFig);
end

%% ===== CURRENT FIGURE CHANGED =====
function CurrentFigureChanged_Callback(hFig) %#ok<DEFNU>
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Dipoles');
    if isempty(ctrl)
        return;
    end
    % Remove JList callback
    bakCallback = java_getcb(ctrl.jListDipoles, 'ValueChangedCallback');
    java_setcb(ctrl.jListDipoles, 'ValueChangedCallback', []);
    % Update dipoles list
    UpdateDipolesList(hFig);
    % Restore callback
    java_setcb(ctrl.jListDipoles, 'ValueChangedCallback', bakCallback);
end

%% ===== CURRENT TIME CHANGED =====
function CurrentTimeChangedCallback(hFig) %#ok<DEFNU>
    PlotSelectedDipoles(hFig);
end

%% ===== LOAD DIPOLES =====
function AddDipoles(hFig, DipolesFile, isDisplay) %#ok<DEFNU>
    global GlobalData;
    if (nargin < 3)
        isDisplay = 1;
    end
    % Make file name relative
    DipolesFile = file_short(DipolesFile);
    % Get dipoles currently shown in figure
    DipolesInfo = getappdata(hFig, 'Dipoles');
    % If dipoles files already displayed: return
    if ~isempty(DipolesInfo) && any(file_compare({DipolesInfo.FileName}, DipolesFile))
        warning('This dipoles file is already displayed. Ignoring... ');
        return
    end
    % Load file
    [iDS, iDipoles] = bst_memory('LoadDipolesFile', DipolesFile);
    if isempty(iDS)
        return
    end
    % Add new dipole structure in DipolesInfo
    iNewDip = length(DipolesInfo) + 1;
    DipolesInfo(iNewDip).FileName = GlobalData.DataSet(iDS).Dipoles(iDipoles).FileName;
    DipolesInfo(iNewDip).Selected = [];
    % Set dipoles in figure application data
    setappdata(hFig, 'Dipoles',     DipolesInfo);
    setappdata(hFig, 'SubjectFile', GlobalData.DataSet(iDS).SubjectFile);
    setappdata(hFig, 'StudyFile',   GlobalData.DataSet(iDS).StudyFile);
    % Update dipoles list
    UpdateDipolesList(hFig, 1);
    % Select "dipoles" tab
    gui_brainstorm('SetSelectedTab', 'Dipoles');
    % Plot selected dipoles
    if isDisplay
        FireUpdateDisplayOptions();
        % Update colormap
        FireUpdateSurfaceColormap();
    end
end

%% ===== GET DIPOLES FOR FIGURE =====
function DipolesInfo = GetDipolesForFigure(hFig)
    global GlobalData;
    % Initialize returned structure
    DipolesInfo.DipoleNames = {};
    DipolesInfo.Dipole = [];
    DipolesInfo.Subset = [];
    DipolesInfo.PreferredTimes = [];
    DipolesInfo.DisplayColorType = 'time';
    DipolesInfo.DisplayAllTime = 0;
    DipolesInfo.DisplayMaxGoodness = 0;
    DipolesInfo.DisplaySelTimes = 0;

    % Get Dipoles description in figure
    DipolesApp = getappdata(hFig, 'Dipoles');
    if isempty(DipolesApp)
        return
    end
    % For each dipoles file displayed in this figure
    for i = 1:length(DipolesApp)
        % Get loaded dipoles
        [iDS, iDipoles] = bst_memory('GetDataSetDipoles', DipolesApp(i).FileName);
        % Add description to dipoles
        DipolesInfo.Dipole = [DipolesInfo.Dipole, GlobalData.DataSet(iDS).Dipoles(iDipoles).Dipole];
        DipolesInfo.DipoleNames = cat(2, DipolesInfo.DipoleNames, GlobalData.DataSet(iDS).Dipoles(iDipoles).DipoleNames);
        % Check for subsets
        DipolesInfo.Subset = cat(2, DipolesInfo.Subset, GlobalData.DataSet(iDS).Dipoles(iDipoles).Subset);
        DipolesInfo.PreferredTimes = cat(2, DipolesInfo.PreferredTimes, GlobalData.DataSet(iDS).Dipoles(iDipoles).PreferredTimes);
    end
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Dipoles');
    if isempty(ctrl)
        return;
    end
    % Goodness
    sScales = GetSliderScales('goodness',hFig);
    if isempty([sScales.min])
        ctrl.jLabelGoodness.setEnabled(0);
        ctrl.jTitleGoodness.setEnabled(0);
        ctrl.jSliderGoodness.setEnabled(0);
        ctrl.jToggleColorGoodness.setEnabled(0);
        ctrl.jToggleMaxGoodness.setEnabled(0);
    end
    % Khi2
    sScales = GetSliderScales('khi2',hFig);
    if isempty([sScales.min])
        ctrl.jLabelKhi2.setEnabled(0);
        ctrl.jTitleKhi2.setEnabled(0);
        ctrl.jSliderKhi2.setEnabled(0);
    end
    % ConfVol
    sScales = GetSliderScales('confvol',hFig);
    if isempty([sScales.min])
        ctrl.jTitleConfVol.setEnabled(0);
        ctrl.jLabelConfVol.setEnabled(0);
        ctrl.jSliderConfVol.setEnabled(0);
    end

    % Make the list visible if more than one set
    nDipoles = length(DipolesInfo.DipoleNames);
    nSubsets = length(DipolesInfo.Subset);
    if (nDipoles > 2) || (nSubsets > 2)
       ctrl.jPanelList.setVisible(1);
    end

    % Set the display color scheme
    if ctrl.jToggleColorTime.isSelected()
        DipolesInfo.DisplayColorType = 'time';
    elseif ctrl.jToggleColorGoodness.isSelected()
        DipolesInfo.DisplayColorType = 'goodness';
    elseif ctrl.jToggleColorGroup.isSelected()
        DipolesInfo.DisplayColorType = 'group';
    end

    % Determine if there is any display filters
    % Show preferred time only
    if ctrl.jToggleSelTimes.isSelected()
        DipolesInfo.DisplaySelTimes = 1;
        ctrl.jToggleMaxGoodness.setSelected(0);
        ctrl.jToggleAllTimes.setSelected(0);
        % Do not allow setting preferred time when this filter is on
        ctrl.jButtonSetSel.setEnabled(0);
    else
        ctrl.jButtonSetSel.setEnabled(1);
    end
    % Show all time
    if ctrl.jToggleAllTimes.isSelected()
        DipolesInfo.DisplayAllTime = 1;
        % Do not allow setting preferred time when all time is displayed
        ctrl.jButtonSetSel.setEnabled(0);
    end
    % Show maximum goodness of selected dipoles
    if ctrl.jToggleMaxGoodness.isSelected()
        DipolesInfo.DisplayMaxGoodness = 1;
    end
end

%% ===== UPDATE DIPOLES LIST =====
function UpdateDipolesList(hFig, isSelectAll)
    % Parse inputs
    if (nargin < 2) || isempty(isSelectAll)
        isSelectAll = 0;
    end
    % Get Dipoles description in figure
    if ~isempty(hFig)
        DipolesInfo = GetDipolesForFigure(hFig);
        if isempty(DipolesInfo)
            return
        end
        nDipoles = length(DipolesInfo.DipoleNames);
        nSubsets = length(DipolesInfo.Subset);
    else
        nDipoles = 0;
        nSubsets = 0;
    end
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Dipoles');
    if isempty(ctrl)
        return;
    end
    % ==== Dipole List
    % Remove JList callback
    bakCallback = java_getcb(ctrl.jListDipoles, 'ValueChangedCallback');
    java_setcb(ctrl.jListDipoles, 'ValueChangedCallback', []);

    % Create a new empty list
    listModel = java_create('javax.swing.DefaultListModel');
    % Add an item in list for each scout found for target figure
    for i = 1:nDipoles
        listModel.addElement(DipolesInfo.DipoleNames{i});
    end
    % Update list model
    ctrl.jListDipoles.setModel(listModel);

    % ==== Subset List
    % Remove JList callback
    sbakCallback = java_getcb(ctrl.jListSubsets, 'ValueChangedCallback');
    java_setcb(ctrl.jListSubsets, 'ValueChangedCallback', []);

    % Create a new empty list
    subsetModel = java_create('javax.swing.DefaultListModel');
    % Add an item in list for each subset
    for i = 1:nSubsets
        subsetModel.addElement(['Subset #' num2str(DipolesInfo.Subset(i))]);
    end

    % Update list model
    ctrl.jListSubsets.setModel(subsetModel);
    drawnow

    % Select all dipoles
    if isSelectAll && (nDipoles > 0)
        iSel = 1:nDipoles;
    % Get dipoles displayed in this figure
    else
        iSel = GetDisplayedDipoles(hFig);
    end
    % Change list selection
    if ~isempty(iSel)
        ctrl.jListDipoles.setSelectedIndices(iSel - 1);
    end

    % Restore callback
    drawnow
    java_setcb(ctrl.jListDipoles, 'ValueChangedCallback', bakCallback);
    java_setcb(ctrl.jListSubsets, 'ValueChangedCallback', sbakCallback);

end

%% ===== UPDATE DISPLAY OPTIONS ======
function FireUpdateDisplayOptions(varargin)
    % Progress bar
    isProgress = bst_progress('IsVisible');
    if ~isProgress
        bst_progress('start', 'Dipoles', 'Updating figure...');
    end
    % Update all figures with selected display options
    hFigs = bst_figures('GetAllFigures');
    for iFig = 1:length(hFigs)
        hFig = hFigs(iFig);
        FigureId = getappdata(hFig, 'FigureId');
        DipolesInfo = GetDipolesForFigure(hFig);
        % Redraw dipoles for current figure
        if ~isempty(FigureId) && ~isempty(DipolesInfo.Dipole)
            switch (FigureId.Type)
                case 'MriViewer'
                    % Nothing specific to update
                case '3DViz'
                    % === 3D: SPECIAL COLORMAPS ===
                    CmapInfo = getappdata(hFig, 'Colormap');
                    switch DipolesInfo.DisplayColorType
                        % Color = Time
                        case 'time'
                            CmapInfo.Type = 'time';
                            CmapInfo.AllTypes{end+1} = 'time';
                            CmapInfo.AllTypes = unique(CmapInfo.AllTypes);
                        % Color = Goodness of fit
                        case 'goodness'
                            CmapInfo.Type = 'percent';
                            CmapInfo.AllTypes{end+1} = 'percent';
                            CmapInfo.AllTypes = unique(CmapInfo.AllTypes);
                        % Color = Group (no colormap)
                        case 'group'
                            CmapInfo.AllTypes = setdiff(CmapInfo.AllTypes, {'time', 'percent', 'stat1'});
                            if ~isempty(CmapInfo.AllTypes)
                                CmapInfo.Type = CmapInfo.AllTypes{1};
                            else
                                CmapInfo.Type = [];
                            end
                    end
                    % Update 3D figure configuration
                    setappdata(hFig, 'Colormap', CmapInfo);
            end
        end
        % Redraw dipoles for current figure
        RedrawDipoles(hFig);
    end
    % Close progress bar
    if ~isProgress
        bst_progress('stop');
    end
end

%% ===== UPDATE COLORMAP =====
function FireUpdateSurfaceColormap()
    [hFigs] = bst_figures('GetAllFigures');
    for iFig = 1:length(hFigs)
        hFig = hFigs(iFig);
        panel_surface('UpdateSurfaceColormap', hFig);
    end
end
%% ===== GET DISPLAYED DIPOLES =====
function iSel = GetDisplayedDipoles(hFig)
    iSel = [];
    if ~isempty(hFig) && ishandle(hFig)
        DipolesInfo = getappdata(hFig, 'Dipoles');
        if isfield(DipolesInfo, 'Selected')
            iSel = DipolesInfo.Selected;
        end
    end
end

%% ===== GET SELECTED DIPOLES =====
function [sDipoles, iDipoles] = GetSelectedDipoles(hFig)
    global GlobalData
    sDipoles = [];
    iDipoles = [];
    % Get Dipoles description in figure
    DipolesInfo = GetDipolesForFigure(hFig);
    if isempty(DipolesInfo) || isempty(DipolesInfo.Dipole)
        return
    end
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Dipoles');
    if isempty(ctrl)
        return;
    end
    % Get the selected dipoles in JList
    iSelList = double(ctrl.jListDipoles.getSelectedIndices() + 1);
    % Save selection list in figure
    DipolesApp = getappdata(hFig, 'Dipoles');
    DipolesApp.Selected = iSelList;
    setappdata(hFig, 'Dipoles', DipolesApp);

    % Select time points
    if DipolesInfo.DisplayAllTime || DipolesInfo.DisplaySelTimes
        % Select all time points when this is asked for explicitly or when
        % using the custom selected times filter
        iSelTime = 1:length(DipolesInfo.Dipole);
    else
        % Select only those time points that are closest to the current time
        curTime = GlobalData.UserTimeWindow.CurrentTime;
        if isempty(curTime)
            iSelTime = find(abs([DipolesInfo.Dipole.Time] - [DipolesInfo.Dipole(1).Time]) < 1e-6);
        else
            iSelTime = find(abs([DipolesInfo.Dipole.Time] - curTime) < 1e-6);
        end
    end

    % Get dipoles that fit the selection
    if isempty(iSelTime)
        return
    end
    iDipoles = iSelTime(ismember([DipolesInfo.Dipole(iSelTime).Index], iSelList));

    % Show only the preferred time point dipoles
    if DipolesInfo.DisplaySelTimes && ~isempty(iDipoles)
        indexNums = [DipolesInfo.Dipole(iDipoles).Index];
        % Find unique index numbers (groups)
        nIndexNums = unique(indexNums);
        % Loop through groups to find the selected time of each group
        iGroupTime = zeros(1,length(nIndexNums));
        for n=1:length(nIndexNums)
            ind = find(indexNums == nIndexNums(n));
            groupTimes = [DipolesInfo.Dipole(ind).Time];
            t = DipolesInfo.PreferredTimes(nIndexNums(n));
            iGroupTime(n) = bst_closest(groupTimes,t) + ind(1) - 1;
        end
        % Get the indices of the maximums
        iDipoles = iDipoles(iGroupTime);
    end

    % Slider thresholds

    % Goodness
    sScales = GetSliderScales('goodness', hFig);
    if ~isempty([sScales.min])
        minGoodness = double(ctrl.jSliderGoodness.getValue()) ./ 100;
        iDipoles = iDipoles([DipolesInfo.Dipole(iDipoles).Goodness] >= minGoodness);
    end
    % Khi2
    sScales = GetSliderScales('khi2', hFig);
    if ~isempty([sScales.min])
        val = double(ctrl.jSliderKhi2.getValue());
        minKhi2 = val/100 * sScales.max;
        iDipoles = iDipoles([DipolesInfo.Dipole(iDipoles).Khi2] <= minKhi2);
        % Update text field
        ctrl.jLabelKhi2.setText(sprintf('%1.1e', minKhi2));
    end
    % ConfVol
    sScales = GetSliderScales('confvol', hFig);
    if ~isempty([sScales.min])
        val = double(ctrl.jSliderConfVol.getValue());
        minConfVol = val/100 * sScales.max;
        iDipoles = iDipoles([DipolesInfo.Dipole(iDipoles).ConfVol] <= minConfVol);
        % Update text field
        ctrl.jLabelConfVol.setText(sprintf('%1.1e', minConfVol));
    end

    % Show only maximum goodness of fit
    if DipolesInfo.DisplayMaxGoodness && ~isempty(iDipoles)
        indexNums = [DipolesInfo.Dipole(iDipoles).Index];
        % Find unique index numbers (groups)
        nIndexNums = unique(indexNums);
        goodVals = [DipolesInfo.Dipole(iDipoles).Goodness];
        % Loop through groups to find the max goodness of each group
        iMaxGoodness = zeros(1,length(nIndexNums));
        for n=1:length(nIndexNums)
            ind = find(indexNums == nIndexNums(n));
            z = zeros(1, length(goodVals));
            % Get the goodness values for the group
            z(ind) = goodVals(ind);
            % Get the maximum goodness for the group
            [maxGoodness, iMaxGoodness(n)] = max(z);
        end
        % Get the indices of the maximums
        iDipoles = iDipoles(iMaxGoodness);
    end

    % Get dipoles selected in time and index at the same time
    sDipoles = DipolesInfo.Dipole(iDipoles);
end

%% ===== REDRAW DIPOLES =====
function RedrawDipoles(hFig)
    % Get figure type
    FigureId = getappdata(hFig, 'FigureId');
    % Redraw dipoles for current figure
    if ~isempty(FigureId)
        switch (FigureId.Type)
            case 'MriViewer'
                panel_surface('UpdateSurfaceData', hFig)
            case '3DViz'
                PlotSelectedDipoles(hFig);
        end
    end
end

%% ===== PLOT DIPOLES =====
function PlotSelectedDipoles(hFig)
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Dipoles');
    if isempty(ctrl)
        return;
    end
    % Get axis
    hAxes = findobj(hFig, '-depth', 1, 'Tag', 'Axes3D');
    if isempty(hAxes)
        return
    end
    % Delete previous dipoles display
    hPoints = findobj(hAxes, '-depth', 1, 'Tag', 'DipolesLoc');
    hLines  = findobj(hAxes, '-depth', 1, 'Tag', 'DipolesOrient');
    if ~isempty(hPoints)
        delete([hPoints(:); hLines(:)]);
    end

    % Get all dipoles
    DipolesInfo = GetDipolesForFigure(hFig);
    % Get selected dipoles
    [sDipoles, iDipoles] = GetSelectedDipoles(hFig);
    if isempty(sDipoles)
        return
    end

    % ===== GROUP DIPOLES =====
    % If color-coded time => group dipoles by time
    ColormapInfo = getappdata(hFig, 'Colormap');
    switch DipolesInfo.DisplayColorType
        case 'time'
            % Get all times for all dipoles
            allTimes = [DipolesInfo.Dipole.Time];
            % Get times for selected dipoles
            selTimes = [sDipoles.Time];
            dipGroupTag = 1:length(sDipoles);
            % Get time colormap
            sColormap = bst_colormaps('GetColormap', ColormapInfo.Type);

            if length(allTimes) < 2
                % Get time limits for this figure
                hAxes = [findobj(hFig, '-depth', 1, 'Tag', 'Axes3D'), ...
                         findobj(hFig, '-depth', 1, 'Tag', 'axc'), ...
                         findobj(hFig, '-depth', 1, 'Tag', 'axa'), ...
                         findobj(hFig, '-depth', 1, 'Tag', 'axs')];
                CLim = get(hAxes(1), 'CLim');
            else
                % Get time limits for the selected group(s)
                CLim = [min(allTimes), max(allTimes)];
            end
            % If start = stop time: cannot create a colormap
            if (CLim(1) == CLim(2))
                CLim(2) = CLim(1) + 0.001;
            end
            % Get color table for dipoles
            if isequal(CLim, [0,1]) && (length(allTimes) < 2)
                iCol = 1;
            else
                iCol = round((selTimes - CLim(1)) / abs(CLim(2) - CLim(1)) * (length(sColormap.CMap) - 1)) + 1;
                iCol = bst_saturate(iCol, [1, length(sColormap.CMap)]);
            end
            ColorTable = sColormap.CMap(iCol, :);

        case 'goodness'
            % Group tag is the Index
            dipGroupTag = 1:length(sDipoles);
            % Get color table for dipoles
            sColormap = bst_colormaps('GetColormap', ColormapInfo.Type);
            selGoodness = bst_saturate([sDipoles.Goodness] * 100, [0, 100]);
            CLim = [0, 100];
            % Get color table for dipoles
            iCol = round((selGoodness - CLim(1)) / abs(CLim(2) - CLim(1)) * (length(sColormap.CMap) - 1)) + 1;
            ColorTable = sColormap.CMap(iCol, :);

        case 'group'
            % Group tag is the Index
            dipGroupTag = [sDipoles.Index];
            % Get color table for dipoles
            ColorTable = GetDipolesColorTable();
    end
    % Group dipoles by index
    dipGroups = unique(dipGroupTag);
    % Get display properties
    pointSize = ctrl.jSliderDipSize.getValue;
    lineWidth = 1.5*(ctrl.jSliderTailWidth.getValue/10);

    % ===== DISPLAY GROUPS =====
    % Get number of different groups
    for iGroup = 1:length(dipGroups)
        % Get all the dipoles of this group
        iDipGroup = find(dipGroupTag == dipGroups(iGroup));
        % Get dipole locations
        Loc = [sDipoles(iDipGroup).Loc];
        % Compute dipole orientations (normalized amplitude)
        Orient = [sDipoles(iDipGroup).Amplitude];
        normAmp = sqrt(sum(Orient.^2,1));
        Orient = Orient ./ repmat(normAmp,3,1) .* 0.02;
        % Get color
        iColor = mod(dipGroups(iGroup)-1, size(ColorTable,1)) + 1;

        % Loop added to be able to track each dipole inpendently
        for i = 1:length(iDipGroup)
            % Plot point
            line(Loc(1,i), Loc(2,i), Loc(3,i), ...
                'LineStyle',       'none', ...
                'MarkerFaceColor', ColorTable(iColor,:), ...
                'Marker',          'o', ...
                'MarkerEdgeColor', [.4 .4 .4], ...
                'MarkerSize',      pointSize, ...
                'Tag',             'DipolesLoc', ...
                'UserData',        iDipoles(iDipGroup(i)), ...
                'Parent',          hAxes);
            % Plot orientation
            if lineWidth > 0.70
                line([Loc(1,i); Loc(1,i) + Orient(1,i)], ...
                     [Loc(2,i); Loc(2,i) + Orient(2,i)], ...
                     [Loc(3,i); Loc(3,i) + Orient(3,i)], ...
                    'Color',      ColorTable(iColor,:), ...
                    'LineStyle',  '-', ...
                    'LineWidth',   lineWidth, ...
                    'Tag',        'DipolesOrient', ...
                    'UserData',        iDipoles(iDipGroup(i)), ...
                    'Parent',     hAxes);
            end
        end
    end
    % Force updating this figure before upd
    drawnow
end

%% ===== GET DIPOLES COLOR TABLE =====
function ColorTable = GetDipolesColorTable()
    ColorTable = [0    1    0   ;
                  1    0    0   ;
                  .4   .4   1   ;
                  1    .694 .392;
                  0    1    1   ;
                  1    0    1   ;
                  .4   0    0  ;
                  0    .5   0];
end

%% ===== COMPUTE DENSITY VOLUMES =====
function Cube = ComputeDensity(sMri, sDipoles) %#ok<DEFNU>
    % Initialize cube
    sizeCube = size(sMri.Cube(:,:,:,1));
    Cube = zeros(sizeCube);

    % Gaussian kernel
    A = zeros(5,5,5);
    A(:,:,1) = [0 0 0 0 0;
                0 0 0 0 0;
                0 0 1 0 0;
                0 0 0 0 0;
                0 0 0 0 0];

    A(:,:,2) = [0 0 0 0 0;
                0 1 2 1 0;
                0 2 3 2 0;
                0 1 2 1 0;
                0 0 0 0 0];

    A(:,:,3) = [0 1 1 1 0;
                1 2 3 2 1;
                1 3 4 3 1;
                1 2 3 2 1;
                0 1 1 1 0];
    A(:,:,4) = A(:,:,2);
    A(:,:,5) = A(:,:,1);

    % Loop on each dipole
    for i = 1:length(sDipoles)
        % Vertices: SCS->Voxels
        mriLoc = round(cs_convert(sMri, 'scs', 'voxel', [sDipoles(i).Loc]'));
        % Get the valid indices in the kernel
        iX = [min(2, mriLoc(:,1) - 1),  min(2, sizeCube(1)-mriLoc(:,1))];
        iY = [min(2, mriLoc(:,2) - 1),  min(2, sizeCube(2)-mriLoc(:,2))];
        iZ = [min(2, mriLoc(:,3) - 1),  min(2, sizeCube(3)-mriLoc(:,3))];
        % Each dipole is displayed by a cube of 5*5*5 voxels (gaussian kernel)
        X = mriLoc(:,1) - iX(1) : mriLoc(:,1) + iX(2);
        Y = mriLoc(:,2) - iY(1) : mriLoc(:,2) + iY(2);
        Z = mriLoc(:,3) - iZ(1) : mriLoc(:,3) + iZ(2);
        Cube(X,Y,Z) = Cube(X,Y,Z) + A(3-iX(1):3+iX(2), 3-iY(1):3+iY(2), 3-iZ(1):3+iZ(2));
    end
    % Convert in percentage
    maxCube = max(Cube(:));
    if (maxCube > 0)
        Cube = Cube ./ maxCube .* 100;
    end
end

%% ===== SLIDER SCALES =====
function sScales = GetSliderScales(varargin)
    global GlobalData;

    if size(varargin,2) < 2
        scale = varargin{1};
        hFig = bst_figures('GetCurrentFigure');
    else
        scale = varargin{1};
        hFig = varargin{2};
    end

    % Initialize returned structures
    sScales.label = {};
    sScales.min = [];
    sScales.max = [];
    DipolesInfo.Dipole = [];

    % Get Dipoles description in figure
    DipolesApp = getappdata(hFig, 'Dipoles');
    if isempty(DipolesApp)
        return
    end
    % For each dipoles file displayed in this figure
    for i = 1:length(DipolesApp)
        % Get loaded dipoles
        [iDS, iDipoles] = bst_memory('GetDataSetDipoles', DipolesApp(i).FileName);
        % Add description to dipoles
        DipolesInfo.Dipole = [DipolesInfo.Dipole, GlobalData.DataSet(iDS).Dipoles(iDipoles).Dipole];
    end

    switch scale
        case 'goodness'
            % find scale of goodness
            sScales.label = {'goodness'};
            sScales.min = min([DipolesInfo.Dipole.Goodness]);
            sScales.max = max([DipolesInfo.Dipole.Goodness]);
        case 'khi2'
            % find scale of khi2 value
            sScales.label = {'khi2'};
            sScales.min = min([DipolesInfo.Dipole.Khi2]);
            sScales.max = max([DipolesInfo.Dipole.Khi2]);
        case 'rkhi2'
            % find scale of Reduced khi2 value
            sScales.label = {'rkhi2'};
            sScales.min = min([DipolesInfo.Dipole.DOF]);
            sScales.max = max([DipolesInfo.Dipole.DOF]);
        case 'confvol'
            % find scale of confidence volume values
            sScales.label = {'confvol'};
            sScales.min = min([DipolesInfo.Dipole.ConfVol]);
            sScales.max = max([DipolesInfo.Dipole.ConfVol]);
    end
end

%% ===== SET GOODNESS THRESHOLD =====
function SetGoodness(Goodness) %#ok<DEFNU>
    % Get panel controls
    ctrl = bst_get('PanelControls', 'Dipoles');
    if isempty(ctrl)
        return;
    end
    % Round goodness of fit
    Goodness = round(Goodness * 100);
    % Set the goodness of fit
    if isfield(ctrl, 'jSliderGoodness') && ctrl.jSliderGoodness.isEnabled() && (Goodness >= 0) && (Goodness < 100)
        % Set slider
        ctrl.jSliderGoodness.setValue(Goodness);
        ctrl.jLabelGoodness.setText(sprintf('%d%%', Goodness));
        % Update figures
        FireUpdateDisplayOptions();
    end
end

author: brainstorm-tools


text:
function [x,fval,exitflag,info] = opti_nl2sol(fun,grad,x0,ydata,lb,ub,opts)
%OPTI_NL2SOL Solve a NLS using NL2SOL (Adaptive Nonlinear Least Squares)
%
%   min sum[ (F(x) - ydata)^2 ]       subject to:   lb <= x <= ub
%    x
%
%   x = opti_nl2sol(fun,grad,x0,ydata) solves a NLS where fun is the fitting
%   function. grad is an optional gradient of the fitting function and x0
%   is a starting guess. ydata is the data to fit the function to. Use
%   ydata = zeros() for solving a SNLE.
%
%   x = opti_nl2sol(fun,grad,x0,ydata,lb,ub) solves subject to decision
%   variables bounds lb <= x <= ub. Infinite bounds are OK.
%
%   x = opti_nl2sol(fun,grad,x0,ydata,lb,ub,opts) uses opts to pass optiset
%   options to the solver.
%
%   [x,fval,exitflag,info] = opti_nl2sol(...) returns the objective value at
%   the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR NL2SOL + NL2SNO
%   See referenced ACM License

%   Copyright (C) 2012 Jonathan Currie (IPL)

if(nargin < 7), opts = optiset; end
if(nargin < 6), ub = []; end
if(nargin < 5), lb = []; end
if(nargin < 4), error('NL2SOL requires at least 4 arguments'); end

%Setup display level
opts.display = dispLevel(opts.display);
opts.optiver = optiver;

%Check we have a valid x0
if(isempty(x0) || any(isnan(x0)))
    error('NL2SOL requires an initial guess, x0!');
end

t = tic;
% Run NL2SOL
[x, fval, exitflag, iter, feval] = nl2sol(fun,grad,x0,ydata,lb,ub,opts);

%Collect Results
info.Iterations = iter;
info.FuncEvals = feval;
info.Time = toc(t);
info.Algorithm = 'NL2SOL: Adaptive Nonlinear Least Squares';

switch(exitflag)
    case 1
        info.Status = 'Optimal';
    case 0
        info.Status = 'Exceeded Iterations';
    case -1
        info.Status = 'Infeasible / Could not Converge';
    case -2
        info.Status = 'NL2SOL Error';
    case -5
        info.Status = 'User Exited';
    otherwise
        info.Status = 'NL2SOL Error';
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [idp,IDP_p,IDP_p0] = p2idp(p,p0)

% P2IDP Point to inverse-depth point conversion.
%   P2IDP(P,P0) is the inverse depth point anchored at P0 corresponding to
%   the Euclidean point P.
%
%   An inverse-depth point (idp) is a 6-vector:
%
%           idp = [x0 y0 z0 el az rho]',
%
%   where:
%       x0, z0, y0: anchor: the 3D point P0 where where distance is referred to.
%       el, az: azimuth and elevation of the ray through P that starts at P0.
%       rho: inverse of the distance from point P to P0.
%
%   [idp,IDP_p,IDP_p0] returns the Jacobians wrt P and P0.
%
%   See also IDP2P.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargout == 1
    m   = p-p0;
    rho = 1/sqrt(dot(m,m));
    py  = vec2py(m);

    idp = [p0;py;rho];
else

    m           = p-p0;
    M_p         = 1;
    M_p0        = -1;
    [mn2,MN2_m] = dotJ(m,m);
    rho         = 1/mn2;
    RHO_mn2     = -1/mn2^2;
    [py,PY_m]   = vec2py(m);

    idp = [p0;py;rho];

    RHO_m = RHO_mn2*MN2_m;

    IDP_p = [zeros(3);PY_m*M_p;RHO_m*M_p];
    IDP_p0 = [eye(3);PY_m*M_p0;RHO_m*M_p0];

end

return

%% Jac

syms x y z x0 y0 z0 real
p = [x;y;z];
p0 = [x0;y0;z0];

idp = p2idp(p,p0);

simplify(IDP_p  - jacobian(idp,p))
simplify(IDP_p0 - jacobian(idp,p0))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function gmt_patch(hPatch, sOutput, sPrefix)

% Get data
mXData = get(hPatch, 'XData');
mYData = get(hPatch, 'YData');

% Transform data and save it
[nRow, nCnt] = size(mXData);
if nRow == 1
  mXData = mXData';
end
[nRow, nCnt] = size(mYData);
if nRow == 1
  mYData = mYData';
end
mXYValues = [mXData mYData];
save([sPrefix '_patch.dat'], 'mXYValues', '-ascii');

% Get attributes
vEdgeColor = get(hPatch, 'EdgeColor');
vFaceColor = get(hPatch, 'FaceColor');
fLineWidth = get(hPatch, 'LineWidth');

% Write GMT
hFile = fopen(sOutput, 'a');
fprintf(hFile, '# Plot patch\n');
sG = sprintf('-G%s/%s/%s', num2str(round(vFaceColor(1)*255)), num2str(round(vFaceColor(2)*255)), num2str(round(vFaceColor(3)*255)));
sW = sprintf('-W%s/%s/%s/%s', num2str(fLineWidth), num2str(vEdgeColor(1)*255), num2str(vEdgeColor(2)*255), num2str(vEdgeColor(3)*255));
fprintf(hFile, 'psxy %s $default %s %s -O -K >> $output\n', [sPrefix '_patch.dat'], sW, sG);
fclose(hFile);

author: CelsoReyes


text:
% Detailed description of the factor data structure and related functions
% -----------------------------------------------------------------------
% We will use structures to implement the factor datatype. The code
%
%   phi = struct('var', [3 1 2], 'card', [2 2 2], 'val', ones(1, 8));
%
% creates a factor over variables X_3, X_1, X_2, which are all binary
% valued, because phi.card(1) (the cardinality of X_3, |Val(X_3)|) is 2,
% and likewise for X_1 and X_2. phi has been initialized so that
% phi(X_3, X_1, X_2) = 1 for any assignment to the variables.
%
% A factor's values are stored in a row vector in the .val field
% using an ordering such that the left-most variables as defined in the
% .var field cycle through their values the fastest. More concretely, for
% the factor phi defined above, we have the following mapping from variable
% assignments to the index of the row vector in the .val field:
%
% -+-----+-----+-----+-------------------+
%  | X_3 | X_1 | X_2 | phi(X_3, X_1, X_2)|
% -+-----+-----+-----+-------------------+
%  |  1  |  1  |  1  |     phi.val(1)    |
% -+-----+-----+-----+-------------------+
%  |  2  |  1  |  1  |     phi.val(2)    |
% -+-----+-----+-----+-------------------+
%  |  1  |  2  |  1  |     phi.val(3)    |
% -+-----+-----+-----+-------------------+
%  |  2  |  2  |  1  |     phi.val(4)    |
% -+-----+-----+-----+-------------------+
%  |  1  |  1  |  2  |     phi.val(5)    |
% -+-----+-----+-----+-------------------+
%  |  2  |  1  |  2  |     phi.val(6)    |
% -+-----+-----+-----+-------------------+
%  |  1  |  2  |  2  |     phi.val(7)    |
% -+-----+-----+-----+-------------------+
%  |  2  |  2  |  2  |     phi.val(8)    |
% -+-----+-----+-----+-------------------+
%
%
% We have provided the AssignmentToIndex and IndexToAssignment functions
% that compute the mapping between the assignments A and the variable indices I,
% given D, the cardinality of the variables. Concretely, given a factor phi, if
% phi.val(I) corresponds to the assignment A, i.e. phi(X = A) = phi.val(I) then
%
%   I = AssignmentToIndex(A, D)
%   A = IndexToAssignment(I, D)
%
% For instance, for the factor phi as defined above, with the assignment
%
%    A = [2 1 2]
%
% to X_3, X_1 and X_2 respectively (as defined by phi.var = [3 1 2]), I = 6
% as phi.val(6) corresponds to the value of phi(X_3 = 2, X_1 = 1, X_2 = 2).
% Thus, AssignmentToIndex([2 1 2], [2 2 2]) returns 6, and conversely,
% IndexToAssignment(6, [2 2 2]) returns the vector [2 1 2]. The second
% argument in the function calls corresponds to the cardinality of the
% sample factor phi, phi.card, which is [2 2 2].
%
% More generally, the assignment vector A is a row vector that corresponds
% to assignments to the variables in a factor, with an understanding that the
% variables for which the assignments refer to are given by the .var field
% of the factor.
%
% Giving AssignmentToIndex a matrix A, one assignment per row, will cause it
% to return a vector of indices I, such that I(k) is the index
% corresponding to the assignment in A(k, :) (row k).
%
% Similarly, giving IndexToAssignment a vector I of indices will yield a
% matrix A of assignments, one per row, such that A(k, :) (the kth row of A)
% corresponds to the assignment mapped to by index I(k).
%
% Getting and setting values to factors
% -------------------------------------
%
% We have provided the convenience functions GetValueOfAssignment and
% SetValueOfAssignment so you do not need to manipulate the .val field
% directly for getting and setting values.
%
% For instance, calling
%
%   GetValueOfAssignment(phi, [1 2 1])
%
% yields the value phi(X_3 = 1, X_1 = 2, X_2 = 1). Again, the variables for
% which the assignments refer to are given by the .var field of the factor.
%
% Similarly, executing
%
%    phi = SetValueOfAssignment(phi, [2 2 1], 6)
%
% causes the value of phi(X_3 = 2, X_1 = 2, X_2 = 1) to be set to 6. Note
% that because MATLAB/Octave passes function arguments by value (not by
% reference), SetValueOfAssignment *does not modify* the factor that you
% passed in. Instead, it returns a new factor with the modified value for
% the specified assignment; this is why we reassigned phi to the result of
% SetValueOfAssignment.
%
% More details about these functions are provided in their respective .m
% files.

% Sample Factors and Outputs
% --------------------------
%
% In the following, we have provided you with some sample input factors, as
% well as the output that you should receive when various operations are
% performed on these factors. You may find these sample inputs and outputs
% helpful in debugging your implementation. For instance, FACTORS.PRODUCT
% is the factor you should get when you execute
%
%   FactorProduct(FACTORS.INPUT(1), FACTORS.INPUT(2))
%
% These sample factors define a simple chain Bayesian network over binary
% variables: X_1 -> X_2 -> X_3
%

% FACTORS.INPUT(1) contains P(X_1)
FACTORS.INPUT(1) = struct('var', [1], 'card', [2], 'val', [0.11, 0.89]);

% FACTORS.INPUT(2) contains P(X_2 | X_1)
FACTORS.INPUT(2) = struct('var', [2, 1], 'card', [2, 2], 'val', [0.59, 0.41, 0.22, 0.78]);

% FACTORS.INPUT(3) contains P(X_3 | X_2)
FACTORS.INPUT(3) = struct('var', [3, 2], 'card', [2, 2], 'val', [0.39, 0.61, 0.06, 0.94]);

% Factor Product
% FACTORS.PRODUCT = FactorProduct(FACTORS.INPUT(1), FACTORS.INPUT(2));
% The factor defined here is correct to 4 decimal places.
FACTORS.PRODUCT = struct('var', [1, 2], 'card', [2, 2], 'val', [0.0649, 0.1958, 0.0451, 0.6942]);

% Factor Marginalization
% FACTORS.MARGINALIZATION = FactorMarginalization(FACTORS.INPUT(2), [2]);
FACTORS.MARGINALIZATION = struct('var', [1], 'card', [2], 'val', [1 1]);

% Observe Evidence
% FACTORS.EVIDENCE = ObserveEvidence(FACTORS.INPUT, [2 1; 3 2]);
FACTORS.EVIDENCE(1) = struct('var', [1], 'card', [2], 'val', [0.11, 0.89]);
FACTORS.EVIDENCE(2) = struct('var', [2, 1], 'card', [2, 2], 'val', [0.59, 0, 0.22, 0]);
FACTORS.EVIDENCE(3) = struct('var', [3, 2], 'card', [2, 2], 'val', [0, 0.61, 0, 0]);

% Compute Joint Distribution
% FACTORS.JOINT = ComputeJointDistribution(FACTORS.INPUT);
FACTORS.JOINT = struct('var', [1, 2, 3], 'card', [2, 2, 2], 'val', [0.025311, 0.076362, 0.002706, 0.041652, 0.039589, 0.119438, 0.042394, 0.652548]);

% Compute Marginal
%FACTORS.MARGINAL = ComputeMarginal([2, 3], FACTORS.INPUT, [1, 2]);
FACTORS.MARGINAL = struct('var', [2, 3], 'card', [2, 2], 'val', [0.0858, 0.0468, 0.1342, 0.7332]);

%p = ComputeJointDistribution(FACTORS.INPUT)

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [g,tfr]=psech(L,p2,p3,p4)
%PSECH  Sampled, periodized hyperbolic secant
%   Usage: g=psech(L);
%          g=psech(L,tfr);
%          g=psech(L,s,'samples);
%          [g,tfr]=psech( ... );
%
%   Input parameters:
%      L   : Length of vector.
%      tfr : ratio between time and frequency support.
%   Output parameters:
%      g   : The periodized hyperbolic cosine.
%
%   `psech(L,tfr)` computes samples of a periodized hyperbolic secant.
%   The function returns a regular sampling of the periodization
%   of the function $sech(\pi\cdot x)
%
%   The returned function has norm equal to 1.
%
%   The parameter *tfr* determines the ratio between the effective support
%   of *g* and the effective support of the DFT of *g*. If $tfr>1$ then *g*
%   has a wider support than the DFT of *g*.
%
%   `psech(L)` does the same setting $tfr=1$.
%
%   `psech(L,s,'samples')` returns a hyperbolic secant with an effective
%   support of *s* samples. This means that approx. 96% of the energy or 74%
%   or the area under the graph is contained within *s* samples. This is
%   equivalent to `psech(L,s^2/L)`.
%
%   `[g,tfr] = psech( ... )` additionally returns the time-to-frequency
%   support ratio. This is useful if you did not specify it (i.e. used
%   the `'samples'` input format).
%
%   The function is whole-point even.  This implies that `fft(psech(L,tfr))`
%   is real for any *L* and *tfr*.
%
%   If this function is used to generate a window for a Gabor frame, then
%   the window giving the smallest frame bound ratio is generated by
%   `psech(L,a*M/L)`.
%
%   Examples:
%   ---------
%
%   This example creates a `psech` function, and demonstrates that it is
%   its own Discrete Fourier Transform:::
%
%     g=psech(128);
%
%     % Test of DFT invariance: Should be close to zero.
%     norm(g-dft(g))
%
%   The next plot shows the `psech` in the time domain compared to the Gaussian:::
%
%     plot((1:128)',fftshift(pgauss(128)),...
%          (1:128)',fftshift(psech(128)));
%     legend('pgauss','psech');
%
%   The next plot shows the `psech` in the frequency domain on a log
%   scale compared to the Gaussian:::
%
%     hold all;
%     magresp(pgauss(128),'dynrange',100);
%     magresp(psech(128),'dynrange',100);
%     legend('pgauss','psech');
%
%   The next plot shows `psech` in the time-frequency plane:::
%
%     sgram(psech(128),'tc','nf','lin');
%
%   See also:  pgauss, pbspline, pherm
%
%   References: jast02-1

complainif_argnonotinrange(nargin,1,4,mfilename);

if nargin==1
  tfr=1;
end;

if size(L,1)>1 || size(L,2)>1
  error('L must be a scalar');
end;

if rem(L,1)~=0
  error('L must be an integer.')
end;

switch(nargin)
 case 1
  tfr=1;
  cent=0;
 case 2
  tfr=p2;
  cent=0;
 case 3
  if ischar(p3)
    switch(lower(p3))
     case {'s','samples'}
      tfr=p2^2/L;
     otherwise
      error('Unknown argument %s',p3);
    end;
    cent=0;
  else
    tfr=p2;
    cent=p3;
  end;
 case 4
  tfr=p2^2/L;
  cent=p4;
end;

safe=12;

g=zeros(L,1);
sqrtl=sqrt(L);

w=tfr;

% Outside the interval [-safe,safe] then sech(pi*x) is numerically zero.
nk=ceil(safe/sqrt(L/sqrt(w)));

lr=(0:L-1).';
for k=-nk:nk
  g=g+sech(pi*(lr/sqrtl-k*sqrtl)/sqrt(w));
end;

% Normalize it.
g=g*sqrt(pi/(2*sqrt(L*w)));

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function ListFacesInHullROI = MakeGeodesicOuterROI (mesh_outer, iV, facesInROI);

FacesByLines=[];
FacesByLines=mesh_outer.faces(facesInROI,:);

[index,value]=find(FacesByLines==iV);
index=index';
t=FacesByLines(index,:);

more off
l=unique(t(:));
l2=l;
oldsize=1;
while size(l,1)>oldsize
        for r=1:size(l2,1)
        [index2,value2]=find(FacesByLines==l2(r));
        index2=index2';
        index=[index index2];
        t2=FacesByLines(index2,:);
        t=[t; t2];
        end
     t=unique(t,'rows');
     oldsize=size(l);
     oldl=l;
     l=unique(t(:));
     l2=setdiff(l,oldl);
     indexListHull=unique(index);
end

ListFacesInHullROI=t;

if size(facesInROI,2) ~= size(indexListHull,2)
  disp([' non geodesic regions removed from the patch '])
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


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

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

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the optimal C and sigma
%               learning parameters found using the cross validation set.
%               You can use svmPredict to predict the labels on the cross
%               validation set. For example,
%                   predictions = svmPredict(model, Xval);
%               will return the predictions on the cross validation set.
%
%  Note: You can compute the prediction error using
%        mean(double(predictions ~= yval))
%
wid=[0.01;0.03;0.1;0.3;1;3;10;30];
errors=0;
for i=1:size(wid,1)
	for j=1:size(wid,1)
		model= svmTrain(X, y, wid(i), @(x1, x2) gaussianKernel(x1, x2, wid(j)));
		predictions=svmPredict(model, Xval);
		errors(i,j)=mean(double(predictions ~= yval));
	end
end
error=min(min(errors));
[i,j]=find(errors==error);
C=wid(i);
sigma=wid(j);

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

end

author: loserChen

repo: Coursera-MachineLearning


text:
% Calculate Pearson cross correlation between two firing rate maps
%
% The size of return variable corrValue depends on input parameter 'output':
% *) 'single', corrValue is a single number. Maps are reshaped in a vector and
%     correlation is calculated.
% *) 'vector', corrValue is a vector. Maps are processed row- or column-wise.
%    corrValue(i) is a correlation coeeficient for row/column i.
% Row/column selection is done by parameter 'processBy'.
%
%  USAGE
%   corrValue = analyses.spatialCrossCorrelation(map1, map2, <options>)
%   map1      2D matrix, rate map
%   map2      2D matrix, rate map
%   <options>   optional list of property-value pairs (see table below)
%
%   =========================================================================
%    Properties    Values
%   -------------------------------------------------------------------------
%   'output'       'single' or 'vector'. Default is 'single'.
%
%   'processBy'     'r' or 'c'. The parameter specifies how maps are processed.
%                   'r' means iteration over y-axis (rows) and calculation of correlation per row.
%                   'c' means iteration over x-axis (columns) and calculation of correlation
%                   per column. Default is 'c'. This only makes difference for 'vector' output.
%                   Single value is independent of processing order.
%   =========================================================================
%
%   corrValue       Correlation value(s). Can be either a single number or a vector.
%
function corrValue = spatialCrossCorrelation(map1, map2, varargin)
    columnWise = 0;
    rowWise = 1;
    single = 'single';

    inp = inputParser;
    defaultProcessBy = 'c';
    defaultReturnFormat = single;

    checkProcessBy = @(x) strcmpi(x, 'c') || strcmpi(x, 'r');
    checkReturnFormat = @(x) strcmpi(x, 'single') || strcmpi(x, 'vector');

    addRequired(inp, 'map1');
    addRequired(inp, 'map2');
    addParamValue(inp, 'processBy', defaultProcessBy, checkProcessBy);
    addParamValue(inp, 'output', defaultReturnFormat, checkReturnFormat);

    parse(inp, map1, map2, varargin{:});

    % get parsed arguments
    if strcmpi(inp.Results.processBy, 'c')
        orientation = columnWise;
    else
        orientation = rowWise;
    end
    returnFormat = inp.Results.output;

    if ~isequal(size(map1), size(map2))
        warning('BNT:mapsSize', 'Maps are of different size. This is suspicious. Correlation won''t be calculated');
        corrValue = nan;
        return;
    end

    if strcmpi(returnFormat, single)
        % Transform the 2-D maps to 1-D arrays by assembling the columns/rows from the maps
        A = reshape(map1, numel(map1), 1);
        B = reshape(map2, numel(map2), 1);

        % Remove the bins with NaN in A from both maps
        nansA = isnan(A);
        A = A(~nansA);
        B = B(~nansA);

        % Remove the bins with NaN in B from both maps
        nansB = isnan(B);
        A = A(~nansB);
        B = B(~nansB);

        % Calculate the correlation for the two rate maps.
        if isempty(A) || isempty(B)
            corrValue = NaN; % Cannot do calculation on empty arrays
        else
            corrValue = corr(A, B);
        end
    else
        nansLeft = isnan(map1);
        map1(nansLeft) = [];
        map2(nansLeft) = [];

        nansRight = isnan(map2);
        map1(nansRight) = [];
        map2(nansRight) = [];

        if orientation == columnWise
            corrMatrix = corr(map1, map2);
            corrValue = diag(corrMatrix)'; % make a row
        else
            corrMatrix = corr(map1', map2');
            corrValue = diag(corrMatrix); % make a column
        end
    end
end

author: kavli-ntnu


text:
function [ev, i] = pne_detected_event(event_log, name, zero)
%PNE_DETECTED_EVENT  Returns the detected event of a particular name
%   EV = PNE_DETECTED_EVENT(EVENT_LOG, NAME)
%   EV = PNE_DETECTED_EVENT(EVENT_LOG, NAME, ZERO)
%   [EV, I] = PNE_DETECTED_EVENT(...)
%
%   Searches through the log of detected events, returning the event with
%   the specified name, if one exists, and optionally only if it is refers
%   to a detected event zero (not an event interval).
%
%   Inputs:
%       EVENT_LOG : struct array of detected events
%       NAME : name of event of interest
%       ZERO : (optional, default = 0)
%           0 - return the event if a zero or an interval are detected for the
%               corresponding event function
%           1 - return the event only if a zero is deteced for the
%               corresponding event function
%
%   Outputs:
%       EV : entry of the EVENT_LOG corresponding to the detected event with
%            the specified NAME, or empty if NAME not in the EVENT_LOG of
%            detected events
%       I : index of EV in EVENT_LOG
%
%   See PNE_DETECT_EVENTS for details on the fields in the EVENT_LOG and EV
%   structs.

%   MP-Opt-Model
%   Copyright (c) 2016-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Shrirang Abhyankar, Argonne National Laboratory
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

if nargin < 3
    zero = 0;
end
ev = [];
i = 0;
for k = 1:length(event_log)
    if strcmp(event_log(k).name, name) && (event_log(k).zero || ~zero)
        ev = event_log(k);
        i = k;
    end
    break
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function result = x2stereo(x)
%X2STEREO Stereographic projection of Euclidean points
%
% result = x2stereo(x);

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

dim = size(x,1)-1;
mask = (x(dim+1,:) == 1);
scale = ones(dim,1)*(1-x(dim+1,~mask));
result = x(1:dim,~mask)./scale;
% end function

repo: gspbox


text:
function [dP] = GradJacobiP(r, alpha, beta, N);

% function [dP] = GradJacobiP(r, alpha, beta, N);
% Purpose: Evaluate the derivative of the Jacobi polynomial of type (alpha,beta)>-1,
%	       at points r for order N and returns dP[1:length(r))]

dP = zeros(length(r), 1);
if(N == 0)
  dP(:,:) = 0.0;
else
  dP = sqrt(N*(N+alpha+beta+1))*JacobiP(r(:),alpha+1,beta+1, N-1);
end;
return

repo: Aero-matlab


text:
% This function do a cross validation test of DNN.
%
% Author: Xiong Xiao, NTU
% Date Created: 10 Oct 2013
% Last Modified: 20 Jul 2015
%
function [cost,cost_pure,subcost,subacc] = CrossValidationTest_tree4(layer, data, para)
para.NET.sentenceMinibatch = 1;
[sentIdxInBlock] = shuffle_data(layer, para, data, 1);  nBlock = length(sentIdxInBlock);
if isfield(para.NET, 'nSequencePerMinibatchCV')
    para.NET.nSequencePerMinibatch = para.NET.nSequencePerMinibatchCV;
end

cost_cv = [];
if para.useGPU; 	cost_cv= gpuArray(cost_cv);        end
cost_cv_pure = cost_cv; subcost = cost_cv;    subacc = cost_cv;

nTestSample = [];
for blk_i = 1:nBlock
    pause(.1);
    if para.useGPU && para.displayGPUstatus==1
        gpuStatus = gpuDevice;
        fprintf('GPU has %2.2E free memory!\n', gpuStatus.FreeMemory);
    end

    [block_data] = BlockDataGeneration(data, para, sentIdxInBlock{blk_i});
    [minibatch] = MinibatchPackaging_tree4(block_data, para);

    clear cost_func
    for batch_i=1:minibatch.nBatch
        batch_data = GetMinibatch2(minibatch, para, batch_i);
        for si=1:length(batch_data)
            batch_nFr(si) = size(batch_data{si},2);
        end

        % Evaluate the cost function and gradient on current batch
        nTestSample(end+1) = max(batch_nFr);
        [cost_func(batch_i)] = DNN_Cost10(layer, batch_data, para, 2);

        if mod(batch_i,para.displayInterval)==0 || (batch_i==minibatch.nBatch && para.displayInterval>batch_i)
            fprintf('CV cost at utterance %d/%d = %.4g / %.4g - %s\n', ...
                batch_i, minibatch.nBatch, mean([cost_func(max(1,batch_i-para.displayInterval+1) : batch_i).cost]), ...
                mean([cost_func(max(1,batch_i-para.displayInterval+1) : batch_i).cost_pure]), datestr(now));
        end
    end
    cost_cv = [cost_cv [cost_func(:).cost]];
    cost_cv_pure = [cost_cv_pure [cost_func(:).cost_pure]];
    subcost = [subcost [cost_func(:).subcost]];
    subacc = [subacc [cost_func(:).subacc]];
    clear minibatch;
    layer = clean_network_layer(layer);
end

subcost =  gather(subcost * nTestSample') / sum(nTestSample);
cost = gather(cost_cv(:)' * nTestSample') / sum(nTestSample);
cost_pure = gather(cost_cv_pure(:)' * nTestSample') / sum(nTestSample);
if isfield(cost_func(1),'subacc')
    subacc =  gather(subacc * nTestSample') / sum(nTestSample);
else
    subacc = [];
end

fprintf('    ----Development cost/cost_pure before training = %f / %f - %s\n', cost, cost_pure, datestr(now));
if length(subcost)>1
    fprintf('        ----subcosts = %f\n', subcost);
end
if length(subacc)>0
    fprintf('    ----(Sub)task accuracy = %2.2f%%\n', 100*subacc);
end

end

author: singaxiong

repo: SignalGraph


text:
function [t1, t2, alpha, TransMat, Displ] = estimate_rigid_from_displ(T, mask)
    % t1 - translation in first dimension
    % t2 - translation in second dimension
    % alpha - angle (in radians). Rotation is performed around (0,0) point.
    %           keep in mind, that upper left pixel has index (1,1).
    % TransMat - transformation matrix for homogeneous coordinates
    % Displ - displacement field for estimated rigid transformation
    %
    % T - displacement field to be approximated
    %
    % Example usage
    % alpha = 0.2;
    % t1=5; t2=-9;
    % TT = displ_from_matrix_2d([cos(alpha), -sin(alpha), t1; sin(alpha), cos(alpha), t2; 0, 0, 1], [256,256], [], true);
    % [et1, et2, ealpha, eTransMat, eDispl] = estimate_rigid_from_displ(TT);
    % (1 - norm(fl(TT-Displ))/norm(fl(TT))) * 100 % -- percent explained by a rigid rotation
    if isempty(mask)
        mask = ones(size(T, 1), size(T, 2));
    end
    mask = logical(mask);

    sz = [size(T, 1), size(T, 2)];
    [n1, n2] = ndgrid(1:sz(1), 1:sz(2));
    T(:,:,1) = T(:,:,1) + n1;
    T(:,:,2) = T(:,:,2) + n2;
    ND = cat(3, n1, n2);

    T = reshape(T, [size(T,1)*size(T, 2), 2]);
    ND = reshape(ND, [size(ND,1)*size(ND, 2), 2]);
    T = T(mask, :)';
    ND = ND(mask, :)';
    T = [T; ones(1, size(T, 2))];
    ND = [ND; ones(1, size(ND, 2))];

    opts = [];
    opts.display = 'off';
    objf = @(x) fgval(x(1), x(2), x(3), T, ND);
%     fl = @(x) x(:);
%     x0 = [mean(fl(T(:,:,1))); mean(fl(T(:,:,2))); 0];
    x0 = [0;0;0];
    xmin = minFunc(objf, x0, opts);
    t1 = xmin(1);
    t2 = xmin(2);
    alpha = xmin(3) * pi / 180;
    TransMat = [cos(alpha), -sin(alpha), t1; ...
                sin(alpha), cos(alpha), t2; ...
                0, 0, 1];
    Displ = displ_from_matrix_2d(TransMat, sz, [], true);
end

function f = fval(t1, t2, alpha, T, ND)
    alpha = alpha * pi / 180; % just to keep t1,t2 and alpha in the same scale

%     pts = reshape(T, [size(T,1)*size(T,2), 2]);
%     pts = [pts, ones(size(pts,1), 1)]';
%     pts0 = reshape(ND, [size(T,1)*size(T,2), 2]);
%     pts0 = [pts0, ones(size(pts0,1), 1)]';

    ca = cos(alpha);
    sa = sin(alpha);
    R = [ca, -sa, t1; sa, ca, t2; 0, 0, 1];
%     ptsm = R * pts0;
%     f = sum(sqrt( sum( mask.*(pts - ptsm).^2, 1) ));
    ptsm = R * ND;
    f = sum(sqrt( sum((T - ptsm).^2, 1) ));
end

function [f, g] = fgval(t1, t2, alpha, T, ND)
    dx = 0.001;
    dx = 1e-5;

    f = fval(t1,t2,alpha,T, ND);
    fp1 = fval(t1 + dx, t2, alpha, T, ND);
    fp2 = fval(t1, t2 + dx, alpha, T, ND);
    fp3 = fval(t1, t2, alpha + dx, T, ND);
    g = [(fp1 - f) / dx; (fp2 - f) / dx; (fp3 - f) / dx];
end

repo: pTVreg


text:
% The COBRAToolbox: testGetDefaultCompartmentSymbols.m
%
% Purpose:
%     - test the getDefaultCompartmentSymbols function
%
% Authors:
%     - Jacek Wachowiak

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('getDefaultCompartmentSymbols'));
cd(fileDir);

% test variables
testCompSymbolList = {'c','m','v','x','e','t','g','r','n','p','l','y'};
testCompNameList = {'Cytoplasm','Mitochondrion','Vacuole','Peroxisome','Extracellular','Pool','Golgi','Endoplasmic_reticulum','Nucleus','Periplasm','Lysosome','Glycosome'};
% function outputs
[defaultCompartmentSymbolList, defaultCompartmentNameList] = getDefaultCompartmentSymbols();

% test
assert(isequal(defaultCompartmentSymbolList, testCompSymbolList))
assert(isequal(defaultCompartmentNameList, testCompNameList))

% change to old directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function [sp_new, scores_new] = change_color_space(color_space, I_rgb, sp, K, im_num, h, w, sph, opts)
% Updates the histograms in the 'sp' variable by recalculating them using
% image transformed into another color space, and color space -specific
% features. This function is not used with default settings.

scores_new = [];

switch color_space
    case 'nrgb'
        I_nrgb = uint8(rgb_to_nrgb(I_rgb));
        [words, k] = compute_features(I_nrgb, I_rgb, 'nrgb', opts, im_num);
        sp_new = compute_histograms(sp, words, k, h, w);
        scores_new = similarity_scores(sp_new, K, opts, 0, sph);
    case 'opp'
        I_opp = uint8(rgb_to_opp(I_rgb));
        [words, k] = compute_features(I_opp, I_rgb, 'opp', opts, im_num);
        sp_new = compute_histograms(sp, words, k, h, w);
        scores_new = similarity_scores(sp_new, K, opts, 0, sph);
    case 'hsv'
        I_hsv = uint8(rgb_to_hsv(I_rgb));
        [words, k] = compute_features(I_hsv, I_rgb, 'hsv', opts, im_num);
        sp_new = compute_histograms(sp, words, k, h, w);
        scores_new = similarity_scores(sp_new, K, opts, 0, sph);
    otherwise
        error('bad image type argument')

end

author: Cloud-CV

repo: object-proposals


text:
clc
clear
close all
addpath(genpath(pwd))

%  basic plotting
x = linspace(-0.1*pi,2*pi, 30);
y = cell(1, 3);
y{1, 1} = 0.4*sinc(x)+0.8;
y{1, 2} = tanh(x);
y{1, 3} = exp(-sinc(x));

figure;
color_ = [0, 114, 189; 126, 47, 142; 162, 20, 47]/255;
ax = axes('Units', 'normalized');
hold(ax, 'on');
box(ax,'on');
set(ax, 'LineWidth', 1.2, 'TickDir', 'in');
for i = 1:3
    plot(x, y{1, i}, 'Parent', ax, 'Color', color_(i, :), 'LineWidth', 3)
end

% add a zoomed zone
zp = BaseZoom();
zp.plot;

author: iqiukp

repo: ZoomPlot-MATLAB


text:
function [fx,dF_dX,dF_dTheta] = f_Lorenz(Xt,Theta,ut,inF)
% Lorenz chaotic evolution function

deltat = inF.deltat;

x       = Xt;
rho     = Theta(1);
sigma   = Theta(2);
beta    = Theta(3);

fx      = zeros(3,1);
fx(1)   = sigma.*(x(2)-x(1));
fx(2)   = x(1).*(rho-x(3))-x(2);
fx(3)   = x(1).*x(2) - beta.*x(3);
fx      = deltat.*fx + x;

J       = zeros(3,3);
J(1,:)  = [-sigma,sigma,0];
J(2,:)  = [rho-x(3),-1,-x(1)];
J(3,:)  = [x(2),x(1),-beta];
dF_dX   = deltat.*J' + eye(3);

dF_dTheta = zeros(3,3);
dF_dTheta(1,:)      = [0,x(1),0];
dF_dTheta(2,:)      = [x(2)-x(1),0,0];
dF_dTheta(3,:)      = [0,0,-x(3)];
dF_dTheta           = deltat.*dF_dTheta;

author: MBB-team

repo: VBA-toolbox


text:
function  Result = Metric(I,V,X)
I=double(I);
V=double(V);
X=double(X);
grey_level=256;
Result.Total = [];
Result.EN=entropy_fusion(X,grey_level);
Result.Total = [Result.Total; Result.EN];

Result.MI=mutural_information(I,V,X,grey_level);
Result.Total = [Result.Total; Result.MI];

Result.Q_G=Qp_ABF(I,V,X);
Result.Total = [Result.Total; Result.Q_G];

%Result.FMI = fmi(I,V,X, 'none', 3);  % it is time consuming
%Result.Total = [Result.Total; Result.FMI];

author: thfylsty


text:
function [xpad, nup, n1, n2] = padsignal(x, padtype)

[nup, n1, n2] = p2up(length(x));

xl = padarray(x(:), n1, padtype, 'pre');
xr = padarray(x(:), n2, padtype, 'post');

xpad = [xl(1:n1); x(:); xr(end-n2+1:end)];

author: HeLiangHIT


text:
function neg_base_exponent(base)
%
%
%
global C
C = 0;

figure;
set(gcf,'position',[350 400 850 350])

subplot(1,2,1);
set(gcf,'windowButtonMotionFcn',@fun_mouseMove)
line([-3 3],[0 0],'color','k','linewidth',3)

function fun_mouseMove (object, eventdata)

C = get(gca, 'CurrentPoint');
title(gca, ['(X,Y) = (', num2str(C(1,1)), ', ',num2str(C(1,2)), ')']);

subplot(1,2,2);
y = base.^C(1,1);
plot(real(y), imag(y),'ro');

xlim([-5 5])
ylim([-5 5])

end

end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function out = dimCheck(f, g, ~)
%DIMCHECK   Check dimension compatability of two CHEBFUN objects.
%   In MATLAB 2016a and below DIMCHECK(F, G) returns:
%       1    if numColumns(F) == numColumns(G)
%     error  otherwise.
%
%   In MATLAB 2016b and above DIMCHECK(F, G) returns:
%       1    if numColumns(F) == numColumns(G)
%      -1    if numColumns(F) == 1 or numColumns(G) == 1
%     error  otherwise.
%
%   DIMCHECK(F, G, 1) will return 0 rather than throw a dimension error.
%
% See also NUMCOLUMNS, SIZE.

% Copyright 2015 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

nf = numColumns(f);
ng = numColumns(g);

% Check dimensions:
if ( nf == ng )
    out = 1;
elseif ( nf == 1 || ( ng == 1 ) )
    out = -1;
else
    out = 0;
end

% Adjust for MATLAB version if out = -1:
if ( (out == -1 ) && verLessThan('matlab', '9.1') )
   out = 0;
end

if ( out == 0 && nargin == 2 )
    % Throw an error if requested (as the caller):
    ME = MException('CHEBFUN:CHEBFUN:dimCheck:dim', ...
        'Matrix dimensions must agree.');
    throwAsCaller(ME)
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function display(SO3F,varargin)
% standard output

if check_option(varargin,'skipHeader')
  disp('  function handle component');
else
  displayClass(SO3F,inputname(1),[],'moreInfo',symChar(SO3F),varargin{:});
  if numel(SO3F) > 1, disp(['  size: ' size2str(SO3F)]); end
end

if SO3F.antipodal, disp('  antipodal: true'); end

% The computation may be to slow , For example:
% F = SO3FunHandle(@(rot) calcTaylor(inv(orientation(rot,cs))*eps,sS.symmetrise)
% try
%   m = mean(SO3F,'resolution',5*degree);
%   disp(['  weight: ' xnum2str(m(1))])
% end
disp(' ')

end

author: mtex-toolbox

repo: mtex

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


text:
function [old_dirs,old_stps,Hdiag,Bcompact] = lbfgsUpdate(y,s,corrections,debug,old_dirs,old_stps,Hdiag)

%B0 = eye(length(y))/Hdiag;
S = old_dirs(:,2:end);
Y = old_stps(:,2:end);
k = size(Y,2);
L = zeros(k);
for j = 1:k
    for i = j+1:k
        L(i,j) = S(:,i)'*Y(:,j);
    end
end
D = diag(diag(S'*Y));
N = [S/Hdiag Y];
M = [S'*S/Hdiag L;L' -D];

ys = y'*s;
Bs = s/Hdiag - N*(M\(N'*s)); % Product B*s
sBs = s'*Bs;

eta = .02;
if ys < eta*sBs
    if debug
        fprintf('Damped Update\n');
    end
    theta = min(max(0,((1-eta)*sBs)/(sBs - ys)),1);
    y = theta*y + (1-theta)*Bs;
end

numCorrections = size(old_dirs,2);
if numCorrections < corrections
    % Full Update
    old_dirs(:,numCorrections+1) = s;
    old_stps(:,numCorrections+1) = y;
else
    % Limited-Memory Update
    old_dirs = [old_dirs(:,2:corrections) s];
    old_stps = [old_stps(:,2:corrections) y];
end

% Update scale of initial Hessian approximation
Hdiag = (y'*s)/(y'*y);

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
addpath data/

patch_dim = 64;
num_patches = 10000;
listing = dir('data/denoise/mit_saliency/*.jpg');

num_file = 20;
for m = 1 : num_file
    fprintf('Extracting patch batch: %d / %d\n', m, num_file);
    % extract random patches
    samples = zeros(patch_dim, patch_dim, 3, num_patches);
    labels = zeros(size(samples));
    for i = 1 : num_patches
        if (mod(i,100) == 0)
            fprintf('Extracting patch: %d / %d\n', i, num_patches);
        end

        r_idx = random('unid', size(listing, 1));
        orig_img = imread(strcat('data/denoise/mit_saliency/', listing(r_idx).name));

        orig_img_size = size(orig_img);
        offset = 30;
        r = random('unid', orig_img_size(1) - patch_dim - 2*offset + 1) + offset;
        c = random('unid', orig_img_size(2) - patch_dim - 2*offset + 1) + offset;

        input_patch = im2double(orig_img(r:r+patch_dim-1, c:c+patch_dim-1, :));

        samples(:,:,:,i) = imnoise(input_patch, 'gaussian', 0, 0.01);
        labels(:,:,:,i) = input_patch;
    end
    samples = single(samples);
    labels = single(labels);
    % save it
    filename = strcat('data\denoise\train\patches_', num2str(m));
    save(filename, '-v7.3', 'samples', 'labels');
end

author: jimmy-ren


text:
function Y = sl2dpca_apply(Mm, PL, PR, data, matsiz, n)
%SL2DPCA_APPLY Applies 2D PCA onto a set of matrices to extract features
%
% $ Syntax $
%   - Y = sl2dpca_apply(Mm, PL, PR, data, matsiz, n)
%
% $ Description $
%   - Mm:           the mean matrix
%   - PL:           the left projection matrix
%   - PR:           the right projection matrix
%   - data:         the matrix samples or the cell array of filenames
%   - matsiz:       the original matrix size
%   - n:            the number of samples
%   - Y:            the extracted 2D features
%
% $ Description $
%   - Y = sl2dpca_apply(data, Mm, PL, PR) extracts 2D features for
%     the matrices given in data, in either a 3D array or a set of
%     array filenames. Suppose the original matrix size is d1 x d2,
%     PL be d1 x k1, PR be d2 x k2, then the feature matrix would be
%     of size k1 x k2. Y is a k1 x k2 x n array.
%
% $ History $
%   - Created by Dahua Lin, on Jul 31st, 2006
%

%% Parse and verify input arguments

if nargin < 6
    raise_lackinput('sl2dpca_apply', 6);
end

matsiz = matsiz(:)';
if length(matsiz) ~= 2
    error('sltoolbox:invalidarg', ...
        'matsiz should be a 2-elem vector');
end

if ~isequal(size(Mm), matsiz)
    error('sltoolbox:sizmismatch', ...
        'the sample size does not match the model');
end
d1 = matsiz(1);
d2 = matsiz(2);

if size(PL, 1) ~= d1 || size(PR, 1) ~= d2
    error('sltoolbox:sizmismatch', ...
        'the size of projection matrices are illegal');
end

%% Compute

if isnumeric(data)

    if size(data, 3) ~= n
        error('sltoolbox:sizmismatch', ...
            'The number of samples is not as specified');
    end

    Y = computeY(data, Mm, PL, PR);

elseif iscell(data)

    Y = zeros(size(PL, 2), size(PR, 2), n);

    nfiles = length(data);
    cf = 0;
    for i = 1 : nfiles
        curdata = slreadarray(data{i});
        curn = size(curdata, 3);
        Y(:,:,cf+1:cf+curn) = computeY(curdata, Mm, PL, PR);
        cf = cf + curn;
    end

else
    error('sltoolbox:invalidarg', ...
        'data should be a numeric array or a cell array of filenames');

end

%% Core compute function

function Y = computeY(X, Mm, PL, PR)

n = size(X, 3);
Y = zeros(size(PL, 2), size(PR, 2), n);
PLT = PL';

for i = 1 : n
    Y(:,:,i) = PLT * (X(:,:,i) - Mm) * PR;
end

author: lmthang


text:
function [COEF] = rmder_inv(P, Vec, Mat)

% rmder_inv returns a vector containing the relative standard deviation of the regression coefficients. The
%                            constant of the regression is not included
%
%
%      Input:
%             P             Property vector
%             Vec           Descriptors vector
%             Mat           Descriptors matrix
%
%
%     Returns:
%
%            COEF           Vector containing the standard deviation for
%                            the regression coefficients without the constant
%
%
%
%
% Andrew G. Mercader, Pablo R. Duchowicz
% INIFTA, La Plata, Argentina
% Created: 30 Jan 2007

if (nargin < 3)
   error('rmder requires at least 3 input variables. Type ''help rmder''.');
end

% Exctracts from descriptor matrix the matrix corresponding to the given
% set of descriptors
% ----------------------------------------------------------------------------
    X = Mat(:,Vec);
% Check that independent (X) and dependent (P) data have compatible dimensions
% ----------------------------------------------------------------------------
     [n_x, k] = size(X);
[n_y,columns] = size(P);

if n_x ~= n_y,
    error('The number of rows in P must equal the number of rows in Mat.');
end

if columns ~= 1,
    error('P must be a vector, not a matrix');
end

n = n_x;

%  Solve for the regression coefficients using ordinary least-squares
%  ------------------------------------------------------------------
     X = [ ones(n,1) X ]   ;

XTXI = X' * X;
if rcond(XTXI) < 1e-25
    S_err = 0.00000000000000000000000000000000001;
    return;
end

 XTXI = inv(XTXI);
 Coefficients = XTXI * X' * P ;

%  Calculate the fitted regression values
%  --------------------------------------

    P_hat = X * Coefficients;

%  Calculate residuals and standard error
%  --------------------------------------

    residuals = P - P_hat;

           S_err = sqrt(residuals' * residuals / (n - k - 1) );

%  Calculate the standard deviation  for the regression coefficients
%  -----------------------------------------------------------------------------

    covariance = XTXI .* S_err^2;

    C = sqrt(diag(covariance, 0));
    Er = abs(C./Coefficients)*100;

COEF=[[Er]'];
COEF(1)=[];

author: Sable

repo: mcbench-benchmarks


text:
function [bndry_lon,bndry_lat] = m_plotbndry(name,varargin)
% M_PLOTBNDRY plots text files of Lat,Lon for political boundaries.
% Text files (derived from the DCW) are obtained from
%   http://www.maproom.psu.edu/cgi-bin/ian/points/index.cgi
%
%     M_PLOTBNDRY(NAME) plots the state or country specified in the
%     string NAME, which may include a path name. The routine will
%     then a) search the specified path for a mat-file of that name
%          b) search the specified path for an ascii *2pts.txt file
%             of that name, and, if found convert it to a mat-file.
%          c) failing that, it will open a file dialog box.
%
%     M_PLOTBNDRY(NAME, ...line properties) will use the specified
%     line properties in drawing the boundary.
%
%     [LON,LAT]=M_PLOTBNDRY(...) returns vectors of the boundary
%     points.
%
% Note: If errors occur when a file is first plotted, check that
% the entire file was downloaded.  It should end with two consecutive
% END lines.

% Original Author: Michael W. Mann
%
% Changes: R. Pawlowicz 21/12/98 - changed interface to read from
%          given directory, allow output, allow various line
%          properties to be specified.
% 6/Nov/00 - eliminate returned stuff if ';' neglected (thx to D Byrne)
% 19/Mar/04 - .mat files not being created because of a bug found by James Connor.

% Set current projection to geographic
Currentmap=m_coord('set');
m_coord('geographic');

targetfile = [name,'2pts.mat']; %try to find binary file
targetpath = dir(targetfile);

if ( length(targetpath)~=0 ) %then load file.
   load(targetfile);
   m_line(bndry_lon,bndry_lat,'tag','m_plotbndry',varargin{:});

else %Can't find binary file, load and process text file.
   %  Create .mat file from data read.

   targetfile = [name,'2pts.txt'];
   targetpath = dir(targetfile);
   if ( length(targetpath)==0 ) %then
      [filename,pathname]=uigetfile('*2pts.txt',['Select ',targetfile,' file']);
      if ( filename == 0 ), error(['Could not find ',targetfile]), end
      targetpath = [pathname,filename];
   end %if
   [fid,message] = fopen(targetfile,'r');
   if ( fid < 0 ) %then
      fclose(fid);
      error(message);
   end %if

   namein = fgetl(fid);
   if ( ~strfind(lower(name),lower(namein)) ) %then
      fclose(fid);
      error(['File contains wrong state! ','Desired: ',name,' Found: ',namein]);
   end %if

   done = 0;
   polynum = 0;
   numpts = [];
   while ( ~done )
      %read in polygon index
      line = fgetl(fid);
      field = sscanf(line,'%3c');
      if ( strcmp(upper(field),'END') ) %then found EOF
         done = 1;
      else %start reading in polygon
         polynum = polynum + 1;
         polynumin = sscanf(line,'%u');
         if ( polynum ~= polynumin ) %then
            fclose(fid);
            error(['Didn''t find polynum ',num2str(polynum),...
                  ', found ',num2str(polynumin),' .']);
         end %if
         line = fgetl(fid); %skip polygon centroid
         line = fgetl(fid); %first point
         ptcount = 0;
         while ( ~strcmp(upper(line(1:3)),'END') )
            ptcount = ptcount + 1;
            line = fgetl(fid);
         end %while
         numpts = [numpts,ptcount];
      end %if
   end %while

   Npolygon = length(numpts);
   fclose(fid); %Close file.

   % Vectors for composite outline
   bndry_lat = zeros(1,sum(numpts)+Npolygon);
   bndry_lon = zeros(1,sum(numpts)+Npolygon);
   i1 = 0; %index to last filled entry in vectors

   fid = fopen(targetfile,'r'); %reopen file
   line = fgetl(fid); %ignore namefield

   for i = 1:Npolygon
      line = fgetl(fid); %ignore polygon index
      line = fgetl(fid); %ignore centroid lat,lon
      [poly,count] = fscanf(fid,'%g',[2,inf]); %read data
      line = fgetl(fid); %ignore END

      % Save data for reuse as .mat file
      i0 = i1 + 1;  i1 = i0 + numpts(i);
      bndry_lon(i0:i1) = [poly(1,:),NaN]; %add NaN to lift pen
      bndry_lat(i0:i1) = [poly(2,:),NaN]; %add NaN to lift pen

   end %for

   fclose(fid);

   m_line(bndry_lon,bndry_lat,'tag','m_plotbndry',varargin{:});
   nchar = length(targetfile);   % Bug fix thanks to James Connor
   matfile = [targetfile(1:(nchar-4)),'.mat'];
   save(matfile,'bndry_lat','bndry_lon');

end %if

m_coord(Currentmap.name);

if nargout==0
 clear bndry_lon bndry_lat
end

author: goGPS-Project


text:
function varargout = circleArcAsCurve(arc, N)
%CIRCLEARCASCURVE Convert a circle arc into a series of points
%
%   deprecated: replaced by circleArcToPolyline
%
%
% ---------
% author : David Legland
% created the 22/05/2006.
% Copyright 2010 INRA - Cepia Software Platform.
%

% HISTORY
% 2011-03-30 use angles in degrees, add default value for N
% 2011-12-09 deprecate

warning('matGeom:deprecated', ...
    'function "circleArcAsCurve" is deprecated, use "circleArcToPolygon" instead');

% format output
if nargout <= 1
    varargout = {circleArcToPolyline(arc, N)};
else
    [x, y] = circleArcToPolyline(arc, N);
    varargout = {x, y};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text: %GE Datafile overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function matrix = genMatEll3D(pde,mesh,fem)
%% Generate global matrices and load vector of FEM for 3D ellipic eq
%     -div(A grad u)  = f,    x\in \Omega
% INPUTS:
% pde --- given data function from equation, e.g.
%         pde.A --- diffusion coefficient
%         pde.f --- right hand side function
%         pde.gD --- Dirichlet boundary value function
%         pde.one --- constant function 1.
% mesh --- mesh structure.
% fem --- global degree of freedom of FEM
%
% OUTPUTS:
% matrix.S --- stiffness matrix (w/o boundary condition)
% matrix.A --- final FEM matrix (after boundary condition)
% matrix.rhsF --- load vector (w/o boundary condition)
% matrix.f --- final RHS matrix (after boundary condition)

% Last Modified: 08/07/2020 by Xu Zhang

%% 1. Stiffness Matrix
S = globMatrix3DStiff(pde.A,mesh,fem,fem);

%% 2. Generate the Right Hand Side Vector
rhsF = globRHS3D(pde.f, mesh, fem, [0,0,0]);

%% 3. Dirichlet Boundary Conditions
Atotal = S;
tu = feval(pde.gD,fem.p(:,1),fem.p(:,2),fem.p(:,3));
ub = tu;
ub(fem.mapper) = 0;
rhsB = Atotal*ub;
A = Atotal(fem.mapper,fem.mapper);
f = rhsF(fem.mapper) - rhsB(fem.mapper);

%% Outputs
matrix = struct('A', A, 'f', f, 'S', S, 'rhsF',rhsF);

repo: ifem


text:
function [data,units] = compute_wing_arear_mm(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = trx(fly).wing_arear ./ trx.pxpermm^2;

end
units = parseunits('mm^2');

author: kristinbranson

repo: JAABA


text:
function D = spm_eeg_inv_mesh_ui(varargin)
% Cortical Mesh user interface
% FORMAT D = spm_eeg_inv_mesh_ui(D, val, sMRI, Msize)
%
% D        - input data struct (optional)
% val      -
% sMRI     -  0 - use template (default), or string with image file name
% Msize    -
%
% D        - same data struct including the meshing files and variables
%__________________________________________________________________________
%
% Invokes spatial normalisation (if required) and the computation of
% the individual mesh.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jeremie Mattout & Christophe Phillips
% $Id: spm_eeg_inv_mesh_ui.m 4027 2010-07-31 12:49:19Z vladimir $

SVNrev = '$Rev: 4027 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','Define head model');

%-Initialisation
%--------------------------------------------------------------------------
[D,val] = spm_eeg_inv_check(varargin{:});

if val == 0
    val = 1;
end

if ~isfield(D, 'inv') || ~isfield(D.inv{val}, 'comment')
    D.inv = {struct('mesh', [])};
    D.inv{val}.date    = strvcat(date,datestr(now,15));
    D.inv{val}.comment = {''};
else
    inv = struct('mesh', []);
    inv.comment = D.inv{val}.comment;
    inv.date    = D.inv{val}.date;
    D.inv{val} = inv;
end

if nargin > 2
    sMRI = varargin{3};
else
    sMRI = [];
end

if isempty(sMRI)
    template = spm_input('Select head  model', '+1','template|individual', [1 0]);
elseif ~ischar(sMRI)
    template = sMRI; % for backward compatibility
else
    template = 0;
end

if template
    sMRI = [];
elseif ~ischar(sMRI)
    % get sMRI file name
    sMRI = spm_select([0 1],'image','Select subject''s structural MRI (Press Done if none)');
    if isempty(sMRI)
        error('No structural MRI selected.');
    end
end

if nargin>3
    Msize = varargin{4};
else
    Msize = spm_input('Cortical mesh', '+1', 'coarse|normal|fine', [1 2 3]);
end

D.inv{val}.mesh = spm_eeg_inv_mesh(sMRI, Msize);

%-Check meshes and display
%--------------------------------------------------------------------------
spm_eeg_inv_checkmeshes(D);

%-Cleanup
%--------------------------------------------------------------------------
spm('FigName','Define head model: done');

author: spm


text:
function y = nanvar(x,w,dim)
%NANVAR Variance, ignoring NaNs.
%   Y = NANVAR(X) returns the sample variance of the values in X, treating
%   NaNs as missing values.  For a vector input, Y is the variance of the
%   non-NaN elements of X.  For a matrix input, Y is a row vector
%   containing the variance of the non-NaN elements in each column of X.
%   For N-D arrays, NANVAR operates along the first non-singleton dimension
%   of X.
%
%   NANVAR normalizes Y by N-1 if N>1, where N is the sample size of the
%   non-NaN elements.  This is an unbiased estimator of the variance of the
%   population from which X is drawn, as long as X consists of independent,
%   identically distributed samples, and data are missing at random.  For
%   N=1, Y is normalized by N.
%
%   Y = NANVAR(X,1) normalizes by N and produces the second moment of the
%   sample about its mean.  NANVAR(X,0) is the same as NANVAR(X).
%
%   Y = NANVAR(X,W) computes the variance using the weight vector W.  The
%   length of W must equal the length of the dimension over which NANVAR
%   operates, and its non-NaN elements must be nonnegative.  Elements of X
%   corresponding to NaN elements of W are ignored.
%
%   Y = NANVAR(X,W,DIM) takes the variance along dimension DIM of X.
%
%   See also VAR, NANSTD, NANMEAN, NANMEDIAN, NANMIN, NANMAX, NANSUM.

%   Copyright 1984-2010 The MathWorks, Inc.
%   $Revision: 1.1.8.2 $  $Date: 2010/10/08 17:25:19 $

if nargin < 2 || isempty(w), w = 0; end

sz = size(x);
if nargin < 3 || isempty(dim)
    % The output size for [] is a special case when DIM is not given.
    if isequal(x,[]), y = NaN(class(x)); return; end

    % Figure out which dimension sum will work along.
    dim = find(sz ~= 1, 1);
    if isempty(dim), dim = 1; end
elseif dim > length(sz)
    sz(end+1:dim) = 1;
end

% Need to tile the mean of X to center it.
tile = ones(size(sz));
tile(dim) = sz(dim);

if isequal(w,0) || isequal(w,1)
    % Count up non-NaNs.
    n = sum(~isnan(x),dim);

    if w == 0
        % The unbiased estimator: divide by (n-1).  Can't do this when
        % n == 0 or 1, so n==1 => we'll return zeros
        denom = max(n-1, 1);
    else
        % The biased estimator: divide by n.
        denom = n; % n==1 => we'll return zeros
    end
    denom(n==0) = NaN; % Make all NaNs return NaN, without a divideByZero warning

    x0 = x - repmat(nanmean(x, dim), tile);
    y = nansum(abs(x0).^2, dim) ./ denom; % abs guarantees a real result

% Weighted variance
elseif numel(w) ~= sz(dim)
    error(message('stats:nanvar:InvalidSizeWgts'));
elseif ~(isvector(w) && all(w(~isnan(w)) >= 0))
    error(message('stats:nanvar:InvalidWgts'));
else
    % Embed W in the right number of dims.  Then replicate it out along the
    % non-working dims to match X's size.
    wresize = ones(size(sz)); wresize(dim) = sz(dim);
    wtile = sz; wtile(dim) = 1;
    w = repmat(reshape(w, wresize), wtile);

    % Count up non-NaNs.
    n = nansum(~isnan(x).*w,dim);

    x0 = x - repmat(nansum(w.*x, dim) ./ n, tile);
    y = nansum(w .* abs(x0).^2, dim) ./ n; % abs guarantees a real result
end


text:
function outputImage = minip(this, applicationDimension)
% Computes minimum intensity projection image along specified dimension,
% uses Matlab min function
%
%
%   Y = MrImage()
%   Y.minip(applicationDimension)
%
% This is a method of class MrImage.
%
% IN
%   applicationDimension    image dimension along which operation is
%                           performed (e.g. 4 = time, 3 = slices)
%                           default: The last dimension with more than one
%                           value is chosen
%                           (i.e. 3 for 3D image, 4 for 4D image)
%
% OUT
%   outputImage
%
% EXAMPLE
%   minip
%
%   See also MrImage MrImage.perform_unary_operation min MrImage.maxip

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

if nargin < 2
    applicationDimension = this.dimInfo.nDims;
else
    applicationDimension = this.dimInfo.convert_application_dimensions(...
        applicationDimension);
end

outputImage = this.perform_unary_operation(@(x) min(x,[],applicationDimension));

author: translationalneuromodeling

repo: tapas


text:
function x = bradford_cdf_inv ( cdf, a, b, c )

%*****************************************************************************80
%
%% BRADFORD_CDF_INV inverts the Bradford CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Input, real A, B, C, the parameters of the PDF.
%    A < B,
%    0.0 < C.
%
%    Output, real X, the corresponding argument of the CDF.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BRADFORD_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'BRADFORD_CDF_INV - Fatal error!' );
  end

  if ( cdf <= 0.0 )
    x = a;
  elseif ( cdf < 1.0 )
    x = a + ( b - a ) * ( ( c + 1.0 )^cdf - 1.0 ) / c;
  elseif ( 1.0 <= cdf )
    x = b;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%{
load('dataset/trafficdb/traffic_patches.mat');
[M,m,n,p] = convert_video3d_to_2d(im2double(imgdb{100}));
out = run_algorithm('MC', 'FPC', M, [])
show_results(M.*out.Omega,out.L,out.S,out.O,p,m,n);
%}

maxiter = 500;
mu_final = .01; tol = 1e-6;
MIdx = M(Idx);

fprintf('\nSolving by FPC...\n');
[U,S,V,numiter] = FPC(size(M),Idx,MIdx,mu_final,maxiter,tol);

L = U*S*V'; % low-rank
S = M - L; % sparse

author: andrewssobral

repo: lrslibrary