text:
function [cfg, artifact] = ft_artifact_tms(cfg, data)

% FT_ARTIFACT_TMS reads the data segments of interest from file and identifies
% artefacts in EEG recordings that were done during TMS stimulation.
%
% Use as
%   [cfg, artifact] = ft_artifact_tms(cfg)
% with the configuration options
%   cfg.dataset     = string with the filename
% or
%   cfg.headerfile  = string with the filename
%   cfg.datafile    = string with the filename
% and optionally
%   cfg.headerformat
%   cfg.dataformat
%
% Alternatively you can use it as
%   [cfg, artifact] = ft_artifact_tms(cfg, data)
% where the input data is a structure as obtained from FT_PREPROCESSING.
%
% In both cases the configuration should also contain
%   cfg.trl         = structure that defines the data segments of interest, see FT_DEFINETRIAL
%   cfg.continuous  = 'yes' or 'no' whether the file contains continuous data (default = 'yes')
% and
%   cfg.method      = 'detect' or 'marker', see below.
%   cfg.prestim     = scalar, time in seconds prior to onset of detected event to mark as artifactual (default = 0.005 seconds)
%   cfg.poststim    = scalar, time in seconds post onset of detected even to mark as artifactual (default = 0.010 seconds)
%
% The different methods are described in detail below.
%
% With cfg.method='detect', TMS-artifact are detected on basis of transient
% high-amplidude gradients that are typical for TMS-pulses. The data is preprocessed
% (again) with the following settings, which are optimal for identifying TMS-pulses.
% Artifacts are identified by means of thresholding the z-transformed value of the
% preprocessed data. This method acts as a wrapper around FT_ARTIFACT_ZVALUE.
%   cfg.artfctdef.tms.derivative  = 'yes'
%   cfg.artfctdef.tms.channel     = Nx1 cell-array with selection of channels, see FT_CHANNELSELECTION for details
%   cfg.artfctdef.tms.cutoff      = z-value at which to threshold (default = 4)
%   cfg.artfctdef.tms.trlpadding  = 0.1
%   cfg.artfctdef.tms.fltpadding  = 0.1
%   cfg.artfctdef.tms.artpadding  = 0.01
% Be aware that if one artifact falls within this specified range of another
% artifact, both artifact will be counted as one. Depending on cfg.prestim and
% cfg.poststim you may not mark enough data as artifactual.
%
% With cfg.method='marker', TMS-artifact onsets and offsets are based on
% markers/triggers that are written into the EEG dataset. This method acts as a
% wrapper around FT_DEFINETRIAL to determine on- and offsets of TMS pulses by reading
% markers in the EEG.
%   cfg.trialfun            = function name, see below (default = 'ft_trialfun_general')
%   cfg.trialdef.eventtype  = 'string'
%   cfg.trialdef.eventvalue = number, string or list with numbers or strings
% The cfg.trialfun option is a string containing the name of a function that you
% wrote yourself and that FT_ARTIFACT_TMS will call. The function should take the
% cfg-structure as input and should give a NxM matrix with M>=3 in the same format as
% "trl" as the output. You can add extra custom fields to the configuration structure
% to pass as arguments to your own trialfun. Furthermore, inside the trialfun you can
% use the FT_READ_EVENT function to get the event information from your data file.
%
% The output argument "artifact" is a Nx2 matrix comparable to the "trl" matrix of
% FT_DEFINETRIAL. The first column of which specifying the beginsamples of an
% artifact period, the second column contains the endsamples of the artifactperiods.
%
% To facilitate data-handling and distributed computing, you can use
%   cfg.inputfile   =  ...
% to read the input data from a *.mat file on disk. This mat files should contain
% only a single variable named 'data', corresponding to the input structure.
%
% See also FT_REJECTARTIFACT, FT_ARTIFACT_CLIP, FT_ARTIFACT_ECG, FT_ARTIFACT_EOG,
% FT_ARTIFACT_JUMP, FT_ARTIFACT_MUSCLE, FT_ARTIFACT_THRESHOLD, FT_ARTIFACT_ZVALUE

% Copyright (C) 2003-2011, Jan-Mathijs Schoffelen & 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 provenance is not needed because just a call to ft_artifact_zvalue
% ft_preamble loadvar data is not needed because ft_artifact_zvalue will do this

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed',    {'datatype', 'continuous'});
cfg = ft_checkconfig(cfg, 'renamedval', {'continuous', 'continuous', 'yes'});
cfg = ft_checkconfig(cfg, 'required', 'method');
cfg = ft_checkconfig(cfg, 'allowedval', {'method', 'detect', 'marker'});

% set default rejection parameters
if ~isfield(cfg, 'artfctdef'),                       cfg.artfctdef                   = [];        end
if ~isfield(cfg.artfctdef, 'tms'),                   cfg.artfctdef.tms               = [];        end
if ~isfield(cfg, 'method'),                          cfg.method                      = 'detect';  end
if ~isfield(cfg, 'prestim'),                         cfg.prestim                     = 0.005;     end
if ~isfield(cfg, 'poststim'),                        cfg.poststim                    = 0.010;     end

if isfield(cfg.artfctdef.tms, 'artifact')
  fprintf('tms artifact detection has already been done, retaining artifacts\n');
  artifact = cfg.artfctdef.tms.artifact;
  return
end

switch cfg.method
  case 'detect'
    % settings for preprocessing
    if ~isfield(cfg.artfctdef.tms, 'derivative'),   cfg.artfctdef.tms.derivative   = 'yes';    end
    % settings for the zvalue subfunction
    if ~isfield(cfg.artfctdef.tms, 'method'),     cfg.artfctdef.tms.method      = 'zvalue';    end
    if ~isfield(cfg.artfctdef.tms, 'channel'),    cfg.artfctdef.tms.channel     = 'all';       end
    if ~isfield(cfg.artfctdef.tms, 'trlpadding'), cfg.artfctdef.tms.trlpadding  = 0.1;         end
    if ~isfield(cfg.artfctdef.tms, 'fltpadding'), cfg.artfctdef.tms.fltpadding  = 0.1;         end
    if ~isfield(cfg.artfctdef.tms, 'artpadding'), cfg.artfctdef.tms.artpadding  = 0.01;        end
    if ~isfield(cfg.artfctdef.tms, 'cutoff'),     cfg.artfctdef.tms.cutoff      = 4;           end

    % the data is either passed into the function by the user or read from file with cfg.inputfile
    hasdata = exist('data', 'var');

    % construct a temporary configuration that can be passed onto artifact_zvalue
    tmpcfg = keepfields(cfg, {'trl', 'continuous', 'dataset', 'datafile', 'headerfile', 'dataformat', 'headerformat'});
    tmpcfg.artfctdef.zvalue = cfg.artfctdef.tms;

    if hasdata
      fsample = data.fsample; % get the sampling rate from the data structure
      [tmpcfg, artifact] = ft_artifact_zvalue(tmpcfg, data);
    else
      [tmpcfg, artifact] = ft_artifact_zvalue(tmpcfg);
      % FT_ARTIFACT_ZVALUE will call FT_CHECKCONFIG with the dataset2files option
      hdr = ft_read_header(tmpcfg.headerfile, 'headerformat', tmpcfg.headerformat);
      fsample = hdr.Fs; % get the sampling rate from the file on disk
    end
    cfg.artfctdef.tms = tmpcfg.artfctdef.zvalue;

    % adjust artifact definition so that Nx2 matrix contains detected TMS
    % events with user-specified pre- and post stimulus period included.
    % The reason for this is that ft_artifact_zvalue centers the period
    % marked as artifactual around the detected event. In the case of a TMS
    % pulse the window you would like to mark as artifactual is not
    % symmetrical around the onset of the pulse.

    % get values and express in samples
    prestim = round(cfg.prestim * fsample);
    poststim = round(cfg.poststim * fsample);

    % adjust Nx2 artifact matrix to be centered non-symmetrically around
    % detected TMS-pulse
    artifact(:,1) = (artifact(:,1)+artifact(:,2))./2 - prestim;
    artifact(:,2) = artifact(:,1) + poststim;
    cfg.artfctdef.tms.artifact = artifact;

  case 'marker'
    % Check if the cfg is correct for this method
    cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
    ft_checkconfig(cfg, 'required', 'trialdef');
    cfg.trialfun = ft_getopt(cfg, 'trialfun', 'ft_trialfun_general');
    trialdef = cfg.trialdef;
    trialdef.prestim = cfg.prestim;
    trialdef.poststim = cfg.poststim;
    cfg.trialdef = ft_checkconfig(trialdef, 'required', {'eventvalue', 'eventtype'});

    % Get the trialfun
    cfg.trialfun = ft_getuserfun(cfg.trialfun, 'trialfun');

    % Evaluate the trialfun
    fprintf('evaluating trialfunction ''%s''\n', func2str(cfg.trialfun));
    trl   = feval(cfg.trialfun, cfg);

    % Prepare the found events for output
    artifact = trl(:,1:2);
    cfg.artfctdef.tms.artifact = artifact;
    fprintf('found %d events\n', size(artifact,1));

  otherwise
    % This should be redundant as ft_checkconfig does not allow other methods than the supported ones.
    ft_error('unsupported method');
end

% FIXME - not removing this causes problems when passing to ft_preprocessing
cfg = rmfield(cfg, 'method');

author: fieldtrip


text:
function all_bbox_gt = remove_difficult_ground_truth_bboxes(all_bbox_gt)
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

num_imgs = numel(all_bbox_gt);
for img_idx = 1:num_imgs
    if size(all_bbox_gt{img_idx},2) >= 6
        is_easy_mask  = all_bbox_gt{img_idx}(:,6) == 0;
        all_bbox_gt{img_idx} = all_bbox_gt{img_idx}(is_easy_mask,:);
    end
end
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
% findInstrument - Finds an instrument by model name, vendor name or VISA resource
%
%    obj = findInstrument('STR') returns an interface object to an
%       instrument whose model name, vendor name or VISA resource matches
%       any part of the string STR. If there is no unique match, obj is []
%       and a list of available instruments is displayed.
%
%       findInstrument currently only works with USB and GPIB instruments
%       and requires a VISA adaptor.
%
%       STR is first compared to the model and vendor names. If none of
%       these match, STR is then compared to the VISA resource names. The
%       matching ignores the case of the strings (i.e., upper/lower case)
%
%    obj = findInstrument('STR', 'ADAPTOR') uses the specified VISA
%       adaptor. This is only applicable if you have multiple VISA
%       installations (by default, findInstrument picks one of the VISA
%       adaptors).
%
%    Examples:
%      % see a list of all the available instruments
%      findInstrument;
%
%      % find an instrument by model number
%      obj = findInstrument('34405');
%      fopen(obj)
%      query(obj,'*IDN?')
%      fclose(obj)
%
%      % find an instrument by model number, using Agilent VISA
%      obj = findInstrument('34405', 'agilent');
%
%      % find an instrument by matching a manufacturer name
%      obj = findInstrument('tek');
%
%      % find an instrument by the model code in the VISA resource
%      obj = findInstrument('0x0618')
%
%    If your instrument doesn't show up with findInstrument, use
%    this Troubleshooting checklist
%
%   Related Functions:
%     instrfindall (Find existing interface objects)
%     instrreset   (Close and delete all interface objects)
%     obj2mfile    (Generate MATLAB code for re-creating an interface object)

% Copyright 2009, The MathWorks, Inc.

function out = findInstrument(searchstr, adaptor)

if exist('instrhwinfo', 'file') ~= 2
    error('This function requires Instrument Control Toolbox');
end

if ~exist('searchstr', 'var')
    searchstr = '';
end
searchstr = strtrim(lower(searchstr));

if strcmp(searchstr, '##troubleshoot')
    helpstr = troubleshootHelp;
    disp(helpstr);
    return;
end

visainfo = instrhwinfo('visa');

if numel(visainfo.InstalledAdaptors) == 0
    error('No VISA adaptors found\n');
end

if ~exist('adaptor', 'var')
    adaptor = visainfo.InstalledAdaptors{end};
end

try
    instr = instrhwinfo('visa', adaptor);
catch %#ok<CTCH>
    error('VISA Adaptor ''%s'' is not installed on this system', adaptor);
end

out = [];

%% Scan over all the instruments returned by the VISA adaptor
% For each instrument: match its VISA resource; try to get its
% ID string; match its vendor or model name

info = []; cnt = 0;
for j=1:numel(instr.ObjectConstructorName)
    constructor = instr.ObjectConstructorName{j};

    res = regexpi(constructor, '[a-z]*\(''[a-z]+''[ ]*,[ ]*''(?<resource>[A-Z0-9:]+)''', 'names');
    if (numel(res) == 0) || strcmpi(res.resource(1:4),'ASRL')
        % serial port, ignore
        continue;
    end

    cnt = cnt + 1;
    info(cnt).resource = res.resource; %#ok<*AGROW>
    info(cnt).constructor = constructor;
    info(cnt).resourceMatched = ~isempty(findstr(searchstr, lower(res.resource)));
    info(cnt).idMatched = false;

    obj = eval(constructor);

    % visa resource strings have "::". For ICT objects, this is replaced
    % by "-". Get the ICT resource name so that we can search for other
    % interface objects for the same resource
    info(cnt).ICT_resource = get(obj, 'RsrcName');

    try
        fopen(obj);
        idString = query(obj, '*IDN?');
        fclose(obj);
        delete(obj);
        info(cnt).connectStatus = true;
    catch %#ok<CTCH>
        % open failed, probably because the connection is already
        % open from within MATLAB or another program
        delete(obj);
        info(cnt).connectStatus = false;
        continue;
    end

    % idString will have the following format
    %   <manufacturer>, <model name>, <serial number>, <firmware revision>
    % we only care about manufacturer and model name
    strs = regexp(idString, '[ ]*,[ ]*', 'split');
    if numel(strs) < 4
        strs{4} = '';
    end
    id.vendor = strs{1};
    id.model = strs{2};

    modelMatched = ~isempty(findstr(searchstr, lower(id.model)));
    vendorMatched = ~isempty(findstr(searchstr, lower(id.vendor)));
    info(cnt).idMatched = modelMatched || vendorMatched;

    info(cnt).id = id;
end % loop over each instrument

if numel(info) == 0
    out = [];
    fprintf('No USB or GPIB instruments found using ''%s'' VISA\n', adaptor);
    return;
end

%% Collate the match information.
% Only create an interface object if there is a unique match

idMatches = find( [info.idMatched] );
numIdMatches = numel(idMatches);

resourceMatches = find( [info.resourceMatched] );
numResourceMatches = numel(resourceMatches);

if numIdMatches == 1
    % A match to vendor or model supersedes any VISA resource matches
    % i.e., from the user's perspective, we'll only examine resource
    % strings if there is no vendor or model match
    idx = idMatches(1);
    if ~isempty(instrfindall('RsrcName', info(idx).ICT_resource))
        warning('findInstrument:MultipleInterfaceObjects', ...
            'Preexisting interface objects for ''%s''', info(idx).resource);
    end
    out = eval(info(idx).constructor);

elseif numIdMatches > 1
    out = [];
    fprintf('%d instruments matched ''%s''\n\n', numIdMatches, searchstr);
    showInstrumentInfo(info, idMatches);
    noConn = find([info.connectStatus] == false);

    if numel(noConn) > 0
        fprintf('%d instrument(s) could not be reached\n\n', numel(noConn));
        showInstrumentInfo(info, noConn);
    end

elseif numIdMatches == 0 && numResourceMatches == 1
    idx = resourceMatches(1);
    if info(idx).connectStatus == false
        out = [];
        fprintf('One instrument matched ''%s'' (VISA resource), but it cannot be reached\n\n', searchstr);
        showInstrumentInfo(info, idx);
        return;
    end
    if ~isempty(instrfindall('RsrcName', info(idx).ICT_resource))
        warning('findInstrument:MultipleInterfaceObjects', ...
            'Preexisting interface objects for ''%s''', info(idx).resource);
    end
    out = eval(info(idx).constructor);

elseif numIdMatches == 0 && numResourceMatches > 0
    out = [];
    fprintf('%d instrument(s) matched ''%s'' (VISA Resource)\n\n', numResourceMatches, searchstr);
    showInstrumentInfo(info, resourceMatches);

elseif numIdMatches == 0 && numResourceMatches == 0
    out = [];
    fprintf('No instruments matched ''%s''. Available instruments:\n\n', searchstr);
    showInstrumentInfo(info, 1:numel(info));

end

end

%%
    function showInstrumentInfo(info, indexList)
        for i=1:numel(indexList)
            idx = indexList(i);
            fprintf('  VISA Resource: %s\n', info(idx).resource);
            if info(idx).connectStatus == false
                fprintf('         Unable to connect (instrument may be in use)\n');
            else
                id = info(idx).id;
                fprintf('         Vendor: %s\n', id.vendor);
                fprintf('          Model: %s\n', id.model);
            end
            fprintf('\n');
        end
    end

 %%
    function out = troubleshootHelp
s = {
'    Troubleshooting:'
'      If your instrument doesn''t show up with findInstrument, use the '
'      following checklist:'
''
'      1) Is the instrument turned on and connected to your computer?'
''
'      2) Check if you have multiple VISAs installed, by typing:'
'            instrhwinfo(''visa'')'
'         If you do have multiple VISAs, explicitly specify the VISA adaptor'
'            obj = findInstrument(''34405'', ''agilent'');'
'            obj = findInstrument(''34405'', ''ni'');'
''
'      3) Does any other program have an open connection to the instrument?'
'         If yes, close that program and try again.'
''
'      4) Do you have an open connection to the instrument from MATLAB?'
'         You can verify this by typing ''instrfindall'' at the command line.'
'         Type ''instrreset'' to close and delete the existing connections'
''
'      5) If the problem persists,'
'         - type ''instrreset'' at the MATLAB command line'
'         - unplug the instrument from the PC,'
'         - wait a few seconds, and'
'         - plug the instrument back in'
};
out = char(s);
    end

% Format of a VISA resource string
% GPIB[board]::primary address[::GPIB secondary address][::INSTR]
% USB[board]::manufacturer-id::model-id::serial #[::USB interface #][::INSTR]
%
% Examples
% USB0::0x0957::0x0588::CN49140544::INSTR
% USB0::0x0957::0x0618::TW46000285::INSTR
% GPIB0::24::INSTR

% Format of an identification string returned by an instrument
% <manufacturer>, <model name>, <serial number>, <firmware revision>
%
% -- Examples
% Agilent Technologies,33210A,MY48000278,1.03-1.02-26-2
% Agilent Technologies,34405A,TW46000285,1.40-3.08
% Agilent Technologies,DSO1014A,CN49140544,00.04.02
% HEWLETT-PACKARD,E3631A,0,2.1-5.0-1.0
% HEWLETT-PACKARD,33120A,0,2.0-2.0-1.0
% TEKTRONIX,TDS 210,0,CF:91.1CT FV:v2.03 TDS2MM:MMV:v1.04
% TEKTRONIX,TDS 3034B,0,CF:91.1CT FV:v3.41 TDS3FFT:v1.00 TDS3TRG:v 1.00
% TEKTRONIX,TDS 3034C,0,CF:91.1CT FV:v4.01 TDS3FFT:v1.00 TDS3TRG:v1.00
% KEITHLEY INSTRUMENTS INC.,MODEL 2701,1216866,A09 /A02
% *IDN LECROY,9354AL,935408135,07.3.0
% *IDN LECROY, WM8300, WAVEMASTER00112,0.5.0 n

author: Sable

repo: mcbench-benchmarks


text:
function plotSelectedFeaturesMap(exp_name, dataOutputDir, featMap, featChnMaps, nFilters)

% figure;
% imagesc(featMap);
% title('Full Feature Map');
% saveas(gcf, [exp_name '_FeatMap.eps'], 'epsc');
% saveas(gcf, [exp_name '_FeatMap.png'], 'png');
%
% for f=1:nFilters
%   figure;
%   for i=1:length(featChnMaps)/nFilters
%     subplot(4,3,i);
%     imagesc(featChnMaps{i*f});
%     title(sprintf('Feature %d', i));
%   end
%   saveas(gcf,sprintf([exp_name '_PerFeatMap_%d.eps'], f), 'epsc');
%   saveas(gcf,sprintf([exp_name '_PerFeatMap_%d.png'], f), 'png');
% end

figure;
imagesc(featMap);
colormap('jet');
title('Full Feature Map');
saveas(gcf, fullfile(dataOutputDir, [exp_name '_FeatMap.eps']), 'epsc');
saveas(gcf, fullfile(dataOutputDir,[exp_name '_FeatMap.png']), 'png');
img = featMap;
img = (img - min(min(img))) ./ (max(max(img)) - min(min(img)));
ime = imresize(img, 100., 'nearest');
img = ind2rgb(gray2ind(img, 255), jet(255));
imwrite(img, fullfile(dataOutputDir, [exp_name '_FeatMap_FULL_IMG.png']), 'png');

figure;
NUM_CHANNELS = 10;
NUM_FILTERS = length(featChnMaps)/NUM_CHANNELS;
images = zeros(size(featMap,1), size(featMap,2)*NUM_CHANNELS);
for i=1:NUM_CHANNELS
  index = 1+(i-1)*size(featMap,2);
  for j=1:NUM_FILTERS
    images(:, index:index+size(featMap,2)-1) = images(:, index:index+size(featMap,2)-1) + featChnMaps{i+(j-1)*NUM_CHANNELS};
  end
end
images = (images - min(min(images))) ./ (max(max(images)) - min(min(images)));
images2 = imresize(images, 20., 'nearest');
images2 = ind2rgb(gray2ind(images2, 255), jet(255));
imshow(images2, 'Border', 'tight');
saveas(gcf, fullfile(dataOutputDir, [exp_name '_PerFeatMap.eps']), 'epsc');
saveas(gcf, fullfile(dataOutputDir, [exp_name '_PerFeatMap.png']), 'png');

for i=1:NUM_CHANNELS
  index = 1+(i-1)*size(featMap,2);
  img = images(:, index:index+size(featMap,2)-1);
  img = imresize(img, 10., 'nearest');
  img = ind2rgb(gray2ind(img, 255), jet(255));
  imwrite(img, fullfile(dataOutputDir, sprintf([exp_name '_PerFeatMap_%d.png'], i)), 'png');
end

figure;
c = 1;
for i=1:NUM_CHANNELS
  for j=1:NUM_FILTERS
    img = featChnMaps{c};
    img = (img - min(min(img))) ./ (max(max(img)) - min(min(img)));
    img = imresize(img, 10., 'nearest');
    img = ind2rgb(gray2ind(img, 255), jet(255));
    imwrite(img, fullfile(dataOutputDir, sprintf([exp_name '_PerFeatMap_filter_%d_channel_%d.png'], i, j)), 'png');
    imshow(img);
    %pause;
    c = c + 1;
  end
end

author: MengyangPu

repo: EDTER

save_path: github-repos/MATLAB/MengyangPu-EDTER


text:
function Show_EWT2D_Filters(fil)

%===========================================================
%
% function Show_EWT2D_Filters(fil)
%
% This function permits to plot the Fourier magnitude of all
% 2D EWT Littlewood-Paley filters.
%
% Input:
%   - fil: cell containing all the filters
%
% Author: J.Gilles
% Institution: UCLA - Department of Mathematics
% email: jegilles@math.ucla.edu
% Date: March, 1st, 2013
%
%===========================================================
p=ceil(length(fil)/2);

figure;
for n=1:length(fil)
    subplot(2,p,n);imshow(fftshift(abs(fil{n})),[]);
end


text:
%SETDATA Set data field in mapping or dataset
%
% 	W = SETDATA(W,DATA,FIELD)
% 	W = SETDATA(W,DATA,N)
%
%   A = SETDATA(A,DATA,FEATLAB)
%
% INPUT
%   W        Mapping
%   DATA     Data to be put in the data field.
%   FIELD    String, name of data field structure to be used for DATA
%            (optional)
%   N        Index of cell array to be used for DATA (optional)
%
%   A        Dataset
%   FEATLAB  Desired feature labels for dataset (optional)
%
% OUTPUT
%   W        Mapping
%   A        Dataset
%
% DESCRIPTION
% This routine can be used to store data in the data fields of either a
% mapping W or a dataset A.
%
% The data field of a mapping consists of a matrix, cell array or a structure.
% In case nor FIELD neither N are given the entire dats field is replaced by
% DATA. In case the data field is a structure DATA is assigned to field FIELD.
% In case the datafield is a cell array, DATA is assigned to cell N.
% DATA cannot be a structure in case of untrained or fixed mappings.
%
% The data field of a dataset A is an array of doubles. It is replaced by
% DATA which may be doubles or another dataset, which will be converted to
% doubles. The numbers of objects in A and DATA should be equal. The feature
% labels of A are replaced by FEATLAB, or, if not supplied and DATA is a
% dataset by the feature labels of DATA. Labels and class probabilities of
% A are preserved.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS

author: marianux

repo: ecg-kit

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


text:
function varargout=mu_util(optn,varargin)
% MU_UTIL Various utility routines
%           This function should not be used directly; instead it is
%           is accessed by other high- and low-level functions.

% MU_UTIL is basically a driver for a number of lower-level routines collected
% together for convenience. They are specified by the first argument:
% 'clip'      - clipping routine
% 'axisticks' - generates a "nice" set of ticks, optimized for degrees/minutes on maps
%               angles of the circle.
% 'x/ygrid'   - generates the lines for axis grids
% 'xylimits'  - finds the x/y limits for a given map
% 'lllimits'  - finds the lat/long limits for a given map (usually for rectboxes)
% 'box'       - returns a line around the boundaries of the map.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 4/April/97
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

switch optn
  case 'clip'
    [varargout{1},varargout{2}]=m_clip(varargin{:});
  case 'axisticks'
    varargout{1}=m_getinterval(varargin{:});
  case {'xgrid','ygrid'}
    [varargout{1},varargout{2},varargout{3},varargout{4}]=m_rectgrid(optn,varargin{:});
  case 'xylimits'
    m_getxylimits;
  case 'lllimits'
    m_getlllimits;
  case 'box'
    [varargout{1},varargout{2}]=m_box(varargin{:});

end

%---------------------------------------------------------
function [Xc,Yc]=m_clip(cliptype,X,Xedge,indx,Y)
% M_CLIP performs clipping of data. Columns of points are
%        assumed to be lines; the first points outside the
%        clip area are recomputed to lie on the edge of the
%        region; others are converted to either NaN or
%        edge points depending on CLIPTYPE.
%
%        indx = 0 for points inside the clip region, 1
%               for those outside
%
%         'on'  - replaces points outside with NaN, but interpolates
%                 to provide points right on the border.
%         'patch' - replaces points outside with nearest border point
%         'point' - does no interpolation (just checks in/out)
%
% In general m_clip will be called 4 times, since it solves a 1-edge problem.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 4/April/97

Xc=X;
Yc=Y;

if ~strcmp(cliptype,'point')

  % Find regions where we suddenly come into the area
  % (indices go from 1 to 0)

  [i,j]=find(diff(indx)==-1);

  if any(i)
    I=i+(j-1)*size(X,1); % 1-d addressing

    % Linearly interpolate to boundary

    bt=(X(I+1)-X(I));
    ibt=abs(bt)<5*eps;
    if any(ibt), bt(ibt)=1*eps; end % In these cases the delta(Y) also = 0, so we just want
                                      % to avoid /0 warnings.
    Yc(I)=Y(I)+(Xedge-X(I)).*(Y(I+1)-Y(I))./bt;
    Yc(I(ibt))=(Y(I(ibt))+Y(I(ibt)+1))/2;
    Xc(I)=Xedge;
    indx(I(isfinite(Yc(I))))=0;
  end

  % Find regions where we suddenly come out of the area
  % (indices go from 0 to 1)

  [i,j]=find(diff(indx)==1);
  if any(i)

    I=i+(j-1)*size(X,1);

    bt=(X(I+1)-X(I));
    ibt=abs(bt)<5*eps;
    if any(ibt), bt(ibt)=eps;  end % In these cases the delta(Y) also = 0, so we just want
                                      % to avoid /0 warnings.

    Yc(I+1)=Y(I)+(Xedge-X(I)).*(Y(I+1)-Y(I))./bt;
    Yc(I(ibt)+1)=(Y(I(ibt))+Y(I(ibt)+1))/2;
    Xc(I+1)=Xedge;
    indx(I(isfinite(Yc(I+1)))+1)=0;
  end

end

switch cliptype
  case {'on','point'}
    Xc(indx)=NaN;
    Yc(indx)=NaN;
  case 'patch'
    Xc(indx)=Xedge;
end

%--------------------------------------------------------------------------
function gval=m_getinterval(gmin,gmax,gtick)
% M_GETINTERVAL picks nice spacing for grid ticks
%        This occurs when the following call is made:
%        TICKS=M_GRID('axisticks',MIN,MAX,APPROX_NUM_TICKS)
%

% Rich Pawlowicz (rich@ocgy.ubc.ca) 2/Apr/1997
%
% 9/Apr/98 - changed things so that max/min limits are not automatically
%            added (this feature made map corners messy sometimes)

% If ticks are specified, we just make sure they are within the limits
% of the map.

if length(gtick)>1
  gval=[gtick(gtick(:)>=gmin & gtick(:)<=gmax)];

% Otherwise, we try to fit approximately gtick ticks in the interval
else

  if gtick>2

    exactint=(gmax-gmin)/(gtick-1)*60; %interval in minutes

    % These are the intervals which we will allow (they are "nice" in the sense
    % that they come to various even multiples of minutes or degrees)
    niceints=[0.1 0.2 0.25 0.5 ...
              1 2 3 4 5 6 10 12 15 20 30 ...
              60*[1 2 3 4 5 6 8 9 10 12 15 18 20 25 30 40 50 60 100 120 180]];

    [dun,I]=min(abs(niceints-exactint));

    gval=niceints(I)/60*[ceil(gmin*60/niceints(I)):fix(gmax*60/niceints(I))];

    gval=[gval(gval>=gmin & gval<=gmax) ];
  else
    gval=[gmin gmax];
  end
end

%--------------------------------------------------------------
function [X,Y,vals,labI]=m_rectgrid(direc,Xlims,Ylims,Nx,Ny,label_pos)
% M_RECTGRID This handles some of the computations involved in creating grids
%            for rectangular maps. Essentially we make our "first guess" using the
%            lat/long limits. Then these curves are clipped to the boundaries, after
%            which we use the clip points as "new" boundaries and recompute the lines.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 4/April/97

global MAP_PROJECTION MAP_VAR_LIST

Ny1=Ny-1;
Ny2=Ny*2;
Ny21=Ny2-1;

% First try some wildly oversampled lines (not including the boundaries)

vals=mu_util('axisticks',Xlims(1),Xlims(2),Nx);

if strcmp(MAP_VAR_LIST.rectbox,'on') |  strcmp(MAP_VAR_LIST.rectbox,'circle')

 % We don't want the end limits here.
 if vals(end) == Xlims(2), vals(end) = []; end
 if vals(1)   == Xlims(1), vals(1)   = []; end

 if direc(1)=='x'
  [lg,lt]=meshgrid(vals,Ylims(1)+diff(Ylims)*[0:1/Ny1:1]);
 else
  [lt,lg]=meshgrid(vals,Ylims(1)+diff(Ylims)*[0:1/Ny1:1]);
 end

 % But sneakily we clip them in transforming, so we end up with isfinite values only
 % inside the axis limits

 [X,Y]=feval(MAP_PROJECTION.routine,'ll2xy',lg,lt,'clip','on');

 % Now we find the first/last unclipped values; these will be our correct starting points.
 % (Note I am converting to one-dimensional addressing).

 istart=sum(cumsum(isfinite(X))==0)+1+[0:size(X,2)-1]*size(X,1);
 iend=size(X,1)-sum(cumsum(isfinite(flipud(X)))==0)+[0:size(X,2)-1]*size(X,1);

 % Now go back and find the lat/longs corresponding to those points; these are our new
 % starting points for the lines (Note that the linear interpolation for clipping at boundaries
 % means that they will not *quite* be the exact longitudes due to curvature, but they should
 % be very close.

 if direc(1)=='x'

  [lgs,lts]=feval(MAP_PROJECTION.routine,'xy2ll',X(istart),Y(istart),'clip','off');
  [lgs,lte]=feval(MAP_PROJECTION.routine,'xy2ll',X(iend),Y(iend),'clip','off');

  % Finally compute the lines within those limits (these *may* include some out-of-bounds points
  % depending on the geometry of the situation; these are converted to NaN as usual.

  [X,Y]=feval(MAP_PROJECTION.routine,'ll2xy',vals(ones(Ny2,1),:),...
              lts(ones(Ny2,1),:)+[0:1/Ny21:1]'*(lte-lts),'clip','on');

 else

  [lgs,lts]=feval(MAP_PROJECTION.routine,'xy2ll',X(istart),Y(istart),'clip','off');
  [lge,lts]=feval(MAP_PROJECTION.routine,'xy2ll',X(iend),Y(iend),'clip','off');

  % Longitudes should be increasing here, but we can run into wrap problems after
  % the tansformations back and forth. It is for this line that I need to make
  % long-lims just a tad less than 360 when really they should be 360 (in m_lllimits)

  lge(lge<=lgs)=lge(lge<=lgs)+360;

  [X,Y]=feval(MAP_PROJECTION.routine,'ll2xy',lgs(ones(Ny2,1),:)+[0:1/Ny21:1]'*(lge-lgs),...
              vals(ones(Ny2,1),:),'clip','on');

 end

else
 if direc(1)=='x'
  [lg,lt]=meshgrid(vals,Ylims(1)+diff(Ylims)*[0:1/Ny1:1]);
 else
  [lt,lg]=meshgrid(vals,Ylims(1)+diff(Ylims)*[0:1/Ny1:1]);
 end
 [X,Y]=feval(MAP_PROJECTION.routine,'ll2xy',lg,lt,'clip','off');

end

switch label_pos
  case {'left','bottom','west','south'}
    labI=1;
  case 'middle'
    labI=round(size(X,1)/2+1/2);
  case {'right','top','east','north'}
    labI=size(X,1);
end

%--------------------------------------------------------------------------------
function m_getxylimits
% M_GET_LIMITS Converts X/Y limits to lat/long limits
%              This is a chunk of code that is needed for most projections.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 4/April/97

global MAP_PROJECTION MAP_VAR_LIST

% Start with the user-specified lat/longs.

MAP_VAR_LIST.lats=MAP_VAR_LIST.ulats;
MAP_VAR_LIST.longs=MAP_VAR_LIST.ulongs;

% Now, let's get the map x/ylims
bX=MAP_VAR_LIST.longs(1)+diff(MAP_VAR_LIST.longs)*[0:1/30:1];
bY=MAP_VAR_LIST.lats(1)+diff(MAP_VAR_LIST.lats)*[0:1/30:1];
bX=[bX MAP_VAR_LIST.longs(2*ones(1,31)) fliplr(bX) MAP_VAR_LIST.longs(ones(1,31)) ];
bY=[MAP_VAR_LIST.lats(ones(1,31)) bY MAP_VAR_LIST.lats(2*ones(1,31)) fliplr(bY) ];
[X,Y]=feval(MAP_PROJECTION.routine,'ll2xy',bX,bY,'clip','off');
MAP_VAR_LIST.xlims=[min(X) max(X)];
MAP_VAR_LIST.ylims=[min(Y) max(Y)];

%-------------------------------------------------------------------------------
function m_getlllimits
% M_GET_LIMITS Converts X/Y limits to lat/long limits
%              This is a chunk of code that is needed for most projections.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 4/April/97

global MAP_PROJECTION MAP_VAR_LIST

[bX,bY]=mu_util('box',31);
% Get its lat/longs.

[lg,lt]=feval(MAP_PROJECTION.routine,'xy2ll',bX,bY,'clip','off');
% Take real part because otherwise funny things might happen if the box is very large

MAP_VAR_LIST.lats=[min(real(lt)) max(real(lt))];

% Are the poles within the axis limits? (Test necessary for oblique mercator and azimuthal)
[px,py]=feval(MAP_PROJECTION.routine,'ll2xy',[0 0],[-90 90],'clip','point');
if isfinite(px(1)), MAP_VAR_LIST.lats(1)=-90; end
if isfinite(px(2)), MAP_VAR_LIST.lats(2)= 90; end

if any(isfinite(px))
  MAP_VAR_LIST.longs=[-179.9 180]+exp(1); % we add a weird number (exp(1)) to get away from
                         % anything that might conceivably be desired as a
                         % boundary - it makes grid generation easier.
                         % Also make the limits just a little less than 180, this
                         % is necessary because I have to have the first and last points
                         % of lines just a little different in order to figure out orientation
                         % in 'm_rectgrid'
else
  MAP_VAR_LIST.longs=[min(lg) max(lg)];
  if all(isnan(px)) & diff(MAP_VAR_LIST.longs)>360*30/31
    ii=lg<mean(MAP_VAR_LIST.longs);
    lg(ii)=lg(ii)+360;
    MAP_VAR_LIST.longs=[min(lg) max(lg)];
  end
end

%------------------------------------------------------------------------
function [X,Y]=m_box(npts)
% M_BOX  Computes coordinates of the map border.
%

% Rich Pawlowicz (rich@ocgy.ubc.ca) 4/April/97

global MAP_PROJECTION MAP_VAR_LIST

n1=npts-1;

switch MAP_VAR_LIST.rectbox
  case 'on'
    X=MAP_VAR_LIST.xlims(1)+diff(MAP_VAR_LIST.xlims)*[0:1/n1:1];
    Y=MAP_VAR_LIST.ylims(1)+diff(MAP_VAR_LIST.ylims)*[0:1/n1:1];
    X=[X MAP_VAR_LIST.xlims(2*ones(1,npts)) fliplr(X) MAP_VAR_LIST.xlims(ones(1,npts))];
    Y=[MAP_VAR_LIST.ylims(ones(1,npts)) Y  MAP_VAR_LIST.ylims(2*ones(1,npts)) fliplr(Y)];
  case 'off'
    lg=MAP_VAR_LIST.longs(1)+diff(MAP_VAR_LIST.longs)*[0:1/n1:1];
    lg=[lg MAP_VAR_LIST.longs(2*ones(1,npts)) fliplr(lg) MAP_VAR_LIST.longs(ones(1,npts))]';
    lt=MAP_VAR_LIST.lats(1)+diff(MAP_VAR_LIST.lats)*[0:1/n1:1];
    lt=[MAP_VAR_LIST.lats(ones(1,npts)) lt  MAP_VAR_LIST.lats(2*ones(1,npts)) fliplr(lt)]';
    [X,Y]=feval(MAP_PROJECTION.routine,'ll2xy',lg,lt,'clip','off');
  case 'circle'
    n1=npts*3-1;
    X=MAP_VAR_LIST.rhomax*cos([0:n1]/n1*pi*2);
    Y=MAP_VAR_LIST.rhomax*sin([0:n1]/n1*pi*2);
end

author: CelsoReyes


text:
function UPO = partial_order_polarized_channels(N)
%N is code length, power of 2. UPO = Univeral partial order
m = log2(N);
max_connection = 32;
UPO = zeros(N, max_connection);
UPO(1, 1) = 1;
% Partial order construction
for layer = 2 : m
    N_tmp = 2^layer;
    UPO_tmp = UPO(1 : N_tmp/2, :);
    for i = 1 : N_tmp/2
        for j = 1 : max_connection
            if UPO_tmp(i, j) == 0
                break;
            else
                UPO_tmp(i, j) = UPO_tmp(i, j) + N_tmp/2;
            end
        end
    end
    UPO(N_tmp/2 + 1 : N_tmp, :) = UPO_tmp;
    for i = N_tmp/4 + 1 : N_tmp/2
        for j = 1 : max_connection
            if UPO(i, j) == 0
                UPO(i, j) = (i - 1) + N_tmp/4;
                break;
            end
        end
    end
end
%Delete redundant all zero columns
while(all(UPO(:, end) == 0))
    UPO(:, end) = [];
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function needed = isneeded(curmethod,allmethods)
% check whether curmethod will need to be calculated

if iscell(curmethod) && length(curmethod)>1
    for i=1:length(curmethod)
        needed(i) = isneeded(curmethod{i},allmethods);
    end
    return;
end

% depedency graph -- each of the measures on the right depend on the
% measure on the left (dependency arrow goes from left to right
% e.g., 'Coh',{'iCoh'} means iCoh depends on Coh (Coh-->iCoh)
% If a new measure is added, make sure to update this depedency table
% (add the method descriptor to the appropriate RHS and make a LHS entry
% for the method)
dependencies = {...
                'dDTF',     {};
                'ffDTF2',   {'dDTF'};
                'ffDTF',    {};
                'nDTF',     {};
                'GGC',      {};
                'iCoh',     {};
                'Coh',      {'iCoh'};
                'S',        {'Coh','GGC'};
                'pCoh',     {'dDTF'};
                'mCoh',     {};
                'RPDC',     {};
                'GPDC',     {};
                'nPDC',     {};
                'DTF',      {'nDTF','ffDTF','ffDTF2'};
                'G',        {'mCoh','pCoh'};
                'PDC',      {'DTF','G'};
                'Rinv',     {'RPDC'}};

if isempty(curmethod)
    needed =0;
    return;
elseif ismember_bc(curmethod,allmethods)
    needed = 1;
    return;
else
    % recursively determine whether any of the other methods the user
    % wants depend on curmethod
    idx=strmatch(curmethod,dependencies(:,1),'exact');
    if isempty(idx) || isempty(dependencies(idx,2)), needed=0; return; end
    children=dependencies(idx,2);
    if any(isneeded(children{1},allmethods))
        needed = 1;
        return;
    end
end

% if we get here, this particular curmethod is not needed
needed = 0;

repo: REST


text:
function [Statistics,significantFeatures] = performStatisticalAnalysis(sampleData,sampleInformation,varargin)
% This function determines if there is a significant difference between
% features computed for two or more groups in a cohort of samples. If the
% cohort contains two groups, the Wilcoxon rank sum test is used. If the
% cohort contains three or more groups, the Kruskal Wallis test is used.
%
% USAGE
% [Statistics,significantFeatures] = performStatisticalAnalysis(sampleData,sampleInformation,stratification)
%
% INPUTS
% sampleData           Table with input data to analyze (e.g., fluxes) with
%                      computed features as rows and sample IDs as columns
% sampleInformation    Table with information on analyzed samples including
%                      group classification with sample IDs as rows
%
% OPTIONAL INPUT
% stratification       Column header containing the desired group
%                      classification in sampleInformation table. If not
%                      provided, the second column will be used.
% groupTest            Decides whether Kruskal-Wallis test(default) or
%                      ANOVA should be used for group comparisons.
%                      Allowed inputs: "Kruskal-Wallis","ANOVA"
%
% OUTPUTS
% Statistics           Table with results of statistical tests for each
%                      computed feature
% significantFeatures  Table with input data reduced to only features that
%                      were statistically significant
%
% AUTHOR
%       - Almut Heinken, 12/2020

parser = inputParser();
parser.addRequired('sampleData', @iscell);
parser.addRequired('sampleInformation', @iscell);
parser.addParameter('stratification', '', @ischar);
parser.addParameter('groupTest', 'Kruskal-Wallis', @ischar);

parser.parse(sampleData, sampleInformation, varargin{:});

sampleData = parser.Results.sampleData;
sampleInformation = parser.Results.sampleInformation;
stratification = parser.Results.stratification;
groupTest = parser.Results.groupTest;

if ~any(strcmp(groupTest,{'Kruskal-Wallis','ANOVA'}))
    error('Wrong input for group test!')
end

% find the column with the sample information to analyze the samples by
if ~isempty(stratification)
stratCol=find(strcmp(sampleInformation(1,:),stratification));
else
    stratCol=2;
end

% delete empty columns in the data
sampleData{1,1}='Averages';
delIndices =cellfun(@isempty, sampleData(1,:));
sampleData(:,delIndices)=[];
% delete columns that are all zeros
cnt=1;
delArray=[];
for i=2:size(sampleData,2)
    if abs(sum(str2double(sampleData(2:end,i))))<0.0000001
        delArray(1,cnt)=i;
        cnt=cnt+1;
    end
end
sampleData(:,delArray)=[];

% delete metadata entries not in the sample data
[C,IA]=setdiff(sampleInformation(:,1),sampleData(1,:),'stable');
sampleInformation(IA(2:end),:)=[];

groups=unique(sampleInformation(2:end,stratCol));

if length(groups) > 1

% Fill out table header
if length(groups)==2
    Statistics={'Feature','Description','p_value_before_FDR_Corr','p_value_after_FDR_Corr','Decision','Rank sum statistic','Z-statistics'};
elseif length(groups)>2
    if strcmp(groupTest,'Kruskal-Wallis')
        Statistics={'Feature','Description','p_value_before_FDR_Corr','p_value_after_FDR_Corr','Decision','Degrees of freedom','Chi-sq'};
    elseif strcmp(groupTest,'ANOVA')
        Statistics={'Feature','Description','p_value_before_FDR_Corr','p_value_after_FDR_Corr','Decision','Degrees of freedom','Sum of squares'};
    end
end
cnt=size(Statistics,2)+1;
for i=1:length(groups)
    Statistics{1,cnt}=['Average_',groups{i}];
    cnt=cnt+1;
    Statistics{1,cnt}=['StandardDeviation_',groups{i}];
    cnt=cnt+1;
end

% test if sample information and sample names match
C=intersect(sampleInformation(:,1),sampleData(1,:));
if isempty(C)
    error('Sample IDs are not present as column headers of the sample data table. Consider transposing sample data table.')
end

for j=2:size(sampleData,2)
    if ~isempty(find(strcmp(sampleInformation(:,1),sampleData{1,j})))
        group{j-1,1}=sampleInformation{find(strcmp(sampleInformation(:,1),sampleData{1,j})),stratCol};
    end
end

%% calculate the statistics
for i=2:size(sampleData,1)
    Statistics{i,1}=sampleData{i,1};
    if contains(version,'(R202') % for Matlab R2020a and newer
        dataAll=cell2mat(sampleData(i,2:end));
    else
        dataAll=str2double(sampleData(i,2:end));
    end

    % separate data by group
    for j=1:length(groups)
        dataGrouped{j}=dataAll';
        delInd=find(~strcmp(group,groups{j}));
        dataGrouped{j}(delInd,:)=[];
    end

    if length(groups)==2
        % use Wilcoxon rank sum test

        [p,h,stats] = ranksum(dataGrouped{1},dataGrouped{2},'method','approximate');
        if isnan(p)
            p=1;
        end
        Statistics{i,3}=p;
        Statistics{i,5}=h;
        Statistics{i,6}=stats.ranksum;
        Statistics{i,7}=stats.zval;

    elseif length(groups)>2
        if strcmp(groupTest,'Kruskal-Wallis')
            % use Kruskal Wallis test
            [p,ANOVATAB] = kruskalwallis(dataAll,group,'off');
            Statistics{i,3}=p;
            if ANOVATAB{2,6} ==0
                Statistics{i,5}='0';
            else
                Statistics{i,5}='1';
            end
            Statistics{i,6}=ANOVATAB{2,3};
            Statistics{i,7}=ANOVATAB{2,5};
        elseif strcmp(groupTest,'ANOVA')
            [p,tbl,stats] = anova1(dataAll,group,'off');
            Statistics{i,3}=p;
            if p<0.05
                Statistics{i,5}='1';
            else
                Statistics{i,5}='0';
            end
            Statistics{i,6}=tbl{2,3};
            Statistics{i,7}=tbl{2,2};
        end
    end

    for j=1:length(groups)
        findCol=find(strcmp(Statistics(1,:),['Average_',groups{j}]));
        Statistics{i,findCol}=mean(dataGrouped{j});
        findCol=find(strcmp(Statistics(1,:),['StandardDeviation_',groups{j}]));
        Statistics{i,findCol}=std(dataGrouped{j});
    end
end

%% correct for false discovery (FDR) rate
pAverages=cell2mat(Statistics(2:end,3));
fdr = mafdr(pAverages,'BHFDR', true);
Statistics(2:end,4)=num2cell(fdr);

for i=2:size(Statistics,1)
    if Statistics{i,4} <0.05
        Statistics{i,5} = 1;
    else
        Statistics{i,5} = 0;
    end
end

%% save only the significant entries as a spreadsheet
nsMets={};
cnt=1;
for i=2:size(Statistics,1)
    if Statistics{i,5}==0
        nsMets{cnt,1}=Statistics{i,1};
        cnt=cnt+1;
    end
end
[C,ia,ib] = intersect(sampleData(:,1),nsMets);
significantFeatures=sampleData;
significantFeatures(ia,:)=[];

%% add reaction/metabolite annotations if possible
database=loadVMHDatabase;

for i=2:size(Statistics,1)
    feat=Statistics{i,1};
    if ~any(~contains(Statistics(:,1),'[fe]'))
        feat=strrep(feat,'EX_','');
        feat=strrep(feat,'[fe]','');
    end
    if ~isempty(find(strcmp(database.metabolites(:,1),feat)))
        Statistics{i,2}=database.metabolites{find(strcmp(database.metabolites(:,1),feat)),2};
    elseif ~isempty(find(strcmp(database.reactions(:,1),feat)))
        Statistics{i,2}=database.reactions{find(strcmp(database.reactions(:,1),feat)),2};
    else
        Statistics{i,2}='NA';
    end
end

else
    Statistics = {};
    significantFeatures = {};
end

end

author: opencobra

repo: cobratoolbox


text:
function wms=links2wms(links);
% LINKS2WMS extract WMS Data URLs from OpenSearch links
% Usage: wms=links2wms(links);
%    [links,struc]=opensearch(q);
%    wms=links2wms(links);
% Example WMS GetCapabilities request:
%    url=[wms{i} '?service=WMS&version=1.3.0&request=GetCapabilities'];
% Rich Signell (rsignell@usgs.gov)
k=0;
wms=[];
for i=1:length(links)
    s=char(links{i});
    if ~isempty(regexpi(s,'wms')),
        k=k+1;
        wms{k}=s;
    end
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function [g_all, vv_all]=vifsub_est_M(org,dist, subbands, M);

% uses convolution for determining the parameters of the distortion channel
% Called by vifvec.m

tol = 1e-15; % tolernace for zero variance. Variance below this is set to zero, and zero is set to this value to avoid numerical issues.

for i=1:length(subbands)
    sub=subbands(i);
    y=org{sub};
    yn=dist{sub};

    % compute the size of the window used in the distortion channel estimation
    lev=ceil((sub-1)/6);
    winsize=2^lev+1; offset=(winsize-1)/2;
    win = ones(winsize);

    % force subband size to be multiple of M
    newsize=floor(size(y)./M)*M;
    y=y(1:newsize(1),1:newsize(2));
    yn=yn(1:newsize(1),1:newsize(2));

    % Correlation with downsampling. This is faster than downsampling after
    % computing full correlation.
    winstep=[M M];
    winstart=[1 1].*floor(M/2)+1;
    winstop=size(y)-ceil(M/2)+1;

    % mean
    mean_x = corrDn(y,win/sum(win(:)),'reflect1',winstep, winstart,winstop);
    mean_y = corrDn(yn,win/sum(win(:)),'reflect1',winstep, winstart,winstop);
    % cov
    cov_xy = corrDn(y.*yn, win, 'reflect1',winstep, winstart,winstop) - sum(win(:)).*mean_x.*mean_y;
    % var
    ss_x = corrDn(y.^2,win, 'reflect1',winstep, winstart,winstop) - sum(win(:)).*mean_x.^2;
    ss_y = corrDn(yn.^2,win, 'reflect1',winstep, winstart,winstop) - sum(win(:)).*mean_y.^2;

    % get rid of numerical problems, very small negative numbers, or very
    % small positive numbers, or other theoretical impossibilities.
    ss_x(ss_x<0)=0;
    ss_y(ss_y<0)=0;

    % Regression
    g = cov_xy./(ss_x+tol);

    % Variance of error in regression
    vv = (ss_y - g.*cov_xy)/(sum(win(:)));

    % get rid of numerical problems, very small negative numbers, or very
    % small positive numbers, or other theoretical impossibilities.
    g (ss_x < tol) = 0;
    vv (ss_x < tol) = ss_y (ss_x < tol);
    ss_x(ss_x<tol)=0;

    g (ss_y < tol) = 0;
    vv (ss_y < tol) = 0;

    % constrain g to be non-negative.
    vv(g<0)=ss_y(g<0);
    g(g<0)=0;

    % take care of numerical errors, vv could be very small negative
    vv( vv <= tol) = tol;

    g_all{i}=g;
    vv_all{i}=vv;

end

repo: SelfExSR


function [output] = conv3d(input, w, b, stride, pad)
% input: Ci x D x H x W
% w: Co x Ci x D x H x W
% b: Co
% stride: int
% pad: int

assert(length(size(input)) == 4);
assert(length(size(w)) == 5 || length(size(w)) == 2);
assert(size(input,2) == size(input,3));
assert(size(input,3) == size(input,4));
assert(size(w,2) == size(input,1));

% padding
input = padarray(input, [0, pad, pad, pad]);

S = stride;
if length(size(w)) == 2 % for cccp kernel
    A = 1;
else
    A = size(w,3);
end
C_out = size(w,1);

% create tensor for output
n = floor((size(input,2)-A) / stride + 1);
output = zeros(size(w,1),n,n,n);

% 3d convolution
for i = 1:n
    for j = 1:n
        for k = 1:n
            patch = input(:,(i-1)*S+1:(i-1)*S+A,...
                    (j-1)*S+1:(j-1)*S+A,...
                    (k-1)*S+1:(k-1)*S+A);
              output(:,i,j,k) = sum(repmat(reshape(patch,1,[]),C_out,1).* ...
                  reshape(w,C_out,[]),2) + reshape(b,[],1);
        end
    end
end

end


text:
classdef CelestialBodyData < matlab.mixin.SetGet & dynamicprops
    %CelestialBodyData Summary of this class goes here
    %   Detailed explanation goes here

    properties(Access=private)
        bodies(1,:) KSPTOT_BodyInfo

        bodyIdCacheArr(1,:) double = [];

        ci CelestialBodyIntegration
    end

    properties(Transient,Access=private)
        topLvlBodyCache KSPTOT_BodyInfo
    end

    properties(Constant)
        emptyBodyInfo KSPTOT_BodyInfo = KSPTOT_BodyInfo.empty(1,0);
        globalBaseFrame GlobalBaseInertialFrame = GlobalBaseInertialFrame();
        numIntBlockSize(1,1) double = 1000*86400;
    end

    methods
        function obj = CelestialBodyData(celBodyData)
            if(isstruct(celBodyData))
                bodyNameFields = fields(celBodyData);

                for(i=1:length(bodyNameFields)) %#ok<*NO4LP>
                    bodyName = bodyNameFields{i};

                    if(not(isprop(obj,bodyName)))
                        addprop(obj,bodyName);
                        bodyInfo = celBodyData.(bodyName);

                        obj.(bodyName) = bodyInfo;
                        obj.bodies(end+1) = bodyInfo;
                    else
                        error('The following body already exists in CelestialBodyData when trying to instantiate the class from the celestial body data structure: %s', bodyName);
                    end
                end

            elseif(isa(celBodyData,'CelestialBodyData'))
                obj = celBodyData;

            else
                error('Input to CelestialBodyData must be the structure representation of celestial body data.');
            end

            for(i=1:length(obj.bodies))
                if(isstruct(obj.bodies(i).celBodyData) || isempty(obj.bodies(i).celBodyData))
                    obj.bodies(i).celBodyData = obj;
                end
            end

            obj.bodyIdCacheArr = [obj.bodies.id];

            obj.ci = CelestialBodyIntegration(obj);
            obj.ci.integrateCelestialBodies(0, obj.numIntBlockSize);
        end

        function resetAllParentNeedsUpdateFlags(obj)
            for(i=1:length(obj.bodies))
                obj.bodies(i).parentBodyInfoNeedsUpdate = true;
            end
        end

        function bodyInfo = getBodyInfoById(obj, bodyId)
%             arr = [obj.bodies.id];
            bodyInfo = obj.bodies(obj.bodyIdCacheArr == bodyId);

            if(numel(bodyInfo) > 1)
                bodyInfo = bodyInfo(1);
            end
        end

        function [minUt, maxUt] = getStateCacheMinUtMaxUt(obj)
            bool = [obj.bodies.propTypeEnum] == BodyPropagationTypeEnum.Numerical;
            numericBodies = obj.bodies(bool);

            if(numel(numericBodies) >= 1)
                numericBody = numericBodies(1);
                times = numericBody.numIntStateCache.times;
                minUt = min(times);
                maxUt = max(times);
            else
                minUt = 0;
                maxUt = 0;
            end
        end

        function updateCelestialBodyStateCache(obj, minUT, maxUT, hWaitbar)
            arguments
                obj(1,1) CelestialBodyData
                minUT(1,1) double
                maxUT(1,1) double
                hWaitbar = [];
            end

            obj.ci.integrateCelestialBodies(minUT, maxUT, hWaitbar);
        end

        %Override the following structure methods for backwards
        %compatibility with structure celBodyData
        function f = fields(obj)
            f = {};

            for(i=1:length(obj.bodies))
                f{end+1} = lower(obj.bodies(i).name); %#ok<AGROW>
            end

            f = f(:);
        end

        function f = fieldnames(obj)
            f = obj.fields();
        end

        function tf = isfield(obj,field)
            tf = isprop(obj,field);
        end

        function allBodyInfo = getAllBodyInfo(obj)
            allBodyInfo = obj.bodies;
        end

        function [listBoxStr, bodyInfos] = getListboxStr(obj)
            listBoxStr = {};

            for(i=1:length(obj.bodies))
                listBoxStr{end+1} = obj.bodies(i).name; %#ok<AGROW>
            end

            bodyInfos = obj.bodies;
        end

        function topLvlBody = getTopLevelBody(obj)
            if(isempty(obj.topLvlBodyCache))
                obj.topLvlBodyCache = getTopLevelCentralBody(obj);
            end

            topLvlBody = obj.topLvlBodyCache;
        end

        function bodies = getAllBodiesWithChildren(obj)
            bodies = KSPTOT_BodyInfo.empty(1,0);
            for(i=1:length(obj.bodies))
                bodyInfo = obj.bodies(i);

                cBodies = bodyInfo.getChildrenBodyInfo(obj);
                if(not(isempty(cBodies)))
                    bodies(end+1) = bodyInfo;
                end
            end
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            obj.bodyIdCacheArr = [obj.bodies.id];

            if(isprop(obj,'sun'))
                obj.topLvlBodyCache = getTopLevelCentralBody(obj);
            end

            if(isempty(obj.ci))
                obj.ci = CelestialBodyIntegration(obj);
                obj.ci.integrateCelestialBodies(0, obj.numIntBlockSize);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function value = r8_factorial ( n )

%*****************************************************************************80
%
%% R8_FACTORIAL returns N!.
%
%  Discussion:
%
%    factorial ( N ) = Product ( 1 <= I <= N ) I
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the function.
%    0 <= N.
%
%    Output, real VALUE, the factorial of N.
%
  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_FACTORIAL - Fatal error!\n' );
    fprintf ( 1, '  N < 0.\n' );
    error ( 'R8_FACTORIAL - Fatal error!' );
  end

  value = 1.0;

  for i = 2 : n
    value = value * i;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x_sh, y_sh]=shift_spectra(spectra, X,C);
% function [x_sh, y_sh]=shift_spectra(spectra, X,C);
%   Function to shift spectra with a 10x better accuracy than original data:
%       -if variable have a step in x-axis 0.1 than the spectra are interpolated
%        to having step of 0.01, the spectra are shifted and then converted back
%        to original number of points
%   This prepares data for multivariate analysis which is sensitive to
%   shift in data
%   Preferred over shift_spectra_or.m function
%
% spectra - matrix of spectra [nxm]
% X-variable (wavelength, mass, etc.)[n x 1]
% C - value at which a maximum of peak  is desired
% created by K.Artyushkova
% kartyush@unm.edu
%

%% converts spectra and X variable into 10X more channels than original
% data
[m,n]=size(spectra);
step=X(1)-X(2);
stepN=step/10*-1;
x = X(1):stepN:(X(m)+stepN);
x=x';
for i=1:n;
y(:,i)=spline(X,spectra(:,i),x);
end

%% finds the maximum for each spectra - its position and shift in number
% of points (S) from the most high value of maximum through the set of
% spectra
for i=1:n
    [k,r(i)]=max(y(:,i));
    xx(i)=x(r(i));
end
[q,s]=max(xx);
for i=1:n
    S(i)=xx(i)-xx(s);
end
stepN=stepN*-1;
S=double(single(1/stepN*abs(S)));
N=max(S);

%% combines all spectra into matrix yy. Now all spectra have the maximum at the
% same position of the variable
[mm,nn]=size(y);
for i=1:n
    Ss=S(i);
    yy(:,i)=y((Ss+1):(mm-(N-Ss)),i);
end

%% shifts X-axis to have a maximum at the desired value of C
X_x=x(1:(mm-(N-1)));
I = size(X_x);
shift=double(single((C-q)));
stepN=stepN*-1;
X_xx=(X_x(1)+shift):stepN:(X_x(I(1))+shift);
X_xx=X_xx';

%% converts data back to original number of points in X
mmm=round(I(1)/10);
for i=1:1:(mmm-1)
    yyy(i,:)=yy(i*10,:);
    X_xxx(i)=X_xx(i*10);
end

y_sh=yyy;
x_sh=X_xxx';
reverplot(x_sh, y_sh)


text:
function b = poisson_rhs ( nrow, ncol, rhs_num )

%*****************************************************************************80
%
%% POISSON_RHS returns the right hand side of a Poisson linear system.
%
%  Discussion:
%
%    The Poisson matrix is associated with an NROW by NCOL rectangular
%    grid of points.
%
%    Assume that the points are numbered from left to right, bottom to top.
%
%    If the K-th point is in row I and column J, set X = I + J.
%
%    This will be the solution to the linear system.
%
%    The right hand side is easily determined from X.  It is 0 for every
%    interior point.
%
%  Example:
%
%    NROW = 3, NCOL = 3
%
%    ^
%    |  7  8  9
%    J  4  5  6
%    |  1  2  3
%    |
%    +-----I---->
%
%    Solution vector X = ( 2, 3, 4, 3, 4, 5, 4, 5, 6 )
%
%    Right hand side B = ( 2, 2, 8, 2, 0, 6, 8, 6, 14 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gene Golub, Charles Van Loan,
%    Matrix Computations, second edition,
%    Johns Hopkins University Press, Baltimore, Maryland, 1989
%    (Section 4.5.4).
%
%  Parameters:
%
%    Input, integer NROW, NCOL, the number of rows and columns
%    in the grid.
%
%    Input, integer RHS_NUM, the number of right hand sides.
%
%    Output, real B(NROW*NCOL,RHS_NUM), the right hand side.
%
  n = nrow * ncol;

  b = zeros ( n, rhs_num );

  k = 0;
  for j = 1 : nrow
    for i = 1 : ncol
      k = k + 1;
      b(k,1) = 0.0;
      if ( i == 1 )
        b(k,1) = b(k,1) + i + j - 1;
      end
      if ( j == 1 )
        b(k,1) = b(k,1) + i + j - 1;
      end
      if ( i == ncol )
        b(k,1) = b(k,1) + i + j + 1;
      end
      if ( j == nrow )
        b(k,1) = b(k,1) + i + j + 1;
      end
    end
  end

  return
end


text:
clear all;
clc;

addpath = './util';

load('300W_test_indoor.mat');

title = '300W-Test-Indoor';

num = 4;

colors = lines(num);

x_limit = 0.08;

linewidth = 3;

fontsize = 12;

% Render curves
plot_ced(x, y, methods, strrep(title,'_',' '), ...
         x_limit, colors, linewidth, fontsize);

ax = gca;
outerpos = ax.OuterPosition;
ti = ax.TightInset;
left = outerpos(1) + ti(1);
bottom = outerpos(2) + ti(2);
ax_width = outerpos(3) - ti(1) - ti(3);
ax_height = outerpos(4) - ti(2) - ti(4);
ax.Position = [left bottom ax_width ax_height];

fig = gcf;
fig.PaperPositionMode = 'auto';
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)]* 1.05;

print(fig,'-dpdf',[title '.pdf']);

author: jiankangdeng

repo: MenpoBenchmark


text:
% std_erpplot() - Command line function to plot STUDY cluster component ERPs. Either
%                 displays grand mean ERPs for all requested clusters in the same figure,
%                 with ERPs for different conditions (if any) plotted in different colors.
%                 Else, displays ERP for each specified cluster in separate figures
%                 (per condition), each containing the cluster component ERPs plus
%                 the grand mean cluster ERP (in bold). ERPs can be plotted only if
%                 component ERPs were computed and saved in the STUDY EEG
%                 datasets.
%                 These can be computed during pre-clustering using the gui-based
%                 function pop_preclust() or the equivalent command line functions
%                 eeg_createdata() and eeg_preclust(). Called by pop_clustedit().
%                 and std_propplot().
% Usage:
%   >> [STUDY] = std_erpplot(STUDY, ALLEEG, key1, val1, key2, val2);
%   >> [STUDY erpdata erptimes pgroup pcond pinter] = std_erpplot(STUDY, ALLEEG, ...);
%
% Inputs:
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - global EEGLAB vector of EEG structures for the datasets included
%                in the STUDY. A STUDY set ALLEEG is typically created by load_ALLEEG().
%
% Optional inputs for component plotting:
%   'clusters' - [numeric vector|'all'] indices of clusters to plot.
%                If no component indices ('comps' below) are given, the average
%                ERPs of the requested clusters are plotted in the same figure,
%                with ERPs for different conditions (and groups if any) plotted
%                in different colors. In 'comps' (below) mode, ERPS for each
%                specified cluster are plotted in separate figures (one per
%                condition), each overplotting cluster component ERPs plus the
%                average cluster ERP in bold. Note this parameter has no effect
%                if the 'comps' option (below) is used. {default: 'all'}
%   'comps'    - [numeric vector|'all'] indices of the cluster components to plot.
%                Note that 'comps', 'all' is equivalent to 'plotsubjects', 'on'.
%
% Optional inputs for channel plotting:
%   'channels' - [numeric vector]  specific channel group to plot. By
%                default, the grand mean channel ERP is plotted (using the
%                same format as for the cluster component means described above)
%   'subject'  - [numeric vector]  In 'changrp' mode (above), index of
%                the subject(s) to plot. Else by default, plot all components
%                in the cluster.
%   'plotsubjects' - ['on'|'off'] When 'on', plot ERP of all subjects.
%
% Other optional inputs:
%   'key','val' - All optional inputs to pop_erpparams() are also accepted here
%                 to plot subset of time, statistics etc. The values used by default
%                 are the ones set using pop_erpparams() and stored in the
%                 STUDY structure.
%
% Outputs:
%   STUDY      - the input STUDY set structure with plotted cluster mean
%                ERPs data to allow quick replotting
%   erpdata    - [cell] ERP data for each condition, group and subjects.
%                size of cell array is [nconds x ngroups]. Size of each element
%                is [times x subjects] for data channels or [times x components]
%                for component clusters. This array may be gicen as input
%                directly to the statcond() function or std_stats function
%                to compute statistics.
%   erptimes   - [array] ERP time point latencies.
%   pgroup     - [array or cell] p-values group statistics. Output of the
%                statcond() function.
%   pcond      - [array or cell] condition statistics. Output of the statcond()
%                function.
%   pinter     - [array or cell] groups x conditions statistics. Output of
%                statcond() function.
%
%   Example:
%            >> [STUDY] = std_erpplot(STUDY,ALLEEG, 'clusters', 2, 'comps', 'all');
%               % Plot cluster-2 component ERPs plus the mean ERP in bold.
%
%  See also  pop_clustedit(), pop_preclust(), eeg_createdata(), eeg_preclust(). std_propplot()
%
% Authors: Arnaud Delorme, CERCO, August, 2006-

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

function [STUDY, erpdata, alltimes, pgroup, pcond, pinter] = std_erpplot(STUDY, ALLEEG, varargin)

if nargin < 2
    help std_erpplot;
    return;
end;
erpdata = []; alltimes = [];
pgroup = []; pcond = []; pinter = [];

% find datatype and default options
% ---------------------------------
dtype = 'erp';
for ind = 1:2:length(varargin)
    if strcmpi(varargin{ind}, 'datatype')
        dtype = varargin{ind+1};
    end;
end;
eval( [ 'STUDY = pop_' dtype 'params(STUDY, ''default'');' ...
        'params = STUDY.etc.' dtype 'params;' ] );
fields     = { 'filter' 'subtractsubjectmean' 'timerange' 'freqrange' 'topotime' 'topofreq' };
defaultval = { [] 'off' [] [] [] [] };
for ind=1:length(fields)
    if ~isfield(params, fields{ind}),
        params = setfield(params, fields{ind}, defaultval{ind});
    end;
end;

opt = finputcheck( varargin, { 'topotime'    'real'    [] params.topotime;
                               'topofreq'    'real'    [] params.topofreq;
                               'filter'      'real'    [] params.filter;
                               'timerange'   'real'    [] params.timerange;
                               'freqrange'   'real'    [] params.freqrange;
                               'ylim'        'real'    [] params.ylim;
                               'caxis'       'real'    [] params.ylim;
                               'statistics'  'string'  [] params.statistics;
                               'groupstats'  'string'  [] params.groupstats;
                               'condstats'   'string'  [] params.condstats;
                               'mcorrect'    'string'  [] params.mcorrect;
                               'plotgroups'  'string'  [] params.plotgroups;
                               'plotconditions' 'string'  [] params.plotconditions;
                               'subtractsubjectmean' 'string' [] params.subtractsubjectmean;
                               'threshold'   'real'    [] params.threshold;
                               'naccu'       'integer' [] params.naccu;
                               'singletrials' 'string' { 'on' 'off' }  params.singletrials;
                               'design'      'integer' []              STUDY.currentdesign;
                               'plotstderr'  'string'  []              'off';
                               'channels'    'cell'    []              {};
                               'clusters'    'integer' []              [];
                               'datatype'    'string'  { 'erp' 'spec' } 'erp';
                               'mode'        'string'  []              ''; % for backward compatibility
                               'comps'       { 'string' 'integer' } [] []; % for backward compatibility
                               'plotmode'    'string' { 'normal' 'condensed' }  'normal';
                               'unitx'       'string' { 'ms' 'Hz' }    'ms';
                               'plotsubjects' 'string' { 'on' 'off' }  'off';
                               'subject'     'string'  []              '';
                               'statmode'    'string'  { 'subjects' 'common' 'trials' } 'subjects'}, 'std_erpplot');
if isstr(opt), error(opt); end;
if isstr(opt.comps), opt.comps = []; opt.plotsubjects = 'on'; end;
if ~isempty(opt.topofreq),  opt.topotime  = opt.topofreq; end;
if ~isempty(opt.freqrange), opt.timerange = opt.freqrange; end;
datatypestr = upper(opt.datatype);
if strcmpi(datatypestr, 'spec'), datatypestr = 'Spectrum'; end;

% =======================================================================
% below this line, all the code should be non-specific to ERP or spectrum
% =======================================================================

allconditions = STUDY.design(opt.design).variable(1).value;
allgroups     = STUDY.design(opt.design).variable(2).value;
paired = { STUDY.design(opt.design).variable(1).pairing ...
           STUDY.design(opt.design).variable(2).pairing };

% for backward compatibility
% --------------------------
if strcmpi(opt.mode, 'comps'), opt.plotsubjects = 'on'; end;
if strcmpi(opt.singletrials, 'off') && ((~isempty(opt.subject) || ~isempty(opt.comps)))
    if strcmpi(opt.condstats, 'on') || strcmpi(opt.groupstats, 'on')
        opt.groupstats = 'off';
        opt.condstats   = 'off';
        disp('No statistics for single subject/component');
    end;
end;

if ~isnan(opt.topotime) & length(opt.channels) < 5
    warndlg2(strvcat('ERP parameters indicate that you wish to plot scalp maps', 'Select at least 5 channels to plot topography'));
    return;
end;

plotcurveopt = {};
if length(opt.clusters) > 1
    plotcurveopt = { 'figure' 'off' };
    opt.plotconditions = 'together';
    opt.plotgroups     = 'together';
    opt.condstats  = 'off';
    opt.groupstats = 'off';
end;
% if length(opt.channels) > 1 && strcmpi(opt.plotconditions, 'together') && strcmpi(opt.plotgroups, 'together')
%     plotcurveopt = { 'figure' 'off' };
%     opt.plotconditions = 'together';
%     opt.plotgroups     = 'together';
%     opt.condstats  = 'off';
%     opt.groupstats = 'off';
% end;
plotcurveopt = { plotcurveopt{:} ...
   'ylim',           opt.ylim, ...
   'threshold',      opt.threshold, ...
   'unitx'           opt.unitx, ...
   'filter',         opt.filter, ...
   'plotgroups',     opt.plotgroups, ...
   'plotconditions', opt.plotconditions };

% channel plotting
% ----------------
if ~isempty(opt.channels)

    chaninds = 1:length(opt.channels);

    if strcmpi(opt.datatype, 'erp')
        [STUDY erpdata alltimes] = std_readerp(STUDY, ALLEEG, 'channels', opt.channels(chaninds), 'timerange', opt.timerange, ...
                'subject', opt.subject, 'singletrials', opt.singletrials, 'design', opt.design);
    else
        [STUDY erpdata alltimes] = std_readspec(STUDY, ALLEEG, 'channels', opt.channels(chaninds), 'freqrange', opt.freqrange, ...
            'rmsubjmean', opt.subtractsubjectmean, 'subject', opt.subject, 'singletrials', opt.singletrials, 'design', opt.design);
    end;
    if isempty(erpdata), return; end;

    % select specific time
    % --------------------
    if ~isempty(opt.topotime) & ~isnan(opt.topotime)
        [tmp ti1] = min(abs(alltimes-opt.topotime(1)));
        [tmp ti2] = min(abs(alltimes-opt.topotime(end)));
        for condind = 1:length(erpdata(:))
            if ~isempty(erpdata{condind})
                erpdata{condind} = mean(erpdata{condind}(ti1:ti2,:,:),1);
            end;
        end;
    end;

    % compute statistics
    % ------------------
    [pcond pgroup pinter] = std_stat(erpdata, 'groupstats', opt.groupstats, 'condstats', opt.condstats, 'paired', paired, ...
                                         'statistics', opt.statistics, 'naccu', opt.naccu, 'threshold', opt.threshold, 'mcorrect', opt.mcorrect);
    if (~isempty(pcond) && length(pcond{1}) == 1) || (~isempty(pgroup) && length(pgroup{1}) == 1), pcond = {}; pgroup = {}; pinter = {}; end; % single subject STUDY
    if length(opt.channels) > 5 && ndims(erpdata{1}) < 3, pcond = {}; pgroup = {}; pinter = {}; end; % topo plotting for single subject

    % plot
    % ----
    locs = eeg_mergelocs(ALLEEG.chanlocs);
    locs = locs(std_chaninds(STUDY, opt.channels(chaninds)));
    [alltitles alllegends ] = std_figtitle('threshold', opt.threshold, 'mcorrect', opt.mcorrect, 'condstat', opt.condstats, 'cond2stat', opt.groupstats, ...
                             'statistics', opt.statistics, 'condnames', allconditions, 'plotsubjects', opt.plotsubjects, 'cond2names', allgroups, 'chanlabels', { locs.labels }, ...
                             'subject', opt.subject, 'valsunit', opt.unitx, 'vals', opt.topotime, 'datatype', datatypestr, 'cond2group', opt.plotgroups, 'condgroup', opt.plotconditions);

    if ~isempty(opt.topotime) && all(~isnan(opt.topotime))
        std_chantopo(erpdata, 'groupstats', pgroup, 'condstats', pcond, 'interstats', pinter, 'caxis', opt.caxis, ...
                                      'chanlocs', locs, 'threshold', opt.threshold, 'titles', alltitles);
    else
        std_plotcurve(alltimes, erpdata, 'groupstats', pgroup, 'legend', alllegends, 'condstats', pcond, 'interstats', pinter, ...
            'chanlocs', locs, 'titles', alltitles, 'plotsubjects', opt.plotsubjects, 'plotstderr', opt.plotstderr, ...
            'condnames', allconditions, 'groupnames', allgroups, 'plottopo', fastif(length(opt.channels) > 1, 'on', 'off'), plotcurveopt{:});
    end;

    set(gcf,'name',['Channel ' datatypestr ]);
    axcopy(gca);
else
    % plot component
    % --------------
    if length(opt.clusters) > 1, figure('color', 'w'); end;
    nc = ceil(sqrt(length(opt.clusters)));
    nr = ceil(length(opt.clusters)/nc);
    comp_names = {};

    for index = 1:length(opt.clusters)

        if length(opt.clusters) > 1, subplot(nr,nc,index); end;
        if strcmpi(opt.datatype, 'erp')
            [STUDY erpdata alltimes] = std_readerp(STUDY, ALLEEG, 'clusters', opt.clusters(index), 'timerange', opt.timerange, ...
                        'component', opt.comps, 'singletrials', opt.singletrials, 'design', opt.design);
        else
            [STUDY erpdata alltimes] = std_readspec(STUDY, ALLEEG, 'clusters', opt.clusters(index), 'freqrange', opt.freqrange, ...
                        'rmsubjmean', opt.subtractsubjectmean, 'component', opt.comps, 'singletrials', opt.singletrials, 'design', opt.design);
        end;
        if isempty(erpdata), return; end;

        % plot specific component
        % -----------------------
        if ~isempty(opt.comps)
            comp_names = { STUDY.cluster(opt.clusters(index)).comps(opt.comps) };
            opt.subject = STUDY.datasetinfo(STUDY.cluster(opt.clusters(index)).sets(1,opt.comps)).subject;
        end;
        [pcond pgroup pinter] = std_stat(erpdata, 'groupstats', opt.groupstats, 'condstats', opt.condstats, 'paired', paired, ...
                                         'statistics', opt.statistics, 'naccu', opt.naccu, 'threshold', opt.threshold, 'mcorrect', opt.mcorrect);

        [alltitles alllegends ] = std_figtitle('threshold', opt.threshold, 'plotsubjects', opt.plotsubjects, 'mcorrect', opt.mcorrect, 'condstat', opt.condstats, 'cond2stat', opt.groupstats, ...
                                 'statistics', opt.statistics, 'condnames', allconditions, 'cond2names', allgroups, 'clustname', STUDY.cluster(opt.clusters(index)).name, 'compnames', comp_names, ...
                                 'subject', opt.subject, 'valsunit', opt.unitx, 'vals', opt.topotime, 'datatype', datatypestr, 'cond2group', opt.plotgroups, 'condgroup', opt.plotconditions);

        if length(opt.clusters) > 1 && index < length(opt.clusters), alllegends = {}; end;
        std_plotcurve(alltimes, erpdata, 'condnames', allconditions, 'legend', alllegends, 'groupnames', allgroups, 'plotstderr', opt.plotstderr, ...
                                          'titles', alltitles, 'groupstats', pgroup, 'condstats', pcond, 'interstats', pinter, ...
                                          'chanlocs', ALLEEG(1).chanlocs, 'plotsubjects', opt.plotsubjects, plotcurveopt{:});
    end;

    set(gcf,'name', ['Component ' datatypestr ] );
    axcopy(gca);
end;

author: PatternRecognition

repo: OpenBMI


text:
function [LMTRad,LMT1,LMT2]=TT2LMT(TT1,TT2,rObsITRS,deltaT,xpyp)
%%TT2LMT Convert from terrestrial time to local mean solar time (LMT). This
%        is a measure of the time at a place on the Earth given the mean
%        location of the Sun, not a high-precision ephemeris.
%
%INPUTS: Jul1,Jul2 Two parts of a Julian date given in TT. The
%                  units of the date are days. The full date is the sum of
%                  both terms. The date is broken into two parts to
%                  provide more bits of precision. It does not matter how
%                  the date is split.
%         rObsITRS The 3X1 location of the observer in the International
%                  Terrestrial Reference System (ITRS). Only the direction
%                  matters, not the magnitude.
%           deltaT An optional parameter specifying the offset between TT
%                  and UT1 in seconds. If this parameter is omitted or if
%                  an empty matrix is passed, then the value of the
%                  function getEOP will be used.
%             xpyp xpyp=[xp;yp] are the polar motion coordinates in radians
%                  including the effects of tides and librations. If this
%                  parameter is omitted or if an empty matrix is passed,
%                  the value from the function getEOP will be used.
%
%OUTPUTS: LMTRad The local apparent solar time in radians. 0 radians is
%                solar midnight. pi radians is solar noon. The mapping is
%                is 2*pi radians for 24 hours. Thus, LATRad*(24/(2*pi))
%                gives the time of day in hours.
%      LMT1,LMT2 Two parts of the local mean solar time in Julian days. The
%                date is split so that LMT1 is the integer part and LMT2 is
%                the fractional part. Like UT1, a zero fractional part of
%                the day corresponds to noon, not midnight.
%
%The local mean solar time is UT1 plus the East longitude of the observer
%in radians in the Terrestrial Intermediate Reference System (TIRS). The
%East longitude is added using the convention that 360 degrees equals 24
%hours of 60 minutes and 60 seconds. For LMT in radians, an additional 12
%hours (pi radians) is added to adhere to the convention that 0 hours/ 0
%radians is midnight, not noon.
%
%The use of UT1 as a definition of Greenwhich mean solar time is given in
%[1]. Since UT1 is defined as a rotation in the TIRS, it only makes sense
%that a local mean solar time is defined by using the longitude offset in
%the TIRS.
%
%REFERENCES:
%[1] D. D. McCarthy, "Evolution of timescales from astronomy to physical
%    metrology," Metrologia, vol. 48, no. 4, pp. S132-S144, Aug. 2011.
%
%April 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Get any Earth orientation parameters that were not provided.
if(nargin<5||isempty(deltaT)||isempty(xpyp))
    [UTC1,UTC2]=TT2UTC(TT1,TT2);
    [xpypNew,~,~,deltaTTUT1]=getEOP(UTC1,UTC2);

    if(nargin<4||isempty(deltaT))
       deltaT=deltaTTUT1;
    end

    if(nargin<5||isempty(xpyp))
        xpyp=xpypNew;
    end
end

[UT11,UT12]=TT2UT1(TT1,TT2,deltaT);
rObsTIRS=ITRS2TIRS(rObsITRS,TT1,TT2,xpyp);
rObsSphere=Cart2Sphere(rObsTIRS);

%2*pi radians per day.
deltaUT1=rObsSphere(2)/(2*pi);

%Add preserving precision.
if(UT11<UT12)
    LMT2=UT11+deltaUT1;
    LMT1=UT12;
else
    LMT2=UT12+deltaUT1;
    LMT1=UT11;
end

%Put the fractional part into LMT2 and the integer part into LMT1.
frac1=LMT1-fix(LMT1);
frac2=LMT2-fix(LMT2);
sumVal=frac1+frac2;
sumInt=fix(sumVal);
sumFrac=sumVal-sumInt;
LMT1=fix(LMT2)+fix(LMT1)+sumInt;
LMT2=sumFrac;

%The local mean solar time in radians. This is equivalent to
%wrapRange(2*pi*(LMT1+LMT2)+pi,-pi,pi); However, since LMT1 is always
%integers, it can just be omitted.
LMTRad=wrapRange(2*pi*LMT2+pi,-pi,pi);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% eeglabexefolder - return the exe folder for EEGLAB. This function is only
%                   relevant for the compiled version of EEGLAB.
%
% Author: Arnaud Delorme, SCCN & CERCO, CNRS, 2008-

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

function str = eeglabexefolder;

str = ctfroot;
inds = find(str == filesep);
str = str(1:inds(end)-1);


text:
%% brute_force_tune
% Code to test the performance of various tuning parameters
% Works sorta like RANSAC I guess?
% Adam Werries 2016, see Apache 2.0 license.

k_max = 50;
% Specify ranges
accel_noise_PSD = logspace(-8, 5, 100);
gyro_noise_PSD = logspace(-8, 5, 100);
% Repeat arrays
accel_noise_PSD = repmat(accel_noise_PSD, [1 k_max]);
gyro_noise_PSD = repmat(gyro_noise_PSD, [1 k_max]);
% Generate random selections of each vector
num_items = length(accel_noise_PSD);
accel_i = randperm(num_items);
gyro_i = randperm(num_items);

rms_error_filter = Inf*ones(1,num_items);
max_error_filter = Inf*ones(1,num_items);
parfor i = 1:num_items
    fprintf('Iteration: %d, ANoise: %08.5e, GNoise: %08.5e\n', i, accel_noise_PSD(accel_i(i)), gyro_noise_PSD(gyro_i(i)));
    temp_conf = LC_KF_config;
    temp_conf.accel_noise_PSD = (accel_noise_PSD(accel_i(i)) * mug_to_mps2)^2;
    temp_conf.gyro_noise_PSD = (gyro_noise_PSD(gyro_i(i)) * deg_to_rad / 60)^2;
    [out_profile,out_IMU_bias_est,out_KF_SD] = Loosely_coupled_INS_GNSS(init_cond, filter_time, epoch, lla, gps, imu, temp_conf, est_IMU_bias);
    xyz = out_profile(:,2:4);
    if ~any(any(isnan(xyz))) && ~any(any(isinf(xyz)))
        llh = ecef2lla(xyz);
        [x,y] = deg2utm(llh(:,1),llh(:,2));
        x = x-min_x;
        y = y-min_y;

        distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
        rms_error_filter(i) = rms(distance);
        max_error_filter(i) = max(distance);
    end
end

[minmax, i] = min(max_error_filter);
fprintf('\nBest max: %08.4f, rms is %08.4f\n', minmax, rms_error_filter(i));
fprintf('Best iteration for max: %d, ANoise: %08.5e, GNoise: %08.5e\n', i, accel_noise_PSD(accel_i(i)), gyro_noise_PSD(gyro_i(i)));
[minrms, i] = min(rms_error_filter);
fprintf('Best rms: %08.4f, max is %08.4f\n', minrms, max_error_filter(i));
fprintf('Best iteration for rms: %d, ANoise: %08.5e, GNoise: %08.5e\n', i, accel_noise_PSD(accel_i(i)), gyro_noise_PSD(gyro_i(i)));
[minrms, i] = min((rms_error_filter+max_error_filter)/2);
fprintf('Best average of RMS and max: %08.4f, rms is  %08.4f, max is %08.4f\n', minrms, rms_error_filter(i), max_error_filter(i));
fprintf('Best iteration for rms: %d, ANoise: %08.5e, GNoise: %08.5e\n', i, accel_noise_PSD(accel_i(i)), gyro_noise_PSD(gyro_i(i)));

author: awerries

repo: kalman-localization


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% Calibrating kl, and ku for a sabr model for Kienitz extrapolation
% the prices are calculated from known sabr parameters

f = 0.03; t =1;                         % forward and time
a = 0.25; b = 0.5; r = -.5;  n = 0.2;   % sabr parameters
mu = 1.5; nu = 3;                       % tail decay

k = 0.001:0.0001:1;                         % strike range

call = sprice(a, b, r, n, f, k, t,1);   % call prices (standard sabr)
put = sprice(a, b, r, n, f, k, t,0);    % put prices (standard sabr)
put(1) = 0;
call(1) = f;                            % assures forward is matched

Nparam = 2;                             % number of calibrated params

% objective function
of = @(x) of_sabr(a, b, r, n, f, k, t, mu, nu, x(1), x(2), call, put);

x0 = [.25*f; 25.5*f];                   % starting values

A = zeros(Nparam,Nparam); bc= zeros(Nparam,1);
Aeq = A; beq = bc;
lb = [.25*f; f];                        % lower bound
ub = [f; 30*f];                         % upper bound

y = fmincon(of,x0,A,bc,Aeq,beq,lb,ub);  % optimization

% verification of results
xval = 0:0.001:.25;                     % x-values
[cl, bl, al, cu, bu, au] = ...
    psabr_param_3(a, b, r, n, f, t,mu,nu,y(1),y(2));
yval = psabr_5(a, b, r, n, f, xval, t, y(1), y(2), ...
           mu, cl, bl,al, nu, cu,bu,au);% y-values calculated

plot(xval,yval);                        % plot the results
Factor = 1000000;                       % used for plotting

yval_call = sprice_5(a, b, r, n, f, xval, t, y(1), y(2), ...
    mu, cl, bl,al, nu, cu,bu,au, 1);    % SABR Call prices
figure; hold on;
    plot(xval,Factor*yval_call,'r');
    plot(k,Factor*call,'g');
hold off;

yval_put = sprice_5(a, b, r, n, f, xval, t, y(1), y(2), ...
    mu, cl, bl,al, nu, cu,bu,au, 0);    % SABR Put prices
figure; hold on;
    plot(xval,Factor*yval_put,'r');
    plot(k,Factor*put,'g');
hold off;

fval = sprice_5(a, b, r, n, f, 0, t, y(1), y(2), ...
    mu, cl, bl,al, nu, cu,bu,au, 1);    % calculate forward value
y                                       % calibrated values
f - fval                                % display difference to forward


text:
classdef VademecumDataLoader < handle

   properties (Access = public)
      cellDataSameSign
      cellDataDifferentSign
   end

   properties (Access = private)
        vademecumPath
        fileNameSameSign
        fileNameDiffSign
        nMx
        nMy
        nPhi
   end

   methods (Access = public)

       function obj = VademecumDataLoader(cParams)
            obj.init(cParams);
            obj.loadVademecumDataBase();
       end

   end

   methods (Access = private)

       function init(obj,cParams)
          obj.nMx = cParams.nMx;
          obj.nMy = cParams.nMy;
          obj.nPhi = cParams.nPhi;
          obj.vademecumPath = '/media/alex/My Passport/Vademecum/';
          obj.fileNameSameSign = 'OptimalSuperEllipseSameStressSign';
          obj.fileNameDiffSign = 'OptimalSuperEllipseDifferentStressSign';
       end

        function loadVademecumDataBase(obj)
            cS = cell(obj.nMx,1);
            cD = cell(obj.nMx,1);
            for iMx = 1:obj.nMx
                [cS{iMx},cD{iMx}] = obj.loadBothVademecum(iMx);
            end
            obj.cellDataSameSign      = cS;
            obj.cellDataDifferentSign = cD;
        end

        function [cS,cD] = loadBothVademecum(obj,iMx)
            cS = obj.loadSameSignVademecum(iMx);
            cD = obj.loadDifferentSignVademecum(iMx);
        end

        function d = loadSameSignVademecum(obj,iMx)
            fName = obj.fileNameSameSign;
            d = obj.loadVademecum(fName,iMx);
        end

        function d = loadDifferentSignVademecum(obj,iMx)
            fName = obj.fileNameDiffSign;
            d = obj.loadVademecum(fName,iMx);
        end

     function cS = loadVademecum(obj,fileName,iMx)
            fName = [obj.vademecumPath,fileName,num2str(iMx),'.mat'];
            d = load(fName);
            cS = cell(obj.nMx,obj.nMy,obj.nPhi);
            for iMy = 1:obj.nMy
                for iphi = 1:obj.nPhi
                    cS{iMx,iMy,iphi}.mesh = d.c{iMx,iMy,iphi}.mesh;
                    cS{iMx,iMy,iphi}.rho = d.c{iMx,iMy,iphi}.rho;
                    cS{iMx,iMy,iphi}.q   = d.c{iMx,iMy,iphi}.q;
                    cS{iMx,iMy,iphi}.xi  = d.c{iMx,iMy,iphi}.xi;
                    cS{iMx,iMy,iphi}.phi = d.c{iMx,iMy,iphi}.phi;
                end
            end
        end

   end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [k, sk, skIC] = heatKernDiagCompute(heatKern, x)

% HEATKERNDIAGCOMPUTE Diagonal of a kernel matrix for a HEAT kernel.
% DESC computes the diagonal of the kernel matrix for the HEAT kernel
% given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
% RETURN sk : unscaled version of the diagonal.
% RETURN sk : unscaled version of the diagonal for the kernel matrix of the
% initial conditions.
%
% SEEALSO : heatKernParamInit, kernDiagCompute, kernCreate, heatKernCompute
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if size(x, 2) ~= 2
    error('Input can only have two columns');
end

option = 'cos';

if strcmp(option, 'cos')
    [k, sk] = heatXheatKernCompute(heatKern, heatKern, x, x);
    k = diag(k);
    sk = diag(sk);
else
    % Split the domain into time domain and spatial domain and account for
    % missing values. If there are no missing values the computation of the
    % kernel is a pointwise prodruct, otherwise it is a kronecker product.
    t = x(x(:,1)~=Inf,1);
    s = x(x(:,2)~=Inf,2);

    if (length(t) == length(s))
        ut = unique(t);
        us = unique(s);
        if (length(ut)*length(us) == length(t))
            t = ut;
            s = us;
            isPointwise = false;
            sk = zeros(length(t)*length(s), 1);
            if heatKern.includeIC
                skIC = zeros(length(t)*length(s), 1);
            end
        else
            isPointwise = true;
            sk = zeros(length(t), 1);
            if heatKern.includeIC
                skIC = zeros(length(t), 1);
            end
        end
    else
        isPointwise = false;
        sk = zeros(length(t)*length(s), 1);
        if heatKern.includeIC
            skIC = zeros(length(t)*length(s), 1);
        end
    end

    % Although this is done in heatKernExpandParam.m, we do it here again as a
    % precaution.

    heatKern.sim.inverseWidth = heatKern.inverseWidthTime;

    sigmax  = sqrt(2/heatKern.inverseWidthSpace);
    lengthX = heatKern.lengthX;
    nterms  = heatKern.nTerms;
    decay   = heatKern.decay;
    diff    = heatKern.diffusion;

    % Precompute some terms
    w = ((1:nterms)*(pi/lengthX))';
    gamma = sqrt(-1)*w;
    beta = decay + diff*(w.^2);
    z1 = sigmax*gamma/2;
    z2 = lengthX/sigmax + z1;
    wz1 = wofzPoppe(sqrt(-1)*z1);
    wz2 = wofzPoppe(sqrt(-1)*z2);
    cK = 4/(lengthX^2);

    simLocal = heatKern.sim;

    if heatKern.includeIC
        sigmah = sqrt(2/heatKern.inverseWidthSpaceIC);
        z1h = sigmah*gamma/2;
        z2h = lengthX/sigmah + z1h;
        wz1h = wofzPoppe(sqrt(-1)*z1h);
        wz2h = wofzPoppe(sqrt(-1)*z2h);
        if isPointwise
            for i=1:nterms
                heatKern.sim.decay = beta(i);
                kt = simKernDiagCompute(heatKern.sim, t);
                ks = sheatKernDiagCompute(sigmax, lengthX, s, w, gamma, wz1, wz2, i, i);
                ktIC = exp(-2*beta(i)*t);
                ksIC = sheatKernDiagCompute(sigmah, lengthX, s, w, gamma, wz1h, wz2h, i, i);
                skIC = skIC + ktIC.*ksIC;
                sk = sk + kt.*ks;
                for j=1:i-1
                    if (mod(i+j,2)==0)
                        heatKern.sim.decay = beta(i);
                        simLocal.decay = beta(j);
                        %kt = diag(simXsimKernCompute(heatKern.sim, simLocal, t, t));
                        kt = simXsimKernDiagCompute(heatKern.sim, simLocal, t);
                        ks = sheatKernDiagCompute(sigmax, lengthX, s, w, gamma, wz1, wz2, i, j);
                        ktIC = exp(-(beta(i)+beta(j))*t);
                        ksIC = sheatKernDiagCompute(sigmah, lengthX, s, w, gamma, wz1h, wz2h, i, j);
                        skIC = skIC + 2*ktIC.*ksIC;
                        sk = sk + 2*kt.*ks;
                    end
                end
            end
        else
            for i=1:nterms
                heatKern.sim.decay = beta(i);
                kt = simKernDiagCompute(heatKern.sim, t);
                ks = sheatKernDiagCompute(sigmax, lengthX, s, w, gamma, wz1, wz2, i, i);
                ktIC = exp(-2*beta(i)*t);
                ksIC = sheatKernDiagCompute(sigmah, lengthX, s, w, gamma, wz1h, wz2h, i, i);
                skIC = skIC + kron(ktIC, ksIC);
                sk = sk + kron(kt, ks);
                for j=1:i-1
                    if (mod(i+j,2)==0)
                        heatKern.sim.decay = beta(i);
                        simLocal.decay = beta(j);
                        %kt = diag(simXsimKernCompute(heatKern.sim, simLocal, t, t));
                        kt = simXsimKernDiagCompute(heatKern.sim, simLocal, t);
                        ks = sheatKernDiagCompute(sigmax, lengthX, s, w, gamma, wz1, wz2, i, j);
                        ktIC = exp(-(beta(i)+beta(j))*t);
                        ksIC = sheatKernDiagCompute(sigmah, lengthX, s, w, gamma, wz1h, wz2h, i, j);
                        skIC = skIC + 2*kron(ktIC, ksIC);
                        sk = sk + 2*kron(kt, ks);
                    end
                end
            end
        end
        sk = cK*sk;
        skIC = cK*skIC;
        k = (heatKern.sensitivity^2)*sk + (heatKern.sensitivityIC^2)*skIC;
    else
        if isPointwise
            for i=1:nterms
                heatKern.sim.decay = beta(i);
                kt = simKernDiagCompute(heatKern.sim, t);
                ks = sheatKernDiagCompute(sigmax, lengthX, s, w, gamma, wz1, wz2, i, i);
                sk = sk + kt.*ks;
                for j=1:i-1
                    if (mod(i+j,2)==0)
                        heatKern.sim.decay = beta(i);
                        simLocal.decay = beta(j);
                        %kt = diag(simXsimKernCompute(heatKern.sim, simLocal, t, t));
                        kt = simXsimKernDiagCompute(heatKern.sim, simLocal, t);
                        ks = sheatKernDiagCompute(sigmax, lengthX, s, w, gamma, wz1, wz2, i, j);
                        sk = sk + 2*kt.*ks;
                    end
                end
            end
        else
            for i=1:nterms
                heatKern.sim.decay = beta(i);
                kt = simKernDiagCompute(heatKern.sim, t);
                ks = sheatKernDiagCompute(sigmax, lengthX, s, w, gamma, wz1, wz2, i, i);
                sk = sk + kron(kt, ks);
                for j=1:i-1
                    if (mod(i+j,2)==0)
                        heatKern.sim.decay = beta(i);
                        simLocal.decay = beta(j);
                        %kt = diag(simXsimKernCompute(heatKern.sim, simLocal, t, t));
                        kt = simXsimKernDiagCompute(heatKern.sim, simLocal, t);
                        ks = sheatKernDiagCompute(sigmax, lengthX, s, w, gamma, wz1, wz2, i, j);
                        sk = sk + kron(2*kt, ks);
                    end
                end
            end
        end
        sk = cK*sk;
        k = (heatKern.sensitivity^2)*sk;
        if nargin > 2
            skIC = 0;
        end
    end
end

k = real(k);
sk = real(sk);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat