text:
% STD_ERSPPLOT - plot STUDY cluster ERSPs. Displays either mean cluster ERSPs,
%                  or else all cluster component ERSPs plus the mean cluster
%                  ERSP in one figure per condition. The ERSPs can be plotted
%                  only if component ERSPs were computed and saved in the
%                  EEG datasets in the STUDY. These may either be computed
%                  during pre-clustering using the gui-based function
%                  POP_PRECLUST, or via the equivalent commandline functions
%                  EEG_CREATEDATA and EEG_PRECLUST. Called by POP_CLUSTEDIT.
% Usage:
%   >> [STUDY] = std_erspplot(STUDY, ALLEEG, key1, val1, key2, val2);
%   >> [STUDY erspdata ersptimes erspfreqs pgroup pcond pinter] = ...
%                std_erspplot(STUDY, ALLEEG ...);
%
% Inputs:
%   STUDY    - STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG   - global vector of EEG structures for the datasets included
%              in the STUDY. ALLEEG for a STUDY set is typically created
%              using LOAD_ALLEEG.
%   either 'channels' or 'cluster' inputs are also mandatory.
%
% Optional inputs for channel plotting:
%   'channels' - [numeric vector]  specific channel group to plot. By
%                default, the grand mean channel ERSP is plotted (using the
%                same format as for the cluster component means described above)
%   'subject'  - [numeric vector]  In 'changrp' mode (above), index of
%                the subject(s) to plot. Else by default, plot all components
%                in the cluster.
%   'plotsubjects' - ['on'|'off'] When 'on', plot ERSP of all subjects.
%   'noplot'   - ['on'|'off'] When 'on', only return output values. Default
%                is 'off'.
%
% Optional inputs:
%   'clusters' - [numeric vector|'all'] indices of clusters to plot.
%                If no component indices ('comps' below) are given, the average
%                ERSPs of the requested clusters are plotted in the same figure,
%                with ERSPs for different conditions (and groups if any) plotted
%                in different colors. In 'comps' (below) mode, ERSP for each
%                specified cluster are plotted in separate figures (one per
%                condition), each overplotting cluster component ERSP plus the
%                average cluster ERSP in bold. Note this parameter has no effect
%                if the 'comps' option (below) is used. {default: 'all'}
%   'comps'    - [numeric vector|'all'] indices of the cluster components to plot.
%                Note that 'comps', 'all' is equivalent to 'plotsubjects', 'on'.
%
% Other optional inputs:
%   'plotmode'  - ['normal'|'condensed'|'none'] 'normal'  -> plot in a new figure;
%                 'condensed' -> plot all curves in the current figure in a
%                 condensed fashion. 'none' toggles off plotting {default: 'normal'}
%   'key','val' - All optional inputs to POP_SPECPARAMS are also accepted here
%                 to plot subset of time, statistics etc. The values used by default
%                 are the ones set using POP_SPECPARAMS and stored in the
%                 STUDY structure.
% Output:
%   STUDY      - the input STUDY set structure with the plotted cluster
%                mean ERSPs stored to allow quick replotting
%   erspdata   - [cell] ERSP data for each condition, group and subjects.
%                size of cell array is [nconds x ngroups]. Size of each element
%                is [freqs x times x subjects] for data channels or
%                [freqs x times x components] for component clusters. This
%                array may be gicen as input  directly to the STATCOND f
%                unction or STD_STATS function to compute statistics.
%   ersptimes  - [array] ERSP time point latencies.
%   erspfreqs  - [array] ERSP point frequency values.
%   pgroup     - [array or cell] p-values group statistics. Output of the
%                STATCOND function.
%   pcond      - [array or cell] condition statistics. Output of the STATCOND
%                function.
%   pinter     - [array or cell] groups x conditions statistics. Output of
%                STATCOND function.
%
% Example:
%        >> [STUDY] = std_erspplot(STUDY,ALLEEG, 'clusters', 'all', ...
%                                       'mode', 'together');
%           % Plot the mean ERSPs of all clusters in STUDY together
%           % on the same figure.
%
% Known limitations: when plotting multiple clusters, the output
%                    contains the last plotted cluster.
%
% See also: POP_CLUSTEDIT, POP_PRECLUST, EEG_CREATEDATA, EEG_PRECLUST, POP_CLUSTEDIT
%
% Authors: Arnaud Delorme, CERCO, August, 2006

% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [STUDY, allersp, alltimes, allfreqs, pgroup, pcond, pinter, events] = std_erspplot(STUDY, ALLEEG, varargin)

if nargin < 2
    help std_erspstatplot;
    return;
end

% find datatype and default options
% ---------------------------------
dtype = 'ersp';
for ind = 1:2:length(varargin)
    if strcmpi(varargin{ind}, 'datatype')
        dtype = varargin{ind+1};
    end
end
if strcmpi(dtype, 'erpim')
     STUDY  = pop_erpimparams(STUDY, varargin{:});
     params = STUDY.etc.erpimparams;
else STUDY  = pop_erspparams( STUDY, varargin{:});
     params = STUDY.etc.erspparams;
end

% get parameters
% --------------
statstruct.etc = STUDY.etc;
statstruct = pop_statparams(statstruct, varargin{:});

% potentially missing fields
% --------------------------
fields     = { 'freqrange'     [];
               'topofreq'      [];
               'topotrial'     [];
               'trialrange'    []
               'concatenate'   'off';
               'colorlimits'   [];
               'ersplim'       [];
               'itclim'        [];
               'maskdata'      'off';
               'subbaseline'   'off' };
for ind=1:size(fields,1)
    if ~isfield(params, fields{ind,1}),
        params = setfield(params, fields{ind,1}, fields{ind,2});
    end
end

% decode input parameters
% -----------------------
options = mystruct(varargin);
options = myrmfield( options, myfieldnames(params));
options = myrmfield( options, myfieldnames(statstruct.etc.statistics));
options = myrmfield( options, { 'threshold' 'statistics' } ); % for backward compatibility
[ opt, moreparams ] = finputcheck( options, { ...
                               'design'      'integer' [] STUDY.currentdesign;
                               'caxis'       'real'    [] [];
                               'statmode'    'string'  [] ''; % deprecated
                               'channels'    'cell'    []              {};
                               'clusters'    'integer' []              [];
                               'datatype'    'string'  { 'itc','ersp','pac' 'erpim' } 'ersp';
                               'plottf'      'real'    []              [];
                               'mode'        'string'  []              ''; % for backward compatibility (now used for statistics)
                               'comps'       {'integer','string'}  []              []; % for backward compatibility
                               'plotsubjects' 'string' { 'on','off' }  'off';
                               'noplot'      'string'  { 'on','off' }  'off';
                               'plotmode'    'string'  { 'normal','condensed','none' }  'normal';
                               'subject'     'string'  []              '' }, ...
                                  'std_erspstatplot', 'ignore');
if ischar(opt), error(opt); end
if strcmpi(opt.noplot, 'on'), opt.plotmode = 'none'; end
if isempty(opt.caxis)
    if strcmpi(opt.datatype, 'ersp')
         opt.caxis = params.ersplim;
    elseif strcmpi(opt.datatype, 'itc') && ~isempty(params.itclim)
        opt.caxis = [-params.itclim(end) params.itclim(end)];
    end
end

allconditions  = {};
allgroups      = {};
condname       = '';
groupname      = '';
if length(STUDY.design(opt.design).variable) > 0, allconditions = STUDY.design(opt.design).variable(1).value; condname  = STUDY.design(opt.design).variable(1).label; end
if length(STUDY.design(opt.design).variable) > 1, allgroups     = STUDY.design(opt.design).variable(2).value; groupname = STUDY.design(opt.design).variable(2).label; end

% for backward compatibility
% --------------------------
if strcmpi(opt.datatype, 'erpim')
    params.topofreq = params.topotrial;
    opt.caxis       = params.colorlimits;
    valunit = 'trials';
else
    valunit = 'Hz';
end
if isempty(opt.plottf) && ~isempty(params.topofreq) && ~isempty(params.topotime) && ~isnan(params.topofreq(1)) && ~isnan(params.topotime(1))
     params.plottf = [ params.topofreq(1) params.topofreq(end) params.topotime(1) params.topotime(end) ];
else params.plottf = opt.plottf;
end
%if strcmpi(opt.mode, 'comps'), opt.plotsubjects = 'on'; end %deprecated
stats = statstruct.etc.statistics;
stats.fieldtrip.channelneighbor = struct([]); % assumes one channel or 1 component
if isempty(STUDY.design(opt.design).variable)
    stats.paired = { };
else
    stats.paired = { STUDY.design(opt.design).variable(:).pairing };
end
if strcmpi(stats.singletrials, 'off') && ((~isempty(opt.subject) || ~isempty(opt.comps)))
    if strcmpi(stats.condstats, 'on') || strcmpi(stats.groupstats, 'on')
        stats.groupstats = 'off';
        stats.condstats  = 'off';
        disp('No statistics for single subject/component');
    end
end

if length(opt.comps) == 1
    stats.condstats = 'off'; stats.groupstats = 'off';
    disp('Statistics cannot be computed for single component');
end

alpha    = fastif(strcmpi(stats.mode, 'eeglab'), stats.eeglab.alpha, stats.fieldtrip.alpha);
mcorrect = fastif(strcmpi(stats.mode, 'eeglab'), stats.eeglab.mcorrect, stats.fieldtrip.mcorrect);
method   = fastif(strcmpi(stats.mode, 'eeglab'), stats.eeglab.method, ['Fieldtrip ' stats.fieldtrip.method ]);
plottfopt = { ...
   'ersplim',     opt.caxis, ...
   'threshold',   alpha, ...
   'effect',      stats.effect, ...
   'maskdata',    params.maskdata ...
   'averagemode'  params.averagemode };
if ~isempty(params.plottf) && length(opt.channels) < 5  && isempty(opt.clusters)
    warndlg2(strvcat('ERSP/ITC parameters indicate that you wish to plot scalp maps', 'Select at least 5 channels to plot topography'));
    allersp = {}; alltimes = []; allfreqs = []; pgroup = []; pcond = []; pinter = []; events  = [];
    return;
end

% plot single scalp map
% ---------------------
if ~isempty(opt.channels)
    if isempty(params.plottf) && length(opt.channels) > 1 && strcmpi(stats.singletrials, 'on')
        error('Cannot plot several channels on the same figure when using single trial statistics');
    end

    [STUDY, allersp, alltimes, allfreqs, events, paramsersp] = std_readdata(STUDY, ALLEEG, 'channels', opt.channels, 'timerange', params.timerange, ...
        'freqrange', params.freqrange, 'subject', opt.subject, 'singletrials', stats.singletrials, 'design', opt.design, 'datatype', opt.datatype, 'subbaseline', params.subbaseline);
    % 'concatenate', params.concatenate NOT TAKEN INTO ACCOUNT
    unitPower = newtimefpowerunit(paramsersp);

    if strcmpi(opt.datatype, 'ersp') && strcmpi(params.subbaseline, 'off')
        if  strcmpi(stats.singletrials, 'off')
            % rational for baseline
            % - no baseline calculation or log transformation at reading time (except single trial baseline if any)
            % - if not single trial and no common baseline, remove baseline and transform data here in each condition (before stats)
            % - otherwise, do so after baseline removal
            paramsersp.singletrials = stats.singletrials;
            paramsersp.commonbase   = params.subbaseline;
            [allersp,basesamples,basevals] = newtimefbaseln(allersp, alltimes, paramsersp);
        else
            opt.subbaseline = 'on';
            disp('Warning: when using single-trial statistics, a common baseline is forced across all conditions');
        end
    end

    %[STUDY allersp alltimes allfreqs tmp events unitPower] = std_readerp(STUDY, ALLEEG, 'channels', opt.channels, 'infotype', opt.datatype, 'subject', opt.subject, ...
    %    'singletrials', stats.singletrials, 'subbaseline', params.subbaseline, 'timerange', params.timerange, 'freqrange', params.freqrange, 'design', opt.design, 'concatenate', params.concatenate);
    %tic
    %[STUDY allersp alltimes allfreqs tmp events unitPower] = std_readersp(STUDY, ALLEEG, 'channels', opt.channels, 'infotype', opt.datatype, 'subject', opt.subject, ...
    %   'singletrials', stats.singletrials, 'subbaseline', params.subbaseline, 'timerange', params.timerange, 'freqrange', params.freqrange, 'design', opt.design, 'concatenate', params.concatenate);
    %toc

    % average single trials
    % ---------------------
    if strcmpi(opt.datatype, 'ersp')
        if  strcmpi(params.subbaseline, 'on')
            disp('Computing common baseline has changed since EEGLAB 14: averaging baselines is now');
            disp('performed before log-transformation of the baseline - in a similar way that baseline');
            disp('is averaged across trials (log transformation is only performed at the end for display)');
            % see above for rational for baseline
            paramsersp.singletrials = stats.singletrials;
            paramsersp.commonbase   = params.subbaseline;
            allersp = newtimefbaseln(allersp, alltimes, paramsersp);
        else
            paramsersp.singletrials = stats.singletrials;
            allersp = cellfun(@(x)newtimefbaseln(x, alltimes, paramsersp), allersp, 'uniformoutput', false);
        end
        % transform to log (except single trials where the transformation
        % is after taking the average - which is after doing stats
        if strcmpi(stats.singletrials, 'off')
            if ~isfield(paramsersp, 'scale') || strcmpi(paramsersp.scale, 'log')
                allersp = cellfun(@(x)10*log10(x), allersp, 'uniformoutput', false);
            end
        end
    end

    % select specific time and freq
    % -----------------------------
    if ~isempty(params.plottf)
        if length(params.plottf) < 3
            params.plottf(3:4) = params.plottf(2);
            params.plottf(2)   = params.plottf(1);
        end
        [~, fi1] = min(abs(allfreqs-params.plottf(1)));
        [~, fi2] = min(abs(allfreqs-params.plottf(2)));
        [~, ti1] = min(abs(alltimes-params.plottf(3)));
        [~, ti2] = min(abs(alltimes-params.plottf(4)));
        for index = 1:length(allersp(:))
            allersp{index} = mean(mean(allersp{index}(fi1:fi2,ti1:ti2,:,:),1),2);
            allersp{index} = reshape(allersp{index}, [1 size(allersp{index},3) size(allersp{index},4) ]);
        end

        % prepare channel neighbor matrix for Fieldtrip
        statstruct = std_prepare_neighbors(statstruct, ALLEEG, 'channels', opt.channels);
        stats.fieldtrip.channelneighbor = statstruct.etc.statistics.fieldtrip.channelneighbor;

        params.plottf = { params.plottf(1:2) params.plottf(3:4) };
        [pcond, pgroup, pinter] = std_stat(allersp, stats);
        if (~isempty(pcond) && length(pcond{1}) == 1) || (~isempty(pgroup) && length(pgroup{1}) == 1), pcond = {}; pgroup = {}; pinter = {}; end % single subject STUDY
    else
        [pcond, pgroup, pinter] = std_stat(allersp, stats);
        if (~isempty(pcond ) && (size( pcond{1},1) == 1 || size( pcond{1},2) == 1)) || ...
           (~isempty(pgroup) && (size(pgroup{1},1) == 1 || size(pgroup{1},2) == 1))
            pcond = {}; pgroup = {}; pinter = {};
            disp('No statistics possible for single subject STUDY');
        end % single subject STUDY
    end

    if strcmpi(stats.singletrials, 'on')
        % For ITC this is optional but it does not change anything
        if strcmpi(opt.datatype, 'ersp')
            if ndims(allersp{1}) == 4, for ind = 1:length(allersp(:)), allersp{ind} = mean(allersp{ind},4); end; end
            if ndims(allersp{1}) == 3, for ind = 1:length(allersp(:)), allersp{ind} = mean(allersp{ind},3); end; end
            if strcmpi(opt.datatype, 'ersp') && (~isfield(paramsersp, 'scale') || strcmpi(paramsersp.scale, 'log'))
                allersp = cellfun(@(x)10*log10(x), allersp, 'uniformoutput', false);
            end
        elseif strcmpi(opt.datatype, 'itc')
            if ~isfield(params, 'itctype'), params.itctype = 'phasecoher'; end
            for iDat = 1:length(allersp(:))
                allersp{iDat} = newtimefitc(allersp{iDat}, params.itctype);
                allersp{iDat} = abs(allersp{iDat});
            end
        end
    end

    % Average channels
    if ~strcmpi(params.averagechan, 'off') && length(opt.channels) > 1
        for index = 1:length(allersp(:))
            if strcmpi(params.averagemode, 'ave')
                allersp{index} = squeeze(mean(allersp{index},3));
            else
                disp('Computing RMS while preserving sign');
                tfsign  = sign(squeeze(mean(allersp{index},3)));
                allersp{index} = squeeze(sqrt(mean(allersp{index}.^2,3))).*tfsign;
            end
        end
    end

    % plot specific channel(s)
    % ------------------------
    if ~strcmpi(opt.plotmode, 'none')
        locsOri = eeg_mergelocs(ALLEEG.chanlocs);
        locs = locsOri(std_chaninds(STUDY, opt.channels));

        % in case channels are being averaged
        if ~strcmpi(params.averagechan, 'off') && length(opt.channels) > 1
            if length(opt.channels) ~= length(locsOri)
                chanlabels = { locs.labels };
                chanlabels(2,:) = {','};
                chanlabels(2,end) = {''};
                locs(1).labels = [ chanlabels{:} ];
            else
                locs(1).labels = 'All channels';
            end
            locs(2:end) = [];
        end

        if ~isempty(params.plottf) % incomtible with averagechan above
            alltitles = std_figtitle('threshold', alpha, 'mcorrect', mcorrect, 'condstat', stats.condstats, 'cond2stat', stats.groupstats, ...
                                     'statistics', method, 'condnames', allconditions, 'cond2names', allgroups, 'chanlabels', { locs.labels }, ...
                                     'subject', opt.subject, 'valsunit', { valunit 'ms' }, 'vals', params.plottf, 'datatype', upper(opt.datatype), ...
                                     'effect', stats.effect, 'factor1', condname, 'factor2', groupname);
            std_chantopo(allersp, 'groupstats', pgroup, 'condstats', pcond, 'interstats', pinter, 'caxis', opt.caxis, ...
                                          'chanlocs', locs, 'threshold', alpha, 'titles', alltitles);
        else
            if length(locs) > 1, opt.plottopo = 'on'; else opt.plottopo = 'off'; end
            if length(locs) == 1 && size(allersp{1},3) > 1
                % channels should be in 3rd dim; reshape data to put subjects in the 4th dim if number of channels is 1
                for index = 1:length(allersp(:))
                    allersp{index} = reshape(allersp{index}, size(allersp{index},1), size(allersp{index},2), 1, size(allersp{index},3));
                end
            end

%             nc = ceil(sqrt(length(opt.channels)));
%             nr = ceil(length(opt.channels)/nc);
%             for index = 1:length(locs)
%                 if length(opt.channels) > 1, try, subplot(nr,nc,index, 'align'); catch, subplot(nr,nc,index); end; end
%                 tmpersp = cell(size(allersp));
%                 for ind = 1:length(allersp(:))
%                     if ~isempty(allersp{ind})
%                         tmpersp{ind} = squeeze(allersp{ind}(:,:,index,:));
%                         tmpersp{ind} = permute(tmpersp{ind}, [2 1 3]); % somehow time/freq are swapped in ntimes = nfreqs
%                     end
%                 end
%                 for ind = 1:length(allersp(:))
%                     if ~isempty(allersp{ind})
%                         allersp{ind} = permute(allersp{ind}, [2 1 3 4]); % somehow time/freq are swapped for tftopo
%                     end
%                 end
                alltitles = std_figtitle('threshold', alpha, 'mcorrect', mcorrect, 'condstat', stats.condstats, 'cond2stat', stats.groupstats, ...
                    'statistics', method, 'condnames', allconditions, 'cond2names', allgroups, ...
                    'subject', opt.subject, 'datatype', upper(opt.datatype), 'plotmode', opt.plotmode, ...
                    'effect', stats.effect, 'factor1', condname, 'factor2', groupname);
                std_plottf(alltimes, allfreqs, allersp, 'datatype', opt.datatype, 'titles', alltitles, ...
                    'groupstats', pgroup, 'condstats', pcond, 'interstats', pinter, 'plottopo', opt.plottopo, 'plotmode', ...
                    opt.plotmode, 'unitcolor', unitPower, 'chanlocs', locs, 'events', events, plottfopt{:});
%             end
        end
    end
else
    if length(opt.clusters) > 1 && strcmpi(stats.singletrials, 'on')
        error('Cannot plot several components on the same figure when using single trial statistics');
    end

    if length(opt.clusters) > 1 && ~strcmpi(opt.plotmode, 'none'), figure; opt.plotmode = 'condensed'; end
    nc = ceil(sqrt(length(opt.clusters)));
    nr = ceil(length(opt.clusters)/nc);
    comp_names = {};

    if length(opt.clusters) > 1 && ( strcmpi(stats.condstats, 'on') || strcmpi(stats.groupstats, 'on'))
        stats.condstats = 'off'; stats.groupstats = 'off';
    end

    for index = 1:length(opt.clusters)

        [STUDY, allersp, alltimes, allfreqs, events, paramsersp] = std_readdata(STUDY, ALLEEG, 'clusters', opt.clusters(index), 'datatype', opt.datatype, ...
            'component', opt.comps, 'singletrials', stats.singletrials, 'subbaseline', params.subbaseline, 'timerange', params.timerange, 'freqrange', params.freqrange, 'design', opt.design, 'concatenate', params.concatenate);
        if length(opt.clusters) > 1, try, subplot(nr,nc,index, 'align'); catch, subplot(nr,nc,index); end; end
        unitPower = newtimefpowerunit(paramsersp);

        % plot specific component
        % -----------------------
        if ~isempty(opt.comps)
            comp_names = { STUDY.cluster(opt.clusters(index)).comps(opt.comps) };
            opt.subject = STUDY.datasetinfo(STUDY.cluster(opt.clusters(index)).sets(1,opt.comps)).subject;
        end

        % average single trials
        % ---------------------
        if strcmpi(opt.datatype, 'ersp')
            if  strcmpi(params.subbaseline, 'on')
                disp('Computing common baseline has changed since EEGLAB 14: averaging baselines is now');
                disp('performed before log-transformation of the baseline - in a similar way that baseline');
                disp('is averaged across trials (log transformation is only performed at the end for display)');
                % see above for rational for baseline
                paramsersp.singletrials = stats.singletrials;
                paramsersp.commonbase   = params.subbaseline;
                allersp = newtimefbaseln(allersp, alltimes, paramsersp);
            else
                paramsersp.singletrials = stats.singletrials;
                allersp = cellfun(@(x)newtimefbaseln(x, alltimes, paramsersp), allersp, 'uniformoutput', false);
            end
            % transform to log (except single trials)
            if strcmpi(stats.singletrials, 'off')
                if ~isfield(paramsersp, 'scale') || strcmpi(paramsersp.scale, 'log')
                    allersp = cellfun(@(x)10*log10(x), allersp, 'uniformoutput', false);
                end
            end
        end

        % statistics
        % ----------
        [pcond pgroup pinter] = std_stat(allersp, stats);
        if (~isempty(pcond ) && (size( pcond{1},1) == 1 || size( pcond{1},2) == 1)) || ...
           (~isempty(pgroup) && (size(pgroup{1},1) == 1 || size(pgroup{1},2) == 1)),
            pcond = {}; pgroup = {}; pinter = {};
            disp('No statistics possible for single subject STUDY');
        end % single subject STUDY

        if strcmpi(stats.singletrials, 'on')
            if strcmpi(opt.datatype, 'ersp')
                if ndims(allersp{1}) == 4, for ind = 1:length(allersp(:)), allersp{ind} = mean(allersp{ind},4); end; end
                if ndims(allersp{1}) == 3, for ind = 1:length(allersp(:)), allersp{ind} = mean(allersp{ind},3); end; end
                if ~isfield(paramsersp, 'scale') || strcmpi(paramsersp.scale, 'log')
                    allersp = cellfun(@(x)10*log10(x), allersp, 'uniformoutput', false);
                end
            elseif strcmpi(opt.datatype, 'itc')
                if ~isfield(params, 'itctype'), params.itctype = 'phasecoher'; end
                for iDat = 1:length(allersp(:))
                    allersp{iDat} = newtimefitc(allersp{iDat}, params.itctype);
                    allersp{iDat} = abs(allersp{iDat});
                end
            end
        end

        % plot specific component
        % -----------------------
        if index == length(opt.clusters), opt.legend = 'on'; end
        if ~strcmpi(opt.plotmode, 'none')
            alltitles = std_figtitle('threshold', alpha, 'mcorrect', mcorrect, 'condstat', stats.condstats, 'cond2stat', stats.groupstats, ...
                                     'statistics', method, 'condnames', allconditions, 'cond2names', allgroups, 'clustname', STUDY.cluster(opt.clusters(index)).name, 'compnames', comp_names, ...
                                     'subject', opt.subject, 'datatype', upper(opt.datatype), 'plotmode', opt.plotmode, ...
                                     'effect', stats.effect, 'factor1', condname, 'factor2', groupname);

            std_plottf(alltimes, allfreqs, allersp, 'datatype', opt.datatype, ...
                                           'groupstats', pgroup, 'condstats', pcond, 'interstats', pinter, 'plotmode', ...
                                           opt.plotmode, 'titles', alltitles, ...
                                          'events', events, 'unitcolor', unitPower, 'chanlocs', ALLEEG(1).chanlocs, plottfopt{:});
        end
    end
end

% remove fields and ignore fields who are absent
% ----------------------------------------------
function s = myrmfield(s, f);

for index = 1:length(f)
    if isfield(s, f{index})
        s = rmfield(s, f{index});
    end
end

% convert to structure (but take into account cells)
% --------------------------------------------------
function s = mystruct(v);

for index=1:length(v)
    if iscell(v{index})
        v{index} = { v{index} };
    end
end
try
    s = struct(v{:});
catch, error('Parameter error'); end

% convert to structure (but take into account cells)
% --------------------------------------------------
function s = myfieldnames(v);

s = fieldnames(v);
if isfield(v, 'eeglab')
    s2 = fieldnames(v.eeglab);
    s = { s{:} s2{:} };
end
if isfield(v, 'fieldtrip')
    s3 = fieldnames(v.fieldtrip);
    for index=1:length(s3)
        s3{index} = [ 'fieldtrip' s3{index} ];
    end
    s = { s{:} s3{:} };
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [y,dy] = affine3D(w,x,varargin)
%
% computes y = Q*w and the derivative wrt. w.
% x = reshape(x,[],3);
% Q = [x(:,1),x(:,2),x(:,3),1, 0      0        0      0  0      0        0
%      0      0      0,     0, x(:,1),x(:,2,1),x(:,3),1, 0      0        0
%      0      0      0,     0, 0      0        0      0, x(:,1),x(:,2,1),x(:,3),1]
% dy = Q:
% if no arguments are given, the parameters for the identity map are returned.
%
% see also transformations/contents.m, trafo.m
%==============================================================================

function [y,dy] = affine3D(w,x,varargin)

% the persitent variable stores the matrix
% Q(x) = kron( I_2 , [x(:,1),x(:,2),1] );
persistent Q

if nargin==0
    help(mfilename)
    runMinimalExample;
    return;
else
  y = mfilename('fullfile');
  dy = reshape(eye(4,3),[],1);         % parameterization of identity
  if ischar(w),  Q  = []; w = []; end; % reset Q
  if isempty(w), return;          end;
end;

if isempty(w) || (size(Q,1) ~= numel(x)),
  n = length(x)/3; x = reshape(x,n,3);
  Q = sparse(kron(speye(3),[x,ones(n,1)]));
  if nargout == 0, return; end;
end;
y  = Q*w;
dy = Q;

%------------------------------------------------------------------------------
function runMinimalExample
fprintf('%s: minimal example\n',mfilename)

omega = [0,10,0,8,0,6]; m = [8,7,6];
w = 22/pi;c = (omega(2:2:end)-omega(1:2:end))'/2;
R = [ cos(w),-sin(w),0;sin(w),cos(w),0;0,0,1];
g = (eye(3)-R)*reshape(c,[],1);
w = reshape([R,g]',[],1)
x = getNodalGrid(omega,m);
z = feval(mfilename,w,x);
FAIRfigure(1); clf;
plotGrid(x,omega,m,'color','r'); axis image; hold on;
plotGrid(z,omega,m,'color','b'); axis image; hold off;

fctn = @(w) feval(mfilename,w,x);
w = w + randn(size(w));
checkDerivative(fctn,w,'fig',2);

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


text:
function figure_size = getnicedialoglocation(figure_size, figure_units)
% adjust the specified figure position to fig nicely over GCBF
% or into the upper 3rd of the screen

%  Copyright 1999-2010 The MathWorks, Inc.

parentHandle = gcbf;
convertData.destinationUnits = figure_units;
if ~isempty(parentHandle)
    % If there is a parent figure
    convertData.hFig = parentHandle;
    convertData.size = get(parentHandle,'Position');
    convertData.sourceUnits = get(parentHandle,'Units');
    c = [];
else
    % If there is no parent figure, use the root's data
    % and create a invisible figure as parent
    convertData.hFig = figure('visible','off');
    convertData.size = get(0,'ScreenSize');
    convertData.sourceUnits = get(0,'Units');
    c = onCleanup(@() close(convertData.hFig));
end

% Get the size of the dialog parent in the dialog units
container_size = hgconvertunits(convertData.hFig, convertData.size ,...
    convertData.sourceUnits, convertData.destinationUnits, get(convertData.hFig,'Parent'));

delete(c);

figure_size(1) = container_size(1)  + 1/2*(container_size(3) - figure_size(3));
figure_size(2) = container_size(2)  + 2/3*(container_size(4) - figure_size(4));

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
% The COBRAToolbox: testTheoretMaxProd.m
%
% Purpose:
%     - test the theoretMaxProd function
%
% Authors:
%     - Jacek Wachowiak
global CBTDIR
% save the current path
currentDir = pwd;

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

% test variables
model = getDistributedModel('ecoli_core_model.mat');

% change solver since qpng is unstable - to be changed after installation of gurobi
% TODO: test this for multiple solvers.
changeCobraSolver('pdco', 'QP');

%Get the molecular weights of the compounds in the model.
[mw,EMatrix,elements] = computeMW(model, [], false);

%lets calc a simple optimization
tempModel = model;
acExPos = ismember(model.rxns,'EX_ac(e)');
tempModel.c = double(acExPos);
FBAsol = optimizeCbModel(tempModel);

% function calls
[ExRxns, MaxTheoOut] = theoretMaxProd(model, 'EX_glc(e)'); % Max 'EX_ac(e)' flux (unique)
allExPos = ismember(model.rxns,ExRxns);
%Determine the Exchangers and their respective exchanged metabolite +
%Weights
exWeights = cellfun(@(x) mw(find(model.S(:,ismember(model.rxns,x)))),ExRxns);
carbEx = cellfun(@(x) EMatrix(find(model.S(:,ismember(model.rxns,x))),ismember(elements,'C')) > 0 ,ExRxns);
carbExPos = ismember(model.rxns, ExRxns(carbEx));
acPosInExch = ismember(ExRxns,'EX_ac(e)');

%Alternative calls
[ExRxns1, MaxTheoOut1] = theoretMaxProd(model, 'EX_glc(e)', '', true, findExcRxns(model,0,0)); % Scaled to 0-1
[ExRxns2, MaxTheoOut2] = theoretMaxProd(model, 'EX_glc(e)', 'pr_mw'); %Compared to other mw
[ExRxns3, MaxTheoOut3] = theoretMaxProd(model, 'EX_glc(e)', 'pr_other_mol');
[ExRxns4, MaxTheoOut4] = theoretMaxProd(model, 'EX_glc(e)', 'pr_other_mw');
[ExRxns5, MaxTheoOut5] = theoretMaxProd(model, 'EX_glc(e)', 'x'); % bad criterion

% We assume, that the solver is deterministic, i.e. the solution found for
% the basic call is the same as for all the other calls.
assert(abs(MaxTheoOut(1) - 20) < 1e-4); % 'EX_ac(e)'
assert(abs(MaxTheoOut1(1) - 2) < 1e-4); % the same but scaled to glucose uptake flux (10)

assert(abs(MaxTheoOut2(1) - FBAsol.x(acExPos)*exWeights(1)) < 1e-4); % 'EX_ac(e)' in weight
%All (molar) exports (>0) without acetate and with
assert(abs(MaxTheoOut3(1) - sum(FBAsol.x(allExPos & ~acExPos & FBAsol.x > 0 & carbExPos) )) < 1e-4); % 'EX_ac(e)'
% Sum of Molecular weights of all Carbon Exporters ( > 0) without Acetate
assert(abs(MaxTheoOut4(1) - sum(FBAsol.x(allExPos).*(FBAsol.x(allExPos)>0).*exWeights.*carbEx.*~acPosInExch)) < 1e-4); % 'EX_ac(e)' in weight yield
assert(isequal(MaxTheoOut5, zeros(20, 1))); % bad criterion returns only 0s

% change to old directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function tracking_env()
addpath('../matconvnet/matlab');
run('vl_setupnn.m') ;
fftw('planner','patient');
end

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet


text:
function [A,B,C,D,E,F,G] = bay_rr(X,Y,gam,level,nb,eigvals,eigvec)
% Bayesian inference of the cost on the three levels of linear ridge regression
%
% >> cost = bay_rr(X, Y, gam, level)
%
% This function implements the cost functions related to the
% Bayesian framework of linear ridge Regression [29]. Optimizing
% this criteria results in optimal model parameters W,b,
% hyperparameters. The criterion can also be used for model
% comparison.
%
% The obtained model parameters w and b are optimal on the first
% level w.r.t J = 0.5*w'*w+gam*0.5*sum(Y-X*w-b).^2.
%
% Full syntax
%
%     * OUTPUTS on the first level: Cost proportional to the posterior of the model parameters.
%
% >> [costL1, Ed, Ew] = bay_rr(X, Y, gam, 1)
%
%         costL1: Cost proportional to the posterior
%         Ed(*) : Cost of the fitting error term
%         Ew(*) : Cost of the regularization parameter
%
%     * OUTPUTS on the second level: Cost proportional to the posterior of gam.
%
% >> [costL2, DcostL2, Deff, mu, ksi, eigval, eigvec] = bay_rr(X, Y, gam, 2)
%
%         costL2    : Cost proportional to the posterior on the second level
%         DcostL2(*): Derivative of the cost proportional to the posterior
%         Deff(*)   : Effective number of parameters
%         mu(*)     : Relative importance of the fitting error term
%         ksi(*)    : Relative importance of the regularization parameter
%         eigval(*) : Eigenvalues of the covariance matrix
%         eigvec(*) : Eigenvectors of the covariance matrix
%
%     * OUTPUTS on the third level: The following commands can be
%     used to compute the level 3 cost function for different
%     models (e.g. models with different selected sets of
%     inputs). The best model can then be chosen as the model with
%     best level 3 cost (CostL3).
%
% >> [costL3, gam_optimal] = bay_rr(X, Y, gam, 3)
%
%         costL3         : Cost proportional to the posterior on the third inference level
%         gam_optimal(*) : Optimal regularization parameter obtained from optimizing the second level
%
%     * INPUTS:
%
% >> cost = bay_rr(X, Y, gam, level)
%
%         X N x d matrix with the inputs of the training data
%         Y N x 1 vector with the outputs of the training data
%         gam Regularization parameter
%         level 1, 2, 3
%
% See also:
%   ridgeregress,bay_lssvm

% -dislaimer-

  if ~exist('fminunc'),
    error('This function needs the optimization function ''fminunc''.');
  end

[N,d] = size(X);

if (level==1 || level==2),
  [W,b] = ridgeregress(X,Y,gam);
  Ew =  .5*W'*W;             % Ew
  Ed = .5*sum((X*W+b-Y).^2); % Ed
  Ewgd = Ew+gam*Ed;          % Ewgd
  A=Ewgd; C=Ed; B=Ew;

  if level==2,
    if nargin>=7,
      if numel(eigvals)>length(eigvals), v = diag(eigvals); else v=eigvals; end
      V = eigvec;
    else
      eval('[V,v] = eigs(X''*X+eye(d)*2,nb);','[V,v] = eig(X''*X+eye(d)*2);v=diag(v);');v=v-2;
      v = v*(N-1)/N;
    end
    Peff = find(v>1000*eps); Neff=length(Peff);
    v = v(Peff); V= V(:,Peff);
    vall = zeros(N-1,1);vall(1:Neff)=v;
    Deff = 1+sum(gam.*v./(1+gam.*v));                 % Deff
    CostL2 = sum(log(vall+1./gam))+(N-1)*log(Ewgd);   % CostL2
    DcostL2 = -sum(1./(gam+vall.*gam^2))+(N-1)*Ed/Ewgd; % DcostL2
    mu = 2*Ed/(N-Deff); ksi = mu*gam;                 % mu and ksi
    A=CostL2; B=DcostL2; C=Deff; D=mu; E=ksi;F=v;G=V;
  end

elseif level==3,

  % check fminunc
  resp = which('fminunc');
  %disp(' ');
  if isempty(resp),
    error(' ''fminunc'' not available');
  end

  eval('nb;','nb=''blabla'';');
  opties=optimset('MaxFunEvals', 2000,'GradObj','on', 'DerivativeCheck', 'off', 'TolFun', .0001, 'TolX', .0001, 'Display','off' );
  [CostL2,DCostL2, Deff, mu,ksi,v, V] = bay_rr(X,Y,gam,2,nb);
  gam_opt = exp(fminunc(@costL2, log(gam), opties,X,Y,nb,v,V));
  [CostL2,DCostL2, Deff, mu,ksi,v, V] = bay_rr(X,Y,gam_opt,2,nb,v,V);
  CostL3 = .5*length(v)*log(mu)+(N-1)*log(ksi) - log(Deff-1)-log(N-Deff) - sum(log(mu+ksi*v));
  A = CostL3; B = gam_opt;

else
  error('level should be ''1'', ''2'' or ''3''.');
end

function [C,Dc] = costL2(log_gam,X,Y,nb,v,V)
%
[C,Dc] = bay_rr(X,Y,exp(log_gam),2,nb,v,V);

author: peterhcharlton

repo: RRest


text:
function Y = invariantMultiply( psfMatData, X, padsize );
%
%           Y = invariantMultiply( psfMatData, X, padsize );
%
%  This function computes the multiplication of a spatially invariant
%  point spread function (PSF) times an image (i.e., convolution):
%                y = A*x
%
%  Here we assume A is made up of a single PSF whose extent may be
%  (much) smaller than the image.
%
%  Input:
%   psfMatData  -  complex array containing the matrix data, usually
%                  computed from onePsfMatrix.m
%            X  -  array containing the image to which the psfMatrix
%                  is to be multiplied.
%
%  Output:
%            Y  -  contains the result after PSF multiplication.
%

%  J. Nagy  1/7/02

imsize = size( X ) - 2*padsize;

%
%  In order for this to be consistent for 2-D and 3-D images, we need to make
%  sure there is a third dimension ...
%
if length(imsize) == 1
  imsize = [imsize, 1, 1];
  padsize = [padsize, 0, 0];
elseif length(imsize) == 2
  imsize = [imsize, 1];
  padsize = [padsize, 0];
end

%
%  partition_info computes number of subregions, and their sizes
%
[nregions, rsize] = partition_info(imsize, padsize);

%
%  Coding the rest of this will be easier if all of the image subregions
%  have the same dimensions.  If it's not, we pad with a few zeros to make
%  it so ...
%
padsize1 = rsize .* nregions - imsize;
if any( padsize1 < 0 )
  error('Something is wrong here ...')
end
X = padarray(X, padsize1, 'post');

Y = zeros( size(X) );

%
%  Now we get information about beginning and ending indices of subregions
%  so we can "put" and "get" subregions correctly ...
%
[RIidx, RJidx, RKidx] = region_indices( nregions, rsize );
[EIidx, EJidx, EKidx] = eregion_indices( RIidx, RJidx, RKidx, 2*padsize );
Tidx = [padsize+1; padsize+rsize];

%
%  Now loop over all the subregions ...
%
for k = 1:nregions(3)
  for j = 1:nregions(2)
    for i = 1:nregions(1)
      Xt = X(EIidx(i,1):EIidx(i,2), EJidx(j,1):EJidx(j,2), EKidx(k,1):EKidx(k,2));

      Yt = multiplyOneRegion( psfMatData, Xt );

      Y(RIidx(i,1):RIidx(i,2), RJidx(j,1):RJidx(j,2), RKidx(k,1):RKidx(k,2)) = ...
           Yt(Tidx(1,1):Tidx(2,1), Tidx(1,2):Tidx(2,2), Tidx(1,3):Tidx(2,3));
    end
  end
end

Y = Y(1:imsize(1), 1:imsize(2), 1:imsize(3));

repo: IRtools


text:
function leap = hammersley_leap_get ( )

%*****************************************************************************80
%
%% HAMMERSLEY_BASE_GET gets the leap vector of the leaped Hammersley subsequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer LEAP(1:DIM_NUM), the Hammersley leap vector.
%
  global hammersley_BASE
  global hammersley_LEAP
  global hammersley_DIM_NUM
  global hammersley_SEED
  global hammersley_STEP

  dim_num = hammersley_DIM_NUM;

  if ( ~halham_dim_num_check ( dim_num ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HAMMERSLEY_LEAP_GET - Fatal error!\n' );
    error ( 'HAMMERSLEY_LEAP_GET - Fatal error!' );
  end

  if ( length ( hammersley_LEAP ) ~= dim_num )

    hammersley_LEAP = [];

    for i = 1 : dim_num
      hammersley_LEAP(i) = 1;
    end

  end

  leap(1:dim_num) = hammersley_LEAP(1:dim_num);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [gk] = tapas_trans_mv2gk(mu, sigma2)
%% Transforms to the scaling of the gamma distribution%
% Input
%   mu          -- Means
%   sigma2      -- Variance
%
% Output
%   gk          -- K parameter of an inverse gamma distribution

% aponteeduardo@gmail.com
% copyright (C) 2015
%

gk = mu.^2 ./ sigma2;

end

author: translationalneuromodeling

repo: tapas


text:
filename = 'ArchSym';
%filename = 'CantileverSquareSYmmetricMesh';
%filename = 'CantileverSquareMedium';
%'CantileverSquareSYmmetricMesh';
%filename = 'CantileverSquare';
%filename = 'CantileverSquareSmall';
%filename = 'CantileverSquareNew';
%'CantileverSquareNewFine';
%filename = 'Cantilever_quad_coarse';
%filename = 'LshapeTriFine';
%filename = 'LshapeTri';
%filename = 'LshapeTriSmall';
%filename = 'Lshape';
%filename = 'LshapeFine';
%filename = 'ArchTriFine';
%'Arch_quad_coarse';
%'Bridge_quad_coarse';
%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';
%filename = 'Bridge';

ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'PDE';
%filterType = 'P1';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-4;
constr_initial = 1e-8;

Vfrac_final = 0.3;
optimality_final = 1e-4;
constr_final = 1e-8;

stressNormExponent_initial = 2;
stressNormExponent_final = 16;
%
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';
optimizerUnconstrained = 'PROJECTED GRADIENT';

%optimizer = 'MMA';
%optimizer = 'IPOPT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
%vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
vademecumFileName = 'SuperEllipseQOptAnalytic';

% designVariable = 'LevelSet';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';
% initial_case = 'full';% optimizerUnconstrained = 'SLERP';

% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';
% initial_case = 'full';
% rho0 = 0.3;

line_search_initiator = 'INCREASING LAST STEP';
incrementFactor = 1.95;
%

%kfrac = 2;
nsteps = 8;%17;

plotting = true;
printing = true;
monitoring = true;
monitoring_interval = 3;
maxiter = 8000;

%

isDirichletPartX = @(x) x > 0.2 & x < 0.4;
isDirichletPartY = @(y) y > 0 & y < 0.04;
isDirichletPart = @(x,y) isDirichletPartX(x) & isDirichletPartY(y);
isNeumannPartX = @(x) x > 0.9 & x < 1;
isNeumannPartY = @(y) y > 0 & y < 0.04;
isNeumannPart = @(x,y) isNeumannPartX(x) & isNeumannPartY(y);

iNotOptimizable = @(coord) isDirichletPart(coord(:,1),coord(:,2)) | isNeumannPart(coord(:,1),coord(:,2));

costDomainNotOptimizable       = {iNotOptimizable};
constraintDomainNotOptimizable = {[]};

isDesignVariableFixed.nodes  = iNotOptimizable;
isDesignVariableFixed.values = @(x) [m1*ones(size(x,1),1);m2*ones(size(x,1),1)];

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% 埃特金方法加速迭代方法求解方程
clear;
format long;
tol = 1e-10;
N = 100;
x0 = 0.5;
phi = @(x) exp(-x);
%phi = @(x) (x+1)^(1/3);

for k = 1 : N
    x1 = phi(x0);
    x2 = phi(x1);
    x2 = x2 - (x2 - x1)^2 / (x2 - 2 * x1 + x0);
    if abs(x2 - x0) < tol
        fprintf('方程的正根: %10.8f\n', x1);
        break;
    end
    x0 = x2;
end
if k == N
    fprintf('迭代方法失败\n');
end
fprintf('迭代次数: %d\n', k);

repo: NumericalAnalysis


text:
function [wX, wY, xdim, ydim, n] = windowizeNARX(X,Y,x_delays,y_delays, steps)
% Re-arrange the data points into a block Hankel matrix for (N)ARX time-series modeling
%
% >> [Xw,Yw] = windowizeNARX(X,Y,xdelays, ydelays, steps)
%
%  Rearrange the points of X and Y in a regressor matrix of the
%  past inputs and outputs (Xw) and the future outputs (Yw).
%
% Full syntax
%
% >> [Xw, Yw, xdim, ydim, n] = windowizeNARX(X, Y, xdelays, ydelays, steps)
%
%       Outputs
%         Xw Matrix of the data used for input including the delays
%         Yw Matrix of the data used for output including the next steps
%         xdim(*) Number of dimensions in new input
%         ydim(*) Number of dimensions in new output
%         n(*) Number of new data points
%       Inputs
%         X       : N x m vector with input data points
%         Y       : N x d vector with output data points
%         xdelays : Number of lags of X in new input
%         ydelays : Number of lags of Y in new input
%         steps(*): Number of future steps of Y in new output (by default 1)
%
% See also:
%   windowize, predict, trainlssvm, simlssvm

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

m=max(x_delays,y_delays);

eval('steps;','steps = 1;');
if steps == 0,
  n = size(X,1)-m;
else
  n = size(X,1)-m -steps+1;
end

wX = zeros(n,size(X,2)*(x_delays+1)+size(Y,2)*(y_delays));
wY = zeros(n,size(Y,2)*steps);
xdim = size(wX,2);
ydim = size(wY,2);

hdx  = (x_delays+1)*size(X,2);

for t=1:n,
  for i=1:x_delays+1,
    wX(t,1+((i-1)*size(X,2):i*size(X,2)-1)) = X(t+m-x_delays+i-1,:);
  end

  for i=1:y_delays,
    wX(t,hdx + (i-1)*size(Y,2) + (1:size(Y,2))) = Y(t+m-y_delays+i-1,:);
  end

  for i=1:steps,
    wY(t,i:i+size(Y,2)-1) = Y(t+m+i-1,:);
  end

end


text:
% compareBoundariesScript

savedir = '~/data/occlusion/labelme/figs/bnd_im';

folder = db(f).annotation.folder;
imfn = db(f).annotation.filename;
bn = strtok(imfn, '.');

im = imread(fullfile(imdir, folder, imfn));
figure(1), hold off, imagesc(im); axis image

load(fullfile(basedir, 'pb2', folder, [bn '_pb2_ucm']), 'ucm');
im = im2double(imresize(im, size(ucm), 'bilinear'));
imwrite(im, fullfile(savedir, imfn), 'Quality', 95);
grayim = rgb2gray(im2double(im));
h = 0*ones(size(grayim));
ucm = im2double(ucm);
pb = ucm / max(ucm(:));
tmpim = ordfilt2(max(pb, 0), 9, ones(3));
dispim = tmpim;%dispim = hsv2rgb(h, tmpim, 0.5*grayim+0.5*tmpim);
figure(2), hold off,  imagesc(dispim), axis image, colormap gray
imwrite(dispim, fullfile(savedir, [bn '_pb2_ucm.jpg']), 'Quality', 95);

load(fullfile(basedir, 'pb', folder, [bn '_pb']), 'pb');
grayim = rgb2gray(im2double(im));
h = 0*ones(size(grayim));
pb = pb / max(pb(:));
tmpim = ordfilt2(max(pb, 0), 9, ones(3));
dispim = tmpim;%dispim = hsv2rgb(h, tmpim, 0.5*grayim+0.5*tmpim);
figure(5), hold off,  imagesc(dispim), axis image
imwrite(dispim, fullfile(savedir, [bn '_pb.jpg']), 'Quality', 95);

load(fullfile(basedir, 'pb2', folder, [bn '_pb2']), 'pb');
grayim = rgb2gray(im2double(im));
h = 0*ones(size(grayim));
pb = pb / max(pb(:));
tmpim = ordfilt2(max(pb, 0), 9, ones(3));
dispim = tmpim;%dispim = hsv2rgb(h, tmpim, 0.5*grayim+0.5*tmpim);
figure(5), hold off,  imagesc(dispim), axis image
imwrite(dispim, fullfile(savedir, [bn '_pb2.jpg']), 'Quality', 95);

load(fullfile(occdir, folder, [bn '_occlusion']), 'bndinfo_all');
occ.po_all = getOcclusionMaps(bndinfo_all);
occim = mean(occ.po_all, 3);
%occim = im2double(imread(fullfile(basedir, 'results', folder, [strtok(imfn, '.') '_occ.jpg'])));
occim = occim/max(occim(:));
%h=0.65*ones(size(grayim));
tmpim = ordfilt2(occim, 9, ones(3));
dispim = tmpim;%dispim = hsv2rgb(h, tmpim, 0.5*grayim+0.5*tmpim);
figure(3), hold off, imagesc(dispim); axis image
imwrite(dispim, fullfile(savedir, [bn '_occave.jpg']), 'Quality', 95);

occim = occ.po_all(:, :, 1);
occim = occim/max(occim(:));
tmpim = ordfilt2(occim, 9, ones(3));
dispim = tmpim;%dispim = hsv2rgb(h, tmpim, 0.5*grayim+0.5*tmpim);
figure(6), hold off, imagesc(dispim); axis image
imwrite(dispim, fullfile(savedir, [bn '_occ1.jpg']), 'Quality', 95);

load(fullfile(basedir, 'labels2', folder, [bn '_labels']));
bmap = seg2bmap(lim, size(pb, 2), size(pb, 1));
bmap(:, [1:10 end-9:end]) = 0;
bmap([1:10 end-9:end], :) = 0;
tmpim = ordfilt2(double(bwmorph(bmap,'thin',inf)), 25, ones(5));
tmpim = cat(3, tmpim>0, zeros([size(tmpim) 2]));
dispim = im2double(im);
dispim(repmat(tmpim(:, :, 1)>0, [1 1 3])) = tmpim(repmat(tmpim(:, :, 1)>0, [1 1 3]));
%dispim = min(max(hsv2rgb(h, tmpim, 0.5*grayim+0.5*tmpim), 0), 1);
figure(4), hold off, imagesc(dispim); axis image,
imwrite(dispim, fullfile(savedir, [bn '_gt.jpg']), 'Quality', 95);

fid = fopen(fullfile(savedir, [bn '_ap.txt']), 'w');
fprintf(fid, 'index = %f   pb = %f   pb2+ucm = %f   occ1 = %f   occave = %f \n', sind(f), pr_pb(f).ap, pr_ucm(f).ap, pr_occ1(f).ap, pr_occave(f).ap);
fclose(fid);

disp([imfn ':  ' num2str([pr_pb2(f).ap pr_occave(f).ap])])

% grayim = imresize(rgb2gray(im2double(im)), size(pb), 'bilinear');
% resultim = repmat(grayim, [1 1 3]);
% tmpim = ordfilt2(max(pb,0), 9, ones(3));
% resultim((tmpim>0.05)) = tmpim(tmpim>0.05);
% tmpim = ordfilt2(occim, 9, ones(3));
% resultim(find(tmpim>0.05)+numel(tmpim)) = tmpim(tmpim>0.05);
% figure(4), hold off, imagesc(resultim), axis image

author: Cloud-CV

repo: object-proposals


text:
classdef nme_load_acc_nln < mp.nme_load_acc & mp.nme_wrapper_ac_nln

%   MATPOWER
%   Copyright (c) 2019-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        nme_class = @mp.nme_load_acc;
    end

    methods
        function obj = nme_load_acc_nln()
            obj@mp.nme_load_acc();
            obj.nme_wrapper_ac_nln_init();
        end

        function build_params(obj, nm, dm)
            build_params@mp.nme_load_acc(obj, nm, dm);
            obj.build_nln_params(nm, dm);
        end

        function nk = count(obj, dm)
            obj.nme.count(dm);
            nk = count@mp.nme_load_acc(obj, dm);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [vals, locs] = minandmax3(f)
%MINANDMAX3   returns the global minimum and maximum value of a CHEBFUN3.
%   VALS = minandmax3(F) returns the global minimum and maximum value of a
%   CHEBFUN3 object F over its domain. VALS is a vector of length 2 such
%   that Y(1) = min3(F(x,y,z)) and Y(2) = max3(F(x,y,z)).
%
%   [VALS, LOCS] = minandmax3(F) also returns the position of the global
%   minimum and maximum.
%
% See also CHEBFUN3/MAX, CHEBFUN3/MAX2, CHEBFUN3/NORM.

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

% check for empty CHEBFUN3.
if ( isempty(f) )
    vals = [];
    locs = [];
    return
end

doNewton = 1;   % Newton polishing?

% Extract low rank representation of f:
[fCore, fCols, fRows, fTubes] = tucker(f);
dom = f.domain;

% Is f the zero function?
if ( iszero(f) )
    locs = [(dom(1) + dom(2))/2  (dom(3) + dom(4))/2  (dom(5) + dom(6))/2];
    locs = [locs; locs];
    vals = [0; 0];
end

[m, n, p] = length(f);
if ( ndf(f) > 5e4 )
    m = min(m, 129);
    n = min(n, 129);
    p = min(p, 129);
end

% We seek a fast initial guess. So we first discretize the object.
xpts = chebpts(m, fCols.domain);
ypts = chebpts(n, fRows.domain);
zpts = chebpts(p, fTubes.domain);
colVals = feval(fCols, xpts);
rowVals = feval(fRows, ypts);
tubeVals = feval(fTubes, zpts);
T = chebfun3.txm(chebfun3.txm(chebfun3.txm(fCore,colVals,1), rowVals,2), ...
    tubeVals,3);

% Minimum entry in the sample:
[ignored, ind] = min(T(:));
[col, row, tub] = ind2sub(size(T), ind);
loc(1,1) = xpts(col);
loc(1,2) = ypts(row);
loc(1,3) = zpts(tub);
vals(1) = feval(f, loc(1, 1), loc(1, 2), loc(1, 3));

% Maximum entry in the sample:
[ignored, ind] = max(T(:));
[col, row, tub] = ind2sub(size(T), ind);
loc(2,1) = xpts(col);
loc(2,2) = ypts(row);
loc(2,3) = zpts(tub);
vals(2) = feval(f, loc(2, 1), loc(2, 2), loc(2, 3));

% Get more digits with optimisation algorithms.
lb = [dom(1); dom(3); dom(5)];
ub = [dom(2); dom(4); dom(6)];

if ( ~isempty(ver('optim')) )
% Matlab's Optimization Toolbox is available. So, use fmincon command.
    options = optimset('Display', 'none', 'TolFun', eps, 'TolX', eps, ...
        'algorithm', 'active-set');
    [minLoc, vals(1)] = fmincon(@(x) feval(f, x(1), x(2), x(3)), ...
        loc(1, :), [], [], [], [], lb, ub, [], options);

    [maxLoc, vals(2)] = fmincon(@(x) -feval(f, x(1), x(2), x(3)), ...
        loc(2,:), [], [], [], [], lb, ub, [], options);
    vals(2) = -vals(2);
    loc(1,:) = minLoc;
    loc(2,:) = maxLoc;
    locs = loc;
else
    try
        % Use core Matlab fminsearch command by converting to an
        % unconstrained optimization problem.
        % Maps from [-1, 1] to dom(1:2), dom(3:4), and dom(5:6)
        % respectively.
        map1 = bndfun.createMap(dom(1:2));
        map2 = bndfun.createMap(dom(3:4));
        map3 = bndfun.createMap(dom(5:6));

        % Unconstrained initial guesses:
        loc(:, 1) = asin(map1.Inv(loc(:, 1)));
        loc(:, 2) = asin(map2.Inv(loc(:, 2)));
        loc(:, 3) = asin(map3.Inv(loc(:, 3)));

        % Maps from R to dom(1:2), dom(3:4), and dom(5:6) respectively.
        map1 = @(x) map1.For(sin(x));
        map2 = @(x) map2.For(sin(x));
        map3 = @(x) map3.For(sin(x));

        % Set options:
        options = optimset('Display', 'off', 'TolFun', eps, 'TolX', eps);
        warnstate = warning;
        warning('off'); % Disable verbose warnings from fminsearch.
        f_mapped = @(x) feval(f, map1(x(1)), map2(x(2)),map3(x(3)));

        [minLoc, vals(1)] = fminsearch(@(x) f_mapped(x), loc(1, :), options);
        [maxLoc, vals(2)] = fminsearch(@(x) -f_mapped(x), loc(2, :), options);

        vals(2) = -vals(2);
        loc(1:2, 1) = map1([minLoc(1); maxLoc(1)]);
        loc(1:2, 2) = map2([minLoc(2); maxLoc(2)]);
        loc(1:2, 3) = map3([minLoc(3); maxLoc(3)]);
        warning(warnstate);
        locs = loc;

    catch
        % Nothing is going to work. So we will have to go with initial
        % guesses.
    end
end

if ( doNewton )
    % Store values before applying any Newton iterations to restore if
    % Newton was diverging.
    locsOld = locs;
    valsOld = vals;

    % If the global max or min is already on the edge or out of domain,
    % do NOT apply Newton iteration:
    if ( isOutOfDomain(locs, dom) )
        return
    end

    % A few steps of Newton for optimization (involves computing Hessian)
    H = Hessian(f);
    gradF = grad(f);
    gradF = @(x,y,z) gradF(x,y,z);

    % Disable verbose warnings from Newton step if e.g., the Hessian
    % matrix is singular which happens e.g., if a 2D function is passed to
    % minandmax3.
    warning_state = warning('off','MATLAB:singularMatrix');
	% Use try-catch to guarantee original warning state is restored.
    try
        lastwarn('')
        k = 1; % k Newton iterations:
        for iter = 1:k
            locs(1,:) = newton_opt(locs(1,:), H, gradF);
            vals(:,1) = feval(f, locs(1, 1), locs(1, 2), locs(1, 3));

            locs(2,:) = newton_opt(locs(2, :), H, gradF);
            vals(:,2) = feval(f,locs(2, 1), locs(2, 2), locs(2, 3));
        end
        [ignored, last_warn] = lastwarn;
        if strcmp(last_warn,'MATLAB:singularMatrix')
            % Return the old value if the Hessian matrix was singular.
            locs = locsOld;
            vals = valsOld;
            return
        end
        warning(warning_state)
    catch err
        warning(warning_state)
        rethrow(err)
    end

    % If Newton's method moved us outside of the domain, then return old
    % value.
    if ( isOutOfDomain(locs, dom) )
        locs = locsOld;
        vals = valsOld;
    end

end

end

%%
function H = Hessian(f)
% Forms the Hessian matrix of a CHEBFUN3 object f

H11 = diff(f, 2, 1);
H12 = diff(diff(f, 1, 2), 1, 1);
H13 = diff(diff(f, 1, 3), 1, 1);

H21 = diff(diff(f, 1, 1), 1, 2);
H22 = diff(f, 2, 2);
H23 = diff(diff(f, 1, 3), 1, 2);

H31 = diff(diff(f, 1, 1), 1, 3);
H32 = diff(diff(f, 1, 2), 1, 3);
H33 = diff(f, 2, 3);
H = @(x,y,z) [feval(H11, x, y, z)  feval(H12, x, y, z)  feval(H13, x, y, z);
              feval(H21, x, y, z)  feval(H22, x, y, z)  feval(H23, x, y, z);
              feval(H31, x, y, z)  feval(H32, x, y, z)  feval(H33, x, y, z)];
end

%%
function sol = newton_opt(init, H, gradF)
% Perform one step of multivariate optimization using Newton's method.
x = init(1);
y = init(2);
z = init(3);
rhs = -gradF(x, y, z);
h = H(x,y,z)\rhs;
sol = init.' + h;
end

%%
function check = isOutOfDomain(locs, dom)
% Check whether both points in the rows of LOCS are inside the cuboid DOM.
if ( locs(1, 1) <= dom(1) || locs(1, 1) >=dom(2) || ...
         locs(1, 2) <= dom(3) || locs(1, 2)>=dom(4) || ...
         locs(1, 3) <= dom(5) || locs(1, 3)>=dom(6) )
     check = true;
elseif ( locs(2, 1) <= dom(1) || locs(2, 1) >= dom(2) || ...
             locs(2, 2) <= dom(3) || locs(2, 2) >= dom(4) || ...
             locs(2, 3) <= dom(5) || locs(2, 3) >= dom(6) )
         check = true;
else
    check = false;
end
end
%%

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function reconstruct(obj, varargin)
% RECONSTRUCT Remove deactivated points from object.
% ------------------------------------------------------------------------------
% OUTPUT
% 1 [obj]
%   Updated object without deactivated points.
% ------------------------------------------------------------------------------
% EXAMPLES
% 1 Thin out of point cloud.
%   pc = pointCloud('Lion.xyz');
%   pc.select('UniformSampling', 2);
%   pc.reconstruct;
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

% Start ------------------------------------------------------------------------

procHierarchy = {'POINTCLOUD' 'RECONSTRUCT'};
msg('S', procHierarchy);
msg('I', procHierarchy, sprintf('Point cloud label = ''%s''', obj.label));

% Create new object with only active points ------------------------------------

% Remove deactivated points from coordinates
obj.X = obj.X(obj.act,:);

% Remove deactivated points from attributes
if isstruct(obj.A)

    % Note: attributes can not be updated field by field due to validation
    % function 'set.A' in pointCloud

    A = obj.A;  % OLD attribute structure
    obj.A = []; % NEW attribute structure

    % Fill NEW attribute structure with values from OLD attribute structure
    att = fields(A);
    for a = 1:numel(att)
        obj.A.(att{a}) = A.(att{a})(obj.act);
        A = rmfield(A, att{a}); % save memory
    end

end

% Update act vector
obj.act = obj.act(obj.act);

% End --------------------------------------------------------------------------

msg('E', procHierarchy);

end

author: pglira


text:
function X = arrange(X,foo)
%ARRANGE Arranges the rank-1 components of a ktensor.
%
%   ARRANGE(X) normalizes the columns of the factor matrices and then sorts
%   the ktensor components by magnitude, greatest to least.
%
%   ARRANGE(X,N) absorbs the weights into the Nth factor matrix instead of
%   lambda.
%
%   ARRANGE(X,P) rearranges the components of X according to the
%   permutation P. P should be a permutation of 1 to NCOMPOMENTS(X).
%
%   See also KTENSOR, NCOMPONENTS, NORMALIZE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Just rearrange and return if second argument is a permutation
if exist('foo','var') && (length(foo) > 1)
    X.lambda = X.lambda(foo);
    for i = 1 : ndims(X)
        X.u{i} = X.u{i}(:,foo);
    end
    return;
end

%% Ensure that matrices are normalized
X = normalize(X);

%% Sort
[X.lambda, idx] = sort(X.lambda, 1, 'descend');
for i = 1 : ndims(X)
    X.u{i} = X.u{i}(:,idx);
end

%% Absorb the weight into one factor, if requested
if exist('foo','var')
    r = length(X.lambda);
    X.u{end} = full(X.u{end} * spdiags(X.lambda,0,r,r));
    X.lambda = ones(size(X.lambda));
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% FRONTIER Rolling efficient frontier.
%    [PORTWGTS,ALLMEAN,ALLCOV] = FRONTIER(UNIVERSE,WINDOW,OFFSET,NPORTS,ACTIVEMAP,CONSET,NUMNONNAN)
%    generates a surface of efficient frontiers given data indicating which
%    stocks are active at each date, portfolio constraints and the number of
%    non missing data points needed for each window showing how asset
%    allocation influences risk and return over time.
%
%    [PW,AM,AC] = FRONTIER(UNIVERSE,WINDOW,OFFSET,NPORTS)
%    generates a surface of efficient frontiers showing how asset allocation
%    influences risk and return over time.
%
%    Inputs:
%
%       UNIVERSE - Portfolio matrix containing the total return data for a
%                  group of securities.  It is an Mx(N+1) matrix where column 1
%                  contains MATLAB date numbers and the remaining columns
%                  are the total return data for each security.  N is the
%                  number of securities.
%       WINDOW - Number of periods of data to use to calculate each frontier.
%       OFFSET - Increment in number of periods between each frontier.
%       NPORTS - The number of portfolios to calculate on each frontier.
%
%    Optional inputs:
%
%
%       ACTIVEMAP - An MxN matrix with boolean elements that correspond to
%                   the UNIVERSE where each element indicates if the asset
%                   is part of the UNIVERSE on the corresponding date.  The
%                   default map is an MxN matrix of ones indicating that all
%                   assets are active at all dates.   Note if constraints,
%                   CONSET, other than the default values are used,
%                   the default ACTIVEMAP must be used.
%       CONSET - Portfolio constraints.  Default constraints are generated by
%                PORTCONS('Default',NumAssets).  This single constraint
%                matrix is applied to each frontier.  Note if an ACTIVEMAP
%                other than the default map, the default constraints must be
%                used.
%       NUMNONNAN - Minimum number of non NaN points for each active asset
%                   in each window of data needed to perform the
%                   optimization.  The default value is WINDOW - N.
%
%    Outputs:
%
%      PORTWGTS - A (Number of Curves)x1 cell array where each element is an
%      NPORTSx(Number of Assets) matrix of weights allocated to each
%      asset.  Number of Assets = length(UNIVERSE).
%
%      ALLMEAN - A (Number of Curves)x1 cell array where each element is an
%      1 x (Number of Assets) vector of the expected asset returns used to
%      generate each curve on the surface.
%
%      ALLCOV - A (Number of Curves)x1 cell array where each element is an
%      (Number of Assets) x (Number of Assets) vector of the covariance
%      matrix used to generate each curve on the surface.
%
%    See also PORTCONS, PORTOPT.
%
%    Reference page in Doc Center
%       doc frontier
%
%    Other functions named frontier
%
%       dsge/frontier
%

author: jmaih


text:
function [K] = ku1u1(x, xp, hyp, i)

logsigmau = hyp(1);
logthetau = hyp(2);
logsigmav = hyp(3);
logthetav = hyp(4);

n_x = size(x,1);
n_xp = size(xp,1);

x = repmat(x,1,n_xp);
xp = repmat(xp',n_x,1);

switch i

case 0

K=exp(1).^(logsigmau+(-1/2).*exp(1).^((-1).*logthetau).*(x+(-1).*xp).^2); ...

case 1 % logsigmau

K=exp(1).^(logsigmau+(-1/2).*exp(1).^((-1).*logthetau).*(x+(-1).*xp).^2); ...

case 2 % logthetau

K=(1/2).*exp(1).^(logsigmau+(-1).*logthetau+(-1/2).*exp(1).^((-1).* ...
  logthetau).*(x+(-1).*xp).^2).*(x+(-1).*xp).^2;

case 3 % logsigmav

K=0;

case 4 % logthetav

K=0;

otherwise

        K = zeros(n_x, n_xp);
end

if K == 0

    K = zeros(n_x, n_xp);

end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function opts = init_LapSRN_opts(scale, depth, gpu)
% -------------------------------------------------------------------------
%   Description:
%       Generate model options for LapSRN
%
%   Input:
%       - scale : SR upsampling scale
%       - depth : number of conv layers in one pyramid level
%       - gpu   : GPU ID, 0 for CPU mode
%
%   Output:
%       - opts  : options for LapSRN
%
%   Citation:
%       Deep Laplacian Pyramid Networks for Fast and Accurate Super-Resolution
%       Wei-Sheng Lai, Jia-Bin Huang, Narendra Ahuja, and Ming-Hsuan Yang
%       IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

    %% network options
    opts.scale              = scale;
    opts.depth              = depth;
    opts.weight_decay       = 0.0001;
    opts.init_sigma         = 0.001;
    opts.conv_f             = 3;
    opts.conv_n             = 64;
    opts.loss               = 'L1';

    %% training options
    opts.gpu                = gpu;
    opts.batch_size         = 64;
    opts.num_train_batch    = 100;     % number of training batch in one epoch
    opts.num_valid_batch    = 100;      % number of validation batch in one epoch
    opts.lr                 = 1e-5;     % initial learning rate
    opts.lr_step            = 50;       % number of epochs to drop learning rate
    opts.lr_drop            = 0.5;      % learning rate drop ratio
    opts.lr_min             = 1e-6;     % minimum learning rate
    opts.patch_size         = 128;
    opts.data_augmentation  = 1;
    opts.scale_augmentation = 1;

    %% dataset options
    opts.train_dataset          = {};
    opts.train_dataset{end+1}   = 'T91';
    opts.train_dataset{end+1}   = 'BSDS200';
    %opts.train_dataset{end+1}   = 'General100';
    opts.valid_dataset          = {};
    %opts.valid_dataset{end+1}   = 'Set5';
    %opts.valid_dataset{end+1}   = 'Set14';
    opts.valid_dataset{end+1}   = 'BSDS100';

    %% setup model name
    opts.data_name = 'train';
    for i = 1:length(opts.train_dataset)
        opts.data_name = sprintf('%s_%s', opts.data_name, opts.train_dataset{i});
    end

    opts.net_name = sprintf('LapSRN_x%d_depth%d_%s', ...
                            opts.scale, opts.depth, opts.loss);

    opts.model_name = sprintf('%s_%s_pw%d_lr%s_step%d_drop%s_min%s_bs%d', ...
                            opts.net_name, ...
                            opts.data_name, opts.patch_size, ...
                            num2str(opts.lr), opts.lr_step, ...
                            num2str(opts.lr_drop), num2str(opts.lr_min), ...
                            opts.batch_size);

    %% setup dagnn training parameters
    if( opts.gpu == 0 )
        opts.train.gpus     = [];
    else
        opts.train.gpus     = [opts.gpu];
    end
    opts.train.batchSize    = opts.batch_size;
    opts.train.numEpochs    = 1000;
    opts.train.continue     = true;
    opts.train.learningRate = learning_rate_policy(opts.lr, opts.lr_step, opts.lr_drop, ...
                                                   opts.lr_min, opts.train.numEpochs);

    opts.train.expDir = fullfile('models', opts.model_name) ; % model output dir
    if( ~exist(opts.train.expDir, 'dir') )
        mkdir(opts.train.expDir);
    end

    opts.train.model_name       = opts.model_name;
    opts.train.num_train_batch  = opts.num_train_batch;
    opts.train.num_valid_batch  = opts.num_valid_batch;

    % setup loss
    opts.level = ceil(log(opts.scale) / log(2));
    opts.train.derOutputs = {};
    for s = opts.level : -1 : 1
        opts.train.derOutputs{end+1} = sprintf('level%d_%s_loss', s, opts.loss);
        opts.train.derOutputs{end+1} = 1;
    end

end

repo: LapSRN


text:
function c = horzcat(varargin)
% function C=horzcat(A,B);
%
% DESCRIPTION
%   Horizontal concatenation of polynomial objects.
%
% INPUTS
%   A,B: polynomials
%
% OUTPUTS
%   C:  horizontal concatenation of input matrices.
%
% SYNTAX
%   [A B]
%       Horizontal concatenation of polynomial matrices A and B.
%       A and B must have the same number of rows.
%   [A1 A2 A3 ...]
%       Horizontal concatenation of several polynomial matrices.
%   C = horzcat(A1,A2,...);
%       Function-call form of horizontal concatenation.
%
% See also vertcat

% 6/8/2002: PJS  Initial Coding

if nargin==1
    c = varargin{1};
else
    % Promote a to polynomial
    a = polynomial(varargin{1});
    [nra,nca] = size(a);

    % Promote b to polynomial
    b = polynomial(varargin{2});
    [nrb,ncb] = size(b);

    if isempty(b);
        c = a;
    elseif isempty(a);
        c = b;
    elseif nra==nrb
        % Get Dimensions
        nta = size(a.degmat,1);
        nva = length(a.varname);
        ntb = size(b.degmat,1);
        nvb = length(b.varname);

        if nva==0 && nvb==0
            % Combine constant terms
            ar = combine(a);
            coef1 = reshape(ar.coefficient,[nra nca]);
            br = combine(b);
            coef2 = reshape(br.coefficient,[nrb ncb]);

            % Stack Coefficients and Form Polynomial
            coefficient = [coef1 coef2];
            c = polynomial(coefficient);
        else
            % Form Degmat, Varname, and Matdim
            adeg = a.degmat;
            bdeg = b.degmat;
            degmat = blkdiag(adeg,bdeg);
            varname = [a.varname(:); b.varname(:)];
            matdim = [nra nca+ncb];

            % Stack up Coefficients
            coef1 = [a.coefficient zeros(nta,nra*ncb)];
            coef2 = [zeros(ntb,nrb*nca) b.coefficient];
            coefficient = [coef1; coef2];

            % Form Polynomial and combine terms
            chkval = 0; % skip validity check
            c = polynomial(coefficient,degmat,varname,matdim,chkval);
            c = combine(c);
        end
    else
        error('All columns must have the same row dimension')
    end

    if nargin>2
        c = horzcat(c,varargin{3:end});
    end
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function s=v_sprintcpx(z,f)
%V_SPRINTCPX  format a complex number for printing S=(Z,F)
%
% Usage: fprintf('%s',v_sprintcpx(z));
%
%  Inputs: z   a complex number to print
%          f   optional formatting string as in fprintf e.g. '0.2f' [default: 'g']
%              may also include 'i' or 'j' [default] to control sqrt(-1) symbol.
%
% Outputs: s   formatted output string

%      Copyright (C) Mike Brookes 2015
%      Version: $Id: v_sprintcpx.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2 || ~numel(f)
    f='g';
end
if any(f=='i')
    ij='i';
else
    ij='j';
end
f((f=='i')|(f=='j'))=[]; % remove i and j specifiers
if ~numel(f)
    f='g';
end
if any(f=='+')
    pl='';
else
    pl='+';
end
f=['%' f];
a=real(z);
b=imag(z);
jx=[1 3 2 4 3 4 1 3 2];
ix=jx(3*sign(a)+sign(b)+5);
switch(ix)
    case 1
        s=sprintf([f f ij],a,b);
    case 2
        s=sprintf([f pl f ij],a,b);
    case 3
        s=sprintf(f,a);
    case 4
        s=sprintf([f ij],b);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function kern = rbfperiodicKernExpandParam(kern, params)

% RBFPERIODICKERNEXPANDPARAM Create kernel structure from RBFPERIODIC kernel's parameters.
% FORMAT
% DESC returns a RBF derived periodic kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : rbfperiodicKernParamInit, rbfperiodicKernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

kern.inverseWidth = params(1);
kern.variance = params(2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function rcnn_exp_cache_features(chunk)

% -------------------- CONFIG --------------------
net_file     = './data/caffe_nets/finetune_voc_2007_trainval_iter_70k';
cache_name   = 'v1_finetune_voc_2007_trainval_iter_70k';
crop_mode    = 'warp';
crop_padding = 16;

% change to point to your VOCdevkit install
VOCdevkit = './datasets/VOCdevkit2007';
% ------------------------------------------------

imdb_train = imdb_from_voc(VOCdevkit, 'train', '2007');
imdb_val   = imdb_from_voc(VOCdevkit, 'val', '2007');
imdb_test  = imdb_from_voc(VOCdevkit, 'test', '2007');
imdb_trainval = imdb_from_voc(VOCdevkit, 'trainval', '2007');

switch chunk
  case 'train'
    rcnn_cache_pool5_features(imdb_train, ...
        'crop_mode', crop_mode, ...
        'crop_padding', crop_padding, ...
        'net_file', net_file, ...
        'cache_name', cache_name);
    link_up_trainval(cache_name, imdb_train, imdb_trainval);
  case 'val'
    rcnn_cache_pool5_features(imdb_val, ...
        'crop_mode', crop_mode, ...
        'crop_padding', crop_padding, ...
        'net_file', net_file, ...
        'cache_name', cache_name);
    link_up_trainval(cache_name, imdb_val, imdb_trainval);
  case 'test_1'
    end_at = ceil(length(imdb_test.image_ids)/2);
    rcnn_cache_pool5_features(imdb_test, ...
        'start', 1, 'end', end_at, ...
        'crop_mode', crop_mode, ...
        'crop_padding', crop_padding, ...
        'net_file', net_file, ...
        'cache_name', cache_name);
  case 'test_2'
    start_at = ceil(length(imdb_test.image_ids)/2)+1;
    rcnn_cache_pool5_features(imdb_test, ...
        'start', start_at, ...
        'crop_mode', crop_mode, ...
        'crop_padding', crop_padding, ...
        'net_file', net_file, ...
        'cache_name', cache_name);
end

% ------------------------------------------------------------------------
function link_up_trainval(cache_name, imdb_split, imdb_trainval)
% ------------------------------------------------------------------------
cmd = {['mkdir -p ./feat_cache/' cache_name '/' imdb_trainval.name '; '], ...
    ['cd ./feat_cache/' cache_name '/' imdb_trainval.name '/; '], ...
    ['for i in `ls -1 ../' imdb_split.name '`; '], ...
    ['do ln -s ../' imdb_split.name '/$i $i; '], ...
    ['done;']};
cmd = [cmd{:}];
fprintf('running:\n%s\n', cmd);
system(cmd);
fprintf('done\n');

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
function z=scb(x,y,varargin)

% Simultaneous Confidence Bands
%
% Example:
%   load ethanol;
%   z = scb(E,NOx,'h',0.5);
%
% result (z) is a matrix with four columns: evaluation points,
% fitted values, lower confidence limit, upper confidence limit.
% Most locfit arguments should work.

fit = locfit(x,y,'ev','grid','mg',20,varargin{:});
kap = kappa0(x,y,varargin{:});
cb = predict(fit,'fitp','band','g','kappa',kap);
z = [fit.fit_points.evaluation_points' cb{1} cb{3}];

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [v,f] = mesh_flip_faces(arg1, arg2, arg3)

if nargin==1
    name = arg1;
    [v,f] = read_mesh(name);
else
    v = arg1;
    f = arg2
end

if size(f,1)~=3
    f = f';
end

f = [f(1,:); f(3,:); f(2,:)];

if nargin==1
    write_mesh(name, v,f);
end

author: gpeyre

repo: matlab-toolboxes


text:
function blockface_correct_frame_shifts(indir, files, t, idxnoprop, outdir)
% blockface_correct_frame_shifts  Correct shift between blockface image
% frames.
%
% When the wax block is sliced with a microtome, and photographs taken of
% the blockface, artifacts invariable happen. Artifacts are typically
% caused by the microtome tray not being pushed all the way back
% (translation artifacts), or by the camera, mirrors or any other elements
% in the light path being moved a little (similarity transform artifacts).
%
% This function assumes that the corrective transforms between consecutive
% frames have been already computed with elastix. The user provides a list
% of frames that need to be corrected, and blockface_correct_frame_shifts
% computes the accumulated transform for each frame, and applies it with
% transformix. Corrected images are saved with the same name as input
% images, but to a different directory.
%
% The transformation of images uses a "parfor" loop, so the function can
% take advantage of multiple cores.
%
% blockface_correct_frame_shifts(INDIR, FILES, T, IDXNOPROP, OUTDIR)
%
%   INDIR is a string with the directory where the input files are kept.
%
%   FILES is the result of a dir() command with a list of the files to
%   correct.
%
%   T is a struct array where T(I) is the correction that registers image I
%   to image I-1. T can be obtained with blockface_intraframe_reg.
%
%   IDXNOPROP is a vector with a list of indices of frames that need to be
%   corrected, but whose correction does not propagate. For example, if the
%   microtome tray is not pushed all the way in in frame I=30, that doesn't
%   mean that that problem will appear in frame I = 31. This is useful if
%   e.g. the camera has been very stable, but now and then the human
%   operator doesn't push the tray all the way back, as it avoids
%   unncessary small transforms and resampling of images.
%
%   OUTDIR is a string with the directory for output images. INDIR and
%   OUTDIR must be different to avoid overwriting the input files.
%
% See also: blockface_intraframe_reg.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2014, 2016 University of Oxford
% Version: 0.2.2
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(5, 5);
nargoutchk(0, 0);

% safety to prevent us from overwriting the source files
if (strcmp(indir, outdir))
    error('Input and output directory are the same, and files would be overwritten')
end

% if destination directory doesn't exist, create it
if (isempty(dir(outdir)))
    ok = mkdir(outdir);
    if (~ok)
        error(['Cannot create output directory: ' outdir])
    end
end

% first frame will be used as a reference for the rest
% copytogreyfile([pathstr filesep files(1).name], ...
%     [outdir filesep files(1).name]);

% initialize array for the accumulated transform for each frame
tAcc = t(1);

% clear up transform for first frame
tAcc(1).TransformParameters = [1 0 0 0];

% for transforms that do not propagate, we need to cancel them with an
% inverse transform in the next slice
for I = idxnoprop

    % if the slice is the last in the stack, there's no next slice to
    % cancel the transform with
    if (I == length(files))
        continue
    end

    % the next slice must have no transformation, otherwise the list of
    % indices is wrong
    if any(t(I+1).TransformParameters ~= [1 0 0 0])
        error(['User says that T(' num2str(I) ') does not propagate, but T(' ...
            num2str(I+1) ') ~= [1 0 0 0]'])
    end

    % cancel the transformation
    tMatrix = elastix_affine_struct2matrix(t(I));
    t(I+1) = elastix_affine_matrix2struct(inv(tMatrix), t(I+1));

end

% accumulate transforms
for I = 2:length(t)

    % update the accumulated transform with the transform for this
    % frame
    tAcc(I) = elastix_compose_afftransf(tAcc(I-1), t(I));

end

% correct frames
parfor I = 1:length(files)

    % is the transform the identity?
    if (all(tAcc(I).TransformParameters == [1 0 0 0]))

        % simply convert to grayscale and copy
        copytogreyfile([indir filesep files(I).name], ...
            [outdir filesep files(I).name]);

    else

        % affine transformation of image to correct frame shifts
        fileout = transformix(tAcc(I), [indir filesep files(I).name]);

        % copy result to file
        copytogreyfile(fileout, ...
            [outdir filesep files(I).name]);

        % delete the temp result file
        delete(fileout);

    end

end

end

function copytogreyfile(filein, fileout)

% check that the image is grayscale. Otherwise, we convert it,
% because elastix uses the first channel and ignores the rest. As
% all transformed images will be converted to grayscale, we also
% convert images that are not shifted, so that they are all
% grayscale
info = imfinfo(filein);

switch (info.ColorType)
    case 'truecolor'

        % load image
        im = imread(filein);

        % convert to grayscale
        im = rgb2gray(im);

        % write image to destination
        imwrite(im, fileout);

    case 'grayscale'

        % if the image is already grayscale, we can directly copy
        % it
        ok = copyfile(filein, fileout);
        if (~ok)
            error(['Cannot copy image ' filein ' to file ' fileout])
        end

    otherwise

        error('Image in input file is neither RGB nor grayscale')
end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
options = optimset('GradObj','on','Hessian','on');
[x,y]=fminunc('fun4',rand(1,2),options)

author: Eurus-Holmes


text:
function z = condEntropy (x, y)
% Compute conditional entropy z=H(x|y) of two discrete variables x and y.
% Input:
%   x, y: two integer vector of the same length
% Output:
%   z: conditional entropy z=H(x|y)
% Written by Mo Chen (sth4nth@gmail.com).
assert(numel(x) == numel(y));
n = numel(x);
x = reshape(x,1,n);
y = reshape(y,1,n);

l = min(min(x),min(y));
x = x-l+1;
y = y-l+1;
k = max(max(x),max(y));

idx = 1:n;
Mx = sparse(idx,x,1,n,k,n);
My = sparse(idx,y,1,n,k,n);
Pxy = nonzeros(Mx'*My/n); %joint distribution of x and y
Hxy = -dot(Pxy,log2(Pxy));

Py = nonzeros(mean(My,1));
Hy = -dot(Py,log2(Py));

% conditional entropy H(x|y)
z = Hxy-Hy;
z = max(0,z);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
%SETSCALE Set scale field (output scaling) in mapping
%
%    W = SETSCALE(W,SCALE)
%
% This sets the SCALE field in W. If W is affine, scaling is is directly
% performed in the weights in W.DATA and W.SCALE = 1.

author: marianux

repo: ecg-kit

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


text:
  function ob = Godwt1(mask, varargin)
%|function ob = Godwt1(mask, varargin)
%| Construct Godwt1 object that computes orthonormal discrete wavelet
%| decomposition of a signal with dimensions [(N)].
%| This is useful for sparsity regularization (aka compressed sensing).
%| (1D or 2D wavelets only)
%|
%| in
%|	'mask'	logical [(Nd)]	image-domain mask, often true(nx,ny)
%|
%| options
%|	'level'	int	decomposition level (default 1)
%|	'wname'	char	wavelet name. default: 'haar'
%|
%| out
%|	ob	[*N np]	fatrix  object, where np = sum(mask(:))
%|
%|
%| Copyright 2012-05-17, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(mask, 'test'), Godwt1_test, return, end

arg.mask = mask;
arg.level = 1;
arg.wname = 'haar';
arg.abs = false;
arg = vararg_pair(arg, varargin);

if isempty(arg.mask), fail 'must provide a mask', end

% transform dimension
idim = size(arg.mask);
if idim(end) == 1
	idim = idim(1:end-1);
end

switch numel(idim)
case 1
	forw = @(arg, x) ir_odwt1(x, ...
		'level', arg.level, 'wname', arg.wname, 'abs', arg.abs);
	back = @(arg, y) fatrix2_maskit(arg.mask, ir_odwt1(y, 'adj', 1, ...
		'level', arg.level, 'wname', arg.wname, 'abs', arg.abs));
	does_many = true;
case 2
	forw = @(arg, x) ir_odwt2(x, ...
		'level', arg.level, 'wname', arg.wname, 'abs', arg.abs);
	back = @(arg, y) fatrix2_maskit(arg.mask, ir_odwt2(y, 'adj', 1, ...
		'level', arg.level, 'wname', arg.wname, 'abs', arg.abs));
	does_many = false;
otherwise
        fail('dim %d unsupported', numel(idim))
end

% build fatrix2 object
arg.idim = idim;
ob = fatrix2('idim', idim, 'odim', idim, 'arg', arg, ...
	'abs', @Godwt1_abs, 'meth', {'codes', @Godwt1_codes, '()'}, ...
	'forw', forw, 'back', back, 'does_many', does_many);

% Godwt1_abs()
function ob = Godwt1_abs(ob)
ob.arg.abs = true;

% Godwt1_codes()
function codes = Godwt1_codes(arg)

switch numel(arg.idim)
case 1
	[dummy codes] = ir_odwt1(zeros([arg.idim 1]), ...
		'level', arg.level, 'wname', arg.wname);
case 2
	[dummy codes] = ir_odwt2(zeros(arg.idim), ...
		'level', arg.level, 'wname', arg.wname);
otherwise
        fail('dim %d unsupported', numel(idim))
end

% Godwt1_test()
function Godwt1_test

if 1 % 1d
	mask = true(8*3,1);
	mask(1:3) = false;
	level = 3;
	U = Godwt1(mask, 'level', level);
%	abs(U) % todo

	if im
		im plc 1 2
		im(1, full(U))
		im(2, U.codes)
		drawnow
	end

	fatrix2_tests(U, 'complex', 1) % check complex data
	test_adjoint(U, 'complex', 1, 'tolre', 1e-9);
end

if 1 % 2d
	mask = true(8*3,16);
	mask(1:3) = false;
	level = 3;
	wname = 'haar';
	wname = 'sym2';
	U = Godwt1(mask, 'level', level, 'wname', wname);

	x = ellipse_im(size(mask));
	if im
		im plc 1 2
		im(1, x)
		im(2, U * x)
		drawnow
	end

	fatrix2_tests(U, 'complex', 1) % check complex data
	test_adjoint(U, 'complex', 1, 'tolre', 1e-9, 'big', 1);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function D = driving_function_mono_nfchoa_fs(x0,xs,f,N,conf)
%DRIVING_FUNCTION_MONO_NFCHOA_FS driving signal for a focused source in NFC-HOA
%
%   Usage: D = driving_function_mono_nfchoa_fs(x0,xs,f,N,conf)
%
%   Input parameters:
%       x0          - position of the secondary sources / m [nx3]
%       xs          - position of focused source / m [nx3]
%       f           - frequency of the monochromatic source / Hz
%       N           - maximum order of spherical harmonics
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function signal [nx1]
%
%   See also: driving_function_mono_nfchoa, driving_function_imp_nfchoa_fs

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 5;
narginchk(nargmin,nargmax);
isargmatrix(x0,xs);
isargpositivescalar(f,N);
isargstruct(conf);

%% ===== Configuration ==================================================
xref = conf.xref;
dimension = conf.dimension;
driving_functions = conf.driving_functions;

%% ===== Computation ====================================================

if strcmp('2D',dimension)

    % === 2-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2D focused source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('2.5D',dimension)

    % === 2.5-Dimensional ================================================

    % Reference point
    xref = repmat(xref,[size(x0,1) 1]);

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2.5D focused source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('3D',dimension)

    % === 3-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 3D focused source.'],upper(mfilename),driving_functions);
    end

else
    error('%s: the dimension %s is unknown.',upper(mfilename),dimension);
end

author: sfstoolbox

repo: sfs-matlab


text:
function uvTformData = sc_src_domain_tform(uvPlaneID, modelPlane, modelReg, srcPos, trgPos, sampleRandReg)

% SC_SRC_DOMAIN_TFORM: The function estimates the source patch domain
% transformation based on the source and target position, and the plane parameters
%
% Input:
%     - uvPlaneID:  numUvPix x 1, the plane label of uv pixels
%     - modelPlane: plane model
%     - modelReg:   regularity model
%     - srcPos:     numUvPix x 2, the center position of source patch
%     - trgPos:     numUvPix x 2, the center position of target patch
%     - sampleRandReg:
%         1: random sampling mode
%         0: regular sampling mode
% Output:
%     - uvTformData: numUvPix x 9

numUvPix    = size(srcPos, 1);

uvTformData = zeros(numUvPix, 9, 'single');
I = eye(3);

for indPlane = 1: modelPlane.numPlane
    % The rectifying transformation for the plane
    rectMat = modelPlane.rectMat{indPlane};
    h7 = rectMat(3,1);
    h8 = rectMat(3,2);

    % Retrieve the uv pixels that have the current plane label
    uvPlaneIndCur = uvPlaneID == indPlane;
    numPlanePixCur = sum(uvPlaneIndCur);

    if(numPlanePixCur)
        % Target patch center position in the rectified domain
        trgPosCur = trgPos(uvPlaneIndCur, :) - 1;
        trgPosCurR = sc_apply_tform_H(trgPosCur, h7, h8);

        % Get dRect
        if(sampleRandReg)    % Random sampling
            % Source patch center position in the rectified domain
            srcPosCur = srcPos(uvPlaneIndCur, :) - 1;
            srcPosCurR = sc_apply_tform_H(srcPosCur, h7, h8);

            % Displacement vector from target to source position
            dRect = srcPosCurR - trgPosCurR;
        else                 % Regularity guided sampling
            dRect = zeros(numPlanePixCur, 2, 'single');

            numDispVecCur = modelReg.numDispVec(indPlane);
            if(numDispVecCur~=0)
                randInd = randi(numDispVecCur, numPlanePixCur, 1);
                dRect = modelReg.dispVec{indPlane}(randInd, :);
            end
        end

        % Compute the transformation that maps from target to source
        % (See Eqn 8 in the paper)
        if(size(dRect,2) ~= 0)
            uvTformCur = zeros(numPlanePixCur, 9, 'single');
            uvTformCur(:,[1,4,7]) = bsxfun(@times, dRect(:,1), [h7, h8, 1]);
            uvTformCur(:,[2,5,8]) = bsxfun(@times, dRect(:,2), [h7, h8, 1]);

            dTemp = dRect*[h7;h8]; % dTemp = dx*h7 + dy*h8
            uvTformCur(:,[3,6,9]) = bsxfun(@times, dTemp, -[h7, h8, 1]);
            uvTformCur = bsxfun(@plus, uvTformCur, I(:)');

            % Apply the offset to cancel out the dependency of the target position
            % (See Eqn 9 in the paper)
            uvTformData(uvPlaneIndCur, :)   = sc_trans_tform(uvTformCur, trgPosCur);
            uvTformData(uvPlaneIndCur, 7:8) = uvTformData(uvPlaneIndCur, 7:8) + 1;
        end
    end
end

end

function y = sc_apply_tform_H(x, h7, h8)
% Apply homography H with third row [h7, h8, 1] to 2D points x

y = x(:,1)*h7 + x(:,2)*h8 + 1;
y = bsxfun(@rdivide, x(:,1:2), y + eps);

end

repo: StructCompletion