text:
function [U,S,output] = lmlra(T,size_core,S0,options)
%LMLRA Low multilinear rank approximation.
%   [U,S,output] = lmlra(T,size_core) computes the factor matrices U{1},
%   ..., U{N} and core tensor S of dimensions size_core belonging to a low
%   multilinear rank approximation of the N-th order tensor T.
%
%   [U,S,output] = lmlra(T,U0,S0) allows the user to provide initial factor
%   matrices U0 and core tensor S0, to be used by the main algorithm.
%   Factor matrices equal to the empty matrix [] will be initialized using
%   the chosen initialization method.
%
%   The structure output contains the output of the selected algorithms:
%
%      output.Initialization - The output of the initialization step.
%      output.Algorithm      - The output of the main algorithm.
%
%   lmlra(T,size_core,options) and lmlra(T,U0,S0,options) allow the user to
%   choose which algorithms will be used in the different steps and also
%   set parameters for those algorithms:
%
%      Use options.Initialization = [{@lmlra_aca}|@lmlra_rnd|@mlsvd] to
%      choose the initialization method. By default, an adaptive cross
%      approximation algorithm is used. The variable
%      options.InitializationOptions will be passed to the chosen
%      initialization method as third argument.
%
%      Use options.Algorithm = [@lmlra_hooi|@lmlra_minf|{@lmlra_nls}| ...
%      lmlra3_dgn|lmlra3_rtr] to choose the main algorithm. The structure
%      options.AlgorithmOptions will be passed to the chosen algorithm.
%
%   Further options are:
%
%      options.Display = false - Set to true to enable printing output
%                                information to the command line.
%
%   See also mlrankest, lmlragen, lmlraerr.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)

% Gather some input information.
T = fmt(T);
N = length(size_core);
if nargin < 3, S0 = []; end
if iscell(size_core) || (nargin >= 3 && isnumeric(S0) && ~isempty(S0))
    U = size_core;
    S = S0;
    size_core = cellfun('size',U,2);
end

% Check the options structure.
isfunc = @(f)isa(f,'function_handle');
xsfunc = @(f)isfunc(f)&&exist(func2str(f),'file');
if nargin < 4
    if nargin == 3 && isstruct(S0), options = S0;
    else options = struct; end
end
if ~isfield(options,'Initialization')
    funcs = {@lmlra_aca,@lmlra_rnd,@mlsvd};
    options.Initialization = funcs{find(cellfun(xsfunc,funcs),1)};
end
if length(size_core) == 1, options.Initialization = @mlsvd; end
if ~isfield(options,'Algorithm')
    funcs = {@lmlra_nls,@lmlra_minf,@lmlra_hooi,@lmlra3_rtr,@lmlra3_dgn};
    options.Algorithm = funcs{find(cellfun(xsfunc,funcs),1)};
end
if ~isfield(options,'Display'), options.Display = false; end
if ~options.Display, print = @(varargin)true; else print = @fprintf; end

% Step 1: initialize the factor matrices unless they were all provided by
% the user.
print('Step 1: Initialization ');
if ~exist('U','var'), U = cell(1,N); end
Uempty = cellfun(@isempty,U);
if any(Uempty) || isempty(S0)
    if isfunc(options.Initialization)
        print('is %s... ',func2str(options.Initialization));
    end
    if ~xsfunc(options.Initialization)
        error('lmlra:Initialization','Not a valid initialization.');
    end
else
    options.Initialization = false;
    print('is manual... ');
end
if isfunc(options.Initialization)
    % Generate initial factor matrices.
    if ~isfield(options,'InitializationOptions')
        options.InitializationOptions = struct;
    end
    [U0,S0,sv] = options.Initialization(T,size_core,...
                 options.InitializationOptions);
    if isstruct(sv)
        output.Initialization = sv;
    else
        output.Initialization.sv = sv;
    end
    % Fill empty factor matrices.
    for n = 1:sum(Uempty), U{n} = U0{n}; end
    if ~isempty(S0), S = S0; end
    output.Initialization.Name = func2str(options.Initialization);
else
    output.Initialization.Name = 'manual';
end
output.Initialization.relerr = frob(lmlrares(T,U,S))/frob(T);
print('relative error = %.6g.\n',output.Initialization.relerr);

% Step 2: run the main LMLRA algorithm.
if xsfunc(options.Algorithm)
    print('Step 2: Algorithm is %s... ',func2str(options.Algorithm));
    if ~isfield(options,'AlgorithmOptions')
        options.AlgorithmOptions = struct;
    end
    [U,S,output.Algorithm] = options.Algorithm(T,U,S,...
                             options.AlgorithmOptions);
    output.Algorithm.Name = func2str(options.Algorithm);
else
    error('lmlra:Algorithm','Not a valid algorithm.');
end
if isfield(output.Algorithm,'iterations')
    print('iterations = %i, ',output.Algorithm.iterations);
end
output.Algorithm.relerr = frob(lmlrares(T,U,S))/frob(T);
print('relative error = %.6g.\n',output.Algorithm.relerr);

% Format output.
fn = fieldnames(output);
for f = 1:length(fn)
    output.(fn{f}) = orderfields(output.(fn{f}));
end

author: andrewssobral

repo: lrslibrary


text:
function [inH2O] = kPa2inH2O(kPa)
% Convert pressure from kilopascals to inches of water.
% Chad Greene 2012
inH2O = kPa*4.01463;

author: Sable

repo: mcbench-benchmarks


text:
function [source] = ft_sourceanalysis(cfg, data, baseline)

% FT_SOURCEANALYSIS performs beamformer dipole analysis on EEG or MEG data
% after preprocessing and a timelocked or frequency analysis
%
% Use as
%   [source] = ft_sourceanalysis(cfg, freq)
% or
%   [source] = ft_sourceanalysis(cfg, timelock)
%
% where the second input argument with the data should be organised in a structure
% as obtained from the FT_FREQANALYSIS or FT_TIMELOCKANALYSIS function. The
% configuration "cfg" is a structure containing information about source positions
% and other options.
%
% The different source reconstruction algorithms that are implemented are
%   cfg.method     = 'lcmv'    linear constrained minimum variance beamformer
%                    'sam'     synthetic aperture magnetometry
%                    'dics'    dynamic imaging of coherent sources
%                    'pcc'     partial cannonical correlation/coherence
%                    'mne'     minimum norm estimation
%                    'rv'      scan residual variance with single dipole
%                    'music'   multiple signal classification
%                    'sloreta' standardized low-resolution electromagnetic tomography
%                    'eloreta' exact low-resolution electromagnetic tomography
% The DICS and PCC methods are for frequency or time-frequency domain data, all other
% methods are for time domain data. ELORETA can be used both for time, frequency and
% time-frequency domain data.
%
% The complete grid with dipole positions and optionally precomputed leadfields is
% constructed using FT_PREPARE_SOURCEMODEL. It can be specified as as a regular 3-D
% grid that is aligned with the axes of the head coordinate system using
%   cfg.xgrid               = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.ygrid               = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.zgrid               = vector (e.g.   0:1:20) or 'auto' (default = 'auto')
%   cfg.resolution          = number (e.g. 1 cm) for automatic grid generation
% If the source model destribes a triangulated cortical sheet, it is described as
%   cfg.sourcemodel.pos     = N*3 matrix with the vertex positions of the cortical sheet
%   cfg.sourcemodel.tri     = M*3 matrix that describes the triangles connecting the vertices
% Alternatively the position of a few dipoles at locations of interest can be
% user-specified, for example obtained from an anatomical or functional MRI
%   cfg.sourcemodel.pos     = N*3 matrix with position of each source
%   cfg.sourcemodel.inside  = N*1 vector with boolean value whether grid point is inside brain (optional)
%   cfg.sourcemodel.dim     = [Nx Ny Nz] vector with dimensions in case of 3-D grid (optional)
%
% Besides the source positions, you may also include previously computed
% spatial filters and/or leadfields using
%   cfg.sourcemodel.filter
%   cfg.sourcemodel.leadfield
%
% The following strategies are supported to obtain statistics for the source parameters using
% multiple trials in the data, either directly or through a resampling-based approach
%   cfg.rawtrial      = 'no' or 'yes'   construct filter from single trials, apply to single trials. Note that you also may want to set cfg.keeptrials='yes' to keep all trial information, especially if using in combination with sourcemodel.filter
%   cfg.jackknife     = 'no' or 'yes'   jackknife resampling of trials
%   cfg.pseudovalue   = 'no' or 'yes'   pseudovalue resampling of trials
%   cfg.bootstrap     = 'no' or 'yes'   bootstrap resampling of trials
%   cfg.numbootstrap  = number of bootstrap replications (e.g. number of original trials)
% If none of these options is specified, the average over the trials will
% be computed prior to computing the source reconstruction.
%
% To obtain statistics over the source parameters between two conditions, you
% can also use a resampling procedure that reshuffles the trials over both
% conditions. In that case, you should call the function with two datasets
% containing single trial data like
%   [source] = ft_sourceanalysis(cfg, freqA, freqB)
%   [source] = ft_sourceanalysis(cfg, timelockA, timelockB)
% and you should specify
%   cfg.randomization      = 'no' or 'yes'
%   cfg.permutation        = 'no' or 'yes'
%   cfg.numrandomization   = number, e.g. 500
%   cfg.numpermutation     = number, e.g. 500 or 'all'
%
% If you have not specified a sourcemodel with pre-computed leadfields,
% the leadfield for each source position will be computed on the fly.
% In that case you can modify the leadfields by reducing the rank
% (i.e.  remove the weakest orientation), or by normalizing each
% column.
%   cfg.reducerank  = 'no', or number (default = 3 for EEG, 2 for MEG)
%   cfg.normalize   = 'no' or 'yes' (default = 'no')
%
% Other configuration options are
%   cfg.channel       = Nx1 cell-array with selection of channels (default = 'all'),
%                       see FT_CHANNELSELECTION for details
%   cfg.frequency     = single number (in Hz)
%   cfg.latency       = single number in seconds, for time-frequency analysis
%   cfg.lambda        = number or empty for automatic default
%   cfg.kappa         = number or empty for automatic default
%   cfg.tol           = number or empty for automatic default
%   cfg.refchan       = reference channel label (for coherence)
%   cfg.refdip        = reference dipole location (for coherence)
%   cfg.supchan       = suppressed channel label(s)
%   cfg.supdip        = suppressed dipole location(s)
%   cfg.keeptrials    = 'no' or 'yes'
%   cfg.keepleadfield = 'no' or 'yes'
%   cfg.projectnoise  = 'no' or 'yes'
%   cfg.keepfilter    = 'no' or 'yes'
%   cfg.keepcsd       = 'no' or 'yes'
%   cfg.keepmom       = 'no' or 'yes'
%   cfg.feedback      = 'no', 'text', 'textbar', 'gui' (default = 'text')
%
% The volume conduction model of the head should be specified as
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% The EEG or MEG sensor positions can be present in the data or can be specified as
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad          = structure with gradiometer definition or filename, see FT_READ_SENS
%
% 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_SOURCEDESCRIPTIVES, FT_SOURCESTATISTICS, FT_PREPARE_LEADFIELD,
% FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

% Undocumented local options:
% cfg.numcomponents
% cfg.refchannel
% cfg.trialweight   = 'equal' or 'proportional'
% cfg.powmethod     = 'lambda1' or 'trace'

% Copyright (c) 2003-2008, F.C. Donders Centre, Robert Oostenveld
%
% 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 data baseline
ft_preamble provenance data baseline
ft_preamble trackconfig

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

% the baseline data can be passed as input argument or can be read from disk
hasbaseline = exist('baseline', 'var');

% check if the input data is valid for this function
data = ft_checkdata(data, 'datatype', {'timelock', 'freq', 'comp'}, 'feedback', 'yes');

if hasbaseline
  baseline = ft_checkdata(baseline, 'datatype', {'timelock', 'freq', 'comp'}, 'feedback', 'yes');
end

% check that the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',    {'toilim', 'latency'});
cfg = ft_checkconfig(cfg, 'renamed',    {'foilim', 'frequency'});
cfg = ft_checkconfig(cfg, 'renamed',    {'jacknife', 'jackknife'});
cfg = ft_checkconfig(cfg, 'renamed',    {'refchannel', 'refchan'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'power',           'dics'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'coh_refchan',     'dics'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'coh_refdip',      'dics'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'dics_cohrefchan', 'dics'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'dics_cohrefdip',  'dics'});
cfg = ft_checkconfig(cfg, 'forbidden',  {'parallel', 'trials'});
cfg = ft_checkconfig(cfg, 'forbidden',  {'foi', 'toi'});
cfg = ft_checkconfig(cfg, 'renamed',    {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'vol',     'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'grid',    'sourcemodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed',    {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed',    {'optofile', 'opto'});

% determine the type of input data
isfreq     = ft_datatype(data, 'freq');
iscomp     = ft_datatype(data, 'comp');
istimelock = ft_datatype(data, 'timelock');
if all(~[isfreq iscomp istimelock])
  ft_error('input data is not recognized');
end

% set the defaults
cfg.keeptrials       = ft_getopt(cfg, 'keeptrials', 'no');
cfg.keepleadfield    = ft_getopt(cfg, 'keepleadfield', 'no');
cfg.trialweight      = ft_getopt(cfg, 'trialweight', 'equal');
cfg.jackknife        = ft_getopt(cfg, 'jackknife',   'no');
cfg.pseudovalue      = ft_getopt(cfg, 'pseudovalue', 'no');
cfg.bootstrap        = ft_getopt(cfg, 'bootstrap',   'no');
cfg.singletrial      = ft_getopt(cfg, 'singletrial', 'no');
cfg.rawtrial         = ft_getopt(cfg, 'rawtrial',    'no');
cfg.randomization    = ft_getopt(cfg, 'randomization', 'no');
cfg.numrandomization = ft_getopt(cfg, 'numrandomization', 100);
cfg.permutation      = ft_getopt(cfg, 'permutation',      'no');
cfg.numpermutation   = ft_getopt(cfg, 'numpermutation',   100);
cfg.wakewulf         = ft_getopt(cfg, 'wakewulf', 'yes');
cfg.killwulf         = ft_getopt(cfg, 'killwulf', 'yes');
cfg.channel          = ft_getopt(cfg, 'channel',  'all');
cfg.supdip           = ft_getopt(cfg, 'supdip',        []);
cfg.latency          = ft_getopt(cfg, 'latency',   'all');
cfg.frequency        = ft_getopt(cfg, 'frequency', 'all');

if istimelock
  cfg.method = ft_getopt(cfg, 'method', 'lcmv');
elseif isfreq
  cfg.method = ft_getopt(cfg, 'method', 'dics');
else
  cfg.method = ft_getopt(cfg, 'method', []);
end

% put the low-level options pertaining to the source reconstruction method in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg',  cfg.method);

% put the low-level options pertaining to the dipole grid in their own field
cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid', 'tight'});  % this is moved to cfg.sourcemodel.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'}); % this is moved to cfg.sourcemodel.unit  by the subsequent createsubcfg

% put the low-level options pertaining to the sourcemodel in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg', {'sourcemodel'});
% move some fields from cfg.sourcemodel back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', {'sourcemodel'});

cfg.(cfg.method).keepfilter    = ft_getopt(cfg.(cfg.method), 'keepfilter',    'no');
cfg.(cfg.method).keepcsd       = ft_getopt(cfg.(cfg.method), 'keepcsd',       'no');
cfg.(cfg.method).keepmom       = ft_getopt(cfg.(cfg.method), 'keepmom',       'yes');
cfg.(cfg.method).projectnoise  = ft_getopt(cfg.(cfg.method), 'projectnoise',  'no');
cfg.(cfg.method).feedback      = ft_getopt(cfg.(cfg.method), 'feedback',      'text');
cfg.(cfg.method).lambda        = ft_getopt(cfg.(cfg.method), 'lambda',        []);
cfg.(cfg.method).kappa         = ft_getopt(cfg.(cfg.method), 'kappa',         []);
cfg.(cfg.method).tol           = ft_getopt(cfg.(cfg.method), 'tol',           []);
cfg.(cfg.method).invmethod     = ft_getopt(cfg.(cfg.method), 'invmethod',     []);
cfg.(cfg.method).powmethod     = ft_getopt(cfg.(cfg.method), 'powmethod',     []);
cfg.(cfg.method).normalize     = ft_getopt(cfg.(cfg.method), 'normalize',     'no');
cfg.(cfg.method).reducerank    = ft_getopt(cfg.(cfg.method), 'reducerank',    []); % the default for this is handled below

if hasbaseline && (strcmp(cfg.randomization, 'no') && strcmp(cfg.permutation, 'no'))
  ft_error('input of two conditions only makes sense if you want to randomize or permute');
elseif ~hasbaseline && (strcmp(cfg.randomization, 'yes') || strcmp(cfg.permutation, 'yes'))
  ft_error('randomization or permutation requires that you give two conditions as input');
end

if isfield(cfg, 'latency') && ischar(cfg.latency) && strcmp(cfg.latency, 'all') && istimelock
  %error('specification of cfg.latency is only required for time-frequency data');
end

if sum([strcmp(cfg.jackknife, 'yes'), strcmp(cfg.bootstrap, 'yes'), strcmp(cfg.pseudovalue, 'yes'), strcmp(cfg.singletrial, 'yes'), strcmp(cfg.rawtrial, 'yes'), strcmp(cfg.randomization, 'yes'), strcmp(cfg.permutation, 'yes')])>1
  ft_error('jackknife, bootstrap, pseudovalue, singletrial, rawtrial, randomization and permutation are mutually exclusive');
end

if strcmp(cfg.rawtrial, 'yes') && isfield(cfg, 'sourcemodel') && ~isfield(cfg.sourcemodel, 'filter')
  ft_warning('Using each trial to compute its own filter is not currently recommended. Use this option only with precomputed filters in cfg.sourcemodel.filter');
end

% start with an empty output structure
source = [];

if istimelock
  % add the time axis to the output
  tmpcfg = keepfields(cfg, {'channel', 'latency', 'showcallinfo'});
  tmpcfg.avgovertime = 'no';
  data = ft_selectdata(tmpcfg, data);
  % restore the provenance information
  [cfg, data] = rollback_provenance(cfg, data);

  % copy the descriptive fields to the output
  source = copyfields(data, source, {'time'});

elseif isfreq
  tmpcfg = keepfields(cfg, {'channel', 'latency', 'frequency', 'refchan', 'nanmean', 'showcallinfo'});

  % ensure that the refchan is kept, if present
  if isfield(tmpcfg, 'refchan') && ~isempty(tmpcfg.refchan) && isempty(match_str(tmpcfg.channel, tmpcfg.refchan))
    hasrefchan = 1;
  else
    hasrefchan = 0;
  end

  if hasrefchan
    if ischar(tmpcfg.refchan), tmpcfg.refchan = {tmpcfg.refchan}; end
    tmpchannel     = ft_channelselection(tmpcfg.channel, data.label); % the channels needed for the spatial filter
    tmpcfg.channel = cat(1, tmpchannel, tmpcfg.refchan);
    tmpcfg         = rmfield(tmpcfg, 'refchan');
  end

  tmpcfg.avgoverfreq = 'yes';
  if isfield(data, 'time')
    tmpcfg.avgovertime = 'yes';
  end
  data = ft_selectdata(tmpcfg, data);
  % restore the provenance information
  [cfg, data] = rollback_provenance(cfg, data);

  if hasrefchan, cfg.channel = match_str(data.label, tmpchannel); end

  % copy the descriptive fields to the output
  source = copyfields(data, source, {'time', 'freq', 'cumtapcnt'});

  % HACK the remainder of the code expects a single number
  cfg.frequency = mean(cfg.frequency);
  if isfield(data, 'time')
    cfg.latency   = mean(cfg.latency);
  end

elseif iscomp
  % FIXME, select the components here
  % FIXME, add the component numbers to the output
  ft_error('the use of component data in ft_sourceanalysis is disabled for the time being: if you encounter this error message and you need this functionality please contact the FieldTrip development team');
end

convertcomp = false;
if iscomp && (strcmp(cfg.method, 'rv') || strcmp(cfg.method, 'music'))
  % these timelock methods are also supported for frequency or component data
  if iscomp
    convertcomp = true;
    % the conversion will be done below, after the latency and channel selection
  end
elseif isfreq && isfield(data, 'labelcmb')
  % ensure that the cross-spectral densities are chan_chan_therest,
  % otherwise the latency and frequency selection could fail, so we don't
  % need to worry about linearly indexed cross-spectral densities below
  % this point, this step may take some time, if multiple trials are
  % present in the data
  fprintf('converting the linearly indexed channelcombinations into a square CSD-matrix\n');
  data = ft_checkdata(data, 'cmbrepresentation', 'full');
end

if isfreq
  % as per the call to ft_checkdata above, the dimord of the freq-data is
  % either (rpt_)chan_chan_otherstuff, or rpttap_chan_otherstuff. The
  % former is with cross-spectra, the latter is with fourierspctrm

  % previously there was some explicit dimord checking here, but I think
  % with the more consistent data handling it is not necessary anymore.
end

% collect and preprocess the electrodes/gradiometer and head model
[headmodel, sens, cfg] = prepare_headmodel(cfg, data);

% set the default for reducing the rank of the leadfields
if isempty(cfg.(cfg.method).reducerank)
  if ft_senstype(sens, 'eeg')
    cfg.(cfg.method).reducerank = 'no';    % for EEG
  elseif ft_senstype(sens, 'meg') && ft_headmodeltype(headmodel, 'infinite')
    cfg.(cfg.method).reducerank = 'no';    % for MEG with a magnetic dipole, e.g. a HPI coil
  elseif ft_senstype(sens, 'meg')
    cfg.(cfg.method).reducerank = 'yes';   % for MEG with a current dipole in a volume conductor
  end
end

% It might be that the number of channels in the data, the number of
% channels in the electrode/gradiometer definition and the number of
% channels in the localspheres volume conduction model are different.
% Hence a subset of the data channels will be used.
Nchans = length(cfg.channel);

if strcmp(cfg.keepleadfield, 'yes') && (~isfield(cfg, 'sourcemodel') || ~isfield(cfg.sourcemodel, 'leadfield'))
  % precompute the leadfields upon the users request
  fprintf('precomputing leadfields\n');
  sourcemodel = ft_prepare_leadfield(cfg, data);
elseif (strcmp(cfg.permutation,   'yes') || ...
    strcmp(cfg.randomization, 'yes') || ...
    strcmp(cfg.bootstrap,     'yes') || ...
    strcmp(cfg.jackknife,     'yes') || ...
    strcmp(cfg.pseudovalue,   'yes') || ...
    strcmp(cfg.singletrial,   'yes') || ...
    strcmp(cfg.rawtrial,      'yes')) && (~isfield(cfg, 'sourcemodel') || ~isfield(cfg.sourcemodel, 'leadfield'))
  % also precompute the leadfields if multiple trials have to be processed
  fprintf('precomputing leadfields for efficient handling of multiple trials\n');
  sourcemodel = ft_prepare_leadfield(cfg, data);
else
  % only prepare the source positions, the leadfield will be computed on the fly if not present

  % copy all options that are potentially used in ft_prepare_sourcemodel
  tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid' 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'showcallinfo'});
  tmpcfg.headmodel = headmodel;
  if ft_senstype(sens, 'eeg')
    tmpcfg.elec = sens;
  elseif ft_senstype(sens, 'meg')
    tmpcfg.grad = sens;
  end
  % construct the dipole positions on which the source reconstruction will be done
  sourcemodel = ft_prepare_sourcemodel(tmpcfg);
  if ischar(cfg.sourcemodel)
    cfg.sourcemodel = sourcemodel;
  end
end

if isfield(cfg.sourcemodel, 'filter')
  if numel(cfg.sourcemodel.filter) == size(sourcemodel.pos, 1)
    sourcemodel.filter = cfg.sourcemodel.filter;
  else
    ft_warning('ignoring predefined filter as it does not match the number of source positions');
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do frequency domain source reconstruction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfreq && any(strcmp(cfg.method, {'dics', 'pcc', 'eloreta', 'mne','harmony', 'rv', 'music'}))

  switch cfg.method
    case 'pcc'
      % this can handle both a csd matrix and a fourier matrix, but needs
      % special handling of the refdip etc.
      cfg.refdip = ft_getopt(cfg, 'refdip', []);
      cfg.supdip = ft_getopt(cfg, 'supdip', []);
      cfg.refchan = ft_getopt(cfg, 'refchan', []);
      cfg.supchan = ft_getopt(cfg, 'supchan', []);
      cfg.refchan = ft_channelselection(cfg.refchan, data.label);
      cfg.supchan = ft_channelselection(cfg.supchan, data.label);

      % HACK: use some experimental code
      if hasbaseline
        ft_error('not supported')
      end

      tmpcfg         = cfg;
      tmpcfg.refchan = ''; % prepare_freq_matrices should not know explicitly about the refchan
      tmpcfg.channel = cfg.channel(:)';
      if isfield(cfg, 'refchan')
        % add the refchan implicitely
        tmpcfg.channel = [tmpcfg.channel cfg.refchan(:)'];
      end
      if isfield(cfg, 'supchan')
        % add the supchan implicitely
        tmpcfg.channel = [tmpcfg.channel cfg.supchan(:)'];
      end

      % select the data in the channels and the frequency of interest
      [Cf, Cr, Pr, Ntrials, tmpcfg] = prepare_freq_matrices(tmpcfg, data);

      if isfield(cfg, 'refchan') && ~isempty(cfg.refchan)
        [dum, refchanindx] = match_str(cfg.refchan, tmpcfg.channel);
      else
        refchanindx = [];
      end
      if isfield(cfg, 'supchan') && ~isempty(cfg.supchan)
        [dum, supchanindx] = match_str(cfg.supchan,tmpcfg.channel);
      else
        supchanindx = [];
      end
      Nchans = length(tmpcfg.channel); % update the number of channels

      % if the input data has a complete fourier spectrum, it can be projected through the filters
      if isfield(data, 'fourierspctrm')
        [dum, datchanindx] = match_str(tmpcfg.channel, data.label);
        fbin = nearest(data.freq, cfg.frequency);
        if strcmp(data.dimord, 'chan_freq')
          avg = data.fourierspctrm(datchanindx, fbin);
        elseif strcmp(data.dimord, 'rpt_chan_freq') || strcmp(data.dimord, 'rpttap_chan_freq')
          avg = transpose(data.fourierspctrm(:, datchanindx, fbin));
        elseif strcmp(data.dimord, 'chan_freq_time')
          tbin = nearest(data.time, cfg.latency);
          avg = data.fourierspctrm(datchanindx, fbin, tbin);
        elseif strcmp(data.dimord, 'rpt_chan_freq_time') || strcmp(data.dimord, 'rpttap_chan_freq_time')
          tbin = nearest(data.time, cfg.latency);
          avg  = transpose(data.fourierspctrm(:, datchanindx, fbin, tbin));
        end
      else
        avg = [];
      end

    case {'eloreta' 'mne' 'rv' 'music' 'harmony'}
      % these can handle both a csd matrix and a fourier matrix
      [Cf, Cr, Pr, Ntrials, cfg] = prepare_freq_matrices(cfg, data);

      % if the input data has a complete fourier spectrum, it can be projected through the filters
      if isfield(data, 'fourierspctrm')
        [dum, datchanindx] = match_str(cfg.channel, data.label);
        fbin = nearest(data.freq, cfg.frequency);
        if strcmp(data.dimord, 'chan_freq')
          avg = data.fourierspctrm(datchanindx, fbin);
        elseif strcmp(data.dimord, 'rpt_chan_freq') || strcmp(data.dimord, 'rpttap_chan_freq')
          avg = transpose(data.fourierspctrm(:, datchanindx, fbin));
        elseif strcmp(data.dimord, 'chan_freq_time')
          tbin = nearest(data.time, cfg.latency);
          avg = data.fourierspctrm(datchanindx, fbin, tbin);
        elseif strcmp(data.dimord, 'rpt_chan_freq_time') || strcmp(data.dimord, 'rpttap_chan_freq_time')
          tbin = nearest(data.time, cfg.latency);
          avg  = transpose(data.fourierspctrm(:, datchanindx, fbin, tbin));
        end
      else % The input data is a CSD matrix, this is enough for computing source power, coherence and residual power.
        avg = Cf;
      end

    case 'dics'

      [Cf, Cr, Pr, Ntrials, cfg] = prepare_freq_matrices(cfg, data);

      % assign a descriptive name to each of the dics sub-methods, the default is power only
      if isfield(cfg, 'refdip') && ~isempty(cfg.refdip);
        submethod = 'dics_refdip';
      elseif isfield(cfg, 'refchan') && ~isempty(cfg.refchan);
        submethod = 'dics_refchan';
      else
        submethod = 'dics_power';
      end

    otherwise
  end

  % This is the place to check for the consistency of the channel order in
  % the pre-computed leadfields/spatial filters, and to correct for it, if
  % necessary. This pertains to bugs 1746 and 3029.
  if isfield(sourcemodel, 'label') && (isfield(sourcemodel, 'leadfield') || isfield(sourcemodel, 'filter'))
    % match the channels in the leadfields/filters with those in the data
    [i1, i2] = match_str(cfg.channel, sourcemodel.label);
    if ~isequal(i2(:), (1:numel(sourcemodel.label))')
      if isfield(sourcemodel, 'leadfield')
        fprintf('\n\nSubselecting/reordering the channels in the precomputed leadfields\n\n');
        inside_indx = find(sourcemodel.inside);
        for k = inside_indx(:)'
          sourcemodel.leadfield{k} = sourcemodel.leadfield{k}(i2, :);
        end
      end
      if isfield(sourcemodel, 'filter')
        fprintf('\n\nSubselecting/reordering the channels in the precomputed filters\n\n');
        inside_indx = find(sourcemodel.inside);
        for k = inside_indx(:)'
          sourcemodel.filter{k} = sourcemodel.filter{k}(:, i2);
        end
      end
      sourcemodel.label = sourcemodel.label(i2);
    end
    if ~isequal(i1(:), (1:numel(cfg.channel))')
      % this is not so easy to deal with, throw an error
      ft_error('There''s a mismatch between the number/order of channels in the data, with respect to the channels in the precomputed leadfield/filter. This is not easy to solve automatically. Please look into this.');
    end
  end

  % fill these with NaNs, so that I dont have to treat them separately
  if isempty(Cr), Cr = nan(Ntrials, Nchans, 1); end
  if isempty(Pr), Pr = nan(Ntrials, 1, 1); end

  if hasbaseline
    % repeat the conversion for the baseline condition
    [bCf, bCr, bPr, Nbaseline, cfg] = prepare_freq_matrices(cfg, baseline);
    % fill these with NaNs, so that I dont have to treat them separately
    if isempty(bCr), bCr = nan(Nbaseline, Nchans, 1); end
    if isempty(bPr), bPr = nan(Nbaseline, 1, 1); end
    % rename the active condition for convenience
    aCf = Cf;
    aCr = Cr;
    aPr = Pr;
    % this is required for averaging 2 conditions using prepare_resampled_data
    cfg2 = [];
    cfg2.numcondition = 2;
    % this is required for randomizing/permuting 2 conditions using prepare_resampled_data
    cfg.numcondition = 2;
  end

  % prepare the resampling of the trials, or average the data if multiple trials are present and no resampling is necessary
  if (Ntrials<=1) && (strcmp(cfg.jackknife, 'yes') || strcmp(cfg.bootstrap, 'yes') || strcmp(cfg.pseudovalue, 'yes') || strcmp(cfg.singletrial, 'yes') || strcmp(cfg.rawtrial, 'yes') || strcmp(cfg.randomization, 'yes') || strcmp(cfg.permutation, 'yes'))
    ft_error('multiple trials required in the data\n');

  elseif strcmp(cfg.permutation, 'yes')
    % compute the cross-spectral density matrix without resampling
    [dum, avg_aCf, avg_aCr, avg_aPr, avg_bCf, avg_bCr, avg_bPr] = prepare_resampled_data(cfg2 , aCf, aCr, aPr, bCf, bCr, bPr);
    % compute the cross-spectral density matrix with random permutation
    [dum, rnd_aCf, rnd_aCr, rnd_aPr, rnd_bCf, rnd_bCr, rnd_bPr] = prepare_resampled_data(cfg, aCf, aCr, aPr, bCf, bCr, bPr);
    % concatenate the different resamplings
    Cf = cat(1, reshape(avg_aCf, [1 Nchans Nchans]), reshape(avg_bCf, [1 Nchans Nchans]), rnd_aCf, rnd_bCf);
    Cr = cat(1, reshape(avg_aCr, [1 Nchans 1     ]), reshape(avg_bCr, [1 Nchans 1     ]), rnd_aCr, rnd_bCr);
    Pr = cat(1, reshape(avg_aPr, [1 1      1     ]), reshape(avg_bPr, [1 1      1     ]), rnd_aPr, rnd_bPr);
    % clear temporary working copies
    clear avg_aCf avg_aCr avg_aPr avg_bCf avg_bCr avg_bPr
    clear rnd_aCf rnd_aCr rnd_aPr rnd_bCf rnd_bCr rnd_bPr
    % the order of the resamplings should be [avgA avgB rndA rndB rndA rndB rndA rndB ....]
    Nrepetitions = 2*cfg.numpermutation + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    Cf = Cf(order,:,:);
    Cr = Cr(order,:,:);
    Pr = Pr(order,:,:);

  elseif strcmp(cfg.randomization, 'yes')
    % compute the cross-spectral density matrix without resampling
    [dum, avg_aCf, avg_aCr, avg_aPr, avg_bCf, avg_bCr, avg_bPr] = prepare_resampled_data(cfg2 , aCf, aCr, aPr, bCf, bCr, bPr);
    % compute the cross-spectral density matrix with random resampling
    [dum, rnd_aCf, rnd_aCr, rnd_aPr, rnd_bCf, rnd_bCr, rnd_bPr] = prepare_resampled_data(cfg, aCf, aCr, aPr, bCf, bCr, bPr);
    % concatenate the different resamplings
    Cf = cat(1, reshape(avg_aCf, [1 Nchans Nchans]), reshape(avg_bCf, [1 Nchans Nchans]), rnd_aCf, rnd_bCf);
    Cr = cat(1, reshape(avg_aCr, [1 Nchans 1     ]), reshape(avg_bCr, [1 Nchans 1     ]), rnd_aCr, rnd_bCr);
    Pr = cat(1, reshape(avg_aPr, [1 1      1     ]), reshape(avg_bPr, [1 1      1     ]), rnd_aPr, rnd_bPr);
    % clear temporary working copies
    clear avg_aCf avg_aCr avg_aPr avg_bCf avg_bCr avg_bPr
    clear rnd_aCf rnd_aCr rnd_aPr rnd_bCf rnd_bCr rnd_bPr
    % the order of the resamplings should be [avgA avgB rndA rndB rndA rndB rndA rndB ....]
    Nrepetitions = 2*cfg.numrandomization + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    Cf = Cf(order,:,:);
    Cr = Cr(order,:,:);
    Pr = Pr(order,:,:);

  elseif strcmp(cfg.jackknife, 'yes')
    % compute the cross-spectral density matrix with jackknife resampling
    [cfg, Cf, Cr, Pr] = prepare_resampled_data(cfg, Cf, Cr, Pr);
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.bootstrap, 'yes')
    % compute the cross-spectral density matrix with bootstrap resampling
    [cfg, Cf, Cr, Pr] = prepare_resampled_data(cfg, Cf, Cr, Pr);
    Nrepetitions = cfg.numbootstrap;

  elseif strcmp(cfg.pseudovalue, 'yes')
    % compute the cross-spectral density matrix with pseudovalue resampling
    [cfg, Cf, Cr, Pr] = prepare_resampled_data(cfg, Cf, Cr, Pr);
    Nrepetitions = Ntrials+1;

  elseif strcmp(cfg.singletrial, 'yes')
    % The idea is that beamformer uses the average covariance to construct the
    % filter and applies it to the single trial covariance/csd. The problem
    % is that beamformer will use the averaged covariance/csd to estimate the
    % power and not the single trial covariance/csd
    ft_error('this option contains a bug, and is therefore not supported at the moment');
    Cf = Cf; % FIXME, should be averaged and repeated for each trial
    Cr = Cr; % FIXME, should be averaged and repeated for each trial
    Pr = Pr; % FIXME, should be averaged and repeated for each trial
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.rawtrial, 'yes')
    % keep all the individual trials, do not average them
    Cf = Cf;
    Cr = Cr;
    Pr = Pr;
    Nrepetitions = Ntrials;

  elseif Ntrials>1
    % compute the average from the individual trials
    Cf = reshape(sum(Cf, 1) / Ntrials, [Nchans Nchans]);
    Cr = reshape(sum(Cr, 1) / Ntrials, [Nchans 1]);
    Pr = reshape(sum(Pr, 1) / Ntrials, [1 1]);
    Nrepetitions = 1;

  elseif Ntrials==1
    % no rearrangement of trials is neccesary, the data already represents the average
    Cf = Cf;
    Cr = Cr;
    Pr = Pr;
    Nrepetitions = 1;
  end

  % reshape so that it also looks like one trial (out of many)
  if Nrepetitions==1
    Cf  = reshape(Cf , [1 Nchans Nchans]);
    Cr  = reshape(Cr , [1 Nchans 1]);
    Pr  = reshape(Pr , [1 1 1]);
  end

  % get the relevant low level options from the cfg and convert into key-value pairs
  tmpcfg = cfg.(cfg.method);
  % disable console feedback for the low-level function in case of multiple
  % repetitions
  if Nrepetitions > 1
    tmpcfg.feedback = 'none';
  end
  optarg = ft_cfg2keyval(tmpcfg);

  if Nrepetitions > 1
    ft_progress('init', cfg.(cfg.method).feedback, 'scanning repetition...');
  end

  for i=1:Nrepetitions
    size_Cf    = size(Cf);
    squeeze_Cf = reshape(Cf(i,:,:), size_Cf(2:end));

    if Nrepetitions > 1
      ft_progress(i/Nrepetitions, 'scanning repetition %d from %d', i, Nrepetitions);
    end

    switch cfg.method
      case 'dics'
        if strcmp(submethod, 'dics_power')
          dip(i) = beamformer_dics(sourcemodel, sens, headmodel, [],  squeeze_Cf, optarg{:});
        elseif strcmp(submethod, 'dics_refchan')
          dip(i) = beamformer_dics(sourcemodel, sens, headmodel, [],  squeeze_Cf, optarg{:}, 'Cr', Cr(i,:), 'Pr', Pr(i));
        elseif strcmp(submethod, 'dics_refdip')
          dip(i) = beamformer_dics(sourcemodel, sens, headmodel, [],  squeeze_Cf, optarg{:}, 'refdip', cfg.refdip);
        end
      case 'pcc'
        if ~isempty(avg) && istrue(cfg.rawtrial)
          % FIXME added by jansch because an appropriate subselection of avg
          % should be done first (i.e. select the tapers that belong to this
          % repetition
          ft_error('rawtrial in combination with pcc has been temporarily disabled');
        else
          dip(i) = beamformer_pcc(sourcemodel, sens, headmodel, avg, squeeze_Cf, optarg{:}, 'refdip', cfg.refdip, 'refchan', refchanindx, 'supdip', cfg.supdip, 'supchan', supchanindx);
        end
      case 'eloreta'
        dip(i) = ft_eloreta(sourcemodel, sens, headmodel, avg, squeeze_Cf, optarg{:});
      case 'mne'
        dip(i) = minimumnormestimate(sourcemodel, sens, headmodel, avg, optarg{:});
      case 'harmony'
        dip(i) = harmony(sourcemodel, sens, headmodel, avg, optarg{:});
        % ft_error(sprintf('method ''%s'' is unsupported for source reconstruction in the frequency domain', cfg.method));
      case {'rv'}
        dip(i) = residualvariance(sourcemodel, sens, headmodel, avg, optarg{:}) ;
      case {'music'}
        ft_error('method ''%s'' is currently unsupported for source reconstruction in the frequency domain', cfg.method);
      otherwise
    end

  end
  if Nrepetitions > 1
    ft_progress('close');
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % do time domain source reconstruction
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif istimelock && any(strcmp(cfg.method, {'lcmv', 'sam', 'mne','harmony', 'rv', 'music', 'pcc', 'mvl', 'sloreta', 'eloreta'}))

  % determine the size of the data
  Nsamples = length(data.time);
  Nchans   = length(data.label);
  if isfield(data, 'cov') && length(size(data.cov))==3
    Ntrials = size(data.cov,1);
  elseif isfield(data, 'trial') && length(size(data.trial))==3
    Ntrials = size(data.trial,1);
  else
    Ntrials = 1;
  end

  if isfield(data, 'cov')
    % use the estimated data covariance matrix
    hascovariance = 1;
  else
    % add a identity covariance matrix, this simplifies the handling of the different source reconstruction methods
    % since the covariance is only used by some reconstruction methods and might not always be present in the data
    if Ntrials==1
      data.cov = eye(Nchans);
    else
      data.cov = zeros(Ntrials,Nchans,Nchans);
      for i=1:Ntrials
        data.cov(i,:,:) = eye(Nchans);
      end
    end
    hascovariance = 0;
    ft_warning('No covariance matrix found - will assume identity covariance matrix (mininum-norm solution)');
  end

  if strcmp(cfg.method, 'pcc')
    % HACK: requires some extra defaults
    if ~isfield(cfg, 'refdip'), cfg.refdip = []; end
    if ~isfield(cfg, 'supdip'), cfg.supdip = []; end

    % HACK: experimental code
    if hasbaseline
      ft_error('not supported')
    end

    tmpcfg = [];
    tmpcfg.channel = cfg.channel(:)';
    if isfield(cfg, 'refchan')
      tmpcfg.channel = [tmpcfg.channel cfg.refchan(:)'];
    end
    if isfield(cfg, 'supchan')
      tmpcfg.channel = [tmpcfg.channel cfg.supchan(:)'];
    end

    % select the data in the channels of interest
    [dum, datchanindx] = match_str(tmpcfg.channel, data.label);
    if Ntrials==1
      data.avg = data.avg(datchanindx,:);
      data.cov = data.cov(datchanindx,datchanindx);
    else
      data.cov   = data.cov(:,datchanindx,datchanindx);
      data.trial = data.trial(:,datchanindx,:);
    end
    data.label = data.label(datchanindx);

    if isfield(cfg, 'refchan') && ~isempty(cfg.refchan)
      [dum, refchanindx] = match_str(cfg.refchan, data.label);
    else
      refchanindx = [];
    end
    if isfield(cfg, 'supchan') && ~isempty(cfg.supchan)
      [dum, supchanindx] = match_str(cfg.supchan, data.label);
    else
      supchanindx = [];
    end
    Nchans = length(tmpcfg.channel); % update the number of channels

  else
    % HACK: use the default code
    % select the channels of interest
    [dum, datchanindx] = match_str(cfg.channel, data.label);
    if strcmp(data.dimord, 'chan_time')
      % It is in principle possible to have timelockanalysis with
      % keeptrial=yes and only a single trial in the raw data.
      % In that case the covariance should be represented as Nchan*Nchan
      data.avg = data.avg(datchanindx,:);
      %data.cov = reshape(data.cov, length(datchanindx), length(datchanindx));
      data.cov = data.cov(datchanindx,datchanindx);
    elseif strcmp(data.dimord, 'rpt_chan_time')
      data.cov   = data.cov(:,datchanindx,datchanindx);
      data.trial = data.trial(:,datchanindx,:);
    else
      ft_error('unexpected dimord');
    end
    data.label = data.label(datchanindx);
    Nchans     = length(data.label);
  end

  if hasbaseline
    % baseline and active are only available together for resampling purposes,
    % hence I assume here that there are multiple trials in both
    baseline.avg   = baseline.avg(datchanindx,:);
    baseline.cov   = baseline.cov(:,datchanindx,datchanindx);
    baseline.trial = baseline.trial(:,datchanindx,:);
    % this is required for averaging 2 conditions using prepare_resampled_data
    cfg2 = [];
    cfg2.numcondition = 2;
  end

  % prepare the resampling of the trials, or average the data if multiple trials are present and no resampling is necessary
  if (strcmp(cfg.jackknife, 'yes') || strcmp(cfg.bootstrap, 'yes') || strcmp(cfg.pseudovalue, 'yes') || strcmp(cfg.singletrial, 'yes') || strcmp(cfg.rawtrial, 'yes') || strcmp(cfg.randomization, 'yes')) && ~strcmp(data.dimord, 'rpt_chan_time')
    ft_error('multiple trials required in the data\n');

  elseif strcmp(cfg.permutation, 'yes')
    % compute the average and covariance without resampling
    [dum, avgA, covA, avgB, covB] = prepare_resampled_data(cfg2 , data.trial, data.cov, baseline.trial, baseline.cov);
    % compute the average and covariance with random permutation
    [cfg, avRA, coRA, avRB, coRB] = prepare_resampled_data(cfg, data.trial, data.cov, baseline.trial, baseline.cov);
    % concatenate the different resamplings
    avg = cat(1, reshape(avgA, [1 Nchans Nsamples]), reshape(avgB, [1 Nchans Nsamples]), avRA, avRB);
    Cy  = cat(1, reshape(covA, [1 Nchans Nchans  ]), reshape(covB, [1 Nchans Nchans  ]), coRA, coRB);
    % clear temporary working copies
    clear avgA avgB covA covB
    clear avRA avRB coRA coRB
    % the order of the resamplings should be [avgA avgB randA randB randA randB randA randB ....]
    Nrepetitions = 2*cfg.numpermutation + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    avg = avg(order,:,:);
    Cy  = Cy (order,:,:);

  elseif strcmp(cfg.randomization, 'yes')
    % compute the average and covariance without resampling
    [dum, avgA, covA, avgB, covB] = prepare_resampled_data(cfg2 , data.trial, data.cov, baseline.trial, baseline.cov);
    % compute the average and covariance with random resampling
    [cfg, avRA, coRA, avRB, coRB] = prepare_resampled_data(cfg, data.trial, data.cov, baseline.trial, baseline.cov);
    % concatenate the different resamplings
    avg = cat(1, reshape(avgA, [1 Nchans Nsamples]), reshape(avgB, [1 Nchans Nsamples]), avRA, avRB);
    Cy  = cat(1, reshape(covA, [1 Nchans Nchans  ]), reshape(covB, [1 Nchans Nchans  ]), coRA, coRB);
    % clear temporary working copies
    clear avgA avgB covA covB
    clear avRA avRB coRA coRB
    % the order of the resamplings should be [avgA avgB randA randB randA randB randA randB ....]
    Nrepetitions = 2*cfg.numrandomization + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    avg = avg(order,:,:);
    Cy  = Cy (order,:,:);

  elseif strcmp(cfg.jackknife, 'yes')
    % compute the jackknife repetitions for the average and covariance
    [cfg, avg, Cy] = prepare_resampled_data(cfg, data.trial, data.cov);
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.bootstrap, 'yes')
    % compute the bootstrap repetitions for the average and covariance
    [cfg, avg, Cy] = prepare_resampled_data(cfg, data.trial, data.cov);
    Nrepetitions = cfg.numbootstrap;

  elseif strcmp(cfg.pseudovalue, 'yes')
    % compute the pseudovalue repetitions for the average and covariance
    [cfg, avg, Cy] = prepare_resampled_data(cfg, data.trial, data.cov);
    Nrepetitions = Ntrials+1;

  elseif strcmp(cfg.singletrial, 'yes')
    % The idea is that beamformer uses the average covariance to construct the
    % filter and applies it to the single trial covariance/csd. The problem
    % is that beamformer will use the averaged covariance/csd to estimate the
    % power and not the single trial covariance/csd
    ft_error('this option contains a bug, and is therefore not supported at the moment');
    % average the single-trial covariance matrices
    Cy = mean(data.cov,1);
    % copy the average covariance matrix for every individual trial
    Cy = repmat(Cy, [Ntrials 1 1]);
    % keep the single-trial ERFs, rename them to avg for convenience
    avg = data.trial;
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.rawtrial, 'yes')
    % do not do any resampling, keep the single-trial covariances
    Cy = data.cov;
    % do not do any resampling, keep the single-trial ERFs (rename them to avg for convenience)
    avg = data.trial;
    Nrepetitions = Ntrials;

  elseif Ntrials>1
    % average the single-trial covariance matrices
    Cy  = reshape(mean(data.cov,1), [Nchans Nchans]);
    % compute the average ERF
    avg = shiftdim(mean(data.trial,1),1);
    Nrepetitions = 1;

  elseif Ntrials==1
    % select the average covariance matrix
    Cy  = data.cov;
    % select the average ERF
    avg = data.avg;
    Nrepetitions = 1;
  end

  % reshape so that it also looks like one trial (out of many)
  if Nrepetitions==1
    Cy  = reshape(Cy , [1 Nchans Nchans]);
    avg = reshape(avg, [1 Nchans Nsamples]);
  end

  % get the relevant low level options from the cfg and convert into key-value pairs
  optarg = ft_cfg2keyval(getfield(cfg, cfg.method));

  % This is the place to check for the consistency of the channel order in
  % the pre-computed leadfields/spatial filters, and to correct for it, if
  % necessary. This pertains to bugs 1746 and 3029.
  if isfield(sourcemodel, 'label') && (isfield(sourcemodel, 'leadfield') || isfield(sourcemodel, 'filter'))
    % match the channels in the leadfields/filters with those in the data
    [i1, i2] = match_str(cfg.channel, sourcemodel.label);
    if ~isequal(i2(:), (1:numel(sourcemodel.label))')
      if isfield(sourcemodel, 'leadfield')
        fprintf('\n\nSubselecting/reordering the channels in the precomputed leadfields\n\n');
        inside_indx = find(sourcemodel.inside);
        for k = inside_indx(:)'
          sourcemodel.leadfield{k} = sourcemodel.leadfield{k}(i2, :);
        end
      end
      if isfield(sourcemodel, 'filter')
        fprintf('\n\nSubselecting/reordering the channels in the precomputed filters\n\n');
        inside_indx = find(sourcemodel.inside);
        for k = inside_indx(:)'
          sourcemodel.filter{k} = sourcemodel.filter{k}(:, i2);
        end
      end
      sourcemodel.label = sourcemodel.label(i2);
    end
    if ~isequal(i1(:), (1:numel(cfg.channel))')
      % this is not so easy to deal with, throw an error
      ft_error('There''s a mismatch between the number/order of channels in the data, with respect to the channels in the precomputed leadfield/filter. This is not easy to solve automatically. Please look into this.');
    end
  end

  size_avg = [size(avg) 1];
  size_Cy  = [size(Cy) 1];
  if strcmp(cfg.method, 'lcmv')% && ~isfield(sourcemodel, 'filter')
    for i = 1:Nrepetitions
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      fprintf('scanning repetition %d\n', i);
      dip(i) = beamformer_lcmv(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, optarg{:});
    end

    % the following has been disabled since it turns out to be wrong (see
    % bugzilla bug 2395)
    %   elseif 0 && strcmp(cfg.method, 'lcmv')
    %     %don't loop over repetitions (slow), but reshape the input data to obtain single trial timecourses efficiently
    %     %in the presence of filters pre-computed on the average (or whatever)
    %     tmpdat = reshape(permute(avg,[2 3 1]),[size_avg(2) size_avg(3)*size_avg(1)]);
    %     tmpdip = beamformer_lcmv(sourcemodel, sens, headmodel, tmpdat, squeeze(mean(Cy,1)), optarg{:});
    %     tmpmom = tmpdip.mom{tmpdip.inside(1)};
    %     sizmom = size(tmpmom);
    %
    %     for i=1:length(tmpdip.inside)
    %       indx = tmpdip.inside(i);
    %       tmpdip.mom{indx} = permute(reshape(tmpdip.mom{indx}, [sizmom(1) size_avg(3) size_avg(1)]), [3 1 2]);
    %     end
    %     try, tmpdip = rmfield(tmpdip, 'pow'); end
    %     try, tmpdip = rmfield(tmpdip, 'cov'); end
    %     try, tmpdip = rmfield(tmpdip, 'noise'); end
    %     for i=1:Nrepetitions
    %       dip(i).pos     = tmpdip.pos;
    %       dip(i).inside  = tmpdip.inside;
    %       dip(i).outside = tmpdip.outside;
    %       dip(i).mom     = cell(1,size(tmpdip.pos,1));
    %       if isfield(tmpdip, 'ori')
    %         dip(i).ori   = cell(1,size(tmpdip.pos,1));
    %       end
    %       dip(i).cov     = cell(1,size(tmpdip.pos,1));
    %       dip(i).pow     = nan(size(tmpdip.pos,1),1);
    %       for ii=1:length(tmpdip.inside)
    %         indx             = tmpdip.inside(ii);
    %         tmpmom           = reshape(tmpdip.mom{indx}(i,:,:),[sizmom(1) size_avg(3)]);
    %         dip(i).mom{indx} = tmpmom;
    %         if isfield(tmpdip, 'ori')
    %           dip(i).ori{indx} = tmpdip.ori{indx};
    %         end
    %
    %         % the following recovers the single trial power and covariance, but
    %         % importantly the latency over which the power is defined is the
    %         % latency of the event-related field in the input and not the
    %         % latency of the covariance window, which can differ from the
    %         % former
    %         dip(i).cov{indx} = (tmpmom*tmpmom')./size_avg(3);
    %         if isempty(cfg.lcmv.powmethod) || strcmp(cfg.lcmv.powmethod, 'trace')
    %           dip(i).pow(indx) = trace(dip(i).cov{indx});
    %         else
    %           [tmpu,tmps,tmpv] = svd(dip(i).cov{indx});
    %           dip(i).pow(indx) = tmps(1);
    %         end
    %       end
    %     end

  elseif strcmp(cfg.method, 'sloreta')
    for i=1:Nrepetitions
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      fprintf('scanning repetition %d\n', i);
      dip(i) = ft_sloreta(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, optarg{:});
    end

  elseif strcmp(cfg.method, 'eloreta')
    for i=1:Nrepetitions
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      fprintf('scanning repetition %d\n', i);
      dip(i) = ft_eloreta(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, optarg{:});
    end
  elseif strcmp(cfg.method, 'sam')
    for i=1:Nrepetitions
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      fprintf('scanning repetition %d\n', i);
      squeeze_avg=reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      dip(i) = beamformer_sam(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, optarg{:});
    end
  elseif strcmp(cfg.method, 'pcc')
    for i=1:Nrepetitions
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      fprintf('scanning repetition %d\n', i);
      squeeze_avg=reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      dip(i) = beamformer_pcc(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, optarg{:}, 'refdip', cfg.refdip, 'refchan', refchanindx, 'supchan', supchanindx);
    end
  elseif strcmp(cfg.method, 'mne')
    for i=1:Nrepetitions
      fprintf('estimating current density distribution for repetition %d\n', i);
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      if hascovariance
        squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
        dip(i) = minimumnormestimate(sourcemodel, sens, headmodel, squeeze_avg, optarg{:}, 'noisecov', squeeze_Cy);
      else
        dip(i) = minimumnormestimate(sourcemodel, sens, headmodel, squeeze_avg, optarg{:});
      end
    end
  elseif strcmp(cfg.method, 'harmony')
    for i=1:Nrepetitions
      fprintf('estimating current density distribution for repetition %d\n', i);
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      if hascovariance
        squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
        dip(i) = harmony(sourcemodel, sens, headmodel, squeeze_avg, optarg{:}, 'noisecov', squeeze_Cy);
      else
        dip(i) = harmony(sourcemodel, sens, headmodel, squeeze_avg, optarg{:});
      end
    end
  elseif strcmp(cfg.method, 'rv')
    for i=1:Nrepetitions
      fprintf('estimating residual variance at each source position for repetition %d\n', i);
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      dip(i) = residualvariance(sourcemodel, sens, headmodel, squeeze_avg,      optarg{:});
    end
  elseif strcmp(cfg.method, 'music')
    for i=1:Nrepetitions
      fprintf('computing multiple signal classification for repetition %d\n', i);
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      if hascovariance
        squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
        dip(i) = music(sourcemodel, sens, headmodel, squeeze_avg, 'cov', squeeze_Cy, optarg{:});
      else
        dip(i) = music(sourcemodel, sens, headmodel, squeeze_avg,                    optarg{:});
      end
    end
  elseif strcmp(cfg.method, 'mvl')
    for i=1:Nrepetitions
      fprintf('estimating current density distribution for repetition %d\n', i);
      fns = fieldnames(cfg);
      optarg = cell(1,length(fns));
      n=1;
      for c=1:length(fns)
        optarg{n} = fns{c};
        optarg{n+1} = cfg.(fns{c});
        n=n+2;
      end
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      dip(i) = mvlestimate(sourcemodel, sens, headmodel, squeeze_avg, optarg{:});
    end
  else
    ft_error('method ''%s'' is unsupported for source reconstruction in the time domain', cfg.method);
  end

elseif iscomp
  ft_error('the use of component data in ft_sourceanalysis is disabled for the time being: if you encounter this error message and you need this functionality please contact the FieldTrip development team');
else
  ft_error('the specified method ''%s'' combined with the input data of type ''%s'' are not supported', cfg.method, ft_datatype(data));
end % if freq or timelock or comp data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clean up and collect the results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

source = copyfields(sourcemodel, source, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label'}); %, 'filter'});

if exist('dip', 'var')
  % the fields in the dip structure might be more recent than those in the sourcemodel structure
  source = copyfields(dip, source, {'pos', 'inside', 'leadfield', 'leadfielddimord', 'label'}); %, 'filter'});

  % prevent duplication of these fields when copying the content of dip into source.avg or source.trial
  dip    = removefields(dip,       {'pos', 'inside', 'leadfield', 'leadfielddimord', 'label'}); %, 'filter'});

  if istrue(cfg.(cfg.method).keepfilter) && isfield(dip(1), 'filter')
    for k = 1:numel(dip)
      dip(k).label        = sens.label;
      dip(k).filterdimord = '{pos}_ori_chan';
    end
  end
end

if ~istrue(cfg.keepleadfield)
  % remove the precomputed leadfields from the output source (if present)
  source = removefields(source, {'leadfield' 'leadfielddimord' 'label'});
end

% remove the precomputed leadfields from the cfg regardless of what keepleadfield is saying
% it should not be kept in cfg, since there it takes up too much space
cfg.sourcemodel = removefields(cfg.sourcemodel, {'leadfield' 'leadfielddimord' 'filter' 'filterdimord' 'label'});

if strcmp(cfg.jackknife, 'yes')
  source.method = 'jackknife';
  source.trial = dip;
  source.df = Ntrials;
elseif strcmp(cfg.bootstrap, 'yes')
  source.method = 'bootstrap';
  source.trial = dip;
  source.df = Ntrials;
elseif strcmp(cfg.pseudovalue, 'yes')
  source.method = 'pseudovalue';
  source.trial = dip;
elseif strcmp(cfg.singletrial, 'yes')
  source.method = 'singletrial';
  source.trial = dip;
  source.df = Ntrials;     % is this correct?
elseif strcmp(cfg.rawtrial, 'yes')
  source.method = 'rawtrial';
  source.trial = dip;
  source.df = Ntrials;     % is this correct?
elseif strcmp(cfg.randomization, 'yes')
  % assign the randomized resamplings to the output, keeping the special meaning of trial 1 and 2 in mind
  source.method = 'randomization';
  source.avgA = dip(1);
  source.avgB = dip(2);
  source.trialA = dip(1+2*(1:cfg.numrandomization));
  source.trialB = dip(2+2*(1:cfg.numrandomization));
elseif strcmp(cfg.permutation, 'yes')
  % assign the randomized resamplings to the output, keeping the special meaning of trial 1 and 2 in mind
  source.method = 'permutation';
  source.avgA = dip(1);
  source.avgB = dip(2);
  source.trialA = dip(1+2*(1:cfg.numpermutation));
  source.trialB = dip(2+2*(1:cfg.numpermutation));
elseif (strcmp(cfg.jackknife, 'yes') || strcmp(cfg.bootstrap, 'yes') || strcmp(cfg.pseudovalue, 'yes') || strcmp(cfg.singletrial, 'yes') || strcmp(cfg.rawtrial, 'yes')) && strcmp(cfg.keeptrials, 'yes')
  % keep the source reconstruction for each repeated or resampled trial
  source.trial = dip;
elseif exist('dip', 'var')
  % it looks like beamformer analysis was done on an average input, keep the average source reconstruction
  source.method = 'average';
  source.avg = dip;
else
  % apparently no computations were performed
end

% remember the trialinfo
if (strcmp(cfg.keeptrials, 'yes') || strcmp(cfg.method, 'pcc')) && isfield(data, 'trialinfo')
  source.trialinfo = data.trialinfo;
end

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

author: spm


text:
classdef MimEditYellow < MimGuiPlugin
    % MimEditYellow. Gui Plugin for setting paint colour
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimEditYellow is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Yellow'
        SelectedText = 'Yellow'
        ToolTip = 'Edit with yellow label'
        Category = 'Segmentation label'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'paint.png'
        IconColour = GemMarkerPoint.DefaultColours{6}
        Location = 36
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.PaintBrushColour = 6;
        end

        function enabled = IsEnabled(gui_app)
             enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Paint');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.PaintBrushColour == 6;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [bound,grad] = svr2w2 (sigma,X,Y,ker,p1,scales,slacks,opt,var2)

 % if ker=='weighted_linear' p1=1; else p1=2; end
 %% NOTE: ASSUMES AT MOMENT THAT IT IS A POLYNOMIAL or rbf, if using 'quadprog'-- not compatible w/ svmlight

 % chirag's edits: rbf kernel (non-weighted) implemented;
 % cannot use with 'svmlight' optimizer as of yet..

  l = size(X,1);
  N=size(X,2);
  Xkeep=X;

  scale=max(scales); % number of scaling factors
  if scale>0
    if length(scales)==1
      D=repmat(sigma(1:N),1,l);
    else
      D=repmat(sigma(scales),1,l);
    end
    X= X .* D';
  end

  slack=max(slacks);
  if slacks==0
    ridge=ones(l,1)*1e-11;
  else
    if length(slacks)==1 ridge=ones(l,1)*sigma(scale+1); end;
    if length(slacks)>1 ridge=sigma(scale+slacks); end;
  end
  switch ker
      case 'weighted_linear'
          K = X*X';
      case 'weighted_poly'
          K = X*X';
          K=(K+1).^p1;
      case 'rbf'
          %% create rbf kernel...
          [n1,n2]=size(X);
          kerneloption=ones(1,n2)*p1;
          metric = diag(1./p1.^2);
          ps = X*metric*X';
          [nps,pps]=size(ps);
          normx = sum(X.^2*metric,2);
          normxsup = sum(X.^2*metric,2);
          ps = -2*ps + repmat(normx,1,pps) + repmat(normxsup',nps,1) ;
          K = exp(-ps/2);
          %%
  end

  K= K + diag(ridge) .^2; %(sum(sigma)/N)*eye(l)*ridge;

  H = K .* (Y*Y');

  % construct the lower/upper bound vector
  lower = [zeros(l,1)];
  C=Inf;
  upper = [C*ones(l,1)];
  A = Y';  b = 0;   c = [-ones(l,1)];

  switch opt
   case {'loqo'}
    [a,b0] = loqo(sparse(H),c,sparse(Y'),0,lower,upper,[],1,0);

   case{'quadprog'}
    options = optimset('LargeScale','off');
    options = optimset(options,'MaxIter',400);
    [a,fval,exit,out,lambda] = quadprog(H,c,[],[],Y',0,lower,upper, ...
					[],options);
    b0=lambda.eqlin(1);
   case{'andre'}
    [a,b0] = quadsolve(H,c,Y',0,10000);

   case{'default'}
    [a,b0]=qp_learn(full(K),Y,1e5);

   case{'svmlight'}
    error('not implemented yet for rbf kernels!');
    %if strcmp(get(al.child,'ker'),'linear') ker=1; p1=1; end;
    %if strcmp(get(al.child,'ker'),'poly') ker=2;
    %p1=get(al.child,'kerparam'); p1=p1{1}; end;
     %if strcmp(get(al.child,'ker'),'gaussian') ker=3; p1= ...
     %get(al.child,'kerparam'); p1=(2*p1{1}^2);
     %end;
     ker='linear';p1=1;
     X=full(X);
     [a b0 ind] = svmlight(X,Y,1e5,ridge,0,ker,p1,0);
     alpha=zeros(size(X,1),1); alpha(ind+1)=a;
     a=abs(alpha);
  end

  % radius calculation

  if var2==0
    switch opt
     case {'loqo'}
      b=loqo(sparse(2*K),-diag(K),sparse(ones(1,l)),1,zeros(l,1),[],[],1,0);
     case {'quadprog'}
      b=quadprog(2*K,-diag(K),[],[],ones(1,l),1,zeros(l,1),[],[],options);
     case {'andre'}
       error(['calc radius: not implemented for this optimizer, use  svmlight or quadprog']);
     case {'default'}
      error(['calc radius: not implemented for this optimizer, use  svmlight or quadprog']);
     case {'svmlight'}   %% use andre anyway for r2w2 calculation
      [b] = quadsolve(2*K,-diag(K),ones(1,l),1,10000);
    end
    r=b'*diag(K)-b'*K*b;
    r=abs(r);
  else
    r = mean(diag(K))-mean(mean(K));
  end

  w = a'*H*a; %w = sum(a); %
  bound=r*w;

  XXtemp=X*X';
  if p1>1
      XXtemp=XXtemp+1;
  end

  grad= zeros(scale+slack,1);

  for i=1:scale                       %% calculate scaling factors
    Xc = Xkeep(:,find(i==scales));
    Kc = Xc*Xc';
    switch ker
        case {'weighted_poly', 'weighted_linear'}
            Kgrad= 2 * sigma(i) *  (Kc * p1) .* (XXtemp).^(p1-1);
        case 'rbf'
            x_k=Xc;
            nbx=size(x_k,1);
            deriv=(x_k*ones(1,nbx)-ones(nbx,1)*x_k').^2;
            Kgrad = (deriv.*K)./(p1^2);
    end
    gradw= -(a.*Y)'*Kgrad*(a .* Y);
    if var2==0
      gradr= b'*diag(Kgrad)- b'*Kgrad*b;
    else
      gradr= mean(diag(Kgrad))-mean(mean(Kgrad));
    end
    grad(i)=  (gradw*r + gradr*w);
  end;

  for i=1:slack                       %% calculate slacking factors
    rs=ridge; rs(find(i~=slacks))=0;
    Kgrad=2*diag(rs);
    gradw= -(a.*Y)'*Kgrad*(a .* Y);
    if var2==0
      gradr= b'*diag(Kgrad)- b'*Kgrad*b;
    else
      gradr= mean(diag(Kgrad))-mean(mean(Kgrad));
    end
    grad(i+scale)=  (gradw*r + gradr*w);
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [y, sigma, p] = rvmRegPred(model, X, t)
% Compute RVM regression model reponse y = w'*X+w0 and likelihood
% Input:
%   model: trained model structure
%   X: d x n testing data
%   t (optional): 1 x n testing response
% Output:
%   y: 1 x n prediction
%   sigma: variance
%   p: 1 x n likelihood of t
% Written by Mo Chen (sth4nth@gmail.com).
index = model.index;
w = model.w;
w0 = model.w0;

X = X(index,:);
y = w'*X+w0;
%% probability prediction
if nargout > 1
    beta = model.beta;
    U = model.U;        % 3.54
    Xo = bsxfun(@minus,X,model.xbar);
    XU = U'\Xo;
    sigma = sqrt((1+dot(XU,XU,1))/beta);   %3.59
end

if nargin == 3 && nargout == 3
    p = exp(-0.5*(((t-y)./sigma).^2+log(2*pi))-log(sigma));
end

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function [label, model, L] = mixGaussVb(X, m, prior)
% Variational Bayesian inference for Gaussian mixture.
% Input:
%   X: d x n data matrix
%   m: k (1 x 1) or label (1 x n, 1<=label(i)<=k) or model structure
% Output:
%   label: 1 x n cluster label
%   model: trained model structure
%   L: variational lower bound
% Reference: Pattern Recognition and Machine Learning by Christopher M. Bishop (P.474)
% Written by Mo Chen (sth4nth@gmail.com).
fprintf('Variational Bayesian Gaussian mixture: running ... \n');
[d,n] = size(X);
if nargin < 3
    prior.alpha = 1;
    prior.kappa = 1;
    prior.m = mean(X,2);
    prior.v = d+1;
    prior.M = eye(d);   % M = inv(W)
end
prior.logW = -2*sum(log(diag(chol(prior.M))));

tol = 1e-8;
maxiter = 2000;
L = -inf(1,maxiter);
model = init(X,m,prior);
for iter = 2:maxiter
    model = expect(X,model);
    model = maximize(X,model,prior);
    L(iter) = bound(X,model,prior);
    if abs(L(iter)-L(iter-1)) < tol*abs(L(iter)); break; end
end
L = L(2:iter);
label = zeros(1,n);
[~,label(:)] = max(model.R,[],2);
[~,~,label(:)] = unique(label);

function model = init(X, m, prior)
n = size(X,2);
if isstruct(m)  % init with a model
    model = m;
elseif numel(m) == 1  % random init k
    k = m;
    label = ceil(k*rand(1,n));
    model.R = full(sparse(1:n,label,1,n,k,n));
elseif all(size(m)==[1,n])  % init with labels
    label = m;
    k = max(label);
    model.R = full(sparse(1:n,label,1,n,k,n));
else
    error('ERROR: init is not valid.');
end
model = maximize(X,model,prior);

% Done
function model = maximize(X, model, prior)
alpha0 = prior.alpha;
kappa0 = prior.kappa;
m0 = prior.m;
v0 = prior.v;
M0 = prior.M;
R = model.R;

nk = sum(R,1); % 10.51
alpha = alpha0+nk; % 10.58
kappa = kappa0+nk; % 10.60
v = v0+nk; % 10.63
m = bsxfun(@plus,kappa0*m0,X*R);
m = bsxfun(@times,m,1./kappa); % 10.61

[d,k] = size(m);
U = zeros(d,d,k);
logW = zeros(1,k);
r = sqrt(R');
for i = 1:k
    Xm = bsxfun(@minus,X,m(:,i));
    Xm = bsxfun(@times,Xm,r(i,:));
    m0m = m0-m(:,i);
    M = M0+Xm*Xm'+kappa0*(m0m*m0m');     % equivalent to 10.62
    U(:,:,i) = chol(M);
    logW(i) = -2*sum(log(diag(U(:,:,i))));
end

model.alpha = alpha;
model.kappa = kappa;
model.m = m;
model.v = v;
model.U = U;
model.logW = logW;

% Done
function model = expect(X, model)
alpha = model.alpha; % Dirichlet
kappa = model.kappa;   % Gaussian
m = model.m;         % Gasusian
v = model.v;         % Whishart
U = model.U;         % Whishart
logW = model.logW;
n = size(X,2);
[d,k] = size(m);

EQ = zeros(n,k);
for i = 1:k
    Q = (U(:,:,i)'\bsxfun(@minus,X,m(:,i)));
    EQ(:,i) = d/kappa(i)+v(i)*dot(Q,Q,1);    % 10.64
end
ElogLambda = sum(psi(0,0.5*bsxfun(@minus,v+1,(1:d)')),1)+d*log(2)+logW; % 10.65
Elogpi = psi(0,alpha)-psi(0,sum(alpha)); % 10.66
logRho = -0.5*bsxfun(@minus,EQ,ElogLambda-d*log(2*pi)); % 10.46
logRho = bsxfun(@plus,logRho,Elogpi);   % 10.46
logR = bsxfun(@minus,logRho,logsumexp(logRho,2)); % 10.49
R = exp(logR);

model.logR = logR;
model.R = R;

% Done
function L = bound(X, model, prior)
alpha0 = prior.alpha;
kappa0 = prior.kappa;
v0 = prior.v;
logW0 = prior.logW;
alpha = model.alpha;
kappa = model.kappa;
v = model.v;
logW = model.logW;
R = model.R;
logR = model.logR;
[d,n] = size(X);
k = size(R,2);

Epz = 0;
Eqz = dot(R(:),logR(:));
logCalpha0 = gammaln(k*alpha0)-k*gammaln(alpha0);
Eppi = logCalpha0;
logCalpha = gammaln(sum(alpha))-sum(gammaln(alpha));
Eqpi = logCalpha;
Epmu = 0.5*d*k*log(kappa0);
Eqmu = 0.5*d*sum(log(kappa));
logB0 = -0.5*v0*(logW0+d*log(2))-logMvGamma(0.5*v0,d);
EpLambda = k*logB0;
logB =  -0.5*v.*(logW+d*log(2))-logMvGamma(0.5*v,d);
EqLambda = sum(logB);
EpX = -0.5*d*n*log(2*pi);
L = Epz-Eqz+Eppi-Eqpi+Epmu-Eqmu+EpLambda-EqLambda+EpX;


text:
function milktruck_demo
% This function is an example on how to integrate GoogleEarth Milktruckdemo in a Matlab figure or GUI and
% control it via Matlab controls.
%
% You need the GoogleEarth plugin to be installed. Therefore just visit a website, where the plugin is needed
% and install it, e.g.:  http://code.google.com/intl/de-DE/apis/earth/
%
% The rest is more simple:
%
% The workaround for integrating GE in Matlab GUI is simple to embed a Microsoft Internet ActiveX-Control in
% the GUI.
%
% Further Readings:
%   http://code.google.com/intl/de-DE/apis/earth/
%   http://code.google.com/intl/de-DE/apis/earth/documentation/reference/index.html
% see in Matlab Help: actx_explore for integrating Shell-Explorer
%
% Programmed by Sven Koerner: koerner(underline)sven(add)gmx.de
% Date: 2011/08/04

% run Demo
MilktruckMA


text:
function [ decision ] = face_check_cnn( img, shape, global_params, cnns )
%FACE_CHECK_CNN Summary of this function goes here
%   Detailed explanation goes here

%%
if(size(img,3) == 3)
    img = rgb2gray(img);
end
% first need to determine the view
centres = cat(1, cnns.centres);

dists = centres*pi/180 - repmat(global_params(2:4)',size(centres,1),1);
[~,view_id] = min(sum(dists.^2,2));

mask_small = cnns(view_id).mask(1:size(cnns(view_id).triX,1), 1:size(cnns(view_id).triX,2));

if(size(cnns(view_id).destination,1) == 66 && size(shape,1) == 68)
    label_inds = [1:60,62:64,66:68];
    shape = shape(label_inds,:);
end
img_crop = Crop(img, shape, cnns(view_id).triangulation,...
    cnns(view_id).triX, mask_small,...
    cnns(view_id).alphas, cnns(view_id).betas,...
    cnns(view_id).nPix, cnns(view_id).minX, ...
    cnns(view_id).minY);

%%
img_crop = reshape(img_crop(logical(cnns(view_id).mask)), 1, cnns(view_id).nPix);
img_crop(isnan(img_crop)) = 0;

% normalisation (local)
img_crop = (img_crop - mean(img_crop));
norms = std(img_crop);
if(norms==0)
    norms = 1;
end
img_crop = img_crop / norms;

% normalisation (global)
img_crop = img_crop - cnns(view_id).mean_ex;
img_crop = img_crop ./ cnns(view_id).std_ex;

mask = cnns(view_id).mask;

img = zeros(size(mask));
img(mask) = img_crop;

cnn = cnns(view_id).cnn;

%%
res = vl_simplenn(cnn, single(img), [], []);
res = gather(res(end).x);
num_bins = numel(res);
[~,res] = sort(res, 3, 'descend') ;
res = squeeze(res);
res = res(1,:);

decision = unQuantizeContinuous(res, 0, 3, num_bins)';

author: TadasBaltrusaitis

repo: OpenFace


text:
function [types, isSecondary] = classifyMoieties(L, S)
% Classifies conserved moieties for a metabolic network
%
% USAGE:
%
%    types = classifyMoieties(L, S)
%
% INPUTS:
%    L:        The `r` x `m` moiety matrix with moiety vectors as columns.
%    S:        The `m` x `n` total stoichiometric matrix.
%
% OUTPUT:
%    types:    an `r` x `1` cell array of with one of the following moiety classifications
%              'Internal' moiety that is also conserved in the open network
%              'Transitive' moiety that is only found in primary metabolites
%              'Integrative' moiety that is not conserved in the open
%               network and found in both primary and secondary metabolites.
%
%
% isSecondary  `m x 1` Boolean vector indicating secondary metabolites
%                      (containing at least one internal conserved moiety')
%
% .. Author: - Hulda S. Haraldsdóttir, June 2015
%              Ronan Fleming, Oct 2020

types = cell(size(L,1),1);

isInternal = ~any(L*S,2); % Internal moieties are conserved in the open network
isSecondary = any(L(isInternal,:),1); % Secondary metabolites contain internal moieties
isTransitive = ~any(L(:,isSecondary),2); % Transitive moieties are only found in primary metabolites
isIntegrative = ~(isTransitive | isInternal); % All other moieties are Integrative

types(isTransitive) = {'Transitive'};
types(isIntegrative) = {'Integrative'};
types(isInternal) = {'Internal'};

author: opencobra

repo: cobratoolbox


text:
%% -------------
function s = ComputeSaliency(img, sigma, alpha)

    % --------- Check the input --------
    if (1 ~= size(img, 3)),
        error('The input image should be GRAY.');
    end

    [dx, dy] = GradientMethod(double(img), 'zhou');
    grad = dx +1j*dy;

    [~, cc] = EigDecBlock(grad, sigma);
    wt = sqrt((sqrt(cc(1,:,:))+sqrt(cc(2,:,:))).^2 + alpha*(sqrt(cc(1,:,:))-sqrt(cc(2,:,:))).^2);
    s = squeeze(wt);
end

%% ------------------------------------------------
% [c11, c12    [dxx, dxy
%  c21, c22] =  dyx, dyy]
% B = -(c11+c22), C = c11*c22-c12*c21
function [postMap, ss] = EigDecBlock(img, sigma)

    winSize = ceil(sigma*6);
    if ~mod(winSize, 2),
        winSize = winSize + 1;
    end

    h = fspecial('gauss', [winSize winSize], sigma);
    [hh, ww] = size(img);
    ss = zeros(2, hh, ww);

    dx = real(img);
    dy = imag(img);
    dxx = imfilter(dx.*dx, h, 'symmetric');
    dxy = imfilter(dx.*dy, h, 'symmetric');
    dyy = imfilter(dy.*dy, h, 'symmetric');

    A = ones(size(img));
    B = -(dxx+dyy);
    C = dxx.*dyy - dxy.*dxy;

    ss(1, :, :) = abs((-B+sqrt(B.^2-4*A.*C))./(2*A));
    ss(2, :, :) = abs((-B-sqrt(B.^2-4*A.*C))./(2*A));

    V12 = (dxx-dyy + sqrt((dxx-dyy).^2+4*dxy.*dxy))./(2*dxy+eps);

    postMap = sqrt(squeeze(ss(1, :, :))).*(V12 + 1i)./sqrt(V12.^2+1+eps);

end

author: thfylsty


text:
%% Analyzing Neural Time Series Data
% Matlab code for Chapter 12
% Mike X Cohen
%
% This code accompanies the book, titled "Analyzing Neural Time Series Data"
% (MIT Press). Using the code without following the book may lead to confusion,
% incorrect data analyses, and misinterpretations of results.
% Mike X Cohen assumes no responsibility for inappropriate or incorrect use of this code.

%% Figure 12.1

load sampleEEGdata

time = -1:1/EEG.srate:1;
f = 4; % frequency of sine wave in Hz

% create sine wave (actually, a cosine wave, for reasons that will become
% clear in Chapter 13)
sine_wave = cos(2*pi*f.*time);

% make a Gaussian
s=4/(2*pi*f);
gaussian_win = exp(-time.^2./(2*s^2));

figure
plot(time,sine_wave.*gaussian_win)

%% Figure 12.2

figure
subplot(511)
plot(squeeze(EEG.data(47,:,1)))
axis tight

subplot(512)
sine_wave = cos(2*pi*12.*time);
plot(sine_wave)
axis tight

subplot(513)
% boxcar envelope
boxcar = zeros(size(sine_wave));
midpoint = (length(time)-1)/2;
boxcar(midpoint-round(EEG.srate/12/5):midpoint+round(EEG.srate/12/1.25)) = 1;
plot(sine_wave.*boxcar)
axis tight

subplot(514)
% boxcar envelope
boxcar = zeros(size(sine_wave));
midpoint = (length(time)-1)/2;
boxcar(midpoint-50:midpoint+50) = 1;
plot(sine_wave.*boxcar)
axis tight

subplot(515)
% redefine gaussian for new sine wave
s=1.5/(2*pi*12); gaussian_win = exp(-time.^2./(2*s^2));
plot(time,sine_wave.*gaussian_win)
axis tight

%% Figure 12.3

srate = 500; % sampling rate in Hz
f = 10; % frequency of the sine wave in Hz
time = -1:1/srate:1; % time, from -1 to 1 second in steps of 1/sampling-rate

sine_wave = exp(2*pi*1i*f.*time); % complex wavelet

% make a Gaussian
s=6/(2*pi*f);
gaussian_win = exp(-time.^2./(2*s^2));

% and together they make a wavelet!
wavelet = sine_wave .* gaussian_win;

% make plots containing each component
figure
subplot(311)
plot(time,real(sine_wave)) % only plot the real component... we�ll learn about this later
title('Sine wave')

subplot(312)
plot(time,gaussian_win) % plots the Gaussian window
title('Gaussian window')

subplot(313)
plot(time,real(wavelet)); % plots the wavelet
title('My first wavelet!')
xlabel('Time (ms)')

%% Figure 12.4

num_wavelets      =  80;  % number of frequency bands
lowest_frequency  =   2;  % in Hz
highest_frequency = 100;  % in Hz

frequencies=linspace(lowest_frequency,highest_frequency,num_wavelets);
% note: the "linspace" function creates linearly spaced numbers between the first and second
% inputs, with the number of steps corresponding to the third input.
figure, plot(frequencies,'-*')
xlabel('Frequency order')
ylabel('Frequency in Hz')

% initialize wavelet family
wavelet_family = zeros(num_wavelets,length(time));

% Loop through frequencies and make a family of wavelets.
for fi=1:num_wavelets

    % create a sine wave at this frequency
    sinewave = exp(2*1i*pi*frequencies(fi).*time); % the "1i" makes it a complex wavelet

    % create a Gaussian window
    gaus_win = exp(-time.^2./(2*(6/(2*pi*frequencies(fi)))^2));

    % create wavelet via element-by-element multiplication of the sinewave and gaussian window
    wavelet_family(fi,:) = sinewave.*gaus_win;

    % note that you can also do this on one line:
    wavelet_family(fi,:) = exp(2*1i*pi*frequencies(fi).*time) .* exp(-time.^2./(2*(6/(2*pi*frequencies(fi)))^2));
end

% Plot a few wavelets
figure
subplot(2,1,1)
plot(time,real(wavelet_family(1:round(rand*30):end,:))')
title('A few wavelets...')

% Note that in the subplot command you don't need commas if you have fewer than 10
% rows/cols and if you are not using any variables.
subplot(212)
plot(time,real(wavelet_family(30,:)))
hold on
plot(time,imag(wavelet_family(30,:)),':')
title('Real and imaginary parts of one wavelet')
legend({'real';'imaginary'})

% finally, image the wavelet family.
figure
imagesc(time,frequencies,real(wavelet_family))
axis xy % equivalent to "set(gca,'ydir','normal')
xlabel('Time (s)')
ylabel('Frequency (Hz)')

%% bonus feature

% Running this cell will generate a movie that shows how a line -- in this
% case, a wavelet -- goes from a line that goes up and down, to a colored
% 'flat' line in a 2-D image. This is the principle that underlies viewing
% time-frequency plots.

figure, set(gcf,'color','k')
surf(repmat(real(wavelet),2,1))
shading interp
axis off
axis([0 length(time) -20 21 -1 1])
view([ -4 4 ])

for i=4:2:90
    view([ -4 i ])
    pause(.1)
end

rotate3d

%% Figure 12.5

% EEG data from one trial (electrode FCz)
eegdata = squeeze(EEG.data(47,:,10));

% create wavelet
time = -1:1/EEG.srate:1;
f = 6; % frequency of sine wave in Hz
sine_wave = exp(1i*2*pi*f.*time);
s = 4.5/(2*pi*f);
gaussian_win = exp(-time.^2./(2*s^2));
wavelet = sine_wave .* gaussian_win;
% half of the wavelet size, useful for chopping off edges after convolution.
halfwaveletsize = ceil(length(wavelet)/2);

% convolve with data
% compute Gaussian
n_conv = length(wavelet) + EEG.pnts - 1;

fft_w = fft(wavelet,n_conv);
fft_e = fft(eegdata,n_conv);
ift   = ifft(fft_e.*fft_w,n_conv)*sqrt(s)/10; % sqrt... is an empirical scaling factor that works here
wavelet_conv_data = real(ift(halfwaveletsize:end-halfwaveletsize+1));

% create filter and apply to data
% (more on how to interpret this code in a few chapters!)
nyquist       = EEG.srate/2;
transition_width = 0.2; % percent
filter_low    = 4; % Hz
filter_high   = 8; % Hz
ffrequencies  = [ 0 filter_low*(1-transition_width) filter_low filter_high filter_high*(1+transition_width) nyquist ]/nyquist;
idealresponse = [ 0 0 1 1 0 0 ];
filterweights = firls(round(3*(EEG.srate/filter_low)),ffrequencies,idealresponse);
eeg_4to8      = filtfilt(filterweights,1,double(eegdata));

% now plot all the pieces
figure

plot(EEG.times,eegdata)
hold on
plot(EEG.times,wavelet_conv_data,'r','linew',2)
plot(EEG.times,eeg_4to8,'m','linew',2)
set(gca,'xlim',[-200 1200],'ydir','r')
xlabel('Time (ms)'), ylabel('Voltage (\muV)')
legend({'Raw data';'wavelet convolved';'band-pass filtered'},0)

%% Figure 12.6

% make a theta-band-centered wavelet
time   = -1:1/EEG.srate:1;
n_conv = EEG.pnts + length(time) - 1;
n2p1   = floor(n_conv/2)+1; % n2p1 = n/2+1

f = 6;
s = 6/(2*pi*f);
wavelet = exp(2*pi*1i*f.*time) .* exp(-time.^2./(2*s^2));
halfwaveletsize = ceil(length(wavelet)/2);

eegdata = squeeze(EEG.data(47,:,10));

figure

subplot(311)
plot(EEG.times,eegdata)
set(gca,'xlim',[-500 1200])

subplot(323)
fft_w = fft(wavelet,n_conv);
hz    = linspace(0,EEG.srate/2,n2p1);
plot(hz,abs(fft_w(1:n2p1))./max(abs(fft_w(1:n2p1))),'k')
hold on

fft_e = fft(eegdata,n_conv);
hz    = linspace(0,EEG.srate/2,n2p1);
plot(hz,abs(fft_e(1:n2p1))./max(abs(fft_e(1:n2p1))),'r')
set(gca,'xlim',[0 40],'ylim',[0 1.05])
title('individual power spectra')

subplot(324)
plot(hz,abs(fft_e(1:n2p1)).*abs(fft_w(1:n2p1)))
set(gca,'xlim',[0 40])

subplot(313)
plot(EEG.times,eegdata)
hold on
ift = ifft(fft_e.*fft_w,n_conv)*sqrt(s)/10;
plot(EEG.times,real(ift(halfwaveletsize:end-halfwaveletsize+1)),'r')
set(gca,'xlim',[-500 1200])

%% Figure 12.7

% create 10 Hz wavelet (kernel)
time = -EEG.pnts/EEG.srate/2 : 1/EEG.srate : EEG.pnts/EEG.srate/2-1/EEG.srate;
f    = 10; % frequency of sine wave in Hz
s    = 4/(2*pi*f);
wavelet = cos(2*pi*f.*time) .* exp(-time.^2./(2*s^2));

% signal is one sine cycle
timeS  = 0:1/EEG.srate:(1/f); % one cycle is 1/f
signal = sin(2*pi*f.*timeS);

% now zero-pad signal
signal = [ zeros(1,EEG.pnts/2-length(timeS)/2) signal zeros(1,EEG.pnts/2-length(timeS)/2) ];

figure

% plot waves
subplot(321)
plot(wavelet)
set(gca,'xlim',[200 length(time)-200])

subplot(323)
plot(signal)
set(gca,'xlim',[200 length(time)-200])

subplot(325)
plot(conv(wavelet,signal,'same'))
set(gca,'xlim',[200 length(time)-200],'ylim',[-12 12])

% now plot dot products at selected phase lags
subplot(322)
plot(wavelet(round(100/f)-2:end),'r')
hold on
plot(signal)
set(gca,'xlim',[200 length(time)-200])
title([ 'dot product: ' num2str( fix(sum(wavelet(round(100/f)-2:end).*signal(1:end-round(100/f)+3))) ) ])

subplot(324)
plot(wavelet(round(2.3*(100/f)-2):end),'r')
hold on
plot(signal)
set(gca,'xlim',[200 length(time)-200])
title([ 'dot product: ' num2str( fix(sum(wavelet(round(2.3*(100/f)-2):end).*signal(1:end-round(2.3*(100/f)-3))) )) ])

subplot(326)
plot(wavelet,'r')
hold on
plot(signal)
set(gca,'xlim',[200 length(time)-200])
title([ 'dot product: ' num2str( fix(sum(wavelet.*signal)) ) ])

%% end

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
% Sample object to test data persistence functions and user controls.

% Copyright 2008-2009 Levente Hunyadi
classdef SampleObject
    properties
        % A scalar of the standard MatLab type.
        RealDoubleScalar = pi;
        % A matrix of the standard MatLab type.
        RealDoubleMatrix = [1.5,2.5;3.5,4.5];
        RealSingleMatrix = single([1.5,2.5;3.5,4.5]);
        ComplexDoubleMatrix = [1.5,2.5i;3.5i,4.5];
        ComplexDoubleScalar = 1.5 + 2.5i;
        IntegerScalar = int32(32);
        IntegerMatrix = int32([1,2;3,4]);
        Logical = true;
        LogicalMatrix = [true,false;false,true];
        String = 'this is a sample string';
        CellArrayOfStrings = {'this','is','a','cell','array','of','strings'};
        Caption = 'SampleObject';
        % An unchangeable property.
        % Whenever the setter method of this property is used to modify the
        % value, an exception is thrown.
        Unchangeable = 'unchangeable';
        NestedObject = SampleNestedObject;
    end
    properties (Constant)
        Re = 6.378e6;  % Earth radius (m)
        RealDoubleConstantMatrix = [1.5,2.5;3.5,4.5];
    end
    properties (Dependent)
        DependentProperty;
        ReadOnlyDependentProperty;
    end
    properties (Hidden)
        HiddenProperty = 'hidden';
    end
    properties (Transient)
        TransientProperty = 'transient';
    end
    properties (Access = private)
        PrivateProperty = 'private';
    end
    properties (Access = protected)
        ProtectedProperty = 'protected';
    end
    methods
        function Method(self)
            disp(self.PrivateProperty);
        end

        function self = set.Caption(self, caption)
            caption = deblank(caption);
            fprintf('Caption has been changed from %s to %s.\n', self.Caption, caption);
            self.Caption = caption;
        end

        function value = get.DependentProperty(self)
            value = self.Caption;
        end

        function self = set.DependentProperty(self, value)
            self.Caption = value;
        end

        function value = get.ReadOnlyDependentProperty(self)
            value = [ self.Caption ' - ' self.Caption ];
        end

        function self = set.Unchangeable(self, value) %#ok<INUSL>
            error('SampleObject:InvalidOperation', 'Cannot change property value to %s', value);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [u,errH1,err,node,elem] = mgfracLap2d(pde,cube,h,s,option)

%% Generate matrices in each level
global mA
level = -log2(h);
mA = cell(level,1);

option.solver = 'none';
mh =2^(level-2)*h;
for i = 2:level
    [u,eqn,info,node,elem,Neumann] = fracLap2d(cube,mh,pde,option);
    mA{i} = eqn.A;
    mh = mh/2;
end
A = mA{level};
b = eqn.b;
freeNode = eqn.freeNode;

%% Iterative method
tol = option.tol;
k = 1;
maxIt = 160;
err = zeros(maxIt,1);
err(1) = 1;
u = zeros(length(b),1);
nb = norm(b(freeNode));
% u = rand(length(b),1);
while err(k) > tol && k < maxIt
    r = b - A*u;
    e = mgVcyclefracLap2d(A,r,pde,cube,h,s,option);
    u(freeNode) = u(freeNode) + e(freeNode);
    err(k+1) = norm(r(freeNode))/nb;
%     fprintf('%d-iteration with error %e \n',k,err(k+1));
    k = k + 1;
end
err = err(1:k);

%% Compute the error
errH1 = getH1error3bd(node,Neumann,pde,u,A,5);
% errH1 = getH1errorbd(node,Neumann,pde,u,eqn.A,5);
% uI = pde.exactu(node); % nodal interpolation
% errH1 = sqrt((u-uI)'*A*(u-uI));

repo: ifem


text:
function bc = specific_bc(xbd,ybd)
%hotwall_bc   Reference problem 3.4  boundary condition
%   bc = specific_bc(xbd,ybd);
%   input
%          xbd          x boundary coordinate vector
%          ybd          y boundary coordinate vector
%
%   specifies hot wall /Morton pp.10/
%   IFISS function: DJS; 5 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
nobd=length(xbd);
bc=zeros(size(xbd));
k=find(xbd==1);
bc(k)=1;
return

author: johannesgerer

repo: jburkardt-m


text:
function OptimizationModel = buildMTAproblemFromModel(model,rxnFBS,Vref,varargin)
% Returns the COBRA Optimization model needed to perform the MTA
%
% USAGE:
%
%       OptimizationModel = buildMTAproblemFromModel(model,rxnFBS,Vref,alpha,epsilon)
%
% INPUT:
%    model:                 Metabolic model (COBRA format)
%    rxnFBS:                Forward, Backward and Unchanged (+1;0;-1) values
%                           corresponding to each reaction.
%    Vref:                  Reference flux of the source state.
%    alpha:                 parameter of the quadratic problem (default = 0.66)
%    epsilon                minimun disturbance for each reaction, (default = 0)
%
% OUTPUTS:
%    OptimizationModel:     COBRA model struct that includes the
%                           stoichiometric contrains, the thermodinamic
%                           constrains and the binary variables.
%
% .. Authors:
%       - Luis V. Valcarcel, 03/06/2015, University of Navarra, CIMA & TECNUN School of Engineering.
%       - Luis V. Valcarcel, 26/10/2018, University of Navarra, CIMA & TECNUN School of Engineering.

p = inputParser; % check the inputs
% check requiered arguments
addRequired(p, 'model');
addRequired(p, 'rxnFBS', @isnumeric);
addRequired(p, 'Vref', @isnumeric);
% Check optional arguments
addOptional(p, 'alpha', 0.66, @isnumeric);
addOptional(p, 'epsilon', zeros(size(model.rxns)), @isnumeric);
% extract variables from parser
parse(p, model, rxnFBS, Vref, varargin{:});
alpha = p.Results.alpha;
epsilon = p.Results.epsilon;

% sometimes epsilon can be given as a single value
if numel(epsilon)==1
    epsilon = epsilon * ones(size(model.rxns));
end

%% --- set the COBRA model ---

% variables
v = 1:length(model.rxns);
y_plus_F = (1:sum(rxnFBS==+1)) + v(end);             % 1 if change in rxnForward, 0 otherwise
y_minus_F = (1:sum(rxnFBS==+1)) + y_plus_F(end);     % 1 if no change in rxnForward, 0 otherwise
y_plus_B = (1:sum(rxnFBS==-1)) + y_minus_F(end);     % 1 if change in rxnBackward, 0 otherwise
y_minus_B = (1:sum(rxnFBS==-1)) + y_plus_B(end);     % 1 if no change in rxnBackward, 0 otherwise
n_var = y_minus_B(end);

% limits of the variables
lb = zeros(n_var,1);
ub = ones (n_var,1);
lb(v) = model.lb;
ub(v) = model.ub;

%type of variables
vartype(1:n_var) = 'B';
vartype(v) = 'C';

% constrains
Eq1 = 1:length(model.mets);                 % Stoichiometric matrix
Eq2 = (1:length(y_plus_F)) + Eq1(end);      % Changes in Forward
Eq3 = (1:length(y_plus_F)) + Eq2(end);      % Change or not change in Forward
Eq4 = (1:length(y_plus_B)) + Eq3(end);      % Changes in Backward
Eq5 = (1:length(y_plus_B)) + Eq4(end);      % Change or not change in Backward
nCon = Eq5(end);

% generate constrain matrix
A = spalloc(nCon, n_var, nnz(model.S) + 5*length(Eq2) + 5*length(Eq4));
b = zeros(nCon,1);
csense = char(zeros(nCon,1));

posF = find(rxnFBS == +1);
posB = find(rxnFBS == -1);
posS = find(rxnFBS == 0);

% First contraint, stoichiometric
A(Eq1,v) = model.S;
b(Eq1) = 0;
csense(Eq1) = 'E';

% Second contraint, Change or not change in Forward
A(Eq2,v(posF)) = eye(length(posF));
A(Eq2,y_plus_F) = - ( Vref(posF) + epsilon(posF) ) .* eye(length(posF));
A(Eq2,y_minus_F) = - model.lb(posF) .* eye(length(posF));
b(Eq2) = 0;
csense(Eq2) = 'G';

% Third contraint, Change or not change in Forward
A(Eq3,y_plus_F) = eye(length(Eq3));
A(Eq3,y_minus_F) = eye(length(Eq3));
b(Eq3) = 1;
csense(Eq3) = 'E';

% Fourth contraint, Backward changes
A(Eq4,posB) = eye(length(posB));
A(Eq4,y_plus_B) = - ( Vref(posB) - epsilon(posB) ) .* eye(length(posB));
A(Eq4,y_minus_B) = - model.ub(posB) .* eye(length(posB));
b(Eq4) = 0;
csense(Eq4) = 'L';

% Fiveth contraint, Change or not change in Backward
A(Eq5,y_plus_B) = eye(length(Eq5));
A(Eq5,y_minus_B) = eye(length(Eq5));
b(Eq5) = 1;
csense(Eq5) = 'E';

% Objective fuction
% linear part
c = zeros(n_var,1);
c(y_minus_F) = alpha/2;
c(y_minus_B) = alpha/2;
c(v(posS)) = -2 * Vref(posS) * (1-alpha);
% quadratic part
F = spalloc(n_var,n_var,length(posS));
F(v(posS),v(posS)) =  2 * (1-alpha) .* eye(length(posS));

% save the resultant model
OptimizationModel = struct();
[OptimizationModel.A, OptimizationModel.lb, OptimizationModel.ub] = deal(A, lb, ub);
[OptimizationModel.b, OptimizationModel.csense] = deal(b, csense);
[OptimizationModel.c, OptimizationModel.F] = deal(c, F);
[OptimizationModel.osense, OptimizationModel.vartype] = deal(+1, vartype); % +1 for minimization

%save the index of the variables
OptimizationModel.idx_variables.v = v;
OptimizationModel.idx_variables.y_plus_F = y_plus_F;
OptimizationModel.idx_variables.y_minus_F = y_minus_F;
OptimizationModel.idx_variables.y_plus_B = y_plus_B;
OptimizationModel.idx_variables.y_minus_B = y_minus_B;

end


text:
function [mi] = yd2mi(yd)
% Convert length from yards to miles.
% Chad Greene 2012
mi = yd*0.0005681818181818;


text:
function yout = poolDataLIST_CONTROL(yin,ahat,nVars,polyorder,usesine)
% Copyright 2015, All Rights Reserved
% Code by Steven L. Brunton
% For Paper, "Discovering Governing Equations from Data:
%        Sparse Identification of Nonlinear Dynamical Systems"
% by S. L. Brunton, J. L. Proctor, and J. N. Kutz

n = size(yin,1);

ind = 1;
% poly order 0
yout{ind,1} = ['1'];
ind = ind+1;

% poly order 1
for i=1:nVars
    yout(ind,1) = yin(i);
    ind = ind+1;
end

if(polyorder>=2)
    % poly order 2
    for i=1:nVars
        for j=i:nVars
            yout{ind,1} = [yin{i},yin{j}];
            ind = ind+1;
        end
    end
end

if(polyorder>=3)
    % poly order 3
    for i=1:nVars
        for j=i:nVars
            for k=j:nVars
                yout{ind,1} = [yin{i},yin{j},yin{k}];
                ind = ind+1;
            end
        end
    end
end

if(polyorder>=4)
    % poly order 4
    for i=1:nVars
        for j=i:nVars
            for k=j:nVars
                for l=k:nVars
                    yout{ind,1} = [yin{i},yin{j},yin{k},yin{l}];
                    ind = ind+1;
                end
            end
        end
    end
end

if(polyorder>=5)
    % poly order 5
    for i=1:nVars
        for j=i:nVars
            for k=j:nVars
                for l=k:nVars
                    for m=l:nVars
                        yout{ind,1} = [yin{i},yin{j},yin{k},yin{l},yin{m}];
                        ind = ind+1;
                    end
                end
            end
        end
    end
end

if(usesine)
    for k=1:10;
        yout{ind,1} = ['sin(',num2str(k),'*yin)'];
        ind = ind + 1;
        yout{ind,1} = ['cos(',num2str(k),'*yin)'];
        ind = ind + 1;
    end
end

output = yout;
newout(1) = {''};
newout{1,2} = 'u';
% for k=1:length(yin)
%     newout{1,1+k} = [yin{k},'dot'];
% end
% newout = {'','xdot','ydot','udot'};
for k=1:size(ahat,1)
    newout(k+1,1) = output(k);
    for j=1:size(ahat,2)
        newout{k+1,1+j} = ahat(k,j);
    end
end
newout

author: eurika-kaiser

repo: SINDY-MPC


text:
% Hopf_Bif - Animation of Hopf bifurcation of a limit cycle from the origin.
% NOTE: Hopf_System.m must be in the same directory as Hopf_Bif.m.
% Copyright Springer 2013.
clear
Max=120;global b;
for j = 1:Max
    F(j) = getframe;
    b=j/40;   % mu goes from 0 to 4.
    options = odeset('RelTol',1e-4,'AbsTol',1e-4);
    x0=1;y0=1;
    [t,x]=ode45(@Hopf_System,[0 100],[x0 y0],options);
    plot(x(:,1),x(:,2),'b');
    axis([0 5 0 5])
    fsize=15;
    set(gca,'xtick',0:1:5,'FontSize',fsize)
    set(gca,'ytick',0:1:5,'FontSize',fsize)
    xlabel('x(t)','FontSize',fsize)
    ylabel('y(t)','FontSize',fsize)
    title('Hopf Bifurcation','FontSize',15);
    F(j) = getframe;
end
movie(F,5)
% End of Program.