text:
% pop_chanplot() - graphic user interface (GUI)-based function with plotting
%                options for visualizing. Only channel measures (e.g., spectra,
%                ERPs, ERSPs, ITCs) that have been computed and saved in the study EEG
%                datasets can be visualized. These can be computed using the GUI-based
%                pop_precomp().
% Usage:
%                >> STUDY = pop_chanplot(STUDY, ALLEEG);
% Inputs:
%   ALLEEG     - Top-level EEGLAB vector of loaded EEG structures for the dataset(s)
%                in the STUDY. ALLEEG for a STUDY set is typically loaded using
%                pop_loadstudy(), or in creating a new STUDY, using pop_createstudy().
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG
%   datasets in ALLEEG.
%
% Outputs:
%   STUDY      - The input STUDY set structure modified according to specified user edits,
%                if any. Plotted channel measure means (maps, ERSPs, etc.) are added to
%                the STUDY structure after they are first plotted to allow quick replotting.
%
% Graphic interface buttons:
%  "Select channel to plot" - [list box] Displays available channels to plot (format is
%                'channel name (number of channels)'). The presented channels depend s
%                on the optional input variable 'channels'. Selecting (clicking on) a
%                channel from the list will display the selected channel channels in the
%                "Select channel(s) to plot" list box. Use the plotting buttons below
%                to plot selected measures of the selected channel. Additional editing
%                options (renaming the channel, rejecting outliers, moving channels to
%                another channel) are also available. The option 'All N channel centroids'
%                at the top of the list displays all the channels in the list except the
%                'Notchannel', 'Outlier' and 'Parentchannel' channels. Selecting this option
%                will plot the channel centroids (i.e. ERP, ERSP, ...) in a single figure.
%  "Select channel(s) to plot" - [list box] Displays the ICA channels of the currently
%                selected channel (in the "Select channel to plot" list box). Each channel
%                has the format: 'subject name, channel index'. Multiple channels can be
%                selected from the list. Use the plotting buttons below to plot different
%                measures of the selected channels on different figures. Selecting the
%                "all subjects" option is  equivalent to using the channel plotting buttons.
%                Additional editing options are reassigning the selected channels to
%                another channel or moving them to the outlier channel.
%  "Plot channel properties" - [button] Displays in one figure all the mean channel measures
%                (e.g., dipole locations, scalp maps, spectra, etc.) that were calculated
%                and saved in the EEG datsets. If there is more than one condition, the ERP
%                and the spectrum will have different colors for each condition. The ERSP
%                and ITC plots will show only the first condition; clicking on the subplot
%                will open a new figure with the different conditions displayed together.
%                Uses the command line function std_propplot().
%  "Plot scalp maps"  - [button] Displays the scalp maps of channel channels.
%                If applied to a channel, scalp maps of the channel channels
%                are plotted along with the channel mean scalp map in one figure.
%                If "All # channel centroids" option is selected, all channel scalp map
%                means are plotted in the same figure. If applied to channels, displays
%                the scalp maps of the specified channel channels in separate figures.
%                Uses the command line functions std_plotmap() and std_plotchanmap().
%  "Plot ERSPs" - [button] Displays the channel channel ERSPs.
%                If applied to a channel, channel ERSPs are plotted in one figure
%                (per condition) with the channel mean ERSP. If "All # channel centroids"
%                option is selected, plots all average ERSPs of the channels in one figure
%                per condition. If applied to channels, display the ERSP images of specified
%                channel channels in separate figures, using one figure for all conditions.
%                Uses the command line functions std_plotersp() and std_plotchannelsp().
%  "Plot ITCs" - [button] Same as  "Plot ERSPs" but with ITC.
%                Uses the command line functions std_plotitc() and std_plotchanitc().
%  "Plot dipoles" - [button] Displays the dipoles of the channel channels.
%                If applied to a channel, plots the channel channel dipoles (in blue)
%                plus the average channel dipole (in red). If "All # channel centroids" option
%                is selected, all channel plots are displayed in one figure each channel in
%                a separate subplot. If applied to channels, displays the ERSP images of the
%                specified channel. For specific channels displays channels dipole (in blue)
%                plus the average channel dipole (in Red) in separate figures.
%                Uses the command line functions std_dipplot() and std_plotchandip().
%  "Plot spectra" - [button] Displays the channel channel spectra.
%                If applied to a channel, displays channel spectra plus the average channel
%                spectrum in bold. For a specific channel, displays the channel channel
%                spectra plus the average channel spectrum (in bold) in one figure per condition.
%                If the "All # channel centroids" option is selected, displays the average
%                spectrum of all channels in the same figure, with spectrum for different
%                conditions (if any) plotted in different colors.
%                If applied to channels, displays the spectrum of specified channel
%                channels in separate figures using one figure for all conditions.
%                Uses the command line functions std_plotspec() and std_plotonechanpec().
%  "Plot ERPs" - [button] Same as "Plot spectra" but for ERPs.
%                Uses the command line functions std_ploterp() and std_plotchannelp().
%  "Create new channel" - [button] Creates a new empty channel.
%                Opens a popup window in which a name for the new channel can be entered.
%                If no name is given the default name is 'Cls #', where '#' is the next
%                available channel number. For changes to take place, press the popup
%                window 'OK' button, else press the 'Cancel' button. After the empty
%                channel is created, channels can be moved into it using,
%                'Reassign selected channel(s)' (see below). Uses the command line
%                function std_createchant().
%  "Rename selected channel" - [button] Renames a channel using the selected (mnemonic) name.
%                Opens a popup window in which a new name for the selected channel can be
%                entered. For changes to take place, press the popup window 'OK' button,
%                else press the 'Cancel' button. Uses the command line function std_renamechant().
%  "Reject outlier channels" - [button] rejects outlier channels to an outlier channel.
%                Opens a popup window to specify the outlier threshold. Move outlier
%                channels that are more than x standard deviations devs from the
%                channel centroid to an outlier channel. For changes to take place,
%                press the popup window 'OK' button, else press the 'Cancel' button.
%                Uses the command line function std_rejectoutliers().
%  "Merge channels" - [button] Merges several channels into one channel.
%                Opens a popup window in which the channels to merge may be specified
%                An optional name can be given to the merged channel. If no name is given,
%                the default name is 'Cls #', where '#' is the next available channel number.
%                For changes to take place, press the popup window 'OK' button, else press
%                the 'Cancel' button. Uses the command line function std_mergechant().
%  "Remove selected outlier channel(s)" - [button] Moves selected channel(s) to the
%                outlier channel. The channels that will be moved are the ones selected
%                in the "Select channel(s) to plot" list box. Opens a popup window in which
%                a list of the selected channel(s) is presented. For changes to take place,
%                press the popup window 'OK' button, else press the 'Cancel' button.
%                Uses the command line function std_moveoutlier().
%  "Reassign selected channel(s)" - [button] Moves selected channel(s) from one channel
%                to another. The channels that will reassign are the ones selected in the
%                "Select channel(s) to plot" list box. Opens a popup window in which
%                a list of possible channels to which to move the selected channel(s) is
%                presented. For changes to take place, press the popup window 'OK' button,
%                else press the 'Cancel' button. Uses the command line function std_movecomp().
%  "Save STUDY set to disk" - [check box] Saves the STUDY set structure modified according
%                to specified user edits to the disk. If no file name is entered will
%                overwrite the current STUDY set file.
%
% See also:  pop_prechant(), pop_chant().
%
% Authors: Arnaud Delorme, Hilit Serby, Scott Makeig, SCCN/INC/UCSD, October 11, 2004

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, October 11, 2004, arno@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% Coding notes: Useful information on functions and global variables used.

function [STUDY, com] = pop_chanplot(varargin)

icadefs;
com = [];
if ~isstr(varargin{1})
    if nargin < 2
        error('pop_chanplot(): You must provide ALLEEG and STUDY structures');
    end
    STUDY  = varargin{1};
    STUDY.etc.erpparams.topotime  = []; % [] for channels and NaN for components
    STUDY.etc.specparams.topofreq = []; % NaN -> GUI disabled
    STUDY.etc.erspparams.topotime = [];
    STUDY.etc.erspparams.topofreq = [];

    % test path
    % ---------
    pathwarn = 'off';
    if ~strcmpi(pwd, STUDY.filepath)
        if length(STUDY.datasetinfo(1).filepath) < 1
            pathwarn = 'on';
        elseif STUDY.datasetinfo(1).filepath(1) == '.'
            pathwarn = 'on';
        end;
    end;
    if strcmpi(pathwarn, 'on')
        warndlg2(strvcat('You have changed your working path and data files are', ...
                         'no longer available; Cancel, and go back to your STUDY folder'), 'warning');
    end;

    oldhistory = STUDY.history;
    STUDY.history = '';
    ALLEEG = varargin{2};
    if ~isfield(STUDY, 'changrp')
        STUDY = std_changroup(STUDY, ALLEEG);
        disp('Warning: history not saved for group creation');
    elseif isempty(STUDY.changrp)
        STUDY = std_changroup(STUDY, ALLEEG);
        disp('Warning: history not saved for group creation');
    end;

    show_chan          = ['pop_chanplot(''showchan'',gcf);'];
    show_onechan       = ['pop_chanplot(''showchanlist'',gcf);'];
	plot_chan_maps     = ['pop_chanplot(''topoplot'',gcf); '];
    plot_onechan_maps  = ['pop_chanplot(''plotchantopo'',gcf); '];
    plot_chan_ersps    = ['pop_chanplot(''erspplot'',gcf); '];
    plot_onechan_ersps = ['pop_chanplot(''plotchanersp'',gcf); '];
    plot_chan_itcs     = ['pop_chanplot(''itcplot'',gcf); '];
    plot_onechan_itcs  = ['pop_chanplot(''plotchanitc'',gcf); '];
    plot_chan_spectra  = ['pop_chanplot(''specplot'',gcf); '];
    plot_onechan_spectra = ['pop_chanplot(''plotchanspec'',gcf); '];
    plot_chan_erp      = ['pop_chanplot(''erpplot'',gcf); '];
    plot_onechan_erp   = ['pop_chanplot(''plotchanerp'',gcf); '];
    plot_chan_dip      = ['pop_chanplot(''dipplot'',gcf); '];
    plot_onechan_dip   = ['pop_chanplot(''plotchandip'',gcf); '];
    plot_chan_sum      = ['pop_chanplot(''plotsum'',gcf); '];
    plot_onechan_sum   = ['pop_chanplot(''plotonechanum'',gcf); '];
    rename_chan        = ['pop_chanplot(''renamechan'',gcf);'];
    move_onechan       = ['pop_chanplot(''movecomp'',gcf);'];
    move_outlier       = ['pop_chanplot(''moveoutlier'',gcf);'];
    create_chan        = ['pop_chanplot(''createchan'',gcf);'];
    reject_outliers    = ['pop_chanplot(''rejectoutliers'',gcf);'];
    merge_channels     = ['pop_chanplot(''mergechannels'',gcf);'];
    erp_opt            = ['pop_chanplot(''erp_opt'',gcf);'];
    spec_opt           = ['pop_chanplot(''spec_opt'',gcf);'];
    ersp_opt           = ['pop_chanplot(''ersp_opt'',gcf);'];
    create_group       = ['pop_chanplot(''create_group'',gcf);'];
    edit_group         = ['pop_chanplot(''edit_group'',gcf);'];
    delete_group       = ['pop_chanplot(''delete_group'',gcf);'];
    saveSTUDY          = [ 'set(findobj(''parent'', gcbf, ''userdata'', ''save''), ''enable'', fastif(get(gcbo, ''value'')==1, ''on'', ''off''));' ];
    browsesave         = [ '[filename, filepath] = uiputfile2(''*.study'', ''Save STUDY with .study extension -- pop_chan()''); ' ...
                           'set(faindobj(''parent'', gcbf, ''tag'', ''studyfile''), ''string'', [filepath filename]);' ];
    sel_all_chans      = ['pop_chanplot(''sel_all_chans'',gcf);'];

    % list of channel groups
    % ----------------------
    show_options = {};
    for index = 1:length(STUDY.changrp)
        show_options{end+1} = [ 'All ' STUDY.changrp(index).name ];
    end;

    % enable buttons
    % --------------
    filename = STUDY.design(STUDY.currentdesign).cell(1).filebase;
    if exist([filename '.datspec']) , spec_enable = 'on'; else spec_enable  = 'off'; end;
    if exist([filename '.daterp'] )  , erp_enable = 'on'; else  erp_enable  = 'off'; end;
    if exist([filename '.datersp']) , ersp_enable = 'on'; else ersp_enable  = 'off'; end;
    if exist([filename '.datitc'])  ,  itc_enable = 'on'; else  itc_enable  = 'off'; end;

    if isfield(ALLEEG(1).dipfit, 'model'), dip_enable   = 'on'; else dip_enable   = 'off'; end;

    % userdata below
    % --------------
    fig_arg{1}{1} = ALLEEG;
    fig_arg{1}{2} = STUDY;
    fig_arg{1}{3} = STUDY.changrp;
    fig_arg{1}{4} = { STUDY.changrp.name };
    fig_arg{2}    = length(STUDY.changrp);

    geometry = { [4] [1] [0.7 0.3 0.3 0.1 0.9] [1 0.3 1] [1 0.3 1] [1 0.3 1] [1 0.3 1] ...
                 [1 0.3 1] [1 0.3 1] };
    str_name       = sprintf('STUDY name ''%s'' - ''%s''', STUDY.name, STUDY.design(STUDY.currentdesign).name);
    if length(str_name) > 80, str_name = [ str_name(1:80) '...''' ]; end;

    uilist   = { ...
        {'style' 'text' 'string' str_name 'FontWeight' 'Bold' 'HorizontalAlignment' 'center'} {} ...
        {'style' 'text'       'string' 'Select channel to plot' 'FontWeight' 'Bold' } ...
        {'style' 'pushbutton' 'string' 'Sel. all' 'callback' sel_all_chans } {} {} ...
        {'style' 'text'       'string' 'Select subject(s) to plot' 'FontWeight' 'Bold'} ...
        {'style' 'listbox'    'string' show_options 'value' 1 'max' 2 'tag' 'chan_list' 'Callback' show_chan } {} ...
        {'style' 'listbox'    'string' '' 'tag' 'chan_onechan' 'max' 2 'min' 1 'callback'    show_onechan } ...
        {'style' 'pushbutton' 'enable'   erp_enable 'string' 'Plot ERPs' 'Callback' plot_chan_erp} ...
        {'style' 'pushbutton' 'enable'   erp_enable 'string' 'Params' 'Callback' erp_opt }  ...
        {'style' 'pushbutton' 'enable'   erp_enable 'string' 'Plot ERP(s)' 'Callback' plot_onechan_erp} ...
        {'style' 'pushbutton' 'enable'  spec_enable 'string' 'Plot spectra' 'Callback' plot_chan_spectra} ...
        {'style' 'pushbutton' 'enable'  spec_enable 'string' 'Params' 'Callback' spec_opt }  ...
        {'style' 'pushbutton' 'enable'  spec_enable 'string' 'Plot spectra' 'Callback' plot_onechan_spectra} ...
        {'style' 'pushbutton' 'enable'  ersp_enable 'string' 'Plot ERSPs' 'Callback' plot_chan_ersps} ...
        {'vertshift' 'style' 'pushbutton' 'enable'  ersp_enable 'string' 'Params' 'Callback' ersp_opt }  ...
        {'style' 'pushbutton' 'enable'  ersp_enable 'string' 'Plot ERSP(s)' 'Callback' plot_onechan_ersps}...
        {'style' 'pushbutton' 'enable'   itc_enable 'string' 'Plot ITCs' 'Callback' plot_chan_itcs} { }  ...
        {'style' 'pushbutton' 'enable'   itc_enable 'string' 'Plot ITC(s)' 'Callback' plot_onechan_itcs}...
        {'style' 'pushbutton' 'string' 'Plot channel properties' 'Callback' plot_chan_sum} {} ...
        {'style' 'pushbutton' 'string' 'Plot channel properties (soon)' 'Callback' plot_onechan_sum 'enable' 'off'} };

   [out_param userdat] = inputgui( 'geometry' , geometry, 'uilist', uilist, ...
                                   'helpcom', 'pophelp(''pop_chanplot'')', ...
                                   'title', 'View and edit current channels -- pop_chanplot()' , 'userdata', fig_arg, ...
                                   'geomvert', [ 1 0.5 1 5 1 1 1 1 1], 'eval', show_chan );

   if ~isempty(userdat)
       ALLEEG = userdat{1}{1};
       STUDY  = userdat{1}{2};
   end

   % history
   % -------
   com = STUDY.history;
   STUDY.history =  sprintf('%s%s', oldhistory, com);

else
    hdl = varargin{2};  %figure handle
    userdat  = get(varargin{2}, 'userdat');
    ALLEEG   = userdat{1}{1};
    STUDY    = userdat{1}{2};
    cls      = userdat{1}{3};
    allchans = userdat{1}{4};

    changrp = get(findobj('parent', hdl, 'tag', 'chan_list')   , 'value');
    onechan = get(findobj('parent', hdl, 'tag', 'chan_onechan'), 'value');

    switch  varargin{1}

        case {'topoplot', 'erspplot','itcplot','specplot', 'erpplot'}
            changrpstr = allchans(changrp);
            plotting_option = varargin{1};
            plotting_option = [ plotting_option(1:end-4) 'plot' ];
            a = ['STUDY = std_' plotting_option '(STUDY,ALLEEG,''channels'','  vararg2str({changrpstr}) ');' ];
             % update Study history
            eval(a); STUDY.history =  sprintf('%s\n%s',  STUDY.history, a);
            userdat{1}{2} = STUDY;
            set(hdl, 'userdat',userdat);

        case {'plotchantopo', 'plotchanersp','plotchanitc','plotchanspec', 'plotchanerp','plotchandip'}
            changrpstr    = allchans(changrp);
            %if length(changrp) > 1
            %    subject = STUDY.subject{onechan-1};
            %else
            %    changrpstruct = STUDY.changrp(changrp);
            %    allsubjects   = unique({ STUDY.datasetinfo([ changrpstruct.setinds{:} ]).subject });
            %    subject = allsubjects{onechan-1};
            %end;
            plotting_option = varargin{1};
            plotting_option = [ plotting_option(9:end) 'plot' ];
            if onechan(1) ~= 1  % check that not all onechan in channel are requested
                 subject = STUDY.design(STUDY.currentdesign).cases.value{onechan-1};
                 a = ['STUDY = std_' plotting_option '(STUDY,ALLEEG,''channels'','  vararg2str({changrpstr}) ', ''subject'', ''' subject ''' );' ];
                 eval(a); STUDY.history =  sprintf('%s\n%s',  STUDY.history, a);
             else
                a = ['STUDY = std_' plotting_option '(STUDY,ALLEEG,''channels'','  vararg2str({changrpstr}) ', ''plotsubjects'', ''on'' );' ];
                eval(a); STUDY.history =  sprintf('%s\n%s',  STUDY.history, a);
             end;
            userdat{1}{2} = STUDY;
            set(hdl, 'userdat',userdat);

        case 'erp_opt' % save the list of selected channels
            [STUDY com] = pop_erpparams(STUDY);
            if ~isempty(com)
                STUDY.history =  sprintf('%s\n%s',  STUDY.history, com);
            end;
            userdat{1}{2} = STUDY;
            set(hdl, 'userdat',userdat); %update information (STUDY)

        case 'spec_opt' % save the list of selected channels
            [STUDY com] = pop_specparams(STUDY);
            if ~isempty(com)
                STUDY.history =  sprintf('%s\n%s',  STUDY.history, com);
            end;
            userdat{1}{2} = STUDY;
            set(hdl, 'userdat',userdat); %update information (STUDY)

        case 'ersp_opt' % save the list of selected channels
            [STUDY com] = pop_erspparams(STUDY);
            if ~isempty(com)
                STUDY.history =  sprintf('%s\n%s',  STUDY.history, com);
            end;
            userdat{1}{2} = STUDY;
            set(hdl, 'userdat',userdat); %update information (STUDY)

        case 'showchanlist' % save the list of selected channels
            if length(changrp) == 1
                STUDY.changrp(changrp).selected = onechan;
            end;
            userdat{1}{2} = STUDY;
            set(hdl, 'userdat',userdat); %update information (STUDY)

       case 'showchan'
            cind     = get(findobj('parent', hdl, 'tag', 'chan_list')   , 'value');
            changrp  = STUDY.changrp(cind);

            % Find datasets availaible
            % ------------------------
            %setind = STUDY.setind .* (changrp.chaninds > 0); % set to 0 the cell not
            %%                                       % containing any electrode
            %allchansets = unique( setind(find(setind(:))) );

            % Generate channel list
            % ---------------------
            chanid{1} = 'All subjects';
            if length(changrp) == 1
                allsubjects = unique({ STUDY.design(STUDY.currentdesign).cell([ changrp.setinds{:} ]).case });
                for l = 1:length(allsubjects)
                    chanid{end+1} = [ allsubjects{l} ' ' changrp.name ];
                end;
            else
                for l = 1:length(STUDY.design(STUDY.currentdesign).cases.value)
                    chanid{end+1} = [ STUDY.design(STUDY.currentdesign).cases.value{l} ];
                end;
            end;

            selected = 1;
            if isfield(changrp, 'selected') & length(cind) == 1
                if ~isempty(STUDY.changrp(cind).selected)
                    selected = min(STUDY.changrp(cind).selected, 1+length(chanid));
                    STUDY.changrp(cind).selected = selected;
                end;
            end;

            set(findobj('parent', hdl, 'tag', 'chan_onechan'), 'value', selected, 'String', chanid);

        case 'sel_all_chans'
            set(findobj('parent', hdl, 'tag', 'chan_list'), 'value', [1:length(STUDY.changrp)]);

            % Generate channel list
            % ---------------------
            chanid{1} = 'All subjects';
            for l = 1:length(STUDY.design(STUDY.currentdesign).cases.value)
                chanid{end+1} = [ STUDY.design(STUDY.currentdesign).cases.value{l} ' All' ];
            end;
            selected = 1;
            set(findobj('parent', hdl, 'tag', 'chan_onechan'), 'value', selected, 'String', chanid);

        case 'plotsum'
            changrpstr = allchans(changrp);
            [STUDY] = std_propplot(STUDY, ALLEEG, allchans(changrp));
            a = ['STUDY = std_propplot(STUDY, ALLEEG, ' vararg2str({ allchans(changrp) }) ' );'  ];
            STUDY.history =  sprintf('%s\n%s',  STUDY.history, a);
            userdat{1}{2} = STUDY;
            set(hdl, 'userdat',userdat);

        case 'create_group'
            channames = { STUDY.changrp(changrp).name };
            for i=1:length(channames), channames{i} = [ ' ' channames{i} ]; end;
            channamestr = strcat(channames{:});
            res = inputdlg2({ 'Name of channel group', 'Channels to group' }, 'Create channel group', 1, { '' channamestr(2:end) });
            if isempty(res), return; end;
            STUDY.changrp(end+1).name = res{1};
            allchans(end+1)         = { res{1} };
            chanlabels = parsetxt(res{2});
            if length(chanlabels) == 1
                warndlg2('Cannot create a channel group with a single channel');
                return;
            end;
            STUDY.changrp(end).channels = chanlabels;
            tmp = std_chanlookup( STUDY, ALLEEG, STUDY.changrp(end));
            STUDY.changrp(end).chaninds = tmp.chaninds;
            userdat{1}{2} = STUDY;
            userdat{1}{4} = allchans;
            set(hdl, 'userdat',userdat);

            % list of channel groups
            % ----------------------
            tmpobj  = findobj('parent', hdl, 'tag', 'chan_list');
            tmptext = get(tmpobj, 'string');
            tmptext{end+1} = [ 'All ' STUDY.changrp(end).name ];
            set(tmpobj, 'string', tmptext, 'value', length(tmptext));

        case 'edit_group'
            if length(changrp) > 1, return; end;
            if length(STUDY.changrp(changrp).channels) < 2, return; end;
            channames = STUDY.changrp(changrp).channels;
            for i=1:length(channames), channames{i} = [ ' ' channames{i} ]; end;
            channamestr = strcat(channames{:});
            res = inputdlg2({ 'Name of channel group', 'Channels to group' }, 'Create channel group', ...
                            1, { STUDY.changrp(changrp).name channamestr(2:end) });
            if isempty(res), return; end;
            STUDY.changrp(end+1).name = '';
            STUDY.changrp(changrp)    = STUDY.changrp(end);
            STUDY.changrp(end)        = [];
            STUDY.changrp(changrp).name = res{1};
            allchans(changrp)         = { res{1} };
            chanlabels = parsetxt(res{2});
            STUDY.changrp(changrp).channels = chanlabels;
            tmp = std_chanlookup( STUDY, ALLEEG, STUDY.changrp(end));
            STUDY.changrp(changrp).chaninds = tmp.chaninds;
            userdat{1}{2} = STUDY;
            userdat{1}{4} = allchans;
            set(hdl, 'userdat',userdat);

            % list of channel groups
            % ----------------------
            show_options = {};
            for index = 1:length(STUDY.changrp)
                show_options{end+1} = [ 'All ' STUDY.changrp(index).name ];
            end;
            tmpobj  = findobj('parent', hdl, 'tag', 'chan_list');
            set(tmpobj, 'string', show_options, 'value', changrp);

        case 'delete_group'
            if length(changrp) > 1, return; end;
            if length(STUDY.changrp(changrp).channels) < 2, return; end;
            STUDY.changrp(changrp)    = [];

            % list of channel groups
            % ----------------------
            show_options = {};
            for index = 1:length(STUDY.changrp)
                show_options{end+1} = [ 'All ' STUDY.changrp(index).name ];
            end;
            tmpobj  = findobj('parent', hdl, 'tag', 'chan_list');
            set(tmpobj, 'string', show_options, 'value', changrp-1);

        case 'renamechan'
            STUDY.saved = 'no';
            chan_name_list = get(findobj('parent', hdl, 'tag', 'chan_list'), 'String');
            chan_num = get(findobj('parent', hdl, 'tag', 'chan_list'), 'Value') -1;
            if chan_num == 0  % 'all subjects' option
                return;
            end
            % Don't rename 'Notchan' and 'Outliers'  channels.
            if strncmpi('Notchan',STUDY.channel(cls(chan_num)).name,8) | strncmpi('Outliers',STUDY.channel(cls(chan_num)).name,8) | ...
                    strncmpi('Parentchannel',STUDY.channel(cls(chan_num)).name,13)
                warndlg2('The Parentchannel, Outliers, and Notchan channels cannot be renamed');
                return;
			end
            old_name = STUDY.channel(cls(chan_num)).name;
            rename_param  = inputgui( { [1] [1] [1]}, ...
                { {'style' 'text' 'string' ['Rename ' old_name] 'FontWeight' 'Bold'} {'style' 'edit' 'string' '' 'tag' 'chan_rename' } {} }, ...
            '', 'Rename channel - from pop_chanplot()' );
            if ~isempty(rename_param) %if not canceled
                new_name = rename_param{1};
                STUDY = std_renamechan(STUDY, ALLEEG, cls(chan_num), new_name);
                % update Study history
                a = ['STUDY = std_renamechan(STUDY, ALLEEG, ' num2str(cls(chan_num)) ', ' STUDY.channel(cls(chan_num)).name  ');'];
                STUDY.history =  sprintf('%s\n%s',  STUDY.history, a);

                new_name = [ STUDY.channel(cls(chan_num)).name ' (' num2str(length(STUDY.channel(cls(chan_num)).onechan))  ' ICs)'];
                chan_name_list{chan_num+1} = renamechan( chan_name_list{chan_num+1}, new_name);
                set(findobj('parent', hdl, 'tag', 'chan_list'), 'String', chan_name_list);
                set(findobj('parent', hdl, 'tag', 'chan_rename'), 'String', '');
                userdat{1}{2} = STUDY;
                set(hdl, 'userdat',userdat); %update STUDY
            end
    end
end

function newname = renamechan(oldname, newname);

    tmpname = deblank(oldname(end:-1:1));
    strpos  = strfind(oldname, tmpname(end:-1:1));

    newname = [ oldname(1:strpos-1) newname ];

author: PatternRecognition

repo: OpenBMI


text:
function error_frobenius = r8mat_is_eigen_right ( n, k, a, x, lambda )

%*****************************************************************************80
%
%% R8MAT_IS_EIGEN_RIGHT determines the error in a right eigensystem.
%
%  Discussion:
%
%    An R8MAT is a matrix of real values.
%
%    This routine computes the Frobenius norm of
%
%      A * X - X * LAMBDA
%
%    where
%
%      A is an N by N matrix,
%      X is an N by K matrix (each of K columns is an eigenvector)
%      LAMBDA is a K by K diagonal matrix of eigenvalues.
%
%    This routine assumes that A, X and LAMBDA are all real.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer K, the number of eigenvectors.
%    K is usually 1 or N.
%
%    Input, real A(N,N), the matrix.
%
%    Input, real X(N,K), the K eigenvectors.
%
%    Input, real LAMBDA(K), the K eigenvalues.
%
%    Output, real ERROR_FROBENIUS, the Frobenius norm
%    of the difference matrix A * X - X * LAMBDA, which would be exactly zero
%    if X and LAMBDA were exact eigenvectors and eigenvalues of A.
%
  c(1:n,1:k) = a(1:n,1:n) * x(1:n,1:k);

  for j = 1 : k
    c(1:n,j) = c(1:n,j) - lambda(j) * x(1:n,j);
  end

  error_frobenius = r8mat_norm_fro ( n, k, c );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Driver script for solving the 2D Poisson equation on curvilinear domains
Globals2D;

% Polynomial order used for approximation
N = 6;

% Read in Mesh
[Nv, VX, VY, K, EToV] = MeshReaderGambitBC2D('circA01.neu');

% Initialize solver and construct grid and metric
StartUp2D;

% hard coded for all Dirichlet boundaries to be on cylinder
[k,f] = find(BCType==Dirichlet);
if(~isempty(k))
  cylfaces = [k,f];
  curved = sort(unique(k)); straight = setdiff(1:K, curved);
  MakeCylinder2D(cylfaces, 1, 0, 0);
end

% choose order to integrate exactly
Nint = ceil(2*N/2);

% build cubature nodes for all elements
CubatureOrder = 2*(Nint+1); cub = CubatureVolumeMesh2D(CubatureOrder);

% build Gauss node data for all element faces
NGauss = (Nint+1); gauss = GaussFaceMesh2D(NGauss);

% build weak Poisson operator matrices
[A, M] = CurvedPoissonIPDG2D();
Abc = CurvedPoissonIPDGbc2D();

% set up right hand side
f = (-2*pi^2-1)*sin(pi*x).*sin(pi*y);

% set up boundary condition
ubc = zeros(gauss.NGauss*Nfaces*K,1);
xbc = gauss.x(gauss.mapD); ybc = gauss.y(gauss.mapD);
ubc(gauss.mapD) = sin(pi*xbc).*sin(pi*ybc);
xbc = gauss.x(gauss.mapN); ybc = gauss.y(gauss.mapN);
ubc(gauss.mapN) = ...
    gauss.nx(gauss.mapN).*(pi*cos(pi*xbc).*sin(pi*ybc)) + ...
    gauss.ny(gauss.mapN).*(pi*sin(pi*xbc).*cos(pi*ybc));

ubc = Abc*ubc;

% solve linear system
solvec = (A+M)\(M*(-f(:)) + ubc);
u = reshape(solvec, Np, K);

repo: Aero-matlab


text:
function initParams(obj)
% INITPARAM  Initialize the paramers of the DagNN
%   OBJ.INITPARAM() uses the INIT() method of each layer to initialize
%   the corresponding parameters (usually randomly).

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

for l = 1:numel(obj.layers)
    p = obj.getParamIndex(obj.layers(l).params) ;
    %params = obj.layers(l).block.initParams() ;
    params=[];
    if(isequal(class(obj.layers(l).block),'dagnn.Conv'))
        sizet = obj.layers(l).block.size;
        h =sizet(1); w =sizet(2); in = sizet(3); out=sizet(4);
        init = obj.layers(l).block.init;
        if( isequal(init,'msra'))
            sc = sqrt(2/(h*w*out)) ;
            params{1,1}= randn(h, w, in, out, 'single')*sc ;
        else
            sc = sqrt(3/(h*w*in)) ;
            params{1,1} = (rand(h, w, in, out, 'single')*2 - 1)*sc;
        end

        if(obj.layers(l).block.hasBias)
            params{1,2} = zeros(out, 1, 'single') ;
        end
    elseif(isequal(class(obj.layers(l).block),'dagnn.ConvTranspose'))
        sizet=obj.layers(l).block.size;
        h =sizet(1); w =sizet(2); in = sizet(3); out=sizet(4);
        %sc = sqrt(3/(h*w*in)) ;
        %params{1,1} = (rand(h, w, in, out, 'single')*2 - 1)*sc;
        sc = sqrt(2/(h*w*in)) ;
        params{1,1}= randn(h, w, in, out, 'single')*sc ;
        %params{1,1} = bilinear_u(h,in,out);
        params{1,2} = zeros(in, 1, 'single') ;
    elseif(isequal(class(obj.layers(l).block),'dagnn.BatchNorm'))
        ss=obj.layers(l-1).block.size;
        %h =size(1); w =size(2); in = size(3); out=size(4);
        %ss = size(obj.vars(obj.layers(l-1).inputIndexes).value);
        in =  ss(4);
        if(in==2)
        end
        params{1,1} = ones(in,1, 'single');  %%zzd
        params{1,2} = zeros(in, 1, 'single') ;
        params{1,3} = zeros(in, 2, 'single') ;
    else
        params = obj.layers(l).block.initParams() ;
    end
    switch obj.device
        case 'cpu'
            params = cellfun(@gather, params, 'UniformOutput', false) ;
        case 'gpu'
            params = cellfun(@gpuArray, params, 'UniformOutput', false) ;
    end
    pppp = isequal(class(obj.layers(l).block),'dagnn.Conv') || ...
        isequal(class(obj.layers(l).block),'dagnn.BatchNorm') ||...
        isequal(class(obj.layers(l).block),'dagnn.ConvTranspose');
    if(pppp&&~isempty(obj.params(p(1)).value))  %for fintune
        continue;
    end
    [obj.params(p).value] = deal(params{:}) ;
    if(isequal(class(obj.layers(l).block),'dagnn.Conv'))
        [obj.params(p(1)).learningRate]=0.1;
        [obj.params(p(1)).trainMethod] = 'gradient';
        if(obj.layers(l).block.hasBias)
            [obj.params(p(2)).learningRate]=2;
            [obj.params(p(2)).trainMethod] = 'gradient';
        end
        if(l<numel(obj.layers) && ~isempty(strfind(class(obj.layers(l+1).block),'Loss')))
            [obj.params(p(1)).learningRate]= 0.01;
            if(obj.layers(l).block.hasBias)
                [obj.params(p(2)).learningRate]= 0.02;
            end
            obj.params(p(1)).value = obj.params(p(1)).value * 0.1;
        end
    end
    if(isequal(class(obj.layers(l).block),'dagnn.BatchNorm'))
        [obj.params(p(1)).learningRate]=0.2;
        [obj.params(p(2)).learningRate]=0.1;
        [obj.params(p(3)).learningRate]=0.05;
        [obj.params(p(1)).weightDecay]=1;
        [obj.params(p(2)).weightDecay]=1;
        [obj.params(p(3)).weightDecay]=1;
    end
end

author: layumi

repo: Image-Text-Embedding


text:
function o = niftistruc(fmt)
% Create a data structure describing NIFTI headers
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: niftistruc.m 7147 2017-08-03 14:07:01Z spm $

if ~nargin, fmt = 'nifti1'; end
switch lower(fmt)
    case {'nifti1','ni1','n+1'}
        o = nifti1struc;
    case {'nifti2','ni2','n+2'}
        o = nifti2struc;
    otherwise
        error('Unknown format.');
end

author: spm


text:
function y = normM(M)
y = sqrt(sum(M.^2, 2));
end

author: xiuyechen

repo: FishExplorer


text:
function [M] = spm_DEM_M(model,varargin)
% Create a [template] model structure
% FORMAT [M] = spm_DEM_M(model,l,n)
% FORMAT [M] = spm_DEM_M(model,X1,X2,...)
%
% model: 'General linear model','GLM'
%        'Factor analysis','FA'
%        'Independent component analysis','ICA'
%        'Sparse coding',"SC'
%        'convolution model'
%        'State space model','SSM',','Double Well'
%        'Lorenz'
%        'Ornstein_Uhlenbeck','OU'
%
% l(i) - number of outputs from level i
% n(i) - number of hidden states in level i
%
% Xi   - deisgn matrix for level i
%
%==========================================================================
% hierarchical generative model
%--------------------------------------------------------------------------
%   M(i).g  = y(t)  = g(x,v,P)    {inline function, string or m-file}
%   M(i).f  = dx/dt = f(x,v,P)    {inline function, string or m-file}
%
%   M(i).pE = prior expectation of p model-parameters
%   M(i).pC = prior covariances of p model-parameters
%   M(i).hE = prior expectation of h hyper-parameters (cause noise)
%   M(i).hC = prior covariances of h hyper-parameters (cause noise)
%   M(i).gE = prior expectation of g hyper-parameters (state noise)
%   M(i).gC = prior covariances of g hyper-parameters (state noise)
%   M(i).Q  = precision components (input noise)
%   M(i).R  = precision components (state noise)
%   M(i).V  = fixed precision (input noise)
%   M(i).W  = fixed precision (state noise)
%
%   M(i).m  = number of inputs v(i + 1);
%   M(i).n  = number of states x(i);
%   M(i).l  = number of output v(i);
%__________________________________________________________________________
% Copyright (C) 2007-2015 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_DEM_M.m 6416 2015-04-21 15:34:10Z guillaume $

switch lower(model)

    % hierarchical linear generative model - the basis for  static models
    % These subsume: Parametric emprical Bayes (PEB) models
    %                Mixed effects (MFX) models
    %======================================================================
    case{'hierarchical linear model','hlm','peb','mfx'}

        % Get design matrices for each level
        %------------------------------------------------------------------
        M(1).E.linear = 1;
        if ~isempty(varargin);
            pE   = varargin;
        else
            error('Please specify design matrices');

            global SPM_DATA;
            fig  = spm_data;
            set(fig,'Name','Please specify design matrices using ''next''')
            waitfor(fig)
            pE   = SPM_DATA;
        end

        for i = 1:length(pE)

            % level i
            %--------------------------------------------------------------
            M(i).m  = size(pE{i},2);
            M(i).l  = size(pE{i},1);
            M(i).g  = inline('P*v','x','v','P');
            M(i).pE = pE{i};
            M(i).hE = 0;
            M(i).hC = 16;

        end

        % non-hierarchical linear generative model (static)
        %==================================================================
    case{'general linear model','glm'}

        % Get design matrix
        %------------------------------------------------------------------
        if ~isempty(varargin);
            pE  = varargin{1};
        else
            global SPM_DATA
            fig = spm_data;
            set(fig,'Name','Please specify a design matrix')
            waitfor(fig)
            pE  = SPM_DATA;
        end

        % This is simply a one-level HLM with flat priors
        %------------------------------------------------------------------
        M       = spm_DEM_M('hlm',pE);

        % Factor analysis model (static) with unknown parameters
        %==================================================================
    case{'factor analysis','fa'}

        % Get orders (this routine will accommodate hierarchical FA)
        %------------------------------------------------------------------
        try
            l = varargin{1};
            l = l(~~l);
        catch
            errordlg('please specify number of inputs and ouputs')
        end

        % This is simply a HLM with unknown parameters
        %------------------------------------------------------------------
        g     = length(l) - 1;
        for i = 1:g
            pE{i} = sparse(l(i),l(i + 1));
            pC{i} = speye(l(i)*l(i + 1));
        end

        % create basic HLM and add prior uncertainty
        %------------------------------------------------------------------
        M     = spm_DEM_M('hlm',pE{:});
        for i = 1:g
            M(i).pC = pC{i};
        end

        % Principal component analysis (static - linear)
        % Equivalent to singular value decomposition (SVD)
        %==================================================================
    case{'principal component analysis','pca','svd'}

        % create basic FA
        %------------------------------------------------------------------
        try
            l    = varargin{1}(1);
        catch
            msgbox('please specify number of outputs')
            error(' ')
        end
        M        = spm_DEM_M('fa',[l l]);
        g        = length(M);

        % assume precisely small error at the first level
        %------------------------------------------------------------------
        M(1).hE  = 8;
        M(1).hC  = 1/32;

        % assume unit causes
        %------------------------------------------------------------------
        M(2).V   = speye(M(2).l,M(2).l);

        % Independent component analysis (static - nonlinear)
        %==================================================================
    case{'independent component analysis','ica'}

        % create basic FA
        %------------------------------------------------------------------
        M        = spm_DEM_M('pca',varargin{:});
        g        = length(M);

        % and add supra-ordinate [nonlinear] level
        %------------------------------------------------------------------
        M(g + 1) = M(g);
        M(g).m   = M(g + 1).l;
        M(g).g   = 'v-tanh(v)*P';
        M(g).pE  = 0;
        M(g).pC  = 0;
        M(g).V   = speye(M(2).l,M(2).l)*exp(8);

        % Sparse coding, pICA (static - nonlinear)
        %==================================================================
    case{'sparse coding','pica','sc'}

        % create basic ica
        %------------------------------------------------------------------
        M          = spm_DEM_M('ica',varargin{:});
        g          = length(M);

        % and add a noise component to the lower levels
        %------------------------------------------------------------------
        for i = 1:(g - 2)
            M(i).hE = 1;
            M(i).V  = sparse(M(i).l,M(i).l);
        end

        % Linear convolution (State-space model)
        %==================================================================
    case{'convolution model'}

        % time-step
        %------------------------------------------------------------------
        try
            dt = varargin{1};
        catch
            dt = 1;
        end

        % smoothness
        %------------------------------------------------------------------
        M(1).E.linear = 1;                          % linear model
        M(1).E.s      = 1/2;                          % smoothness

        % level 1
        %------------------------------------------------------------------
        pE.f    = [-1  4;                           % the Jacobian for the
                   -2 -1]/(4*dt);                   % hidden sates
        pE.g    = [spm_dctmtx(4,2)]/4;              % the mixing parameters
        pE.h    = [1 0]';                           % input parameter
        M(1).n  = 2;
        M(1).f  = inline('P.f*x + P.h*v','x','v','P');
        M(1).g  = inline('P.g*x','x','v','P');
        M(1).pE = pE;                               % prior expectation

        % level 2
        %------------------------------------------------------------------
        M(2).l  = 1;                                % inputs
        M(2).V  = 1;                                % with shrinkage priors

        % Nonlinear [Polynomial] factor analysis (static - nonlinear)
        %==================================================================
    case{'nonlinear model','nlfa'}

        % non-hierarchical nonlinear generative model (dynamic)
        %==================================================================
    case{'state space model','ssm','double well'}

        % temporal correlations
        %------------------------------------------------------------------
        M(1).E.s  = 1/2;

        % model specification - 1st level
        %------------------------------------------------------------------
        f      = '(16*x./(1 + x.^2) - x/2 + v*2)/8';
        g      = '(x.^2)/16';
        M(1).x = 1;
        M(1).f = inline(f,'x','v','P');
        M(1).g = inline(g,'x','v','P');
        M(1).V = exp(2);
        M(1).W = exp(16);

        % model specification - 2nd level
        %------------------------------------------------------------------
        M(2).v = 0;
        M(2).V = 1/8;

                % non-hierarchical nonlinear generative model (dynamic)
        %==================================================================
    case{'lorenz'}

        % correlations
        %--------------------------------------------------------------------------
        M(1).E.linear = 3;
        M(1).E.s      = 1/8;

        % level 1
        %--------------------------------------------------------------------------
        P       = [18.0; -4.0; 46.92];
        x       = [0.9; 0.8; 30];
        f       = '[-P(1) P(1) 0; (P(3)-x(3)) -1 -x(1); x(2) x(1) P(2)]*x/32;';
        M(1).f  = f;
        M(1).g  = inline('sum(x)','x','v','P');
        M(1).x  = x;
        M(1).pE = P;
        M(1).V  = exp(0);
        M(1).W  = exp(16);

        % level 2
        %--------------------------------------------------------------------------
        M(2).v  = 0;
        M(2).V  = exp(16);

        % Ornstein_Uhlenbeck linear generative model (dynamic)
        %==================================================================
    case{'Ornstein_Uhlenbeck','ou'}

        % temporal correlations
        %------------------------------------------------------------------
        M(1).E.linear = 1;                          % linear model
        M(1).E.s  = 1/2;

        % model specification - 1st level
        %------------------------------------------------------------------
        f      = '-1/16*x + v';
        g      = 'x';
        M(1).x = 0;
        M(1).f = inline(f,'x','v','P');
        M(1).g = inline(g,'x','v','P');
        M(1).V = exp(8);
        M(1).W = exp(32);

        % model specification - 2nd level
        %------------------------------------------------------------------
        M(2).v = 0;
        M(2).V = exp(-8);

    otherwise

        errordlg('unknown model; please add to spm_DEM_M')
        return

end

% check and conplete model specification
%---------------------------------------------------------------------------
M       = spm_DEM_M_set(M);


text:
function [x, infos] = ns_nmf(V, rank, in_options)
% Nonsmooth nonnegative matrix factorization (Nonsmooth-NMF)
%
% The problem of interest is defined as
%
%      min || V - W*S*H ||_F^2,
%
%      or
%
%      min  D(V||W*S*H),
%
%      where
%      {V, W, S, H} > 0.
%
% Given a non-negative matrix V, factorized non-negative matrices {W, S, H} are calculated.
%
%
% Inputs:
%       V           : (m x n) non-negative matrix to factorize
%       rank        : rank
%       in_options
%
%
% Output:
%       x           : non-negative matrix solution, i.e., x.W: (m x rank), x.H: (rank x n)
%       infos       : log information
%           epoch   : iteration nuber
%           cost    : objective function value
%           optgap  : optimality gap
%           time    : elapsed time
%           grad_calc_count : number of sampled data elements (gradient calculations)
%
% Reference:
%       A. Pascual-Montano, J. M. Carazo, K. Kochi, D. Lehmann, and R. D. Pascual-Marqui,
%       "Nonsmooth nonnegative matrix factorization (Nonsmooth-NMF),"
%       IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), vol.28, no.3, pp.403-415, 2006.
%
%       Z. Yang, Y. Zhang, W. Yan, Y. Xiang, and S. Xie,
%       "A fast non-smooth nonnegative matrix factorization for learning sparse representation,"
%       IEEE Access, vol.4, pp.5161-5168, 2016.
%
%
% This file is part of NMFLibrary.
%
% This file is originally created by Graham Grindlay.
%
% 2010-01-14 Graham Grindlay (grindlay@ee.columbia.edu)
%
% Copyright (C) 2008-2010 Graham Grindlay (grindlay@ee.columbia.edu)
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%
% This file is partially created by Silja Polvi-Huttunen, University of Helsinki, Finland, 2014
%
%
% Created by modifiying the original code by H.Kasai on Jul. 23, 2018
%
% Change log:
%
%       Jul. 29, 2019 (Hiroyuki Kasai): Modified.
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jun. 24, 2022 (Hiroyuki Kasai): Fixed a bug of W normalization.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options.theta         = 0.5; % decides the degree in [0,1] of nonsmoothing (use 0 for standard NMF)
    local_options.metric_type   = 'euc'; % 'euc' (default) or 'kl-div'
    local_options.update_alg    = 'mu';  % 'mu' or 'apg'
    local_options.apg_maxiter   = 100;
    local_options.myeps         = 1e-16;
    local_options.norm_w        = 1;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = 'Nonsmooth-NMF';
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % initialize for this algorithm
    I = eye(rank);
    S = (1-options.theta) * I + (options.theta/rank) * ones(rank);

    % store initial info
    clear infos;
    WS = W * S;
    [infos, f_val, optgap] = store_nmf_info(V, WS, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        if strcmp(options.update_alg, 'mu')

            % update H
            WS = W * S;
            if strcmp(options.metric_type, 'euc')
                H = H .* (WS' * V) ./ ((WS' * WS) * H + 1e-9);
            elseif strcmp(options.metric_type, 'kl-div')
                H = H .* (WS' * (V./(WS * H + 1e-9))) ./ (sum(WS, 1)' * ones(1,n));
            else
                error('Invalid metric type')
            end

            % normalize rows in H
            [W, H] = normalize_WH(V, W, H, rank, 'type2');

            % update W
            SH = S * H;
            if strcmp(options.metric_type, 'euc')
                W = W .* (V * SH') ./ (W * (SH * SH') + 1e-9);
            elseif strcmp(options.metric_type, 'kl-div')
                W = W .* ((V ./ (W * SH + 1e-9)) * SH') ./ (ones(m, 1) * sum(SH,2)');
            end

        else % support only 'euc' metric.

            if 0
                % update H
                WS = W*S;
                [H, ~, ~] = nesterov_mnls_general(V, WS, [], H, 1, options.apg_maxiter, 'basic');

                % normalize rows in H
                [W, H] = normalize_WH(V, W, H, rank, 'type2');

                % update W
                SH = S*H;
                [W, ~, ~] = nesterov_mnls_general(V, [], SH', W, 1, options.apg_maxiter, 'basic');

            else

                % update W
                SH = S * H;
                [W, ~, ~] = nesterov_mnls_general(V, [], SH', W, 1, options.apg_maxiter, 'basic');
                %W_prev = W;
                W = W + (W<options.myeps) .* options.myeps;

                % normalize W
                if options.norm_w
                    %W11 = bsxfun(@rdivide,W,sqrt(sum(W.^2,1)));
                    W = normalize_W(W, 2);
                end

                % update H
                WS = W*S;
                [H, ~, ~] = nesterov_mnls_general(V, WS, [], H, 1, options.apg_maxiter, 'basic');

            end

        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % update epoch
        epoch = epoch + 1;

        % store info
        WS = W * S;
        infos = store_nmf_info(V, WS, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;
    x.S = S;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function h=m_windbarb(long,lat,u,v,varargin)
%  M_WINDBARB Project wind barbs onto map axes
%
%  M_WINDBARB(LONG,LAT,U,V) projects two dimensional wind barbs onto the
%  current map axes. The vector components (U,V) are assumed to be in units
%  of knots and are specified at the points (LON,LAT). It handles winds up
%  to 130 knots. Winds exceeding 130 knots will appear as 130 knots.
%
%  If (u,v) are both positive, barbs appear to the SW of the vector head.
%  (in quiver, the arrowhead would be in the NE)
%
%  M_WINDBARB(...,S) uses the input S to scale the vectors after
%  they have been automatically scaled to fit within the grid. If omitted,
%  S = 0.9 is assumed.
%
%  M_WINDBARB(...,'PropertyName',PropertyValue,...) allows you to specify
%  additional LINE properties.
%
%  An additional parameter/value pair allows for u/v vectors in units
%  other than knots:
%          'units' : 'knots' | 'm/s' | 'kmh' | 'mph'
%

% Original code:
%  MFILE:   m_windbarb.m
%  MATLAB:  9.0.0 (R2016a)
%  VERSION: 1.0 (19 March 2017)
%  AUTHOR:  Erye
%  CONTACT: tfoterye@gmail.com
%
%  Oct/2017 - code "improved" for m_map style

global MAP_PROJECTION MAP_VAR_LIST

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

scale=0.9;
if nargin>4 && ~ischar(varargin{1})
    scale=varargin{1};
    varargin(1)=[];
end

if scale==0
    error(['map:' mfilename ':invalidScale'], ...
            'Invalid scale factor - must be greater than zero.')
end
%      1 knot = 0.5144 m/s = 1.852 kmh = 1.151 mph

scf=1;  % assume knots
k=1;
while k<=length(varargin)
    switch lower(varargin{k}(1:3))
        case 'uni'
            switch lower(varargin{k+1})
                case 'knots'
                    scf=1;
                case {'m/s','meters/sec','meters/second'}
                    scf=0.5144;
                case 'kmh'
                    scf=1.852;
                case 'mph'
                    scf=1.151;
                otherwise
                    error(['map:' mfilename ':invalidUnits'], ...
                        'Units specified not recognized');
            end
            varargin([k k+1])=[];
        otherwise
            k=k+2;
    end
end

[X,Y]=m_ll2xy(long,lat,'clip','point');

[XN ,YN ]=m_ll2xy([long(:) long(:)]',[lat(:) lat(:)+.001]','clip','off');
[XE ,YE ]=m_ll2xy([long(:) long(:)+(.001)./cos(lat(:)*pi/180)]',[lat(:) lat(:)]','clip','off');
mU=u.*reshape(diff(XE),size(lat))*1000 + v.*reshape(diff(XN),size(lat))*1000;
mV=u.*reshape(diff(YE),size(lat))*1000 + v.*reshape(diff(YN),size(lat))*1000;

umag = sqrt(u.^2+v.^2)/scf ; %wind speed (should be in knots as input)
theta = atan2(mV,mU);

%create 18 logical matrices for 18 possible barbs. Non-zero when the barb
%is called for at that gridpoint.
g{1} =  umag > 7.5  & umag <= 47.5;
g{2} =  umag > 17.5 & umag <= 47.5;
g{3} =  umag > 27.5;
g{4} = (umag > 37.5 & umag <= 47.5) | (umag > 57.5 & umag <= 97.5);
g{5} =  umag > 67.5;
g{6} = (umag > 77.5 & umag <  97.5) | umag > 107.5;
g{7} =  umag > 87.5 & umag < 97.5   | umag > 117.5;
g{8} =  umag > 127.5;
g{9} = (umag > 2.5  & umag <= 7.5 ) | (umag > 12.5 & umag <= 17.5);
g{10} = umag > 22.5 & umag <= 27.5;
g{11} =(umag > 32.5 & umag <= 37.5) | (umag > 52.5 & umag <= 57.5);
g{12} =(umag > 42.5 & umag <= 47.5) | (umag > 62.5 & umag <= 67.5);
g{13} =(umag > 72.5 & umag <= 77.5) | (umag > 102.5 & umag <= 107.5);
g{14} =(umag > 82.5 & umag <= 87.5) | (umag > 112.5 & umag <= 117.5);
g{15} =(umag > 92.5 & umag <= 97.5) | (umag > 122.5 & umag <= 127.5);
g{16} = umag > 47.5;
g{17} = umag > 97.5;
g{18} = umag>0;    % All have background line

%position of each barb relative to grid point: [x0 y0; x1 y1]
c{1} = [-1    0; -1.125 .325];  % Full barb 1
c{2} = [-.875 0; -1     .325];  % Full barb 2
c{3} = [-.75  0; -.875  .325];  % Full barb 3
c{4} = [-.625 0; -.75   .325];  % Full barb 4
c{5} = [-.5   0; -.625  .325];  % Full barb 5
c{6} = [-.375 0; -.5    .325];  % Full barb 6
c{7} = [-.25  0; -.375  .325];  % Full barb 7
c{8} = [-.125 0; -.25   .325];  % Full barb 8
c{9} = [-.875 0; -.9375 .1625]; % Half barb 2
c{10} =[-.75  0; -.8125 .1625]; % Half barb 3
c{11} =[-.625 0; -.6875 .1625]; % Half barb 4
c{12} =[-.5   0; -.5625 .1625]; % Half barb 5
c{13} =[-.375 0; -.4375 .1625]; % Half barb 6
c{14} =[-.25  0; -.3125 .1625]; % Half barb 7
c{15} =[-.125 0; -.1875 .1625]; % Half barb 8
c{16} =[-1    0; -.875  .325];  % Slant other way for triangle 1
c{17} =[-.75  0; -.625  .325];  % Slant other way for taingle
c{18} =[0     0; -1      0  ];  % Base

%set scale based on average latitude spacing
[m,n]=size(X);
scale2 = scale*(max(max(X))-min(min(X)))/n;

%draw the barbs
Ax=[];Ay=[];
for nn = 1:18

   ivals=g{nn}(:);
   if sum(ivals)>0  % number of barbs to draw

      %rotation operations
      cthet=cos(theta(ivals));
      sthet=sin(theta(ivals));

      x1=c{nn}(1,1)*cthet-c{nn}(1,2)*sthet;
      y1=c{nn}(1,1)*sthet+c{nn}(1,2)*cthet;
      x2=c{nn}(2,1)*cthet-c{nn}(2,2)*sthet;
      y2=c{nn}(2,1)*sthet+c{nn}(2,2)*cthet;

      x1 = x1*scale2+X(ivals);
      x2 = x2*scale2+X(ivals);
      y1 = y1*scale2+Y(ivals);
      y2 = y2*scale2+Y(ivals);
      x = [x1 x2 NaN(size(x1))]';   % Speed up by vectorizing this call with
      y = [y1 y2 NaN(size(y1))]';   %    Nans between line segments.

      Ax=[Ax;x(:)];
      Ay=[Ay;y(:)];
   end
end
h=line(Ax,Ay,varargin{:});

set(h,'tag','m_windbarb');

if nargout==0
 clear h
end

author: goGPS-Project


text:
function [rval_space,rval_time,ind_space,ind_time] = classify_comp_corr(Yr,A,C,b,f,options)

% Classify components by looking at how well the spatial and temporal components
% correlate with the raw data at point of high activation

% INPUTS
% Yr:       reshaped data in 3D matrix or memmaped file
% A,C,b,f:  identified components
% options:  options structure
%    space_thresh:  r-value threshold for spatial components (default: 0.4)
%    time_thresh:  r-value threshold for temporal components (default: 0.4)
%    Athresh:  threshold for determining spatial overlap (default: 0.1)
%    Np:       number of high activity intervals for each component (default: 20)
%    peak_int: Interval around each local peak to be considered (default: -2:6)
%    MinPeakDist:   minimum peak distance for finding points of high activity  (default: 10)

% OUTPUTS:
% rval_space:     r-values between spatial components and data patches
% rval_time:      r-values between temporal components and trace averages
% ind_space:      components with rval_space > space_thresh
% ind_time:       components with rval_time > time_thresh
% MinPeakDist:

% Written by Eftychios A. Pnevmatikakis, Simons Foundation, 2016
% based on discussions with Matt Kaufman and Farzaneh Najafi, CSHL

defoptions = CNMFSetParms;
if nargin < 6 || isempty(options)
    options = defoptions;
end

if ~isfield(options,'peak_int') || isempty(options.peak_int); int = defoptions.peak_int; else int = options.peak_int; end
if ~isfield(options,'Npeaks') || isempty(options.peak_int); Np = defoptions.Npeaks; else Np = options.Npeaks; end
if ~isfield(options,'A_thresh') || isempty(options.A_thresh); Athresh = defoptions.A_thresh; else Athresh = options.A_thresh; end
if ~isfield(options,'space_thresh') || isempty(options.space_thresh); options.space_thresh = defoptions.space_thresh; end
if ~isfield(options,'time_thresh') || isempty(options.time_thresh); options.time_thresh = defoptions.time_thresh; end
if ~isfield(options,'MinPeakDist') || isempty(options.MinPeakDist); options.MinPeakDist = defoptions.MinPeakDist; end

memmaped = isobject(Yr);

[K_m,T] = size(C);
pk = zeros(K_m,Np);

parfor i = 1:K_m
    [~,pk_temp] = findpeaks(C(i,:),'SortStr','descend','Npeaks',Np,'MinPeakDistance',options.MinPeakDist);
    if ~isempty(pk_temp)
        if length(pk_temp) < Np
            pk_temp(length(pk_temp)+1:Np) = pk_temp(1);
        end
        pk(i,:) = pk_temp;
    end
end

%% expand intervals
lint = length(int);

locs = repmat(pk,1,lint) + kron(int,ones(size(pk)));
locs = sort(locs,2,'ascend');
locs(locs<1)=1;
locs(locs>T)=T;
LOCS = mat2cell(locs,ones(K_m,1),Np*lint);
for i = 1:K_m
    LOCS{i} = unique(LOCS{i});
end

%% compute r-values
nA = full(sqrt(sum(A.^2)))';
tAA = (A'*A)./(nA*nA') > Athresh;
tAA(1:K_m+1:K_m^2) = 0;

rval_space = zeros(K_m,1);
rval_time  = zeros(K_m,1);

if options.d3 == 1; b_rs = reshape(b,options.d1,options.d2,[]); else b_rs = reshape(b,options.d1,options.d2,options.d3,[]); end
if memmaped
    parfor i = 1:K_m
        ovlp_cmp = find(tAA(:,i));
        indeces = LOCS{i};
        for j = 1:length(ovlp_cmp)
            indeces = setdiff(indeces,LOCS{ovlp_cmp(j)});
        end
        a_temp = reshape(A(:,i),options.d1,options.d2*options.d3);
        [rows,temp] = find(a_temp>0);
        [cols,plns] = ind2sub([options.d2,options.d3],temp);
        if options.d3 > 1
            a_temp = reshape(full(A(:,i)),options.d1,options.d2,options.d3);
            a_temp = a_temp(min(rows):max(rows),min(cols):max(cols),min(plns):max(plns));
        else
            a_temp = a_temp(min(rows):max(rows),min(cols):max(cols));
        end

        if options.d3 == 1
            b_temp = reshape(b_rs(min(rows):max(rows),min(cols):max(cols),:),numel(a_temp),[]);
        else
            b_temp = reshape(b_rs(min(rows):max(rows),min(cols):max(cols),min(plns):max(plns),:),numel(a_temp),[]);
        end
        if ~isempty(indeces)
            if options.d3 == 1
                yytemp = Yr.Y(min(rows):max(rows),min(cols):max(cols),min(indeces):max(indeces));
                y_temp = yytemp(:,:,indeces-min(indeces)+1);
            else
                yytemp = Yr.Y(min(rows):max(rows),min(cols):max(cols),min(plns):max(plns),min(indeces):max(indeces));
                y_temp = yytemp(:,:,:,indeces-min(indeces)+1);
            end
            y_temp = reshape(y_temp,[],length(indeces));
            mY_space = double(mean(y_temp,2) - b_temp*mean(f(:,indeces),2));
            mY_time = double(mean(y_temp,1)-mean(b_temp,1)*f(:,indeces));
            rval_space(i) = corr(full(a_temp(:)),mY_space);
            rval_time(i) = corr(C(i,indeces)',mY_time');
        else
            rval_space(i) = NaN;
            rval_time(i) = NaN;
        end
    end
else
    parfor i = 1:K_m
        ovlp_cmp = find(tAA(:,i));
        indeces = LOCS{i};
        for j = 1:length(ovlp_cmp)
            indeces = setdiff(indeces,LOCS{ovlp_cmp(j)});
        end
        a_temp = reshape(A(:,i),options.d1,options.d2*options.d3);
        [rows,temp] = find(a_temp>0);
        [cols,plns] = ind2sub([options.d2,options.d3],temp);
        if options.d3 > 1
            a_temp = reshape(full(A(:,i)),options.d1,options.d2,options.d3);
            a_temp = a_temp(min(rows):max(rows),min(cols):max(cols),min(plns):max(plns));
        else
            a_temp = a_temp(min(rows):max(rows),min(cols):max(cols));
        end

        if options.d3 == 1
            b_temp = reshape(b_rs(min(rows):max(rows),min(cols):max(cols),:),numel(a_temp),[]);
        else
            b_temp = reshape(b_rs(min(rows):max(rows),min(cols):max(cols),min(plns):max(plns),:),numel(a_temp),[]);
        end
        if ~isempty(indeces)
            if options.d3 == 1
                y_temp = Yr(min(rows):max(rows),min(cols):max(cols),indeces);
            else
                y_temp = Yr(min(rows):max(rows),min(cols):max(cols),min(plns):max(plns),indeces);
            end
            y_temp = reshape(y_temp,[],length(indeces));
            mY_space = double(mean(y_temp,2) - b_temp*mean(f(:,indeces),2));
            mY_time = double(mean(y_temp,1)-mean(b_temp,1)*f(:,indeces));
            rval_space(i) = corr(full(a_temp(:)),mY_space);
            rval_time(i) = corr(C(i,indeces)',mY_time');
        else
            rval_space(i) = NaN;
            rval_time(i) = NaN;
        end
    end
end

ind_space = rval_space > options.space_thresh;
ind_time = rval_time > options.time_thresh;

author: flatironinstitute

repo: CaImAn-MATLAB


text:
clear all; close all;
I=imread('rice.png');
J=im2bw(I, graythresh(I));
K=bwlabel(J);
RGB=label2rgb(K);
figure;
subplot(121);  imshow(J);
subplot(122);  imshow(RGB);

author: HuangCongQing


text:
here=pwd;
cd /home/bernard/Data/ABC/animals
attrClasses=dlmread('all-perclass.attributes');
indices=dlmread('split0.mask');
classes=dlmread('all.classid');
cd (here);

trainInstancesIndices=find(indices==1);
trainInstancesLabels=classes(indices==1);
trainClassesIndices=unique(trainInstancesLabels,'stable');

auxLabels=zeros(size(trainInstancesLabels));
for i=1:length(unique(trainInstancesLabels))
    c=find(trainInstancesLabels>0, 1);
    class=trainInstancesLabels(c);
    auxLabels(trainInstancesLabels==class)=i;
    trainInstancesLabels(trainInstancesLabels==class)=0;
end
trainInstancesLabels=auxLabels';

testInstancesIndices=find(indices==0);
testInstancesLabels=classes(indices==0);
testClassesIndices=unique(testInstancesLabels,'stable');

auxLabels=zeros(size(testInstancesLabels));
for i=1:length(unique(testInstancesLabels))
    c=find(testInstancesLabels>0, 1);
    class=testInstancesLabels(c);
    auxLabels(testInstancesLabels==class)=i;
    testInstancesLabels(testInstancesLabels==class)=0;
end
testInstancesLabels=auxLabels';

attrClasses=unique(attrClasses, 'rows', 'stable');
uniqueClasses=unique(classes, 'stable');
newAttrClasses=zeros(size(attrClasses));
for i=1:length(uniqueClasses)
    newAttrClasses(uniqueClasses(i),:)=attrClasses(i,:);
end
attrClasses=newAttrClasses;

repo: Embarrassingly-simple-ZSL


text:
% invfirwsord() - Estimate windowed sinc FIR filter transition band width
%                 depending on filter order and window type
%
% Usage:
%   >> [df, dev] = invfirwsord(wtype, fs, m);
%   >> df = invfirwsord('kaiser', fs, m, dev);
%
% Inputs:
%   wtype - char array window type. 'rectangular', 'bartlett', 'hann',
%           'hamming', 'blackman', or 'kaiser'
%   fs    - scalar sampling frequency}
%   m     - scalar filter order
%   dev   - scalar maximum passband deviation/ripple (Kaiser window
%           only)
%
% Output:
%   df    - scalar estimated transition band width
%   dev   - scalar maximum passband deviation/ripple
%
% References:
%   [1] Smith, S. W. (1999). The scientist and engineer's guide to
%       digital signal processing (2nd ed.). San Diego, CA: California
%       Technical Publishing.
%   [2] Proakis, J. G., & Manolakis, D. G. (1996). Digital Signal
%       Processing: Principles, Algorithms, and Applications (3rd ed.).
%       Englewood Cliffs, NJ: Prentice-Hall
%   [3] Ifeachor E. C., & Jervis B. W. (1993). Digital Signal
%       Processing: A Practical Approach. Wokingham, UK: Addison-Wesley
%
% Author: Andreas Widmann, University of Leipzig, 2005
%
% See also:
%   firws, firwsord

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2005-2014 Andreas Widmann, University of Leipzig, widmann@uni-leipzig.de
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
% $Id$

function [ df, dev ] = invfirwsord(wintype, fs, m, dev)

winTypeArray = {'rectangular', 'bartlett', 'hann', 'hamming', 'blackman', 'kaiser'};
winDfArray = [0.9 2.9 3.1 3.3 5.5];
winDevArray = [0.089 0.056 0.0063 0.0022 0.0002];

% Check arguments
if nargin < 3 || isempty(fs) || isempty(m) || isempty(wintype)
    ft_error('Not enough input arguments.')
end

% Window type
wintype = find(strcmp(wintype, winTypeArray));
if isempty(wintype)
    ft_error('Unknown window type.')
end

if wintype == 6 % Kaiser window
    if nargin < 4 || isempty(dev)
        ft_error('Not enough input arguments.')
    end
    devdb = -20 * log10(dev);
    df = (devdb - 8) / (2.285 * 2 * pi * (m - 1));
else
    df = winDfArray(wintype) / m;
    dev = winDevArray(wintype);
end

% df is normalized
df = df * fs;

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function psi = data2psiX(data,srate,freqbins,permtest)
% calculates phase slope index (PSI) as formulated in the paper:
%    Nolte G, Ziehe A, Nikulin VV, Schlogl A, Kramer N, Brismar T, Muller KR.
%    Robustly estimating the flow direction of information in complex physical systems.
%    Physical Review Letters. To appear.
%    (for further information:    http://doc.ml.tu-berlin.de/causality/ )
%
% Usage:
%   psi = data2psiX(data,srate,freqbins,permtest);
%
% Input:
%     data:  MxNxT matrix for M channels, N timepoints, and T trials
%    srate:  data sampling rate in Hz
% freqbins:  KxQ matrix containing frequency boundaries in rows K.
% permtest:  Permutation test (boolean). If true, psi output values are in standardizes Z values
%               relative to a null hypothesis distribution
%
% Output:
%      psi:  phase-slope-index values. For M channels PSI is an MxM matrix if
%               one frequency bin, or MxMxK if freqbins has K rows (with K>1).
%               psi(i,j) is the directed connectivity from channel i to
%               channel j, (e.g., channel i is the sender if psi(i,j) is
%               positive; channel i is the receiver if negative)
%
%
%  (This function was modified from the original by Mike X Cohen)

% License
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see http://www.gnu.org/licenses/.

% modified from the original dat2psi.m by mikexcohen@gmail.com

%% check inputs and initialize

if nargin<2
    help data2psiX
    error('Read help file.');
elseif nargin==2
    permtest = 0;
    freqbins = [];
elseif nargin==3
    permtest = 0;
end

n_permutes = 1000;

% data dimensions
[nchan,npnts,ntrials]=size(data);

% define frequencies
hz  = linspace(0,srate/2,floor(npnts/2)+1);
nhz = length(hz);

% initialize
cs  = zeros(nchan,nchan,nhz);
psi = zeros(nchan,nchan,size(freqbins,1));

%% compute cross-spectral density

datafft = fft(bsxfun(@times,data,hanning(npnts)'),[],2);

% trial-average cross-spectral density
for triali=1:ntrials
    for freqi=1:nhz

        tempfftdat = squeeze(datafft(:,freqi,triali))';
        cs(:,:,freqi) = cs(:,:,freqi) + tempfftdat'*conj(tempfftdat);
    end
end
cs = cs./triali;

%% compute PSI

for freqbini=1:size(freqbins,1)

    % find FFT indices of requested frequency bands
    freqidx = dsearchn(hz',freqbins(freqbini,1)) : dsearchn(hz',freqbins(freqbini,2));
    nfidx   = length(freqidx);

    if nfidx<4
        warning('There are fewer than four frequency bins. Consider using longer time windows or wider frequency bands.')
    end

    % temporary phase-frequency matrix from this frequency band
    pp = zeros(nchan,nchan,nfidx);

    for fi=1:nfidx
        pp(:,:,fi) = cs(:,:,freqidx(fi))./sqrt(diag(cs(:,:,freqidx(fi)))*diag(cs(:,:,freqidx(fi)))');
    end

    % average phase slope for each frequency band
    psi(:,:,freqbini) = sum(imag(conj(pp(:,:,1:end-1)).*pp(:,:,2:end)),3);

    %% optional permutation testing

    if permtest

        nulldist = zeros(nchan,nchan,n_permutes);
        for permi=1:n_permutes
            nulldist(:,:,permi) = sum(imag(conj(pp(:,:,randperm(nfidx))).*pp(:,:,randperm(nfidx))),3);
        end

        psi(:,:,freqbini) = ( psi(:,:,freqbini)-mean(nulldist,3) ) ./ std(nulldist,[],3);
    end

    % zero-out diagonals
    tmp = squeeze(psi(:,:,freqbini));
    tmp(logical(eye(nchan))) = 0;
    psi(:,:,freqbini) = tmp;

end

%% end

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
%This program localizes single molecules and determines their positions
%by fitting a gaussian to their intensity profile.  A position dependent
%background subtraction is used.

function einzelreader(handles, h)
    global xpix ypix wbox psf_w02;

    version = '6/23/08'; %record version (date) of gui being used

    rbox=str2double(get(handles.rbox_edit,'String'));
    wbox=2*rbox+1;
    [xpix,ypix] = meshgrid(-rbox:rbox,-rbox:rbox);

    q       = str2double(get(handles.cam_pix_edit,'String'));
    wvlnth  = str2double(get(handles.wvlnth,'String'))/1000; %convert wavelength from nm to um
    NA      = str2double(get(handles.NA,'String'));
    psf_scale = str2double(get(handles.psf_scale_edit,'String'));

    psf_w0 = psf_scale*0.55*wvlnth/NA/1.17; % 1/e2 radius of PSF in um, use 1/e2 = FWHM/1.17 from Pawley
                                      % with scale factor (20% for "real objective" due to
                                      % measured PSF from Hess and Webb
                                      % 2002)
    psf_std=psf_w0/2; %standard deviation of psf
    psf_w02=(psf_w0/q)*(psf_w0/q); %square of 1/e^2 radius in pixels

    yfit_psf=exp(-2*((xpix).*(xpix)+(ypix).*(ypix))/psf_w02);
    npix=sum(sum(yfit_psf));     % area of molecule in square pixels

    total_molecules=0;
    xcm_all=zeros(1,1);
    ycm_all=zeros(1,1);
    framenum_all=zeros(1,1);
    n_fail_a0=0;
    n_fail_outbox=0;

    imagefile = get(handles.im_file_edit, 'String');
    outpath = get(handles.out_dir_edit, 'String');
    if outpath(length(outpath))~='/'
        outpath = [outpath,'/'];
    end
    preview = get(handles.show_preview,'Value');

    fs0 = str2double(get(handles.fs,'String')); %Setup for preview frame skipping
    fs1 = fs0;

    if get(handles.rolling_ball_cb,'Value') %if using rolling ball subtraction
        %Obtain info on image type, etc
        file_info = get_file_info2(imagefile, handles);

        n_end   = file_info.stop;
        n_start = file_info.start;

        if get(handles.custom_roi,'Value')
            x_offset   = str2double(get(handles.x_off_edit,'String'));
            y_offset   = str2double(get(handles.y_off_edit,'String'));
            x_size     = str2double(get(handles.x_size_edit,'String'));
            y_size     = str2double(get(handles.y_size_edit,'String'));
        else
            x_offset = 1;
            y_offset = 1;
            x_size   = file_info.width;
            y_size   = file_info.height;
        end

        rball=str2double(get(handles.rb_radius_edit,'String')); %radius of rolling ball
        se = strel('ball',rball,rball,0); %structural element, i.e. rolling ball

        FWHM=str2double(get(handles.FWHM_edit,'String'));; %FWHM of gaussian smoothing in pixels
        rk=(FWHM)/sqrt(2*log(2)); %1/e^2 smoothing radius in pixels
        kw=20; %kernal width of smoothing function
        [X,Y]=meshgrid(-kw/2:kw/2,-kw/2:kw/2);
        kd=sqrt(X.*X+Y.*Y);
        gs=exp(-2*kd.*kd/(rk*rk));
        gs=gs/sum(sum(gs)); %smoothing function normalized to have area = 1

    else %else use integrated widefield sum
        %Load wf sum & mean background mat file
        meanbkg_mat_file = getappdata(0,'last_mbkg_save'); %strcat(dir,'analysis\',date,'/',base_name,'sum_wf_meanbkg.mat');
        load(meanbkg_mat_file,'meanbkg_all','image_sum','x_size','y_size','x_offset','y_offset','field0','n_field','file_info');%,'zero_level'); %mat file from sum_wf_meanbkg.m

        n_end   = file_info.stop;
        n_start = file_info.start;
        %give warning if starting/ending frames are not compatible with wf sum & mean background mat file
        if(n_start<field0)
            warndlg(sprintf('Warning: Starting frame is less than corresponding frame used in Compute WF & Mean Background, need to change starting frame or re-run Compute WF & Mean Background'));
            return
        end
        if(n_end>n_field)
            warndlg(sprintf('Warning: Ending frame is greater than corresponding frame used in Compute WF & Mean Background, need to change ending frame or re-run Compute WF & Mean Background'));
            return
        end

        bkg_pc = str2double(get(handles.bkg_percent,'String'))/100;
        zero_level = str2double(get(handles.zero_lvl,'String'));

        %normalize so that ave pixel value = 1
        image_sum_mean=mean(mean(image_sum(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1)));
        image_sum_norm = (image_sum(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1))/image_sum_mean;
    end

    %Load variables from the GUI edit boxes
    %--------------------------------------
    if get(handles.pixel_photon,'Value')
        pix_to_pho = str2double(get(handles.pix_to_pho,'String'));
    else
        pix_to_pho = 1;
    end

    iprod_thresh    = str2double(get(handles.iprod_thresh_edit,'String'))*pix_to_pho; %initial threshold for a "bright object"
    threshold       = str2double(get(handles.threshold_edit,'String'))*pix_to_pho; %this value to be included
    upper_threshold = str2double(get(handles.upperthresh_edit,'String'))*pix_to_pho; %this value to not be excluded

    n_bright_pixel_threshold   = str2double(get(handles.min_above_edit,'String')); %min # of pixels above
    box_overlap_factor         = str2double(get(handles.box_overlap_edit,'String')); %if center to center closer, don't include either
    max_pixels_above_threshold = str2double(get(handles.max_above_edit,'String')); %max # of pixels above

    %Initialize temporary runtime file for passing variables
    %-------------
    for x=1:99
        temp_runtime = [outpath,'temp_runtime_',num2str(x,'%2.2d'),'.mat'];
        if ~exist(temp_runtime,'file')
            break
        end
    end
    setappdata(0,'temp_runtime',temp_runtime);
    %-------------

    %Initialize waitbar and pause component
    w = waitbarxmod(0,'Executing "einzelreader.m" ...','CreateCancelBtn','delete(gcf)');
    set(w,'Name','Progress Bar');
    uicontrol('Style','pushbutton','Parent',w,'String','Pause','Position',[210,10,60,23], ...
              'UserData',1,'Callback',@pause_gui);
    drawnow;   %Draw the extra button immediately
    pause(.1); %Pause to ensure window completes drawing
    %keepontop('Progress Bar');

    wb_norm = n_end-n_start;
    if wb_norm == 0 %Avoid divide by zero errors
        wb_norm = 1;
    end

    if h~=0                             %If preview window is active, store the handle of the progress bar
        setappdata(h,'PrBar_Handle',w); %in the appdata of the figure
    end

    for fileloop=n_start:n_end
        if strcmp(file_info.type,'singlepage')
            index = num2str(fileloop,file_info.prec);
%             infile = [file_info.path,'/',file_info.part_name,index,file_info.ext];
            infile = [file_info.path,'/',file_info.part_name,index,file_info.ext];
            i1=double(imread(infile));
        elseif strcmp(file_info.type,'multipage')
            i1=double(imread(imagefile,fileloop));
        end

        compl = (fileloop-n_start)/wb_norm;
        drawnow;
        waitbarxmod(compl,w,'Executing "einzelreader.m" ...'); %Rename the waitbar

        if get(handles.rolling_ball_cb,'Value') %if using rolling ball subtraction
            i1_gs = uint16(conv2(i1,gs,'same')); %smoothed original
            rb_im = imopen(i1_gs,se);
            bkg   = zeros(y_size,x_size);
            bkg   = double(rb_im(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1));
        else
            bkg    = zeros(y_size,x_size);
            bkg    = image_sum_norm*bkg_pc*meanbkg_all(fileloop)+zero_level;
        end

        iprod = zeros(y_size,x_size);
        iprod = i1(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1)-bkg;
        iprod=iprod.*(iprod>0); % set any negative pixel values to zero

        high_pixel_mask = zeros(y_size,x_size);
        high_pixel_xy   = zeros(10000,3);  % x,y,active
        n_high_pixels   = 0;
        pix_val=iprod_thresh;

        while(pix_val >= iprod_thresh) %continue until reach minimum threshold
            pix_val=0;
            for i=rbox+1:y_size-rbox-1
                for j=rbox+1:x_size-rbox-1
                    if high_pixel_mask(i,j)==0 && iprod(i,j)>pix_val
                        pix_val=iprod(i,j); %find brightest pixel
                        high_pixel_y=i;
                        high_pixel_x=j;
                    end
                end
            end
            if(pix_val < iprod_thresh)
                break
            end
            x0box=high_pixel_x-rbox;
            y0box=high_pixel_y-rbox;
            x1box=high_pixel_x+rbox;
            y1box=high_pixel_y+rbox;

            high_pixel_mask(y0box:y1box,x0box:x1box)=1;
            n_high_pixels=n_high_pixels+1;
            high_pixel_xy(n_high_pixels,1)=high_pixel_x;
            high_pixel_xy(n_high_pixels,2)=high_pixel_y;
            high_pixel_xy(n_high_pixels,3)=1; % active
        end

        drawnow;              %**These added to keep waitbar responsive...
                              % Without these 'drawnow's, the queue gets large and
                              % ignores button presses during execution...
                              %(this function is fast when there is no
                              % actual drawing to do)

        waitbarxmod(compl,w); % <-check waitbar (fast,see waitbarxmod.m)

        for i=1:n_high_pixels
          for j=1:i
            if (i~=j & high_pixel_xy(i,3)==1 & high_pixel_xy(j,3)==1)
              dx=high_pixel_xy(i,1)-high_pixel_xy(j,1);
              dy=high_pixel_xy(i,2)-high_pixel_xy(j,2);
              dmin_nearest_box=sqrt(dx*dx+dy*dy);

              if (dmin_nearest_box<box_overlap_factor*rbox) % these boxes overlap
                if j>i
                  high_pixel_xy(j,3)=-3;       % make one of the boxes inactive
                else
                  high_pixel_xy(i,3)=-3;       % make one of the boxes inactive
                end
              end
            end
          end
        end

        drawnow;
        waitbarxmod(compl,w); %check waitbar

        n_boxes=0;
        boxes_xy=zeros(10000,3);

        for i=1:n_high_pixels
           if high_pixel_xy(i,3)==1
               n_boxes=n_boxes+1;
               boxes_xy(n_boxes,1)=high_pixel_xy(i,1);
               boxes_xy(n_boxes,2)=high_pixel_xy(i,2);
               boxes_xy(n_boxes,3)=1;
           end
        end

        drawnow;
        waitbarxmod(compl,w); %check waitbar

        failed=zeros(n_boxes,1);
        for i=1:n_boxes
           x0_box=boxes_xy(i,1)-rbox;
           y0_box=boxes_xy(i,2)-rbox;
           x1_box=boxes_xy(i,1)+rbox;
           y1_box=boxes_xy(i,2)+rbox;

           grab=iprod(y0_box:y1_box,x0_box:x1_box);
           minbkg=min(min(grab));
           grab=grab-minbkg;

           xm_sum=0;
           ym_sum=0;
           m_sum=0;
           for x=x0_box:x1_box
             for y=y0_box:y1_box
                xind=floor(x);
                yind=floor(y);
                intens=iprod(yind,xind);
                xm_sum=xm_sum+xind*intens;
                ym_sum=ym_sum+yind*intens;
                m_sum=m_sum+intens;
             end
           end

           x_cm(i)=xm_sum/m_sum;
           y_cm(i)=ym_sum/m_sum;

           %re-center box around center of mass
           if round(x_cm(i)) > 1+rbox && round(x_cm(i)) < x_size-rbox
               boxes_xy(i,1)=round(x_cm(i));
           end
           if round(y_cm(i)) > 1+rbox && round(y_cm(i)) < y_size-rbox
               boxes_xy(i,2)=round(y_cm(i));
           end

           x0_box=boxes_xy(i,1)-rbox;
           y0_box=boxes_xy(i,2)-rbox;
           x1_box=boxes_xy(i,1)+rbox;
           y1_box=boxes_xy(i,2)+rbox;

           grab=zeros(wbox,wbox);
           grab=iprod(y0_box:y1_box,x0_box:x1_box);
           minbkg=min(min(grab));
           grab=grab-minbkg;
           grab_sum(i)=sum(sum(grab));
           image_bright_in_box(:,:,i)=grab;

           xc_box=(x0_box+x1_box)*0.5;
           yc_box=(y0_box+y1_box)*0.5;

           xguess=x_cm(i)-xc_box;
           yguess=y_cm(i)-yc_box;

           for n=1:wbox
                for p=1:wbox
                    k=(n-1)*wbox+p;
                    xymerge(k)=0;
                    zmerge(k)=grab(n,p);
                end
           end

           beta=[xguess,yguess,50];
           [betafit,resid,J,COVB,mse] = nlinfit(xymerge,zmerge,@gaussian_merge,beta);
           ci = nlparci(betafit,resid,'covar',COVB); %calculate error estimates on parameters
           ci_err=(ci(:,2)-ci(:,1))/2;
           xf_err(i)=ci_err(1);
           yf_err(i)=ci_err(2);
           a0_err(i)=ci_err(3);

           zfitmerge=gaussian_merge(betafit,xymerge);

           for n=1:wbox
                for p=1:wbox
                    k=(n-1)*wbox+p;
                    zfit(n,p)=zfitmerge(k);
                end
           end

           yf_box(i)=betafit(2)+yc_box;
           xf_box(i)=betafit(1)+xc_box;
           a0_box(i)=betafit(3);

           if(a0_box(i) < 0)
                n_fail_a0=n_fail_a0+1;
                failed(i)=1;
           end
           if xf_box(i) > x1_box || xf_box(i) < x0_box || yf_box(i) > y1_box || yf_box(i) < y0_box
                n_fail_outbox=n_fail_outbox+1;
                failed(i)=1;
           end

           drawnow   %Flush event queue to keep cancel/pause button responsive on every pass.
                     %This is the slowest function by far and drawnow barely changes its exec time

           waitbarxmod(compl,w); %check waitbar
        end

        n_pixels_above_threshold=zeros(n_boxes,1);
        n_pixels_above_upper_threshold=zeros(n_boxes,1);

        for i=1:n_boxes
          if (boxes_xy(i,3)==1)
            n_pixels_above_threshold(i)=0;
            for j=1:2*rbox+1
              for k=1:2*rbox+1
                intens=image_bright_in_box(j,k,i);
                if intens>=threshold
                  n_pixels_above_threshold(i)=n_pixels_above_threshold(i)+1;
                end
                if intens>=upper_threshold
                  n_pixels_above_upper_threshold(i)=n_pixels_above_upper_threshold(i)+1;
                end
              end
            end
          end
        end

        drawnow;
        waitbarxmod(compl,w); %check waitbar

        n_molecules_found=0;

        for i=1:n_boxes
            if boxes_xy(i,3)==1
                if (n_pixels_above_threshold(i)>=n_bright_pixel_threshold && n_pixels_above_upper_threshold(i)<=max_pixels_above_threshold)
                    n_molecules_found=n_molecules_found+1;
                else
                    if (n_pixels_above_threshold(i)<n_bright_pixel_threshold)
                      boxes_xy(i,3)=-2; %draw red box
                    end
                    if (n_pixels_above_threshold(i)>max_pixels_above_threshold)
                      boxes_xy(i,3)=-1; %draw yellow box
                    end
                end
            end
        end

        drawnow;
        waitbarxmod(compl,w); %check waitbar

        if preview  %If preview active, update images on screen
            set(0,'CurrentFigure',h);
            if ~exist('plot2','var')
                position2 = [.01,.37,.25,.25];
                plot2 = axes('Position', position2);
                setappdata(plot2,'pos',position2);
            else
                set(h,'CurrentAxes',plot2);
            end
            i1display = i1(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1);
            imagesc(i1display,'HitTest','off');
            axis image;
            colormap gray;
            title(['Last Frame Processed: ', num2str(fileloop)])
            set(plot2, 'ButtonDownFcn',{@focus_swap,guidata(gcbo)});
            clear i1display;

            if ~exist('plot1','var')
                position1 = [.01,.05,.25,.25];
                plot1 = axes('Position', position1,'Color','k','YTickLabel','','XTickLabel','');
                axis image;
                set(plot1,'Position',position1);
                setappdata(plot1,'pos', position1);
            end

            if get(handles.rolling_ball_cb,'Value')
                if ~exist('plot3','var')
                    position3 = [.01,.70,.25,.25];
                    plot3 = axes('Position', position3);
                    setappdata(plot3,'pos',position3);
                else
                    set(h,'CurrentAxes',plot3);
             end
                imagesc(rb_im,'HitTest','off')
                colormap gray, axis image;
                title(['Rolling Ball Profile: ',num2str(fileloop)])
                set(plot3,'ButtonDownFcn',{@focus_swap,guidata(gcbo)});
            end

            if ~exist('plot4','var')
                if x_size < y_size
                    positionf = [.18,.065,.87,.87];
                else
                    positionf = [.32,.3,.6,.65];
                end
                plot4 = axes('Position', positionf);
                setappdata(0,'focused',plot4);
            else
                set(h,'CurrentAxes',plot4);
            end

            if get(handles.mnl_scale_check,'Value')
                clim = str2double(get(handles.mnl_scale_edit,'String'));
                imagesc(iprod,'HitTest','off',[0 clim]);
            else
                imagesc(iprod,'HitTest','off');
            end

            axis image;
            colormap gray;
            title(['Einzel Reader: Frame: ',num2str(fileloop)]);
            set(plot4, 'ButtonDownFcn',{@focus_swap,guidata(gcbo)});

            if (get(handles.show_colorbar,'Value') == 1)
                focused = getappdata(0, 'focused');
                set(h,'CurrentAxes', focused);
                pos = get(focused, 'Position');
                colorbar('Position',[.947, pos(2), .015, pos(4)],'DrawMode','fast')
            end

            set(h,'CurrentAxes', plot4);
        end

        if n_boxes>0

            if preview %Only draw if preview is active
                hold on
                draw_boxes(n_boxes,boxes_xy,rbox);
            end

            for i=1:n_boxes
                if boxes_xy(i,3)==1 && failed(i) ~= 1
                    if preview %Only plot if preview is active
                        %plot(x_cm(i),y_cm(i),'.m','HitTest','off'); %centroid
                        %plot(xf_box(i),yf_box(i),'.b','HitTest','off'); %gaussian fit
                    end
                    total_molecules=total_molecules+1;
                    xcm_all(total_molecules)=x_cm(i);
                    ycm_all(total_molecules)=y_cm(i);
                    xf_all(total_molecules)=xf_box(i);
                    yf_all(total_molecules)=yf_box(i);
                    a0_all(total_molecules)=a0_box(i);
                    grab_sum_all(total_molecules)=grab_sum(i);
                    framenum_all(total_molecules)=fileloop;
                    xf_err_all(total_molecules)=xf_err(i);
                    yf_err_all(total_molecules)=yf_err(i);
                    a0_err_all(total_molecules)=a0_err(i);
                end
            end
            nmol_all(fileloop)=n_boxes;%%

            if preview %Remove the hold
                hold off
            end

            drawnow;              %Draw any boxes and update einzelreader image (flush queue and check waitbar if preview not active)
            waitbarxmod(compl,w); %Update the waitbar
        end

        if get(handles.frame_delay_check,'Value')
            delay = str2double(get(handles.frame_delay_edit,'String'))/1000;
            if delay ~= inf
                pause(delay);
            end
        end

        if preview  %If preview active, compute wf & merge on first, last, and nth frames (always)
            if fileloop == n_end
                setappdata(0,'Save',1); %Trigger image save on last frame
                save(temp_runtime);     %Save variables to a temporary file for function passing
%                 fpalm_render_einzelreader(handles,h);
            elseif fileloop == n_start+fs1-1
                fs1=fs1+fs0;
                save(temp_runtime);
%                 fpalm_render_einzelreader(handles,h);
            elseif fileloop == n_start
                save(temp_runtime);
%                 fpalm_render_einzelreader(handles,h);
            end
        end
    end

    if ~preview %If preview is off, continue to wf & merge at end and save the result
        plot1 = 0;
        save(temp_runtime);
        setappdata(0,'Save',1);
        fpalm_render_einzelreader(handles,h);
    else        %Otherwise delete progress bar and handle stored in appdata
        delete(w);
        rmappdata(h,'PrBar_Handle');
    end

    if get(handles.oc_wf,'Value') %save widefield image if selected
        if (file_info.part_name(length(file_info.part_name)) == '_') %Avoid saving with repeated seperators
            wf_file = [outpath,'/',file_info.part_name,'wf.tif'];
        else
            wf_file = [outpath,'/',file_info.part_name,'_wf.tif'];
        end
        imwrite(uint8((image_sum/max(image_sum(:)))*255),wf_file,'Compression','none');
    end

%Save Variables
%----------------
    if (file_info.part_name(length(file_info.part_name)) == '_') %Avoid saving with repeated seperators
        out_file = [outpath,file_info.part_name,num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'.mat'];
    else
        out_file = [outpath,file_info.part_name,'_',num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'.mat'];
    end

    if ~exist(outpath,'dir') %Make sure the directory exists
        answer = questdlg('Output directory could not be accessed. Choose a new save location?','Warning!','Yes','Discard','Yes');

        if strcmp(answer,'Yes')
            outpath = uigetdir();
            if ~outpath
                return
            end
            if outpath(length(outpath))~='/'
                outpath = [outpath,'/'];
            end

            if (file_info.part_name(length(file_info.part_name)) == '_')
                out_file = [outpath,file_info.part_name,num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'.mat'];
            else
                out_file = [outpath,file_info.part_name,'_',num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'.mat'];
            end
        else
            return
        end
    end

    prec_edit     = get(handles.prec_edit,'String');                %<- These to be saved with the
    bkgn          = str2double(get(handles.bkgn_noise,'String'));   %<- mat file for loading purposes
    ppp           = str2double(get(handles.pix_to_pho,'String'));   %<-
    jshift        = str2double(get(handles.jshift,'String'));       %<-
    ishift        = str2double(get(handles.ishift,'String'));       %<-
    exf           = str2double(get(handles.exf_edit,'String'));     %<-   expansion factor

    a0_phot=a0_all/ppp;    % peak amplitude of each molecule in photons
    N=npix*a0_phot;        % number of photons for each molecule

    %localization precision in um
    lp2=((psf_std^2)+(q^2)/12)*1./N+8*pi*(psf_std^4)*(bkgn^2)/(q^2)*1./(N.*N);
    lp=sqrt(lp2);

    list = who; %Use regular expressions to ensure handles are not saved
    match = regexp(list,'^plot|handles|\<h\>|\<w\>');
    for i=length(match):-1:1
        if match{i}
            list(i)=[];
        end
    end
    save(out_file,list{:}); %Save variables used in einzelreader.m
    setappdata(0,'last_einzel_save',out_file);  %Save location of einzelreader mat file for other functions to use
%----------------
    %Clean up temp files and appdata
    if exist(temp_runtime,'file')
        delete(temp_runtime)
    end
    if isappdata(0,'temp_runtime')
        rmappdata(0,'temp_runtime')
    end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB