text:
classdef VariationalRefinement < handle
    %VARIATIONALREFINEMENT  Variational optical flow refinement
    %
    % This class implements variational refinement of the input flow field,
    % i.e. it uses input flow to initialize the minimization of the following
    % functional:
    %
    %     E(U) = integral_{Omega}(delta * Psi(E_I) + gamma * Psi(E_G) + alpha * Psi(E_S))
    %
    % where `E_I`, `E_G`, `E_S` are color constancy, gradient constancy and
    % smoothness terms respectively. `Psi(s^2) = sqrt(s^2 + epsilon^2)` is a
    % robust penalizer to limit the influence of outliers. A complete
    % formulation and a description of the minimization procedure can be found
    % in [Brox2004].
    %
    % ## References
    % [Brox2004]:
    % > Thomas Brox, Andres Bruhn, Nils Papenberg, and Joachim Weickert.
    % > "High accuracy optical flow estimation based on a theory for warping".
    % > In Computer Vision-ECCV 2004, pages 25-36. Springer, 2004.
    %
    % See also: cv.VariationalRefinement.calc
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % Number of outer (fixed-point) iterations in the minimization
        % procedure. default 5
        FixedPointIterations
        % Number of inner successive over-relaxation (SOR) iterations in the
        % minimization procedure to solve the respective linear system.
        % default 5
        SorIterations
        % Relaxation factor in SOR. default 1.6
        Omega
        % Weight of the smoothness term. default 20.0
        Alpha
        % Weight of the color constancy term. default 5.0
        Delta
        % Weight of the gradient constancy term. default 10.0
        Gamma
    end

    %% Constructor/destructor
    methods
        function this = VariationalRefinement()
            %VARIATIONALREFINEMENT  Creates an instance of VariationalRefinement
            %
            %     obj = cv.VariationalRefinement()
            %
            % See also: cv.VariationalRefinement.calc
            %
            this.id = VariationalRefinement_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.VariationalRefinement
            %
            if isempty(this.id), return; end
            VariationalRefinement_(this.id, 'delete');
        end
    end

    %% DenseOpticalFlow
    methods
        function flow = calc(this, I0, I1, varargin)
            %CALC  Calculates an optical flow
            %
            %     flow = obj.calc(I0, I1)
            %     flow = obj.calc(I0, I1, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __I0__ first 8-bit single-channel input image.
            % * __I1__ second input image of the same size and the same type
            %   as `I0`.
            %
            % ## Output
            % * __flow__ computed flow image that has the same size as `I0`
            %   and type `single` (2-channels).
            %
            % ## Options
            % * __InitialFlow__ specify the initial flow. Not set by default.
            %
            % See also: cv.VariationalRefinement
            %
            flow = VariationalRefinement_(this.id, 'calc', I0, I1, varargin{:});
        end

        function collectGarbage(this)
            %COLLECTGARBAGE  Releases all inner buffers
            %
            %     obj.collectGarbage()
            %
            VariationalRefinement_(this.id, 'collectGarbage');
        end
    end

    %% VariationalRefinement
    methods
        function [flow_u, flow_v] = calcUV(this, I0, I1, varargin)
            %CALCUV  calc function overload to handle separate horizontal (u) and vertical (v) flow components (to avoid extra splits/merges)
            %
            %     [flow_u, flow_v] = obj.calcUV(I0, I1)
            %     [flow_u, flow_v] = obj.calcUV(I0, I1, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __I0__ first 8-bit single-channel input image.
            % * __I1__ second input image of the same size and the same type
            %   as `I0`.
            %
            % ## Output
            % * **flow_u** computed horizontal flow image that has the same
            %   size as `I0` and type `single` (1-channel1).
            % * **flow_v** computed vertical flow image that has the same size
            %   as `I0` and type `single` (1-channel1).
            %
            % ## Options
            % * __InitialFlowU__ specify initial U-flow. Not set by default.
            % * __InitialFlowV__ specify initial V-flow. Not set by default.
            %
            % Note that `flow(:,:,1)==flow_u`, and `flow(:,:,2)==flow_v`.
            %
            % See also: cv.VariationalRefinement.calc
            %
            [flow_u, flow_v] = VariationalRefinement_(this.id, 'calcUV', I0, I1, varargin{:});
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.VariationalRefinement.empty
            %
            VariationalRefinement_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the algorithm is empty (e.g. in the very
            %   beginning or after unsuccessful read).
            %
            % See also: cv.VariationalRefinement.clear
            %
            b = VariationalRefinement_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.VariationalRefinement.save, cv.VariationalRefinement.load
            %
            name = VariationalRefinement_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in a file storage.
            %
            % See also: cv.VariationalRefinement.load
            %
            VariationalRefinement_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from a file storage.
            % The previous model state is discarded.
            %
            % See also: cv.VariationalRefinement.save
            %
            VariationalRefinement_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.FixedPointIterations(this)
            value = VariationalRefinement_(this.id, 'get', 'FixedPointIterations');
        end
        function set.FixedPointIterations(this, value)
            VariationalRefinement_(this.id, 'set', 'FixedPointIterations', value);
        end

        function value = get.SorIterations(this)
            value = VariationalRefinement_(this.id, 'get', 'SorIterations');
        end
        function set.SorIterations(this, value)
            VariationalRefinement_(this.id, 'set', 'SorIterations', value);
        end

        function value = get.Omega(this)
            value = VariationalRefinement_(this.id, 'get', 'Omega');
        end
        function set.Omega(this, value)
            VariationalRefinement_(this.id, 'set', 'Omega', value);
        end

        function value = get.Alpha(this)
            value = VariationalRefinement_(this.id, 'get', 'Alpha');
        end
        function set.Alpha(this, value)
            VariationalRefinement_(this.id, 'set', 'Alpha', value);
        end

        function value = get.Delta(this)
            value = VariationalRefinement_(this.id, 'get', 'Delta');
        end
        function set.Delta(this, value)
            VariationalRefinement_(this.id, 'set', 'Delta', value);
        end

        function value = get.Gamma(this)
            value = VariationalRefinement_(this.id, 'get', 'Gamma');
        end
        function set.Gamma(this, value)
            VariationalRefinement_(this.id, 'set', 'Gamma', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function[phi] = orth_proj(P)
%%% Funciton to compute the projection of a vector orthogonal to the
%%% subspace spanned by a matrix

%%%         Inputs
%%%     P - basis matrix
%%%     x - input vector

%%%         Outputs
%%%     phi- orthogonal projection matrix
%%%     y - output vector

[nr, ~] = size(P);

phi = eye(nr) - P * P';
%y = x - P * (P' * x);

end

author: andrewssobral

repo: lrslibrary


text:
function varargout=myidentification(varargin)

[varargout{1:nargout}]=dsge_tools.identification(varargin{:});

end

author: jmaih


text:
function precisions = precision_plot(positions, ground_truth, video_name, show)
%PRECISION_PLOT
%   Calculates precision for a series of distance thresholds (percentage of
%   frames where the distance to the ground truth is within the threshold).
%   The results are shown in a new figure if SHOW is true.
%
%   Accepts positions and ground truth as Nx2 matrices (for N frames), and
%   a title string.
%
%   Joao F. Henriques, 2014
%   http://www.isr.uc.pt/~henriques/

	max_threshold = 50;  %used for graphs in the paper

	precisions = zeros(max_threshold, 1);

	if size(positions,1) ~= size(ground_truth,1)
% 		fprintf('%12s - Number of ground truth frames does not match number of tracked frames.\n', title)

		%just ignore any extra frames, in either results or ground truth
		n = min(size(positions,1), size(ground_truth,1));
		positions(n+1:end,:) = [];
		ground_truth(n+1:end,:) = [];
	end

	%calculate distances to ground truth over all frames
    gt_center = [ground_truth(:,1)+(ground_truth(:,3)-1)/2 ground_truth(:,2)+(ground_truth(:,4)-1)/2];
    positions_center = [positions(:,1)+(positions(:,3)-1)/2 positions(:,2)+(positions(:,4)-1)/2];

    distances = sqrt(sum((positions_center-gt_center).^2,2));

    index = ground_truth>0;
    ind = (sum(index,2)==4);

    distances(~ind) = -1;

    %compute precisions
    for p = 1:max_threshold
        precisions(p) = nnz(distances <= p) / numel(distances);
    end

    %get annotation of plot
    score = precisions(20);
    tmp = sprintf('%.3f', score);
    tmpName = ['TrackerName' ' [' tmp ']'];
    titleName = 'Precision plots of error';

	%plot the precisions
	if show == 1
		figure('NumberTitle','off', 'Name',['Precisions - ' video_name])
		plot(precisions, 'r-', 'LineWidth',2)
        legend1 = legend(tmpName);
		xlabel('Threshold'), ylabel('Precision')
        title(titleName)
        saveas(gca,'PrecisionPlot','jpg')
	end

end

author: Daikenan

repo: ASRCF

save_path: github-repos/MATLAB/Daikenan-ASRCF


text:
% demos for ch07

%% sparse signal recovery demo
clear; close all;

d = 512; % signal length
k = 20;  % number of spikes
n = 100; % number of measurements
%
% random +/- 1 signal
x = zeros(d,1);
q = randperm(d);
x(q(1:k)) = sign(randn(k,1));

% projection matrix
A = unitize(randn(d,n),1);
% noisy observations
sigma = 0.005;
e = sigma*randn(1,n);
y = x'*A + e;

[model,llh] = rvmRegSeq(A,y);
plot(llh);

m = zeros(d,1);
m(model.index) = model.w;

h = max(abs(x))+0.2;
x_range = [1,d];
y_range = [-h,+h];
figure;
subplot(2,1,1);plot(x); axis([x_range,y_range]); title('Original Signal');
subplot(2,1,2);plot(m); axis([x_range,y_range]); title('Recovery Signal');

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function [kls,xx1,xx2] = vbmc_kldiv(vp1,vp2,Ns,gaussflag)
%VBMC_KLDIV Kullback-Leibler divergence between two variational posteriors.
%   KLS = VBMC_KLDIV(VP1,VP2) returns an estimate of the (asymmetric)
%   Kullback-Leibler (KL) divergence between two variational posterior
%   distributions VP1 and VP2. KLS is a 2-element vector whose first element
%   is KL(VP1||VP2) and the second element is KL(VP2||VP1). The symmetrized
%   KL divergence can be computed as mean(KLS).
%
%   KLS = VBMC_KLDIV(VP1,VP2,NS) uses NS random draws to estimate each
%   KL divergence (default NS=1e5).
%
%   KLS = VBMC_KLDIV(VP1,VP2,NS,GAUSSFLAG) computes the "Gaussianized"
%   KL-divergence if GAUSSFLAG=1, that is the KL divergence between two
%   multivariate normal distibutions with the same moments as the variational
%   posteriors given as inputs. Otherwise, the standard KL-divergence is
%   returned for GAUSSFLAG=0 (default).
%
%   [KLS,XX1,XX2] = VBMC_KLDIV(...) returns NS samples from the variational
%   posteriors VP1 and VP2 as, respectively, NS-by-D matrices XX1 and XX2,
%   where D is the dimensionality of the problem.
%
%   If GAUSSFLAG is 1, VP1 and/or VP2 can be N-by-D matrices of samples
%   from variational posteriors (they do not need have the same number
%   of samples).
%
%   See also VBMC, VBMC_MTV, VBMC_PDF, VBMC_RND, VBMC_DIAGNOSTICS.

if nargin < 3 || isempty(Ns); Ns = 1e5; end
if nargin < 4 || isempty(gaussflag); gaussflag = false; end

% This was removed because the comparison *has* to be in original space,
% given that the transform might change for distinct variational posteriors
% if nargin < 5 || isempty(origflag); origflag = true; end
origflag = true;

kls = NaN(1,2);

if ~gaussflag && (~vbmc_isavp(vp1) || ~vbmc_isavp(vp2))
    error('vbmc_kldiv:WrongInputs', ...
        'Unless the KL divergence is Gaussianized, VP1 and VP2 need to be variational posteriors.');
end

%try
    if gaussflag
        if Ns == 0  % Analytical calculation
            if origflag
                error('vbmc_kldiv:NoAnalyticalMoments', ...
                    'Analytical moments are available only for the transformed space.')
            end
            [q1mu,q1sigma] = vbmc_moments(vp1,0);
            [q2mu,q2sigma] = vbmc_moments(vp2,0);
            xx1 = []; xx2 = [];
        else        % Numerical moments
            if vbmc_isavp(vp1)
                [q1mu,q1sigma] = vbmc_moments(vp1,origflag,Ns);
            else
                q1mu = mean(vp1,1);
                q1sigma = cov(vp1);
            end
            if vbmc_isavp(vp2)
                [q2mu,q2sigma] = vbmc_moments(vp2,origflag,Ns);
            else
                q2mu = mean(vp2,1);
                q2sigma = cov(vp2);
            end
        end
        [kls(1),kls(2)] = mvnkl(q1mu,q1sigma,q2mu,q2sigma);

    else
        MINP = realmin;

        xx1 = vbmc_rnd(vp1,Ns,origflag,1);
        q1 = vbmc_pdf(vp1,xx1,origflag);
        q2 = vbmc_pdf(vp2,xx1,origflag);
        q1(q1 == 0 | ~isfinite(q1)) = 1;    % Ignore these points
        q2(q2 == 0 | ~isfinite(q2)) = MINP;
        kls(1) = -mean(log(q2) - log(q1));

        xx2 = vbmc_rnd(vp2,Ns,origflag,1);
        q1 = vbmc_pdf(vp1,xx2,origflag);
        q2 = vbmc_pdf(vp2,xx2,origflag);
        q1(q1 == 0 | ~isfinite(q1)) = MINP;
        q2(q2 == 0 | ~isfinite(q2)) = 1;    % Ignore these points
        kls(2) = -mean(log(q1) - log(q2));

    end

    kls = max(kls,0); % Correct for numerical errors

%catch

    % Could not compute KL divs

%end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [source] = ft_sourcedescriptives(cfg, source)

% FT_SOURCEDESCRIPTIVES computes descriptive parameters of the source
% analysis results.
%
% Use as
%   [source] = ft_sourcedescriptives(cfg, source)
%
% where cfg is a structure with the configuration details and source is the
% result from a beamformer source estimation. The configuration can contain
%   cfg.cohmethod        = 'regular', 'lambda1', 'canonical'
%   cfg.powmethod        = 'regular', 'lambda1', 'trace', 'none'
%   cfg.supmethod        = 'chan_dip', 'chan', 'dip', 'none' (default)
%   cfg.projectmom       = 'yes' or 'no' (default = 'no')
%   cfg.eta              = 'yes' or 'no' (default = 'no')
%   cfg.kurtosis         = 'yes' or 'no' (default = 'no')
%   cfg.keeptrials       = 'yes' or 'no' (default = 'no')
%   cfg.keepcsd          = 'yes' or 'no' (default = 'no')
%   cfg.keepnoisecsd     = 'yes' or 'no' (default = 'no')
%   cfg.keepmom          = 'yes' or 'no' (default = 'yes')
%   cfg.keepnoisemom     = 'yes' or 'no' (default = 'yes')
%   cfg.resolutionmatrix = 'yes' or 'no' (default = 'no')
%   cfg.feedback         = 'no', 'text' (default), 'textbar', 'gui'
%
% The following option only applies to timecourses.
%   cfg.flipori          = 'yes' or 'no' (default = 'no')
%
% The following option only applies to single-trial timecourses.
%   cfg.fixedori         = 'within_trials' or 'over_trials' (default = 'over_trials')
%
% If repeated trials are present that have undergone some sort of
% resampling (i.e. jackknife, bootstrap, singletrial or rawtrial), the mean,
% variance and standard error of mean will be computed for all source
% parameters. This is done after applying the optional transformation
% on the power and projected noise.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_SOURCEANALYSIS, FT_SOURCESTATISTICS, FT_MATH

% Copyright (C) 2004-2015, Robert Oostenveld & Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar source
ft_preamble provenance source

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
% source = ft_checkdata(source, 'datatype', 'source', 'feedback', 'yes');

% cfg = ft_checkconfig(cfg, 'forbidden',   {'trials'});    % trial selection is not implented here, you may want to consider ft_selectdata

% DEPRECATED by roboos on 13 June 2013
% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2199 for more details
% support for this functionality can be removed at the end of 2013
cfg = ft_checkconfig(cfg, 'deprecated',  {'transform'}); % please use ft_math instead

% set the defaults
cfg.transform        = ft_getopt(cfg, 'transform',        []);
cfg.projectmom       = ft_getopt(cfg, 'projectmom',       'no'); % if yes -> svdfft
cfg.numcomp          = ft_getopt(cfg, 'numcomp',          1);
cfg.powmethod        = ft_getopt(cfg, 'powmethod',        []); % see below
cfg.cohmethod        = ft_getopt(cfg, 'cohmethod',        []); % see below
cfg.feedback         = ft_getopt(cfg, 'feedback',         'textbar');
cfg.supmethod        = ft_getopt(cfg, 'supmethod',        'none');
cfg.resolutionmatrix = ft_getopt(cfg, 'resolutionmatrix', 'no');
cfg.eta              = ft_getopt(cfg, 'eta',              'no');
cfg.fa               = ft_getopt(cfg, 'fa',               'no');
cfg.kurtosis         = ft_getopt(cfg, 'kurtosis',         'no');
cfg.keeptrials       = ft_getopt(cfg, 'keeptrials',       'no');
cfg.trials           = ft_getopt(cfg, 'trials',           'all');
cfg.keepcsd          = ft_getopt(cfg, 'keepcsd',          'no');
cfg.keepmom          = ft_getopt(cfg, 'keepmom',          'yes');
cfg.keepnoisecsd     = ft_getopt(cfg, 'keepnoisecsd',     'no');
cfg.keepnoisemom     = ft_getopt(cfg, 'keepnoisemom',     'yes');
cfg.fwhm             = ft_getopt(cfg, 'fwhm',             'no');
cfg.fwhmremovecenter = ft_getopt(cfg, 'fwhmremovecenter', 0);
cfg.fwhmmethod       = ft_getopt(cfg, 'fwhmmethod',       'barnes');
cfg.fwhmmaxdist      = ft_getopt(cfg, 'fwhmmaxdist',      []);
cfg.fixedori         = ft_getopt(cfg, 'fixedori',         'over_trials');
cfg.flipori          = ft_getopt(cfg, 'flipori',          'no');

% only works for mne
cfg.demean         = ft_getopt(cfg, 'demean',         'yes');
cfg.baselinewindow = ft_getopt(cfg, 'baselinewindow', [-inf 0]);
cfg.zscore         = ft_getopt(cfg, 'zscore',         'yes');

zscore = strcmp(cfg.zscore, 'yes');
demean = strcmp(cfg.demean, 'yes');

if ischar(cfg.trials) && strcmp(cfg.trials,'all')
  % do nothing
elseif ischar(cfg.trials)
  ft_error('only ''all'' is allowed for string input for cfg.trials');
else
  % check whether there's a trial field in the source structure, and
  % subselect, otherwise error
  if isfield(source, 'trial')
    source.trial = source.trial(cfg.trials);
    if isfield(source, 'cumtapcnt'), source.cumtapcnt = source.cumtapcnt(cfg.trials,:); end
  else
    ft_error('subselecting trials in ft_sourcedescriptives is currently only possible with a ''trial'' field');
  end
end

% get desired method from source structure
source.method = ft_getopt(source,'method',[]);

% this is required for backward compatibility with the old sourceanalysis
if isfield(source, 'method') && strcmp(source.method, 'randomized')
  source.method = 'randomization';
elseif isfield(source, 'method') && strcmp(source.method, 'permuted')
  source.method = 'permutation';
elseif isfield(source, 'method') && strcmp(source.method, 'jacknife')
  source.method = 'jackknife';
end

% determine the type of data, this is only relevant for a few specific types
ispccdata = isfield(source, 'avg')   && isfield(source.avg, 'csdlabel');
islcmvavg = isfield(source, 'avg')   && isfield(source, 'time') && isfield(source.avg,   'mom') && any(size(source.avg.pow)==1);
islcmvtrl = isfield(source, 'trial') && isfield(source, 'time') && isfield(source.trial, 'mom');
ismneavg  = isfield(source, 'avg')   && isfield(source, 'time') && isfield(source.avg,   'mom') && size(source.avg.pow, 2)==numel(source.time);

% check the consistency of the defaults
if strcmp(cfg.projectmom, 'yes')
  if isempty(cfg.powmethod)
    cfg.powmethod = 'regular'; % set the default
  elseif ~strcmp(cfg.powmethod, 'regular')
    ft_error('unsupported powmethod in combination with projectmom');
  end
  if isempty(cfg.cohmethod)
    cfg.cohmethod = 'regular'; % set the default
  elseif ~strcmp(cfg.cohmethod, 'regular')
    ft_error('unsupported cohmethod in combination with projectmom');
  end
else
  if isempty(cfg.powmethod)
    cfg.powmethod = 'lambda1'; % set the default
  end
  if isempty(cfg.cohmethod)
    cfg.cohmethod = 'lambda1'; % set the default
  end
end

% this is required for backward compatibility with an old version of sourcedescriptives
if isfield(cfg, 'singletrial'), cfg.keeptrials = cfg.singletrial;  end

% do a validity check on the input data and specified options
if strcmp(cfg.resolutionmatrix, 'yes')
  if ~isfield(source.avg, 'filter')
    ft_error('The computation of the resolution matrix requires keepfilter=''yes'' in sourceanalysis.');
  elseif ~isfield(source, 'leadfield')
    ft_error('The computation of the resolution matrix requires keepleadfield=''yes'' in sourceanalysis.');
  end
end

if istrue(cfg.fwhm) && ~isfield(source.avg, 'filter')
  ft_error('computation of the fwhm requires keepfilter=''yes'' in sourceanalysis.');
end

if istrue(cfg.eta) && strcmp(cfg.cohmethod, 'svdfft')
  ft_error('eta cannot be computed in combination with the application of svdfft');
end

if istrue(cfg.keeptrials) && ~strcmp(cfg.supmethod, 'none')
  ft_error('you cannot keep trials when you want to partialize something');
end

if istrue(cfg.flipori) && ~istrue(cfg.projectmom)
  ft_error('flipori requires projectmom=''yes''');
end

% set some flags for convenience
isnoise    = isfield(source, 'avg') && isfield(source.avg, 'noisecsd');
keeptrials = strcmp(cfg.keeptrials, 'yes');
projectmom = strcmp(cfg.projectmom, 'yes');
flipori    = strcmp(cfg.flipori, 'yes');

% determine the subfunction used for computing power
switch cfg.powmethod
  case 'regular'
    powmethodfun = @powmethod_regular;
  case 'lambda1'
    powmethodfun = @powmethod_lambda1;
  case 'trace'
    powmethodfun = @powmethod_trace;
  case 'none'
    powmethodfun = [];
  otherwise
    ft_error('unsupported powmethod');
end

% represent the selection of sources in the brain as a row-vector with indices
insideindx = find(source.inside(:)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ispccdata
  % the source reconstruction was computed using the pcc beamformer
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  Ndipole = size(source.pos,1);

  if ischar(source.avg.csdlabel{1}), source.avg.csdlabel = {source.avg.csdlabel}; end
  if numel(source.avg.csdlabel)==1
    source.avg.csdlabel = repmat(source.avg.csdlabel, [Ndipole 1]);
  end

  dipsel     = find(strcmp(source.avg.csdlabel{1}, 'scandip'));
  refchansel = find(strcmp(source.avg.csdlabel{1}, 'refchan'));
  refdipsel  = find(strcmp(source.avg.csdlabel{1}, 'refdip'));
  supchansel = find(strcmp(source.avg.csdlabel{1}, 'supchan'));
  supdipsel  = find(strcmp(source.avg.csdlabel{1}, 'supdip'));

  % cannot handle reference channels and reference dipoles simultaneously
  if numel(refchansel)>0 && numel(refdipsel)>0
    ft_error('cannot simultaneously handle reference channels and reference dipole');
  end

  % these are only used to count the number of reference/suppression dipoles and channels
  refsel = [refdipsel refchansel];
  supsel = [supdipsel supchansel];

  % first do the projection of the moment, if requested
  if projectmom
    source.avg.ori = cell(1, Ndipole);
    ft_progress('init', cfg.feedback, 'projecting dipole moment');
    for i=insideindx
      ft_progress(i/length(insideindx), 'projecting dipole moment %d/%d\n', i, length(insideindx));

      if numel(source.avg.csdlabel)>1
        dipsel     = find(strcmp(source.avg.csdlabel{i}, 'scandip'));
        refchansel = find(strcmp(source.avg.csdlabel{i}, 'refchan'));
        refdipsel  = find(strcmp(source.avg.csdlabel{i}, 'refdip'));
        supchansel = find(strcmp(source.avg.csdlabel{i}, 'supchan'));
        supdipsel  = find(strcmp(source.avg.csdlabel{i}, 'supdip'));

        % these are only used to count the number of reference/suppression dipoles and channels
        refsel = [refdipsel refchansel];
        supsel = [supdipsel supchansel];
      end

      mom     = source.avg.mom{i}(dipsel,     :);
      ref     = source.avg.mom{i}(refdipsel,  :);
      sup     = source.avg.mom{i}(supdipsel,  :);
      refchan = source.avg.mom{i}(refchansel, :);
      supchan = source.avg.mom{i}(supchansel, :);
      % compute the projection of the scanning dipole along the direction of the dominant amplitude
      if length(dipsel)>1, [mom, rmom]  = svdfft(mom, cfg.numcomp, source.cumtapcnt); else rmom = []; end
      source.avg.ori{i} = rmom;
      % compute the projection of the reference dipole along the direction of the dominant amplitude
      if length(refdipsel)>1, [ref, rref] = svdfft(ref, 1, source.cumtapcnt); else rref = []; end
      % compute the projection of the supression dipole along the direction of the dominant amplitude
      if length(supdipsel)>1, [sup, rsup] = svdfft(sup, 1, source.cumtapcnt); else rsup = []; end

      % compute voxel-level fourier-matrix
      source.avg.mom{i} = cat(1, mom, ref, sup, refchan, supchan);

      % create rotation-matrix
      rotmat = zeros(0, length(source.avg.csdlabel{i}));
      if ~isempty(rmom)
        rotmat = [rotmat; rmom zeros(1,numel(refsel)+numel(supsel))];
      end
      if ~isempty(rref)
        rotmat = [rotmat; zeros(1, numel(dipsel)), rref, zeros(1,numel(refchansel)+numel(supsel))];
      end
      if ~isempty(rsup)
        rotmat = [rotmat; zeros(1, numel(dipsel)+numel(refdipsel)), rsup, zeros(1,numel(refchansel)+numel(supchansel))];
      end
      for j=1:length(supchansel)
        rotmat(end+1,:) = 0;
        rotmat(end,numel(dipsel)+numel(refdipsel)+numel(supdipsel)+j) = 1;
      end
      for j=1:length(refchansel)
        rotmat(end+1,:) = 0;
        rotmat(end,numel(dipsel)+numel(refdipsel)+numel(supdipsel)+numel(supchansel)+j) = 1;
      end

      % compute voxel-level csd-matrix
      if isfield(source.avg, 'csd'), source.avg.csd{i}           = rotmat * source.avg.csd{i} * rotmat'; end
      % compute voxel-level noisecsd-matrix
      if isfield(source.avg, 'noisecsd'), source.avg.noisecsd{i} = rotmat * source.avg.noisecsd{i} * rotmat'; end
      % compute rotated filter
      if isfield(source.avg, 'filter'),   source.avg.filter{i}   = rotmat * source.avg.filter{i}; end
      if isfield(source.avg, 'csdlabel')
        % remember what the interpretation is of all CSD output components
        scandiplabel = repmat({'scandip'}, 1, cfg.numcomp);          % only one dipole orientation remains
        refdiplabel  = repmat({'refdip'},  1, length(refdipsel)>0);  % for svdfft at max. only one dipole orientation remains
        supdiplabel  = repmat({'supdip'},  1, length(supdipsel)>0);  % for svdfft at max. only one dipole orientation remains
        refchanlabel = repmat({'refchan'}, 1, length(refchansel));
        supchanlabel = repmat({'supchan'}, 1, length(supchansel));
        % concatenate all the labels
        source.avg.csdlabel{i} = cat(2, scandiplabel, refdiplabel, supdiplabel, refchanlabel, supchanlabel);
      end

      % compute rotated leadfield
      % FIXME in the presence of a refdip and/or supdip, this does not work; leadfield is Nx3
      if isfield(source,  'leadfield')
        % FIXME this is a proposed dirty fix
        n1 = size(source.leadfield{i},2);
        %n2 = size(rotmat,2) - n1;
        n2 = size(rotmat,2) - n1 +1; %added 1 JM
        source.leadfield{i}    = source.leadfield{i} * rotmat(1:n2, 1:n1)';
      end
    end % for i=insideindx
    ft_progress('close');

    % update the indices
    dipsel     = find(strcmp(source.avg.csdlabel, 'scandip'));
    refchansel = find(strcmp(source.avg.csdlabel, 'refchan'));
    refdipsel  = find(strcmp(source.avg.csdlabel, 'refdip'));
    supchansel = find(strcmp(source.avg.csdlabel, 'supchan'));
    supdipsel  = find(strcmp(source.avg.csdlabel, 'supdip'));
    refsel     = [refdipsel refchansel];
    supsel     = [supdipsel supchansel];
  end % if projectmom

  if keeptrials
    cumtapcnt = source.cumtapcnt(:);
    sumtapcnt = cumsum([0;cumtapcnt]);
    Ntrial = length(cumtapcnt);

    ft_progress('init', cfg.feedback, 'computing singletrial voxel-level cross-spectral densities');
    for triallop = 1:Ntrial
      source.trial(triallop).csd = cell(Ndipole, 1);  % allocate memory for this trial
      source.trial(triallop).mom = cell(Ndipole, 1);  % allocate memory for this trial

      ft_progress(triallop/Ntrial, 'computing singletrial voxel-level cross-spectral densities %d%d\n', triallop, Ntrial);
      for i=insideindx
        dat = source.avg.mom{i};
        tmpmom = dat(:, sumtapcnt(triallop)+1:sumtapcnt(triallop+1));
        tmpcsd = (tmpmom * tmpmom') ./cumtapcnt(triallop);
        source.trial(triallop).mom{i} = tmpmom;
        source.trial(triallop).csd{i} = tmpcsd;
      end % for i=insideindx
    end % for triallop
    ft_progress('close');
    % remove the average, continue with separate trials, but keep track of
    % the csdlabel
    csdlabel = source.avg.csdlabel;
    source   = rmfield(source, 'avg');
  else
    fprintf('using average voxel-level cross-spectral densities\n');
    csdlabel = source.avg.csdlabel;
  end % if keeptrials

  % process the csdlabel for each of the dipoles
  hasrefdip  = true;
  hasrefchan = true;
  hassupdip  = true;
  hassupchan = true;

  dipselcell     = cell(Ndipole,1);
  refdipselcell  = cell(Ndipole,1);
  refchanselcell = cell(Ndipole,1);
  supdipselcell  = cell(Ndipole,1);
  supchanselcell = cell(Ndipole,1);

  for i = insideindx
    dipsel     = find(strcmp(csdlabel{i}, 'scandip'));
    refchansel = find(strcmp(csdlabel{i}, 'refchan'));
    refdipsel  = find(strcmp(csdlabel{i}, 'refdip'));
    supchansel = find(strcmp(csdlabel{i}, 'supchan'));
    supdipsel  = find(strcmp(csdlabel{i}, 'supdip'));

    hasrefdip  = ~isempty(refdipsel)  && hasrefdip; %NOTE: it has to be true for all dipoles!
    hasrefchan = ~isempty(refchansel) && hasrefchan;
    hassupdip  = ~isempty(supdipsel)  && hassupdip;
    hassupchan = ~isempty(supchansel) && hassupchan;

    dipselcell{i}     = dipsel;
    refdipselcell{i}  = refdipsel;
    refchanselcell{i} = refchansel;
    supdipselcell{i}  = supdipsel;
    supchanselcell{i} = supchansel;
  end

  if keeptrials
    % do the processing of the CSD matrices for each trial
    if ~strcmp(cfg.supmethod, 'none')
      ft_error('suppression is only supported for average CSD');
    end
    %dipselcell = mat2cell(repmat(dipsel(:)', [Ndipole 1]), ones(Ndipole,1), length(dipsel));
    %if hasrefdip,  refdipselcell  = mat2cell(repmat(refdipsel(:)',  [Ndipole 1]), ones(Ndipole,1), length(refdipsel));  end
    %if hasrefchan, refchanselcell = mat2cell(repmat(refchansel(:)', [Ndipole 1]), ones(Ndipole,1), length(refchansel)); end
    %if hassupdip,  supdipselcell  = mat2cell(repmat(supdipsel(:)',  [Ndipole 1]), ones(Ndipole,1), length(supdipsel));  end
    %if hassupchan, supchanselcell = mat2cell(repmat(supchansel(:)', [Ndipole 1]), ones(Ndipole,1), length(supchansel)); end

    ft_progress('init', cfg.feedback, 'computing singletrial voxel-level power');
    for triallop = 1:Ntrial
      %initialize the variables
      source.trial(triallop).pow = zeros(Ndipole, 1);
      if hasrefdip,  source.trial(triallop).refdippow     = zeros(Ndipole, 1); end
      if hasrefchan, source.trial(triallop).refchanpow    = zeros(Ndipole, 1); end
      if hassupdip,  source.trial(triallop).supdippow     = zeros(Ndipole, 1); end
      if hassupchan, source.trial(triallop).supchanpow    = zeros(Ndipole, 1); end

      ft_progress(triallop/Ntrial, 'computing singletrial voxel-level power %d%d\n', triallop, Ntrial);
      source.trial(triallop).pow(source.inside) = cellfun(powmethodfun, source.trial(triallop).csd(source.inside), dipselcell(source.inside));
      if hasrefdip,  source.trial(triallop).refdippow(source.inside)  = cellfun(powmethodfun,source.trial(triallop).csd(source.inside), refdipselcell(source.inside));  end
      if hassupdip,  source.trial(triallop).supdippow(source.inside)  = cellfun(powmethodfun,source.trial(triallop).csd(source.inside), supdipselcell(source.inside));  end
      if hasrefchan, source.trial(triallop).refchanpow(source.inside) = cellfun(powmethodfun,source.trial(triallop).csd(source.inside), refchanselcell(source.inside)); end
      if hassupchan, source.trial(triallop).supchanpow(source.inside) = cellfun(powmethodfun,source.trial(triallop).csd(source.inside), supchanselcell(source.inside)); end
      %FIXME kan volgens mij niet
      if isnoise && isfield(source.trial(triallop), 'noisecsd')
        % compute the power of the noise projected on each source component
        source.trial(triallop).noise = cellfun(powmethodfun,source.trial(triallop).csd, dipselcell);
        if hasrefdip,  source.trial(triallop).refdipnoise  = cellfun(powmethodfun,source.trial(triallop).noisecsd, refdipselcell);  end
        if hassupdip,  source.trial(triallop).supdipnoise  = cellfun(powmethodfun,source.trial(triallop).noisecsd, supdipselcell);  end
        if hasrefchan, source.trial(triallop).refchannoise = cellfun(powmethodfun,source.trial(triallop).noisecsd, refchanselcell); end
        if hassupchan, source.trial(triallop).supchannoise = cellfun(powmethodfun,source.trial(triallop).noisecsd, supchanselcell); end
      end % if isnoise
    end % for triallop
    ft_progress('close');

    if strcmp(cfg.keepcsd, 'no')
      source.trial = rmfield(source.trial, 'csd');
    end

  else
    % do the processing of the average CSD matrix
    for i=insideindx
      switch cfg.supmethod
        case 'chan_dip'
          supindx = [supdipsel supchansel];
          if i==insideindx(1), refsel  = refsel - length(supdipsel); end % adjust index only once
          refchanselcell{i} = refchanselcell{i} - length(supdipsel);
          refdipselcell{i}  = refdipselcell{i} - length(supdipsel);
        case 'chan'
          supindx = supchansel;
        case 'dip'
          supindx = supdipsel;
          if i==insideindx(1), refsel  = refsel - length(supdipsel); end
          refchanselcell{i} = refchanselcell{i} - length(supdipsel);
          refdipselcell{i}  = refdipselcell{i} - length(supdipsel);
        case 'none'
          % do nothing
          supindx = [];
      end
      tmpcsd  = source.avg.csd{i};
      scnindx = setdiff(1:size(tmpcsd,1), supindx);
      tmpcsd  = tmpcsd(scnindx, scnindx) - tmpcsd(scnindx, supindx)*pinv(tmpcsd(supindx, supindx))*tmpcsd(supindx, scnindx);
      source.avg.csd{i}   = tmpcsd;
    end % for insideindx
    % source.avg.csdlabel = source.avg.csdlabel(scnindx);

    if isnoise && ~strcmp(cfg.supmethod, 'none')
      source.avg = rmfield(source.avg, 'noisecsd');
    end

    % initialize the variables
    source.avg.pow           = nan(Ndipole, 1);
    if hasrefdip,  source.avg.refdippow     = nan(Ndipole, 1); end
    if hasrefchan, source.avg.refchanpow    = nan(Ndipole, 1); end
    if hassupdip,  source.avg.supdippow     = nan(Ndipole, 1); end
    if hassupchan, source.avg.supchanpow    = nan(Ndipole, 1); end
    if isnoise
      source.avg.noise         = nan(Ndipole, 1);
      if hasrefdip,  source.avg.refdipnoise     = nan(Ndipole, 1); end
      if hasrefchan, source.avg.refchannoise    = nan(Ndipole, 1); end
      if hassupdip,  source.avg.supdipnoise     = nan(Ndipole, 1); end
      if hassupchan, source.avg.supchannoise    = nan(Ndipole, 1); end
    end % if isnoise
    if hasrefdip||hasrefchan, source.avg.coh    = nan(Ndipole, 1); end
    if strcmp(cfg.eta, 'yes')
      source.avg.eta           = nan(Ndipole, 1);
      source.avg.ori             = cell(1, Ndipole);
    end
    if strcmp(cfg.eta, 'yes') && ~isempty(refsel)
      source.avg.etacsd = nan(Ndipole, 1);
      source.avg.ucsd   = cell(1, Ndipole);
    end
    if strcmp(cfg.fa, 'yes')
      source.avg.fa = nan(Ndipole, 1);
    end

    for i=insideindx
      dipsel = dipselcell{i};
      refchansel = refchanselcell{i};
      refdipsel  = refdipselcell{i};
      refsel     = [refchansel refdipsel];
      supchansel = supchanselcell{i};

      % compute the power of each source component
      if strcmp(cfg.projectmom, 'yes') && cfg.numcomp>1
        source.avg.pow(i) = powmethodfun(source.avg.csd{i}(dipselcell{i},dipselcell{i}), 1);
      else
        source.avg.pow(i) = powmethodfun(source.avg.csd{i}(dipselcell{i},dipselcell{i}));
      end

      if hasrefdip,  source.avg.refdippow(i)  = powmethodfun(source.avg.csd{i}(refdipsel,refdipsel));   end
      %if hassupdip,  source.avg.supdippow(i)  = powmethodfun(source.avg.csd{i}(supdipsel,supdipsel));   end
      if hasrefchan, source.avg.refchanpow(i) = powmethodfun(source.avg.csd{i}(refchansel,refchansel)); end
      %if hassupchan, source.avg.supchanpow(i) = powmethodfun(source.avg.csd{i}(supchansel,supchansel)); end
      if isnoise
        % compute the power of the noise projected on each source component
        if strcmp(cfg.projectmom, 'yes') && cfg.numcomp>1
          source.avg.noise(i) = powmethodfun(source.avg.noisecsd{i}(dipselcell{i},dipselcell{i}), 1);
        else
          source.avg.noise(i) = powmethodfun(source.avg.noisecsd{i}(dipselcell{i},dipselcell{i}));
        end
        if hasrefdip,  source.avg.refdipnoise(i)  = powmethodfun(source.avg.noisecsd{i}(refdipsel,refdipsel));   end
        if hassupdip,  source.avg.supdipnoise(i)  = powmethodfun(source.avg.noisecsd{i}(supdipsel,supdipsel));   end
        if hasrefchan, source.avg.refchannoise(i) = powmethodfun(source.avg.noisecsd{i}(refchansel,refchansel)); end
        if hassupchan, source.avg.supchannoise(i) = powmethodfun(source.avg.noisecsd{i}(supchansel,supchansel)); end
      end % if isnoise

      if ~isempty(refsel)
        % compute coherence
        csd = source.avg.csd{i};
        switch cfg.cohmethod
          case 'regular'
            % assume that all dipoles have been projected along the direction of maximum power
            Pd                = abs(csd(dipsel, dipsel));
            Pr                = abs(csd(refsel, refsel));
            Cdr               = csd(dipsel, refsel);
            source.avg.coh(i) = (Cdr.^2) ./ (Pd*Pr);
          case 'lambda1'
            % compute coherence the Joachim Gross' way
            Pd                = lambda1(csd(dipsel, dipsel));
            Pr                = lambda1(csd(refsel, refsel));
            Cdr               = lambda1(csd(dipsel, refsel));
            source.avg.coh(i) = abs(Cdr).^2 ./ (Pd*Pr);
          case 'canonical'
            % compute canonical coherence

            ccoh = ft_connectivity_cancorr(csd([dipsel refsel],[dipsel refsel]), 'indices', [ones(1,numel(dipsel)) ones(1,numel(refsel))*2]);
            source.avg.coh(i)  = ccoh(1,2);
          otherwise
            ft_error('unsupported cohmethod');
        end % cohmethod
      end

      % compute eta
      if strcmp(cfg.eta, 'yes')
        [source.avg.eta(i), source.avg.ori{i}] = csd2eta(source.avg.csd{i}(dipselcell{i},dipselcell{i}));
        if ~isempty(refsel)
          %FIXME this only makes sense when only a reference signal OR a dipole is selected
          [source.avg.etacsd(i), source.avg.ucsd{i}] = csd2eta(source.avg.csd{i}(dipsel,refsel));
        end
      end

      %compute fa
      if strcmp(cfg.fa, 'yes')
        source.avg.fa(i) = csd2fa(source.avg.csd{i}(dipsel,dipsel));
      end
    end % for diplop

    if strcmp(cfg.keepcsd, 'no')
      source.avg = rmfield(source.avg, 'csd');
    end
    if strcmp(cfg.keepnoisecsd, 'no') && isnoise
      source.avg = rmfield(source.avg, 'noisecsd');
    end

  end % if keeptrials

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ismneavg
  % the source reconstruction was computed using mne and contains an average timecourse
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if demean
    begsmp = nearest(source.time, cfg.baselinewindow(1));
    endsmp = nearest(source.time, cfg.baselinewindow(2));
    ft_progress('init', cfg.feedback, 'baseline correcting dipole moments');
    for diplop=1:length(insideindx)
      ft_progress(diplop/length(insideindx), 'baseline correcting dipole moments %d/%d\n', diplop, length(insideindx));
      mom = source.avg.mom{insideindx(diplop)};
      mom = ft_preproc_baselinecorrect(mom, begsmp, endsmp);
      source.avg.mom{insideindx(diplop)} = mom;
    end
    ft_progress('close');
  end

  if projectmom
    if isfield(source, 'tri')
      nrm = surface_normals(source.pos, source.tri, 'vertex');
      source.avg.phi = zeros(size(source.pos,1),1);
    end
    ft_progress('init', cfg.feedback, 'projecting dipole moment');
    for diplop=1:length(insideindx)
      ft_progress(diplop/length(insideindx), 'projecting dipole moment %d/%d\n', diplop, length(insideindx));
      mom = source.avg.mom{insideindx(diplop)};
      [mom, rmom] = svdfft(mom, 1);
      source.avg.mom{insideindx(diplop)} = mom;
      source.avg.ori{insideindx(diplop)} = rmom;
    end
    if isfield(source, 'tri')
      for diplop=insideindx
        source.avg.phi(diplop) = source.avg.ori{diplop}*nrm(diplop,:)';
      end
    end
    if isfield(source.avg, 'noisecov')
      source.avg.noise = nan+zeros(size(source.pos,1),1);
      for diplop=insideindx
        rmom = source.avg.ori{diplop};
        source.avg.noise(diplop) = rmom*source.avg.noisecov{diplop}*rmom';
      end
    end
    ft_progress('close');
  end % if projectmom

  if flipori
    tmpmom = cat(1, source.avg.mom{source.inside});
    [u, s, v] = svd(tmpmom, 'econ');
    flip( source.inside) = sign(u(:,1));
    flip(~source.inside) = nan;
    for i=1:numel(source.inside)
      if source.inside(i)
        source.avg.mom{i} = flip(i) * source.avg.mom{i};
        source.avg.ori{i} = flip(i) * source.avg.ori{i};
      end
    end
  end % if flipori

  if zscore
    begsmp = nearest(source.time, cfg.baselinewindow(1));
    endsmp = nearest(source.time, cfg.baselinewindow(2));
    % zscore using baselinewindow for power
    ft_progress('init', cfg.feedback, 'computing power');
    %source.avg.absmom = source.avg.pow;
    for diplop=1:length(insideindx)
      ft_progress(diplop/length(insideindx), 'computing power %d/%d\n', diplop, length(insideindx));
      mom = source.avg.mom{insideindx(diplop)};
      mmom = mean(mom(:,begsmp:endsmp),2);
      smom = std(mom(:,begsmp:endsmp),[],2);
      pow  = sum(((mom-mmom(:,ones(size(mom,2),1)))./smom(:,ones(size(mom,2),1))).^2,1);
      source.avg.pow(insideindx(diplop),:) = pow;
      source.avg.mom{insideindx(diplop)}   = diag(1./smom)*mom;
    end
    ft_progress('close');
  else
    % just square for power
    ft_progress('init', cfg.feedback, 'computing power');
    %source.avg.absmom = source.avg.pow;
    for diplop=1:length(insideindx)
      ft_progress(diplop/length(insideindx), 'computing power %d/%d\n', diplop, length(insideindx));
      mom = source.avg.mom{insideindx(diplop)};
      pow = sum(mom.^2,1);
      source.avg.pow(insideindx(diplop),:) = pow;
      %source.avg.absmom(insideindx(diplop),:) = sum(mom,1);
    end
    ft_progress('close');
  end % if zscore

  if strcmp(cfg.kurtosis, 'yes')
    fprintf('computing kurtosis based on dipole timecourse\n');
    source.avg.k2 = nan(size(source.pos,1),1);
    for diplop=1:length(insideindx)
      mom = source.avg.mom{insideindx(diplop)};
      if length(mom)~=prod(size(mom))
        ft_error('kurtosis can only be computed for projected dipole moment');
      end
      source.avg.k2(insideindx(diplop)) = kurtosis(mom);
    end
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif islcmvavg
  % the source reconstruction was computed using the lcmv beamformer and contains an average timecourse
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if projectmom
    ft_progress('init', cfg.feedback, 'projecting dipole moment');
    for diplop=1:length(insideindx)
      ft_progress(diplop/length(insideindx), 'projecting dipole moment %d/%d\n', diplop, length(insideindx));
      mom = source.avg.mom{insideindx(diplop)};
      [mom, rmom] = svdfft(mom, 1);
      source.avg.mom{insideindx(diplop)} = mom;
      source.avg.ori{insideindx(diplop)} = rmom;
    end
    ft_progress('close');
  end

  if flipori
    tmpmom = cat(1, source.avg.mom{source.inside});
    [u, s, v] = svd(tmpmom, 'econ');
    flip( source.inside) = sign(u(:,1));
    flip(~source.inside) = nan;
    for i=1:numel(source.inside)
      if source.inside(i)
        source.avg.mom{i} = flip(i) * source.avg.mom{i};
        source.avg.ori{i} = flip(i) * source.avg.ori{i};
      end
    end
  end

  if ~strcmp(cfg.powmethod, 'none')
    fprintf('recomputing power based on dipole timecourse\n')
    source.avg.pow = nan(size(source.pos,1),1);
    for diplop=1:length(insideindx)
      mom = source.avg.mom{insideindx(diplop)};
      cov = mom * mom';
      source.avg.pow(insideindx(diplop)) = powmethodfun(cov);
    end
  end

  if strcmp(cfg.kurtosis, 'yes')
    fprintf('computing kurtosis based on dipole timecourse\n');
    source.avg.k2 = nan(size(source.pos,1),1);
    for diplop=1:length(insideindx)
      mom = source.avg.mom{insideindx(diplop)};
      if length(mom)~=prod(size(mom))
        ft_error('kurtosis can only be computed for projected dipole moment');
      end
      source.avg.k2(insideindx(diplop)) = kurtosis(mom);
    end
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif islcmvtrl
  % the source reconstruction was computed using the lcmv beamformer and contains a single-trial timecourse
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  ntrial = length(source.trial);

  if projectmom && strcmp(cfg.fixedori, 'within_trials')
    % the dipole orientation is re-determined for each trial
    ft_progress('init', cfg.feedback, 'projecting dipole moment');
    for trllop=1:ntrial
      ft_progress(trllop/ntrial, 'projecting dipole moment %d/%d\n', trllop, ntrial);
      for diplop=1:length(insideindx)
        mom = source.trial(trllop).mom{insideindx(diplop)};
        [mom, rmom] = svdfft(mom, 1);
        source.trial(trllop).mom{insideindx(diplop)} = mom;
        source.trial(trllop).ori{insideindx(diplop)} = rmom;  % remember the orientation
      end
    end
    ft_progress('close');
  elseif projectmom && strcmp(cfg.fixedori, 'over_trials')
    ft_progress('init', cfg.feedback, 'projecting dipole moment');
    % compute average covariance over all trials
    for trllop=1:ntrial
      for diplop=1:length(insideindx)
        mom = source.trial(trllop).mom{insideindx(diplop)};
        if trllop==1
          cov{diplop} = mom*mom'./size(mom,2);
        else
          cov{diplop} = mom*mom'./size(mom,2) + cov{diplop};
        end
      end
    end
    % compute source orientation over all trials
    for diplop=1:length(insideindx)
      [dum, ori{diplop}] = svdfft(cov{diplop}, 1);
    end
    % project the data in each trial
    for trllop=1:ntrial
      ft_progress(trllop/ntrial, 'projecting dipole moment %d/%d\n', trllop, ntrial);
      for diplop=1:length(insideindx)
        mom = source.trial(trllop).mom{insideindx(diplop)};
        mom = ori{diplop}*mom;
        source.trial(trllop).mom{insideindx(diplop)} = mom;
        source.trial(trllop).ori{insideindx(diplop)} = ori{diplop};
      end
    end
    ft_progress('close');
  end

  if ~strcmp(cfg.powmethod, 'none')
    fprintf('recomputing power based on dipole timecourse\n')
    for trllop=1:ntrial
      for diplop=1:length(insideindx)
        mom = source.trial(trllop).mom{insideindx(diplop)};
        cov = mom * mom';
        source.trial(trllop).pow(insideindx(diplop)) = powmethodfun(cov);
      end
    end
  end

  if strcmp(cfg.kurtosis, 'yes')
    fprintf('computing kurtosis based on dipole timecourse\n');
    for trllop=1:ntrial
      source.trial(trllop).k2 = nan(size(source.pos,1),1);
      for diplop=1:length(insideindx)
        mom = source.trial(trllop).mom{insideindx(diplop)};
        if length(mom)~=numel(mom)
          ft_error('kurtosis can only be computed for projected dipole moment');
        end
        source.trial(trllop).k2(insideindx(diplop)) = kurtosis(mom);
      end
    end
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % dealing with pcc, lcmv, dics or mne input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isfield(source, 'avg') && isfield(source.avg, 'pow') && isfield(source.avg, 'noise') && ~ismneavg
  % compute the neural activity index for the average
  source.avg.nai = source.avg.pow(:) ./ source.avg.noise(:);
end

if isfield(source, 'trial') && isfield(source.trial, 'pow') && isfield(source.trial, 'noise')
  % compute the neural activity index for the trials
  ntrials = length(source.trial);
  for trlop=1:ntrials
    source.trial(trlop).nai = source.trial(trlop).pow ./ source.trial(trlop).noise;
  end
end

if strcmp(source.method, 'randomization') || strcmp(source.method, 'permutation')
  % compute the neural activity index for the two randomized conditions
  source.avgA.nai = source.avgA.pow ./ source.avgA.noise;
  source.avgB.nai = source.avgB.pow ./ source.avgB.noise;
  for trlop=1:length(source.trialA)
    source.trialA(trlop).nai = source.trialA(trlop).pow ./ source.trialA(trlop).noise;
  end
  for trlop=1:length(source.trialB)
    source.trialB(trlop).nai = source.trialB(trlop).pow ./ source.trialB(trlop).noise;
  end
end

if ~isempty(cfg.transform)
  fprintf('applying %s transformation on the power and projected noise\n', cfg.transform);
  % apply the specified transformation on the power
  if isfield(source, 'avg'   ) && isfield(source.avg   , 'pow'), source.avg .pow = feval(cfg.transform, source.avg .pow); end
  if isfield(source, 'avgA'  ) && isfield(source.avgA  , 'pow'), source.avgA.pow = feval(cfg.transform, source.avgA.pow); end
  if isfield(source, 'avgB'  ) && isfield(source.avgB  , 'pow'), source.avgB.pow = feval(cfg.transform, source.avgB.pow); end
  if isfield(source, 'trial' ) && isfield(source.trial , 'pow'), for i=1:length(source.trial ), source.trial (i).pow = feval(cfg.transform, source.trial (i).pow); end; end
  if isfield(source, 'trialA') && isfield(source.trialA, 'pow'), for i=1:length(source.trialA), source.trialA(i).pow = feval(cfg.transform, source.trialA(i).pow); end; end
  if isfield(source, 'trialB') && isfield(source.trialB, 'pow'), for i=1:length(source.trialB), source.trialB(i).pow = feval(cfg.transform, source.trialB(i).pow); end; end
  % apply the specified transformation on the projected noise
  if isfield(source, 'avg'   ) && isfield(source.avg   , 'noise'), source.avg .noise = feval(cfg.transform, source.avg .noise); end
  if isfield(source, 'avgA'  ) && isfield(source.avgA  , 'noise'), source.avgA.noise = feval(cfg.transform, source.avgA.noise); end
  if isfield(source, 'avgB'  ) && isfield(source.avgB  , 'noise'), source.avgB.noise = feval(cfg.transform, source.avgB.noise); end
  if isfield(source, 'trial' ) && isfield(source.trial , 'noise'), for i=1:length(source.trial ), source.trial (i).noise = feval(cfg.transform, source.trial (i).noise); end; end
  if isfield(source, 'trialA') && isfield(source.trialA, 'noise'), for i=1:length(source.trialA), source.trialA(i).noise = feval(cfg.transform, source.trialA(i).noise); end; end
  if isfield(source, 'trialB') && isfield(source.trialB, 'noise'), for i=1:length(source.trialB), source.trialB(i).noise = feval(cfg.transform, source.trialB(i).noise); end; end
end

if strcmp(source.method, 'pseudovalue')
  % compute the pseudovalues for the beamformer output
  avg = source.trial(1);        % the first is the complete average
  Ntrials = length(source.trial)-1; % the remaining are the leave-one-out averages
  pseudoval = [];
  if isfield(source.trial, 'pow')
    allavg = getfield(avg, 'pow');
    for i=1:Ntrials
      thisavg = getfield(source.trial(i+1), 'pow');
      thisval = Ntrials*allavg - (Ntrials-1)*thisavg;
      pseudoval(i).pow = thisval;
    end
  end
  if isfield(source.trial, 'coh')
    allavg = getfield(avg, 'coh');
    for i=1:Ntrials
      thisavg = getfield(source.trial(i+1), 'coh');
      thisval = Ntrials*allavg - (Ntrials-1)*thisavg;
      pseudoval(i).coh = thisval;
    end
  end
  if isfield(source.trial, 'nai')
    allavg = getfield(avg, 'nai');
    for i=1:Ntrials
      thisavg = getfield(source.trial(i+1), 'nai');
      thisval = Ntrials*allavg - (Ntrials-1)*thisavg;
      pseudoval(i).nai = thisval;
    end
  end
  if isfield(source.trial, 'noise')
    allavg = getfield(avg, 'noise');
    for i=1:Ntrials
      thisavg = getfield(source.trial(i+1), 'noise');
      thisval = Ntrials*allavg - (Ntrials-1)*thisavg;
      pseudoval(i).noise = thisval;
    end
  end
  % store the pseudovalues instead of the original values
  source.trial = pseudoval;
end

if strcmp(source.method, 'jackknife') || strcmp(source.method, 'bootstrap') || strcmp(source.method, 'pseudovalue') || strcmp(source.method, 'singletrial') || strcmp(source.method, 'rawtrial')
  % compute descriptive statistics (mean, var, sem) for multiple trial data
  % compute these for as many source parameters as possible

  % for convenience copy the trials out of the source structure
  dip = source.trial;

  % determine the (original) number of trials in the data
  if strcmp(source.method, 'bootstrap') %VERANDERD ER ZAT GEEN .RESAMPLE IN SOURCE
    Ntrials = size(source.trial,2); % WAS size(source.resample, 2);
  else
    Ntrials = length(source.trial);
  end
  fprintf('original data contained %d trials\n', Ntrials);

  % allocate memory for all elements in the dipole structure
  sumdip = [];
  if isfield(dip(1), 'var'),   sumdip.var    = zeros(size(dip(1).var  )); sumdip.var(~source.inside) = nan; end
  if isfield(dip(1), 'pow'),   sumdip.pow    = zeros(size(dip(1).pow  )); sumdip.pow(~source.inside) = nan; end
  if isfield(dip(1), 'coh'),   sumdip.coh    = zeros(size(dip(1).coh  )); sumdip.coh(~source.inside) = nan; end
  if isfield(dip(1), 'rv'),    sumdip.rv     = zeros(size(dip(1).rv   )); sumdip.rv(~source.inside) = nan; end
  if isfield(dip(1), 'noise'), sumdip.noise  = zeros(size(dip(1).noise)); sumdip.noise(~source.inside) = nan; end
  if isfield(dip(1), 'nai'),   sumdip.nai    = zeros(size(dip(1).nai  )); sumdip.nai(~source.inside) = nan; end
  sqrdip = [];
  if isfield(dip(1), 'var'),   sqrdip.var    = zeros(size(dip(1).var  )); sqrdip.var(~source.inside) = nan; end
  if isfield(dip(1), 'pow'),   sqrdip.pow    = zeros(size(dip(1).pow  )); sqrdip.pow(~source.inside) = nan; end
  if isfield(dip(1), 'coh'),   sqrdip.coh    = zeros(size(dip(1).coh  )); sqrdip.coh(~source.inside) = nan; end
  if isfield(dip(1), 'rv'),    sqrdip.rv     = zeros(size(dip(1).rv   )); sqrdip.rv(~source.inside) = nan; end
  if isfield(dip(1), 'noise'), sqrdip.noise  = zeros(size(dip(1).noise)); sqrdip.noise(~source.inside) = nan; end
  if isfield(dip(1), 'nai'),   sqrdip.nai    = zeros(size(dip(1).nai  )); sqrdip.nai(~source.inside) = nan; end
  if isfield(dip(1), 'mom')
    sumdip.mom = cell(size(dip(1).mom));
    sqrdip.mom = cell(size(dip(1).mom));
    for i=1:length(dip(1).mom)
      sumdip.mom{i} = zeros(size(dip(1).mom{i}));
      sqrdip.mom{i} = zeros(size(dip(1).mom{i}));
    end
  end
  if isfield(dip(1), 'csd')
    sumdip.csd = cell(size(dip(1).csd));
    sqrdip.csd = cell(size(dip(1).csd));
    for i=1:length(dip(1).csd)
      sumdip.csd{i} = zeros(size(dip(1).csd{i}));
      sqrdip.csd{i} = zeros(size(dip(1).csd{i}));
    end
  end

  for trial=1:length(dip)
    % compute the sum of all values
    if isfield(dip(trial), 'var'),    sumdip.var   = sumdip.var    + dip(trial).var;    end
    if isfield(dip(trial), 'pow'),    sumdip.pow   = sumdip.pow    + dip(trial).pow;    end
    if isfield(dip(trial), 'coh'),    sumdip.coh   = sumdip.coh    + dip(trial).coh;    end
    if isfield(dip(trial), 'rv'),     sumdip.rv    = sumdip.rv     + dip(trial).rv;     end
    if isfield(dip(trial), 'noise'),  sumdip.noise = sumdip.noise  + dip(trial).noise;  end
    if isfield(dip(trial), 'nai'),    sumdip.nai   = sumdip.nai    + dip(trial).nai;    end
    % compute the sum of squared values
    if isfield(dip(trial), 'var'),    sqrdip.var    = sqrdip.var   + (dip(trial).var  ).^2; end
    if isfield(dip(trial), 'pow'),    sqrdip.pow    = sqrdip.pow   + (dip(trial).pow  ).^2; end
    if isfield(dip(trial), 'coh'),    sqrdip.coh    = sqrdip.coh   + (dip(trial).coh  ).^2; end
    if isfield(dip(trial), 'rv'),     sqrdip.rv     = sqrdip.rv    + (dip(trial).rv   ).^2; end
    if isfield(dip(trial), 'noise'),  sqrdip.noise  = sqrdip.noise + (dip(trial).noise).^2; end
    if isfield(dip(trial), 'nai'),    sqrdip.nai    = sqrdip.nai   + (dip(trial).nai  ).^2; end
    % do the same for the cell-array with mom
    if isfield(dip(trial), 'mom')
      for i=1:length(dip(1).mom)
        sumdip.mom{i} = sumdip.mom{i} +  dip(trial).mom{i};
        sqrdip.mom{i} = sqrdip.mom{i} + (dip(trial).mom{i}).^2;
      end
    end
    % do the same for the cell-array with csd
    if isfield(dip(trial), 'csd')
      for i=1:length(dip(1).csd)
        sumdip.csd{i} = sumdip.csd{i} +  dip(trial).csd{i};
        sqrdip.csd{i} = sqrdip.csd{i} + (dip(trial).csd{i}).^2;
      end
    end
  end

  % compute the mean over all repetitions
  if isfield(sumdip, 'var'),    dipmean.var    = sumdip.var   / length(dip); end
  if isfield(sumdip, 'pow'),    dipmean.pow    = sumdip.pow   / length(dip); end
  if isfield(sumdip, 'coh'),    dipmean.coh    = sumdip.coh   / length(dip); end
  if isfield(sumdip, 'rv'),     dipmean.rv     = sumdip.rv    / length(dip); end
  if isfield(sumdip, 'noise'),  dipmean.noise  = sumdip.noise / length(dip); end
  if isfield(sumdip, 'nai'),    dipmean.nai    = sumdip.nai   / length(dip); end
  % for the cell-array with mom, this is done further below
  % for the cell-array with csd, this is done further below

  % the estimates for variance and SEM are biased if we are working with the jackknife/bootstrap
  % determine the proper variance scaling that corrects for this bias
  % note that Ntrials is not always the same as the length of dip, especially in case of the bootstrap
  if strcmp(source.method, 'singletrial')
    bias = 1;
  elseif strcmp(source.method, 'rawtrial')
    bias = 1;
  elseif strcmp(source.method, 'jackknife')
    % Effron gives SEM estimate for the jackknife method in equation 11.5 (paragraph 11.2)
    % to get the variance instead of SEM, we also have to multiply with the number of trials
    bias = (Ntrials - 1)^2;
  elseif strcmp(source.method, 'bootstrap')
    % Effron gives SEM estimate for the bootstrap method in algorithm 6.1 (equation 6.6)
    % to get the variance instead of SEM, we also have to multiply with the number of trials
    bias = Ntrials;
  elseif strcmp(source.method, 'pseudovalue')
    % note that I have not put any thought in this aspect yet
    ft_warning('don''t know how to compute bias for pseudovalue resampling');
    bias = 1;
  end

  % compute the variance over all repetitions
  if isfield(sumdip, 'var'),    dipvar.var    = bias*(sqrdip.var    - (sumdip.var   .^2)/length(dip))/(length(dip)-1); end
  if isfield(sumdip, 'pow'),    dipvar.pow    = bias*(sqrdip.pow    - (sumdip.pow   .^2)/length(dip))/(length(dip)-1); end
  if isfield(sumdip, 'coh'),    dipvar.coh    = bias*(sqrdip.coh    - (sumdip.coh   .^2)/length(dip))/(length(dip)-1); end
  if isfield(sumdip, 'rv' ),    dipvar.rv     = bias*(sqrdip.rv     - (sumdip.rv    .^2)/length(dip))/(length(dip)-1); end
  if isfield(sumdip, 'noise' ), dipvar.noise  = bias*(sqrdip.noise  - (sumdip.noise .^2)/length(dip))/(length(dip)-1); end
  if isfield(sumdip, 'nai' ),   dipvar.nai    = bias*(sqrdip.nai    - (sumdip.nai   .^2)/length(dip))/(length(dip)-1); end

  % compute the SEM over all repetitions
  if isfield(sumdip, 'var'),    dipsem.var    = (dipvar.var   /Ntrials).^0.5; end
  if isfield(sumdip, 'pow'),    dipsem.pow    = (dipvar.pow   /Ntrials).^0.5; end
  if isfield(sumdip, 'coh'),    dipsem.coh    = (dipvar.coh   /Ntrials).^0.5; end
  if isfield(sumdip, 'rv' ),    dipsem.rv     = (dipvar.rv    /Ntrials).^0.5; end
  if isfield(sumdip, 'noise' ), dipsem.noise  = (dipvar.noise /Ntrials).^0.5; end
  if isfield(sumdip, 'nai' ),   dipsem.nai    = (dipvar.nai   /Ntrials).^0.5; end

  % compute the mean and SEM over all repetitions for the cell-array with mom
  if isfield(dip(trial), 'mom')
    for i=1:length(dip(1).mom)
      dipmean.mom{i} = sumdip.mom{i}/length(dip);
      dipvar.mom{i} = bias*(sqrdip.mom{i} - (sumdip.mom{i}.^2)/length(dip))/(length(dip)-1);
      dipsem.mom{i} = (dipvar.mom{i}/Ntrials).^0.5;
    end
  end

  % compute the mean and SEM over all repetitions for the cell-array with csd
  if isfield(dip(trial), 'csd')
    for i=1:length(dip(1).csd)
      dipmean.csd{i} = sumdip.csd{i}/length(dip);
      dipvar.csd{i} = bias*(sqrdip.csd{i} - (sumdip.csd{i}.^2)/length(dip))/(length(dip)-1);
      dipsem.csd{i} = (dipvar.csd{i}/Ntrials).^0.5;
    end
  end

  if strcmp(source.method, 'pseudovalue')
    % keep the trials, since they have been converted to pseudovalues
    % and hence the trials contain the interesting data
  elseif keeptrials
    % keep the trials upon request
  else
    % remove the original trials
    source = rmfield(source, 'trial');
    % assign the descriptive statistics to the output source structure
    source.avg = dipmean;
    source.var = dipvar;
    source.sem = dipsem;
  end
end

if strcmp(cfg.resolutionmatrix, 'yes')
  % this is only implemented for pcc and no refdips/chans at the moment
  Nchan        = size(source.leadfield{insideindx(1)}, 1);
  Ninside      = length(insideindx);
  allfilter    = zeros(Ninside,Nchan);
  allleadfield = zeros(Nchan,Ninside);
  dipsel       = match_str(source.avg.csdlabel, 'scandip');
  ft_progress('init', cfg.feedback, 'computing resolution matrix');
  for diplop=1:length(insideindx)
    ft_progress(diplop/length(insideindx), 'computing resolution matrix %d/%d\n', diplop, length(insideindx));
    % concatenate all filters
    allfilter(diplop,:)    = source.avg.filter{insideindx(diplop)}(dipsel,:);
    % concatenate all leadfields
    allleadfield(:,diplop) = source.leadfield{insideindx(diplop)};
  end
  ft_progress('close');
  % multiply the filters and leadfields to obtain the resolution matrix
  % see equation 1 and 2 in De Peralta-Menendez RG, Gonzalez-Andino SL: A critical analysis of linear inverse solutions to the neuroelectromagnetic inverse problem. IEEE Transactions on Biomedical Engineering 45: 440-448, 1998.
  source.resolution = nan(Ndipole, Ndipole);
  source.resolution(insideindx, insideindx) = allfilter*allleadfield;
end

% compute fwhm
if strcmp(cfg.fwhm, 'yes')
  switch cfg.fwhmmethod
    case 'barnes'
      if ~isfield(source, 'dim')
        ft_error('computation of fwhm is not possible with method ''barnes'' is not possible when the dipoles are not defined on a regular 3D grid');
      end
      fprintf('computing fwhm of spatial filters using method ''barnes''\n');
      source = estimate_fwhm1(source, cfg.fwhmremovecenter);
    case 'gaussfit'
      fprintf('computing fwhm of spatial filters using method ''gaussfit''\n');
      source = estimate_fwhm2(source, cfg.fwhmmaxdist);
    otherwise
      ft_error('unknown method for fwhm estimation');
  end
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   source
ft_postamble provenance source
ft_postamble history    source
ft_postamble savevar    source

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute eta from a csd-matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [eta, u] = csd2eta(csd)
[u,s,v] = svd(real(csd));
eta     = s(2,2)./s(1,1);
u       = u'; %orientation is defined in the rows

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute fa from a csd-matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fa] = csd2fa(csd)
s  = svd(real(csd));
ns = rank(real(csd));
s  = s(1:ns);
ms = mean(s);
fa = sqrt( (ns./(ns-1)) .* (sum((s-ms).^2))./(sum(s.^2)) );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute power
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = powmethod_lambda1(x, ind)

if nargin==1
  ind = 1:size(x,1);
end
s = svd(x(ind,ind));
p = s(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute power
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = powmethod_trace(x, ind)

if nargin==1
  ind = 1:size(x,1);
end
p = trace(x(ind,ind));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute power
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = powmethod_regular(x, ind)

if nargin==1
  ind = 1:size(x,1);
end
p = abs(x(ind,ind));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to obtain the largest singular value or trace of the
% source CSD matrices resulting from DICS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = lambda1(x)
s = svd(x);
s = s(1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x, c, funVal, ValueL]=fusedLogisticR(A, y, lambda, opts)
%
%%
% Function fusedLogisticR
%      Logistic Loss with the Fused Lasso Penalty
%
%% Problem
%
%  min  f(x,c) = - weight_i * log (p_i) + lambda * ||x||_1 +   lambda_2 ||Rx||_1
%                               + 1/2 rsL2 * ||x||_2^2
%
%  By default, rsL2=0.
%  When rsL2 is nonzero, this correspons the well-know elastic net.
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  y_i (either 1 or -1) is the response
%
%  p_i= 1/ (1+ exp(-y_i (x' * a_i + c) ) ) denotes the probability
%
%  weight_i denotes the weight for the i-th sample
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  lambda -        Lq/L1 norm regularization parameter (lambda >=0)
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         The obtained weight of size n x 1
%  c-         The obtained intercept (scalar)
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on February 18, 2010.
%
%% Related papers
%
% [1]  Jun Liu, Lei Yuan, and Jieping Ye, An Efficient Algorithm for
%      a Class of Fused Lasso Problems, KDD, 2010.
%
%% Related functions
%
%  sll_opts, initFactor, pathSolutionLeast
%  LogisticC, nnLogisticR, nnLogisticC
%
%%

%% Verify and initialize the parameters
%%
if (nargin <3)
    error('\n Inputs: A, y, and lambda should be specified!\n');
end

[m,n]=size(A);

if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

if (lambda<=0)
    error('\n lambda should be positive!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & others

% The parameter 'weight' contains the weight for each training sample.
% See the definition of the problem above.
% The summation of the weights for all the samples equals to 1.
if (isfield(opts,'sWeight'))
    sWeight=opts.sWeight;

    if ( length(sWeight)~=2 || sWeight(1) <=0 || sWeight(2) <= 0)
        error('\n Check opts.sWeight, which contains two positive values');
    end

    % we process the weight, so that the summation of the weights for all
    % the samples is 1.

    p_flag=(y==1);                  % the indices of the postive samples
    m1=sum(p_flag) * sWeight(1);    % the total weight for the positive samples
    m2=sum(~p_flag) * sWeight(2);   % the total weight for the positive samples

    weight(p_flag,1)=sWeight(1)/(m1+m2);
    weight(~p_flag,1)=sWeight(2)/(m1+m2);
else
    weight=ones(m,1)/m;             % if not specified, we apply equal weight
end

%% Starting point initialization

% process the regularization parameter

% L2 norm regularization
if isfield(opts,'rsL2')
    rsL2=opts.rsL2;
    if (rsL2<0)
        error('\n opts.rsL2 should be nonnegative!');
    end
else
    rsL2=0;
end

p_flag=(y==1);                  % the indices of the postive samples
m1=sum(weight(p_flag));         % the total weight for the positive samples
m2=1-m1;                        % the total weight for the positive samples

% L1 norm regularization
if (opts.rFlag==0)
    lambda=lambda;
else % lambda here is the scaling factor lying in [0,1]
    if (lambda<0 || lambda>1)
        error('\n opts.rFlag=1, and lambda should be in [0,1]');
    end

    % we compute ATb for computing lambda_max, when the input lambda is a ratio

    b(p_flag,1)=m2;  b(~p_flag,1)=-m1;
    b=b.*weight;

    % compute AT b
    if (opts.nFlag==0)
        ATb =A'*b;
    elseif (opts.nFlag==1)
        ATb= A'*b - sum(b) * mu';  ATb=ATb./nu;
    else
        invNu=b./nu;               ATb=A'*invNu-sum(invNu)*mu';
    end

    lambda_max=max(abs(ATb));
    lambda=lambda*lambda_max;

    if ~isfield(opts,'fusedPenalty')
        lambda2=0;
    else
        lambda2=lambda_max * opts.fusedPenalty;
    end

    rsL2=rsL2*lambda_max; % the input rsL2 is a ratio of lambda_max
end

% initialize a starting point
if opts.init==2
    x=zeros(n,1); c=log(m1/m2);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        x=zeros(n,1);
    end

    if isfield(opts,'c0')
        c=opts.c0;
    else
        c=log(m1/m2);
    end
end

% compute A x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./nu; mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);    Ax=Ax./nu;
end

% restart the program for better efficiency
%  this is a newly added function
if (~isfield(opts,'rStartNum'))
    opts.rStartNum=opts.maxIter;
else
    if (opts.rStartNum<=0)
        opts.rStartNum=opts.maxIter;
    end
end

%% The main program

z0=zeros(n-1,1);
% z0 is the starting point for flsa

%% The Armijo Goldstein line search scheme + accelearted gradient descent

if (opts.mFlag==0 && opts.lFlag==0)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m + rsL2; % the intial guess of the Lipschitz continuous gradient

    weighty=weight.*y;
    % the product between weight and y

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);
    cp=c;         ccp=0;

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;  sc=c + beta* ccp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute As=A*s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        aa=- y.*(As+ sc);

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
            rsL2/2 * s'*s;

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y.* weight) * (1 - prob)
        b= -weighty.*(1-prob);

        gc=sum(b); % the gradient of c

        % compute g= AT b, the gradient of x
        if (opts.nFlag==0)
            g=A'*b;
        elseif (opts.nFlag==1)
            g=A'*b - sum(b) * mu';  g=g./nu;
        else
            invNu=b./nu;              g=A'*invNu-sum(invNu)*mu';
        end

        g=g+ rsL2 * s;  % add the squared L2 norm regularization term

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the Lq/L1-norm regularized projection
            v=s-g/L;  c= sc- gc/L;

            [x, z, infor]=flsa(v, z0,...
                lambda / L, lambda2 / L, n,...
                1000, 1e-8, 1, 6);
            z0=z;

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Ax+ c);

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb ) +...
                rsL2/2 * x'*x;

            r_sum= (v'*v + (c-sc)^2) / 2;
            l_sum=fun_x - fun_s - v'* g - (c-sc)* gc;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + v'* g + c * gc
            %                           + L/2 * (v'*v + (c-sc)^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%e, r_sum=%e',L, r_sum);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;
        % store values for L

        xxp=x-xp;    ccp=c-cp;
        funVal(iterStep)=fun_x + lambda* sum(abs(x)) +...
            lambda2 * sum( abs( x(2:n) -x(1:(n-1)) ));

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end

        % restart the program every opts.rStartNum
        if (~mod(iterStep, opts.rStartNum))
            alphap=0; alpha=1;
            xp=x; Axp=Ax; xxp=zeros(n,1); L =1;
        end
    end
else
    error('\n This function only supports opts.mFlag=0 & opts.lFlag=0');
end

author: thomaskuestner


text:
function [rx ry rz]= GetEulerAngles(R)

% This function return the rotation along x,y and z direction from a
% Rotation Matrix

%Inputs:
    % R= 3x3 Rotation Matrix
%Outputs:
    % rx= Rotation along x direction in radians
    % ry= Rotation along y direction in radians
    % rz= Rotation along z direction in radians

%     R =
%
% [                           cos(ry)*cos(rz),                          -cos(ry)*sin(rz),          sin(ry)]
% [ cos(rx)*sin(rz) + cos(rz)*sin(rx)*sin(ry), cos(rx)*cos(rz) - sin(rx)*sin(ry)*sin(rz), -cos(ry)*sin(rx)]
% [ sin(rx)*sin(rz) - cos(rx)*cos(rz)*sin(ry), cos(rz)*sin(rx) + cos(rx)*sin(ry)*sin(rz),  cos(rx)*cos(ry)]

% Author : Sandeep Sasidharan
% http://sandeepsasidharan.webs.com

ry=asin(R(1,3));
rz=acos(R(1,1)/cos(ry));
rx=acos(R(3,3)/cos(ry));

author: Sable

repo: mcbench-benchmarks