text:
classdef GemProgressPanel < CoreProgressInterface
    % GemProgressPanel. A panel used to report progress informaton
    %
    %     GemProgressPanel implements the CoreProgressInterface, which is an
    %     interface used by applications that use CoreMat to report the progress of
    %     operations. This panel is displayed over the centre of the panel
    %     object whose handle is supplied in the constructor.
    %
    %
    %     Licence
    %     -------
    %     Part of GEM. https://github.com/tomdoel/gem
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ShowDebuggingControls = false
    end

    properties (Access = private)

        % Handles to the GUI components
        Parent
        PanelHandle
        UiControlTitle
        UiControlText
        UiControlCancel
        UiControlQuit
        UiControlHide
        ProgressBarHandle

        IncrementThreshold = 5

        % Flags
        Hold = false;
        ShowProgressBar = false
        UserClickedCancel = false

        % The current values
        DialogText
        DialogTitle = ''
        ProgressValue = 0
        PanelVisibility = false;

        % The previous values
        CurrentlyDisplayedDialogTitle
        CurrentlyDisplayedDialogText
        CurrentlyDisplayedProgressValue
        CurrentlyDisplayedVisibility = false
        ProgressBarCurrentlyVisible = false;

        TimerRef
        MaxTimeBetweenUpdates = 0.25
        ChangePending
    end

    properties (GetAccess = private, Constant)
        OpacityLevel = 0.85
        Width = 500
        Height = 250
    end

    methods
        function obj = GemProgressPanel(parent)
            obj.ChangePending = false;
            obj.CurrentlyDisplayedDialogTitle = '';
            obj.CurrentlyDisplayedDialogText = '';
            obj.CurrentlyDisplayedProgressValue = [];

            obj.Parent = parent.GraphicalComponentHandle;

            set(obj.Parent, 'Units', 'Pixels');

            % Create the panel
            panel_background_colour = parent.StyleSheet.BackgroundColour;
            text_color = parent.StyleSheet.TextPrimaryColour;

            progress_position = obj.GetPanelPosition;
            obj.PanelHandle = uipanel('Parent', obj.Parent, 'Title', '', 'BorderType', 'etchedin', 'ForegroundColor', text_color, ...
                'BackgroundColor', panel_background_colour, 'Units', 'pixels', 'Position', progress_position, 'Visible', 'off' ...
            );

            title_position = [20, 165, 460, 55];
            text_position = [20, 115, 460, 50];
            cancel_position = [175, 20, 150, 40];
            quit_position = [370, 20, 80, 40];
            hide_position = [50, 20, 80, 40];
            progress_bar_position = [50, 90, 400, 18];

            obj.UiControlTitle = uicontrol('parent', obj.PanelHandle, 'style', 'text', 'units', 'pixel', 'Position', title_position, ...
                'string', 'Please wait', 'FontUnits', 'pixels', 'FontSize', 24, 'FontWeight', 'bold', 'Fore', text_color, 'Back', panel_background_colour, 'Visible', 'off');
            obj.UiControlText = uicontrol('parent', obj.PanelHandle, 'style', 'text', 'units', 'pixel', 'Position', text_position, ...
                'string', 'Please wait', 'FontUnits', 'pixels', 'FontSize', 16, 'Fore', text_color, 'Back', panel_background_colour, 'Visible', 'off');
            obj.UiControlCancel = uicontrol('parent', obj.PanelHandle, 'string', 'Cancel', ...
                'Position', cancel_position, 'FontUnits', 'pixels', 'Callback', @obj.CancelButton, 'Visible', 'off');
            obj.UiControlQuit = uicontrol('parent', obj.PanelHandle, 'string', 'Force Quit', ...
                'Position', quit_position, 'FontUnits', 'pixels', 'Callback', @obj.QuitButton, 'Visible', 'off');
            obj.UiControlHide = uicontrol('parent', obj.PanelHandle, 'string', 'Hide Dialog', ...
                'Position', hide_position, 'FontUnits', 'pixels', 'Callback', @obj.HideButton, 'Visible', 'off');

            [obj.ProgressBarHandle, ~] = javacomponent('javax.swing.JProgressBar', ...
                progress_bar_position, obj.PanelHandle);
            obj.ProgressBarHandle.setBackground(java.awt.Color(panel_background_colour(1), panel_background_colour(2), panel_background_colour(3)));
            set(obj.ProgressBarHandle, 'visible', 0);
        end

        function Resize(obj)
            set(obj.PanelHandle, 'Position', obj.GetPanelPosition);
        end

        function ShowAndHold(obj, text)
            if nargin < 2
                text = 'Please wait';
            end
            obj.DialogTitle = CoreTextUtilities.RemoveHtml(text);
            obj.DialogText = '';
            obj.ProgressValue = 0;
            obj.Hold = true;
            obj.UserClickedCancel = false;
            obj.PanelVisibility = true;
            obj.ShowProgressBar = false;
            obj.TimerRef = [];
            obj.Update();
        end

        function Hide(obj)
            obj.DialogTitle = 'Please wait';
            obj.Hold = false;
            obj.PanelVisibility = false;
            obj.ShowProgressBar = false;
            obj.Update();
        end

        function Complete(obj)
            % Call to complete a progress operaton, which will also hide the dialog
            % unless the dialog is being held

            obj.ShowProgressBar = false;
            obj.ProgressValue = 100;
            if ~obj.Hold
                obj.PanelVisibility = false;
            else
                obj.SetProgressText('');
            end
            obj.Update();
        end

        function SetProgressText(obj, text)
            if nargin < 2
               text = 'Please wait';
            end
            obj.DialogText = CoreTextUtilities.RemoveHtml(text);
            obj.PanelVisibility = true;
            obj.Update();
        end

        function SetProgressValue(obj, progress_value)
            obj.ShowProgressBar = true;
            obj.ProgressValue = progress_value;
            obj.PanelVisibility = true;
            obj.Update();
        end

        function SetProgressAndMessage(obj, progress_value, text)
            obj.ShowProgressBar = true;
            obj.DialogText = CoreTextUtilities.RemoveHtml(text);
            obj.ProgressValue = progress_value;
            obj.PanelVisibility = true;
            obj.Update();
        end

        function cancelled = CancelClicked(obj)
            cancelled = obj.UserClickedCancel;
            obj.UserClickedCancel = false;
        end
    end

    methods (Access = private)
        function changed = Update(obj)
            changed = false;
            changed = obj.UpdateDialogTitle || changed;
            changed = obj.UpdateDialogText || changed;
            changed = obj.UpdateProgressValue || changed;
            changed = obj.UpdateVisibility || changed;
            changed = obj.UpdateProgressBarVisibility || changed;

            if changed
                obj.ChangePending = true;
            end

            if obj.ChangePending
                if isempty(obj.TimerRef) || toc(obj.TimerRef) > obj.MaxTimeBetweenUpdates
                    obj.TimerRef = tic;
                    drawnow;
                    obj.ChangePending = false;
                end
            end

        end

        function changed = UpdateDialogTitle(obj)
            changed = false;
            if ~strcmp(obj.DialogTitle, obj.CurrentlyDisplayedDialogTitle)
                set(obj.UiControlTitle, 'String', obj.DialogTitle);
                obj.CurrentlyDisplayedDialogTitle = obj.DialogTitle;
                changed = true;
            end
        end

        function changed = UpdateDialogText(obj)
            changed = false;
            if ~strcmp(obj.DialogText, obj.CurrentlyDisplayedDialogText)
                set(obj.UiControlText, 'String', obj.DialogText);
                obj.CurrentlyDisplayedDialogText = obj.DialogText;
                changed = true;
            end
        end

        function changed = UpdateProgressValue(obj)
            changed = false;
            if isempty(obj.CurrentlyDisplayedProgressValue) || (obj.ProgressValue ~= obj.CurrentlyDisplayedProgressValue)
                obj.ProgressBarHandle.setValue(obj.ProgressValue);
                obj.CurrentlyDisplayedProgressValue = obj.ProgressValue;
                changed = true;
            end
        end

        function changed = UpdateProgressBarVisibility(obj)
            changed = false;
            if ~obj.ProgressBarCurrentlyVisible && obj.ShowProgressBar && obj.PanelVisibility
                set(obj.ProgressBarHandle, 'visible', 1);
                obj.ProgressBarCurrentlyVisible = true;
                changed = true;
            elseif obj.ProgressBarCurrentlyVisible && (~obj.ShowProgressBar || ~obj.PanelVisibility)
                set(obj.ProgressBarHandle, 'visible', 0);
                obj.ProgressBarCurrentlyVisible = false;
                changed = true;
            end
        end

        function changed = UpdateVisibility(obj)
            changed = false;
            if obj.PanelVisibility && ~obj.CurrentlyDisplayedVisibility
                set(obj.PanelHandle, 'Visible', 'on');
                set(obj.UiControlText, 'Visible', 'on');
                set(obj.UiControlTitle, 'Visible', 'on');
                if obj.ShowDebuggingControls
                    set(obj.UiControlQuit, 'Visible', 'on');
                    set(obj.UiControlHide, 'Visible', 'on');
                end
                set(obj.UiControlCancel, 'Visible', 'on');
                obj.CurrentlyDisplayedVisibility = true;
                changed = true;
            elseif ~obj.PanelVisibility && obj.CurrentlyDisplayedVisibility
                set(obj.UiControlText, 'Visible', 'off');
                set(obj.UiControlTitle, 'Visible', 'off');
                set(obj.UiControlQuit, 'Visible', 'off');
                set(obj.UiControlHide, 'Visible', 'off');
                set(obj.UiControlCancel, 'Visible', 'off');
                set(obj.ProgressBarHandle, 'visible', 0);
                set(obj.PanelHandle, 'Visible', 'off');
                obj.CurrentlyDisplayedVisibility = false;
                obj.ProgressBarCurrentlyVisible = false;
                changed = true;
            end
        end

        function CancelButton(obj, ~, ~)
            obj.UserClickedCancel = true;
            obj.PanelVisibility = false;
            obj.Hold = false;
            obj.Update();
        end

        function QuitButton(obj, ~, ~)
            obj.PanelVisibility = false;
            obj.Hold = false;
            obj.UpdateVisibility();
            throw(MException('GemProgressPanel:UserForceQuit', 'User forced plugin to quit'));
        end

        function HideButton(obj, ~, ~)
            obj.PanelVisibility = false;
            obj.Hold = false;
            obj.UpdateVisibility();
        end

        function progress_position = GetPanelPosition(obj)
            panel_position = get(obj.Parent, 'Position');

            window_width = panel_position(3);
            window_height = panel_position(4);
            progress_width = min(obj.Width, window_width);
            progress_height = min(obj.Height, window_height);
            progress_x = max(1, round((window_width - progress_width)/2));
            progress_y = max(1, (2/3)*round((window_height - progress_height)));
            progress_position = [progress_x, progress_y, progress_width, progress_height];
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function []=plot3DDICPPresults(varargin)
%% function for plotting 3D-DIC results in STEP3.
% Plotting the 3D reconstruction of points correlated with Ncorr
% The function opens a selection window for all the possible measures to plot
% After selection, the animation figures are plotted
%
% INPUT options:
% DIC3DPPresults

%%
switch nargin
    case 0 % in case no results were entered
        % ask user to load results from 1 or more camera pairs and turn
        % into a cell array
        [file,path] = uigetfile(pwd,'Select a DIC3DPPresults structure');
        result=load([path file]);
        DIC3DPPresults=result.DIC3DPPresults;
        optStruct=struct;
    case 1
        % use given struct
        DIC3DPPresults=varargin{1};
        optStruct=struct;
    case 2
        % use given struct
        DIC3DPPresults=varargin{1};
        optStruct=varargin{2};
    otherwise
        error('wrong number of input arguments');
end

%% select what to plot
Prompt={'\bf{Select which parameters to plot}';... % 1
    'Surfaces with color as pair index';... % 2
    'Points with color as pair index';... % 3
    'Surfaces with color as combined correlation coefficient';... % 4
    'Points with color as correlation coefficient';... % 5
    'Surfaces with color as dispMgn (Displacement magnitude)';... % 6
    'Points with color as dispMgn (Displacement magnitude)';... % 7
    'Surfaces with color as dispX (X Displacement)';... % 8
    'Points with color as dispX (X Displacement)';... % 9
    'Surfaces with color as dispY (Y Displacement)';... % 10
    'Points with color as dispY (Y Displacement)';... % 11
    'Surfaces with color as dispZ (Z Displacement)';... % 12
    'Points with color as dispZ (Z Displacement)';... % 13
    'Surfaces with color as J (surface area change)';... % 14
    'Surfaces with color as lambda1 (1st principal stretch)';... % 15
    'Surfaces with color as Lambda2 (2nd principal stretch)';... % 16
    'Surfaces with color as Epc1 (1st principal Lagrangian strain)';... % 17
    'Surfaces with color as Epc2 (2nd principal Lagrangian strain)';... % 18
    'Surfaces with color as epc1 (1st principal Almansi strain)';... % 19
    'Surfaces with color as epc2 (2nd principal Almansi strain)';... % 20
    'Surfaces with color as Emgn (Lagrangian strain tensor magnitude)';... % 21
    'Surfaces with color as emgn (Almansi strain tensor magnitude)';... % 22
    'Surfaces with color as Lagrangian equivalent strain';... % 23
    'Surfaces with color as Eulerian equivalent strain';... % 24
    'Surfaces with color as Max Lagrangian shear strain';... % 25
    'Surfaces with color as Max Eulerian shear strain';... % 26
    'Surfaces with color as Lamda1+direction (1st principal stretch value and direction)';... % 27
    'Surfaces with color as Lamda2+direction (2nd principal stretch value and direction)';... % 28
    'Surfaces with color as Epc1+direction (1st principal Lagrangian strain value and direction)';... % 29
    'Surfaces with color as Epc2+direction (2nd principal Lagrangian strain value and direction)';... % 30
    'Surfaces with color as epc1+direction (1st principal Almansi strain value and direction)';... % 31
    'Surfaces with color as epc2+direction (2nd principal Almansi strain value and direction)';... % 32
    'Surfaces with color as Epc1+Epc2+direction (1st and 2nd principal Lagrangian strain values and directions)';... % 33
    'Surfaces with color as epc1+epc2+direction (1st and 2nd principal Almansi strain values and directions)'; % 34
    'Surfaces with color as Max Lagrangian shear strain + directions';... % 35
    'Surfaces with color as Max Eulerian shear strain + directions';... % 36
    'Surfaces with color as Lamda1+Lamda2+direction (1st and 2nd principal stretch values and directions)';... % 37
    'Surfaces with color as FaceIsoInd (triangular face isotropy index)';... % 38
    'Select to remove rigid body motion';... % 39
    'Surfaces with color as FaceColors (grayscale from images)';... % 40
    'Select All'; }; % 41

Title='Select which parameters to plot';

Formats=struct;
Formats(1,1).type='text'; %1
Formats(1,2).type='none';
Formats(2,1).type='check'; %2
Formats(2,2).type='check'; %3
Formats(3,1).type='check'; %4
Formats(3,2).type='check'; %5
Formats(4,1).type='check'; %6
Formats(4,2).type='check'; %7
Formats(5,1).type='check'; %8
Formats(5,2).type='check'; %9
Formats(6,1).type='check'; %10
Formats(6,2).type='check'; %11
Formats(7,1).type='check'; %12
Formats(7,2).type='check'; %13
Formats(8,1).type='check'; %14
Formats(8,2).type='none';
Formats(9,1).type='check'; %15
Formats(9,2).type='check'; %16
Formats(10,1).type='check'; %17
Formats(10,2).type='check'; %18
Formats(11,1).type='check'; %19
Formats(11,2).type='check'; %20
Formats(12,1).type='check'; %21
Formats(12,2).type='check'; %22
Formats(13,1).type='check'; %23
Formats(13,2).type='check'; %24
Formats(14,1).type='check'; %25
Formats(14,2).type='check'; %26
Formats(15,1).type='check'; %27
Formats(15,2).type='check'; %28
Formats(16,1).type='check'; %29
Formats(16,2).type='check'; %30
Formats(17,1).type='check'; %31
Formats(17,2).type='check'; %32
Formats(18,1).type='check'; %33
Formats(18,2).type='check'; %34
Formats(19,1).type='check'; %35
Formats(19,2).type='check'; %36
Formats(20,1).type='check'; %37
Formats(20,2).type='none';
Formats(21,1).type='check'; %38
Formats(21,2).type='check'; %39
Formats(22,1).type='check'; %40
Formats(22,2).type='check'; %41

DefAns=cell(numel(Prompt),1);
DefAns{1}=[];
for ii=2:numel(Prompt)
    DefAns{ii}=false;
end

Options.Resize='on';
Options.FontSize=10;

[Answer,Canceled] = inputsdlg(Prompt, Title, Formats, DefAns, Options);
if Answer{41}
    for ii=[2:38 40]
        Answer{ii}=true;
    end
end
if Answer{39}
    RBMlogic=true;
else
    RBMlogic=false;
end
if Canceled
    return
end
%% Select plotting options
if ~Canceled && sum(cell2mat(Answer(2:end)))>0
    answer = inputdlg({'Enter maximum correlation coefficient to keep points (leave blank for keeping all points)'},'Input',[1,50]);
    CorCoeffCutOff=str2double(answer{1}); % maximal correlation coefficient for display (use [] for default which is max)
    if isnan(CorCoeffCutOff)
        CorCoeffCutOff=[];
    end
else
    return
end

%% create option struct for plotting
% complete the struct fields
if ~isfield(optStruct,'zDirection')
    optStruct.zDirection=1;
end
optStruct.maxCorrCoeff=CorCoeffCutOff;

% optStruct.lineColor='k';
% optStruct.smoothLogic=0;
% optStruct.colorMap=cMap;
% optStruct.supTitleString=pointMeasureString;

%% plot according to answer

if Canceled
    return
end

if Answer{2}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'FacePairInds',RBMlogic,optStruct);
end
if Answer{3}
    anim8_DIC3DPP_pointMeasure(DIC3DPPresults,'pairInd',RBMlogic,optStruct);
end
if Answer{4}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'FaceCorrComb',RBMlogic,optStruct);
end
if Answer{5}
    anim8_DIC3DPP_pointMeasure(DIC3DPPresults,'corrComb',RBMlogic,optStruct);
end
if Answer{6}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'DispMgn',RBMlogic,optStruct);
end
if Answer{7}
    anim8_DIC3DPP_pointMeasure(DIC3DPPresults,'DispMgn',RBMlogic,optStruct);
end
if Answer{8}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'DispX',RBMlogic,optStruct);
end
if Answer{9}
    anim8_DIC3DPP_pointMeasure(DIC3DPPresults,'DispX',RBMlogic,optStruct);
end
if Answer{10}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'DispY',RBMlogic,optStruct);
end
if Answer{11}
    anim8_DIC3DPP_pointMeasure(DIC3DPPresults,'DispY',RBMlogic,optStruct);
end
if Answer{12}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'DispZ',RBMlogic,optStruct);
end
if Answer{13}
    anim8_DIC3DPP_pointMeasure(DIC3DPPresults,'DispZ',RBMlogic,optStruct);
end
if Answer{14}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'J',RBMlogic,optStruct);
end
if Answer{15}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'Lamda1',RBMlogic,optStruct);
end
if Answer{16}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'Lamda2',RBMlogic,optStruct);
end
if Answer{17}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'Epc1',RBMlogic,optStruct);
end
if Answer{18}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'Epc2',RBMlogic,optStruct);
end
if Answer{19}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'epc1',RBMlogic,optStruct);
end
if Answer{20}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'epc2',RBMlogic,optStruct);
end
if Answer{21}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'Emgn',RBMlogic,optStruct);
end
if Answer{22}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'emgn',RBMlogic,optStruct);
end
if Answer{23}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'Eeq',RBMlogic,optStruct);
end
if Answer{24}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'eeq',RBMlogic,optStruct);
end
if Answer{25}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'EShearMax',RBMlogic,optStruct);
end
if Answer{26}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'eShearMax',RBMlogic,optStruct);
end
if Answer{27}
    anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,'Lamda1',RBMlogic,optStruct);
end
if Answer{28}
    anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,'Lamda2',RBMlogic,optStruct);
end
if Answer{29}
    anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,'Epc1',RBMlogic,optStruct);
end
if Answer{30}
    anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,'Epc2',RBMlogic,optStruct);
end
if Answer{31}
    anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,'epc1',RBMlogic,optStruct);
end
if Answer{32}
    anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,'epc2',RBMlogic,optStruct);
end
if Answer{33}
    anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,{'Epc1','Epc2'},RBMlogic,optStruct);
end
if Answer{34}
    anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,{'epc1','epc2'},RBMlogic,optStruct);
end
if Answer{35}
    % develop shear direction plot
    msgbox('Shear directions cannot be plotted. It is under construction...');
end
if Answer{36}
    % develop shear direction plot
    msgbox('Shear directions cannot be plotted. It is under construction...');
end
if Answer{37}
    anim8_DIC3DPP_faceMeasureDirection(DIC3DPPresults,{'Lamda1','Lamda2'},RBMlogic,optStruct);
end
if Answer{38}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'FaceIsoInd',RBMlogic,optStruct);
end
if Answer{40}
    anim8_DIC3DPP_faceMeasure(DIC3DPPresults,'FaceColors',RBMlogic,optStruct);
end

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function p = multivariateGaussian(X, mu, Sigma2)
%MULTIVARIATEGAUSSIAN Computes the probability density function of the
%multivariate gaussian distribution.
%    p = MULTIVARIATEGAUSSIAN(X, mu, Sigma2) Computes the probability
%    density function of the examples X under the multivariate gaussian
%    distribution with parameters mu and Sigma2. If Sigma2 is a matrix, it is
%    treated as the covariance matrix. If Sigma2 is a vector, it is treated
%    as the \sigma^2 values of the variances in each dimension (a diagonal
%    covariance matrix)
%

k = length(mu);

if (size(Sigma2, 2) == 1) || (size(Sigma2, 1) == 1)
    Sigma2 = diag(Sigma2);
end

X = bsxfun(@minus, X, mu(:)');
p = (2 * pi) ^ (- k / 2) * det(Sigma2) ^ (-0.5) * ...
    exp(-0.5 * sum(bsxfun(@times, X * pinv(Sigma2), X), 2));

end

author: Ayatans

repo: Machine-Learning-homework


text:
function y = mtrapezlogpdf(x,a,u,v,b)
%MTRAPEZLOGPDF Multivariate trapezoidal probability log pdf.
%   Y = MTRAPEZLOGPDF(X,A,U,V,B) returns the logarithm of the pdf of the
%   multivariate trapezoidal distribution with external bounds A and B and
%   internal points U and V, evaluated at the values in X. The multivariate
%   trapezoidal pdf is the product of univariate trapezoidal pdfs in each
%   dimension.
%
%   For each dimension i, the univariate trapezoidal pdf is defined as:
%
%                 |       __________
%                 |      /|        |\
%         p(X(i)) |     / |        | \
%                 |    /  |        |  \
%                 |___/___|________|___\____
%                    A(i) U(i)     V(i) B(i)
%                             X(i)
%
%   X can be a matrix, where each row is a separate point and each column
%   is a different dimension. Similarly, A, B, C, and D can also be
%   matrices of the same size as X.
%
%   The log pdf is typically preferred in numerical computations involving
%   probabilities, as it is more stable.
%
%   See also MTRAPEZPDF, MTRAPEZRND.

% Luigi Acerbi 2022

[N,D] = size(x);

if D > 1
    if isscalar(a); a = a*ones(1,D); end
    if isscalar(u); u = u*ones(1,D); end
    if isscalar(v); v = v*ones(1,D); end
    if isscalar(b); b = b*ones(1,D); end
end

if size(a,2) ~= D || size(u,2) ~= D || size(v,2) ~= D || size(b,2) ~= D
    error('mtrapezpdf:SizeError', ...
        'A, B, C, D should be scalars or have the same number of columns as X.');
end

if size(a,1) == 1; a = repmat(a,[N,1]); end
if size(u,1) == 1; u = repmat(u,[N,1]); end
if size(v,1) == 1; v = repmat(v,[N,1]); end
if size(b,1) == 1; b = repmat(b,[N,1]); end

y = -inf(size(x));
lnf = log(0.5) + log(b - a + v - u) + log(u - a);

for ii = 1:D
    idx = x(:,ii) >= a(:,ii) & x(:,ii) < u(:,ii);
    y(idx,ii) = log(x(idx,ii) - a(idx,ii)) - lnf(idx,ii);

    idx = x(:,ii) >= u(:,ii) & x(:,ii) < v(:,ii);
    y(idx,ii) = log(u(idx,ii)-a(idx,ii)) - lnf(idx,ii);

    idx = x(:,ii) >= v(:,ii) & x(:,ii) < b(:,ii);
    y(idx,ii) = log(b(idx,ii) - x(idx,ii)) - log(b(idx,ii) - v(idx,ii)) + log(u(idx,ii)-a(idx,ii)) - lnf(idx,ii);
end

y = sum(y,2);

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function val=calcGAE(xTrue,xEst,is3D)
%%CALCGAE Compute the scalar geometric average error (GAE). Unlike the
%         RMSE, it is not dominated by large individual terms.
%
%INPUTS: xTrue The truth data. This is either an xDimXNumSamples matrix or
%              an xDimXNXnumSamples matrix. The latter formulation is
%              useful when the MSE over multiple Monte Carlo runs of an
%              entire length-N track is desired. In the first formulation,
%              we take N=1. Alternatively, if the same true value is used
%              for all numSamples, then xTrue can just be an xDimXN matrix.
%              Alternatively, if xTrue is the same for all numSamples and
%              all N, then just an xDimX1 matrix can be passed. N and
%              numSamples are inferred from xEst.
%         xEst An xDimXnumSamples set of estimates or an xDimXNXnumSamples
%              set of estimates (if values at N times are desired).
%         is3D An optional indicating that xEst is 3D. This is only used if
%              xEst is a matrix. In such an instance, there is an ambiguity
%              whether xEst is truly 2D, so N=1, or whether numSamples=1
%              and xEst is 3D with the third dimension being 1. If this
%              parameter is omitted or an empty matrix is passed, it is
%              assumed that N=1.
%
%OUTPUTS: val The 1XN set of scalar GAE values.
%
%The GAE is given in Equation 3 in [1].
%
%EXAMPLE:
%For a Gaussian random vector, the root-trace of the covariance matrix is
%the RMSE. The RMSE is larger than the average Euclidean error, which is
%also larger than the  geometric average error
% R=[28,   4, 10;
%     4,  22, 16;
%     10, 16, 16];%The covariance matrix.
% xTrue=[10;-20;30];
% numRuns=100000;
% xEst=GaussianD.rand(numRuns,xTrue,R);
% valRMSE=rootTrace=sqrt(trace(R))
% valAEE=calcAEE(xTrue,xEst)
% valGAE=calcGAE(xTrue,xEst)
%
%REFERENCES:
%[1] X. R. Li and Z. Zhao, "Measures of performance for evaluation of
%    estimators and filters," in Proceedings of SPIE: Conference on Signal
%    and Data processing of Small Targets, vol. 4473, San Diego, CA, 29
%    Jul. 2001, pp. 530-541.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(is3D))
    is3D=false;
end

xDim=size(xEst,1);
if(ismatrix(xEst)&&is3D==false)
    N=1;
    numSamples=size(xEst,2);
    xEst=reshape(xEst,[xDim,1,numSamples]);

    if(size(xTrue,2)==1)
        %If the true values are the same for all samples.
        xTrue=repmat(xTrue,[1,1,numSamples]);
    else
        xTrue=reshape(xTrue,[xDim,1,numSamples]);
    end
else
    N=size(xEst,2);
    numSamples=size(xEst,3);

    if(ismatrix(xTrue))
        if(size(xTrue,2)==1)
            %If the true values are the same for all samples and for all N.
            xTrue=repmat(xTrue,[1,N,numSamples]);
        else
            %If the true values are the same for all samples.
            xTrue=repmat(xTrue,[1,1,numSamples]);
        end
    end
end

val=zeros(1,N);
for k=1:N
    val(k)=exp((1/(2*numSamples))*sum(log(sum((xEst(:,k,:)-xTrue(:,k,:)).^2,1))));
end
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:
hydro = struct();

hydro = readNEMOH(hydro,'../Sphere/');
% hydro = readWAMIT(hydro,'../../WAMIT/Sphere/sphere.out',[]);
% hydro = combineBEM(hydro); % Compare to WAMIT
hydro = radiationIRF(hydro,15,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,15,[],[],[],[]);
writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function varargout=gausprod(varargin)
%
% For calling details please see v_gausprod.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_gausprod.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: gausprod.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_gausprod(varargin{:});
else
    v_gausprod(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function f = calccod(x,y,dim,wantgain,wantmeansub,wantsafe)
%
% USAGE::
%
%   f = calccod(x,y,dim,wantgain,wantmeansub,wantsafe)
%
% <x>,<y> are matrices with the same dimensions
%
% <dim> (optional) is the dimension of interest.
% default to 2 if <x> is a (horizontal) vector and to 1 if not.
% special case is 0 which means to calculate globally.
%
% <wantgain> (optional) is
%
%   - 0 means normal
%   - 1 means allow a gain to be applied to each case of <x>
%     to minimize the squared error with respect to <y>.
%     in this case, there cannot be any NaNs in <x> or <y>.
%   - 2 is like 1 except that gains are restricted to be non-negative.
%     so, if the gain that minimizes the squared error is negative,
%     we simply set the gain to be applied to be 0.
%     default: 0.
%
% <wantmeansub> (optional) is
%
%   - 0 means do not subtract any mean.  this makes it such that
%     the variance quantification is relative to 0.
%   - 1 means subtract the mean of each case of <y> from both
%     <x> and <y> before performing the calculation.  this makes
%     it such that the variance quantification
%     is relative to the mean of each case of <y>.
%     note that <wantgain> occurs before <wantmeansub>.
%     default: 1.
%
% <wantsafe> (optional) is whether to protect against NaNs. Default: 1.
%
% calculate the coefficient of determination (R^2) indicating
% the percent variance in <y> that is explained by <x>.  this is achieved
% by calculating 100*(1 - sum((y-x).^2) / sum(y.^2)).  note that
% by default, we subtract the mean of each case of <y> from both <x>
% and <y> before proceeding with the calculation.
%
% the quantity is at most 100 but can be 0 or negative (unbounded).
% note that this metric is sensitive to DC and scale and is not symmetric
% (i.e. if you swap <x> and <y>, you may obtain different results).
% it is therefore fundamentally different than Pearson's correlation
% coefficient (see calccorrelation.m).
%
% if <wantsafe> is true, NaNs are handled gracefully (a NaN causes
% that data point to be ignored). if you are sure there are no NaNs,
% you can turn off <wantsafe> to increase execution speed.
%
% if there are no valid data points (i.e. all data points are
% ignored because of NaNs), we return NaN for that case.
%
% Note some weird cases:
%
% calccod([],[]) is []
%
% Example::
%
%   x = randn(1,100);
%   calccod(x,x+0.1*randn(1,100))
%

% history:
% 2013/08/18 - fix pernicious case where <x> is all zeros and <wantgain> is 1 or 2.
% 2010/11/28 - add <wantgain>==2 case
% 2010/11/23 - changed the output range to percentages.  thus, the range is (-Inf,100].
%              also, we removed the <wantr> input since it was dumb.

% input
if ~exist('dim','var') || isempty(dim)
  if isvector(x) && size(x,1)==1
    dim = 2;
  else
    dim = 1;
  end
end
if ~exist('wantgain','var') || isempty(wantgain)
  wantgain = 0;
end
if ~exist('wantmeansub','var') || isempty(wantmeansub)
  wantmeansub = 1;
end
if ~exist('wantsafe','var') || isempty(wantsafe)
  wantsafe = 1;
end

% handle weird case up front
if isempty(x)
  f = [];
  return;
end

% handle 0 case
if dim==0
  x = x(:);
  y = y(:);
  dim = 1;
end

% handle gain
if wantgain
  % to get the residuals, we want to do something like y-x*inv(x'*x)*x'*y
  temp = 1./dot(x,x,dim) .* dot(x,y,dim);
  if wantgain==2
    temp(temp < 0) = 0;  % if the gain was going to be negative, rectify it to 0.
  end
  x = bsxfun(@times,x,temp);
end

% propagate NaNs (i.e. ignore invalid data points)
if wantsafe
  x(isnan(y)) = NaN;
  y(isnan(x)) = NaN;
end

% handle mean subtraction
if wantmeansub
  mn = nanmean(y,dim);
  y = bsxfun(@minus,y,mn);
  x = bsxfun(@minus,x,mn);
end

% finally, compute it
f = 100*(1 - zerodiv(nansum((y-x).^2,dim),nansum(y.^2,dim),NaN,0));

% JUNK:
%
% % <wantr> (optional) is whether to apply signedarraypower(f,0.5)
% %   at the very end, giving something like a correlation coefficient (r).
% %   default: 1.
%
% if ~exist('wantr','var') || isempty(wantr)
%   wantr = 1;
% end
%
% if wantr
%   f = signedarraypower(f,0.5);
% end

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function [f,h] = extractPole(f)
%EXTRACTPOLE  Removes from F the term accounting for non-zero pole/origin
%
% [G,H] = EXTRACTPOLE(F) if F is non-zero at the origin, then this function
% returns DISKFUNs G and H, such that H is rank 1, G = F - H, and G is
% zero at the origin. If F is zero at the origin then G is an empty
% DISKFUN.

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

if ( ~isempty(f) && f.nonZeroPoles )
    h = f;
    % The col and row for the pole is always assumed to be at the first
    % index
    h.cols = f.cols(:,1);
    h.rows = f.rows(:,1);
    h.pivotValues = f.pivotValues(1);
    h.idxPlus = 1;
    h.idxMinus = [];
    h.pivotLocations = f.pivotLocations(1,:);

    % Now remove these rows and columns from f
    if size(f.cols,2) > 1
        f.cols = f.cols(:,2:end);
        f.rows = f.rows(:,2:end);
    else     % Bug in chebfun requires checking for these cases.
        f.cols = [];
        f.rows = [];
    end
    f.pivotValues = f.pivotValues(2:end);
    f.idxPlus = f.idxPlus(2:end)-1;
    f.idxMinus = f.idxMinus - 1;
    f.pivotLocations = f.pivotLocations(2:end,:);
    f.nonZeroPoles = 0;
else
    h = diskfun([]);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [m,n] = size(a,dim)
%SIZE         Implements  size(a)  for slope
%
%   [m,n] = size(a,dim)
%
% functionality as Matlab function size, second argument dim optional
%

% written  12/06/98     S.M. Rump
% modified 09/28/01     S.M. Rump  matrices and multi-dimensional arrays
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargout==2
    [m,n] = size(zeros(a.size));
  elseif nargin==1
    m = a.size;
  else
    m = size(zeros(a.size),dim);
  end

repo: OpenMAS


text:
function m=ccwarpf(f,n,s)
%CCWARPF  Warp cepstral coefficients M=(F,N,S)
% f(1) is the original sample freq, f(2) is the new sample freq
% n(1) is the original number of coefficients, n(2) is the new number
% s is a string: s(1),s(2) =l for linear, m for mel frequency, use capitals if c0 included

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: ccwarpf.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3
   s='ll';
end
if length(f)<2
   f(2)=1;
end
if length(n)<2
   n(2)=n(1);
end
z=s<'a';
s=s+32*z;
if all(s=='l')
   k=1:n(2)-z(2);
   ff=((1:n(1)).'-z(1))*f(2)/f(1);
   fa=2*sin(ff*pi).*ff/pi;
   fb=ff.^2;
   ka=1-2*rem(k,2);
   kb=k.^2;
   r1=ones(n(1),1);
   c1=ones(1,n(2)-z(2));
   a=fa(:,c1).*ka(r1,:);
   b=fb(:,c1)-kb(r1,:);
   f0=find(fix(ff)==ff);
   if length(f0)
      a(f0,:)=ff(f0,c1)==k(ones(length(f0),1),:);
      b(f0,:)=1;
   end
   m=a./b;
   if z(2)
      m=[[1; 0.5*fa(2:n(1))./fb(2:n(1))] m];
   end
end

author: decouples


text:
function [Results, OPTIONS, HeadModel] = bst_wmne_mosher(HeadModel,OPTIONS)
% BST_WMNE: Compute the whitened and weighted minimum-norm operator (wMNE imaging kernel)
%
% USAGE:  [Results,OPTIONS] = bst_wmne(HeadModel, OPTIONS) : Compute mininum operator
%                   OPTIONS = bst_wmne()                   : Return default options
%
% DESCRIPTION:
%     This program computes the whitened and weighted minimum-norm operator,
%     (the wMNE imaging kernel), which is used to compute whitened
%     and weighted minimum-norm estimates (MNE).
%         (e.g., J=wMNEoperator*B; where B is the unwhitened data).
%     It can also compute the whitened and noise-normalized dynamic
%     statistical parametric mapping (dSPM) inverse operator, and/or the
%     whitened standardized low resolution brain electromagnetic tomography
%     (sLORETA) inverse operator, which are used to compute whitened source
%     activity dSPM and sLORETA maps.
%         (e.g., S_dSPM=dSPMoperator*B; where B is the unwhitened data).
%         (e.g., S_sLORETA=sLORETAoperator*B; where B is the unwhitened data).
%
%     The function was written with the goal of providing some of the same
%     functionality of the MNE software written by Matti Hamalainen, but no
%     guarantees are made that it performs all computations in the same
%     exact way. It also provides some functionalities not available in the
%     MNE software.
%
% INPUTS:
%    - HeadModel: Array of Brainstorm head model structures
%         |- Gain       : Forward field matrix for all the channels (unconstrained source orientations)
%         |- GridOrient : Dipole orientation matrix
%         |- area       : Vector with the areas (or possibly volumes) associated with the vertices of the source space.
%    - OPTIONS: structure
%         |- NoiseCov      : NoiseCov is the noise covariance matrix.
%         |- ChannelTypes  : Type of each channel (for each row of the Leadfield and the NoiseCov matrix)
%         |- InverseMethod : {'wmne', 'dspm', 'sloreta'}
%         |- SourceOrient  : String or a cell array of strings specifying the type of orientation constraints for each HeadModel (default: 'fixed')
%         |- SNR        : Signal-to noise ratio defined as in MNE (default: 3).
%         |- diagnoise  : Flag to discard off-diagonal elements of NoiseCov (assuming heteroscedastic uncorrelated noise) (default: 0)
%         |- loose      : Value that weights the source variances of the dipole components defining the tangent space of the cortical surfaces (default: []).
%         |- depth      : Flag to do depth weighting (default: 1).
%         |- weightexp  : Order of the depth weighting. {0=no, 1=full normalization, default=0.8}
%         |- weightlimit: Maximal amount depth weighting (default: 10).
%         |- magreg     : Amount of regularization of the magnetometer noise covariance matrix
%         |- gradreg    : Amount of regularization of the gradiometer noise covariance matrix.
%         |- eegreg     : Amount of regularization of the EEG noise covariance matrix.
%         |- ecogreg    : Amount of regularization of the ECOG noise covariance matrix.
%         |- seegreg    : Amount of regularization of the SEEG noise covariance matrix.
%         |- fMRI       : Vector of fMRI values are the source points.
%         |- fMRIthresh : fMRI threshold. The source variances of source points with OPTIONS.fMRI smaller
%         |               than fMRIthresh will be multiplied by OPTIONS.fMRIoff.
%         |- fMRIoff    : Weight assigned to non-active source points according to fMRI and fMRIthresh.
%
% OUTPUTS:
%    - Results : structure with the wMNE inverse operator and possibly the dSPM
%                and/or sLORETA inverse operators, and other information:

% NOTES:
%     - More leadfield matrices can be used: the solution will combine all
%       leadfield matrices appropriately.
%
%     - This leadfield structure allows to combine surface and volume
%       source spaces with and without dipole orientation constraints,and
%       with and without area or volumetric current density computations.
%       If using a single sphere headmodel, the silent radial
%       component could be eliminated using the SVD (e.g., use
%       bst_remove_silent.m).
%
%     - NoisCov: This should be computed from the pre-stimulus period for
%       averaged ERF data (e.g., using MNE), or from an empty room recording
%       for unaveraged spontaneous or induced data.
%
%     - Orientation constrains for dipoles (.SourceOrient field)
%         - "fixed"     : Dipoles constrained to point normal to cortical surfaces
%         - "free"      : No constraints, dipoles point in x,y,z directions
%         - "loose"     : Source variances of dipole components pointing tangentially to the cortical surfaces are multipled by OPTIONS.loose
%         - "truncated" : An SVD of the gain matrix for each source point is
%                         used to remove the dipole component with least variance, which for
%                         the Single Sphere Head Model, corresponds to the radialsilent component).
%       => For dealing with multiple source spaces with different types of orientation constraints use for example,
%          OPTION.SourceOrient{1}='fixed';
%          OPTION.SourceOrient{2}='free';
%          This has to correspond with the HeadModel Structure
%
%     - The output HeadModel param is used here in return to save LOTS of memory in the bst_wmne function,
%       event if it seems to be absolutely useless. Having a parameter in both input and output have the
%       effect in Matlab of passing them "by referece". So please, do NOT remove it from the function description
%
%     - sLORETA: Output values are multiplied by 1e12 for display in Brainstorm (time series and cortical maps).

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Copyright (C) 2010 - Rey Rene Ramirez
%
% Authors:  Rey Rene Ramirez, Ph.D.   e-mail: rrramirez at mcw.edu
%           Francois Tadel, 2010-2013
%           John Mosher, 2013

%% ===== DEFINE DEFAULT OPTIONS =====
Def_OPTIONS.NoiseCov    = [];
Def_OPTIONS.InverseMethod = 'wmne';
Def_OPTIONS.SNR         = 3;
Def_OPTIONS.diagnoise   = 0;
%Def_OPTIONS.SourceOrient= {'free'};
Def_OPTIONS.SourceOrient= {'fixed'};
Def_OPTIONS.loose       = 0.2;
Def_OPTIONS.depth       = 1;
Def_OPTIONS.weightexp   = 0.5;
Def_OPTIONS.weightlimit = 10;
Def_OPTIONS.regnoise    = 1;
Def_OPTIONS.magreg      = .1;
Def_OPTIONS.gradreg     = .1;
Def_OPTIONS.eegreg      = .1;
Def_OPTIONS.ecogreg     = .1;
Def_OPTIONS.seegreg     = .1;
Def_OPTIONS.fMRI        = [];
Def_OPTIONS.fMRIthresh  = [];
Def_OPTIONS.fMRIoff     = 0.1;
Def_OPTIONS.pca         = 1;
% Return the default options
if (nargin == 0)
    Results = Def_OPTIONS;
    return
end
% Make the default for all the leadfields
numL = size(HeadModel,2);
Def_OPTIONS.SourceOrient = repmat(Def_OPTIONS.SourceOrient, [1 numL]);
% Copy default options to OPTIONS structure (do not replace defined values)
OPTIONS = struct_copy_fields(OPTIONS, Def_OPTIONS, 0);

%% ===== CHECK FOR INVALID VALUES =====
disp(' ');
% Detect if the input noise covariance matrix is or should be diagonal
C_noise = OPTIONS.NoiseCov;
variances = diag(C_noise);
if isequal(C_noise, diag(variances))
    OPTIONS.diagnoise = 1;
    disp('wMNE> Detected diagonal noise covariance: setting diagnoise to 1');
end
% If OPTIONS.diagnoise is 1, then OPTIONS.pca=0
if OPTIONS.diagnoise
    OPTIONS.pca=0;
    disp('wMNE> If using diagonal noise covariance, PCA option should be off. Setting PCA option off.')
end
if isempty(OPTIONS.NoiseCov)
    error('You need to input the noise covariance in the NoiseCov field of OPTIONS.')
end
if (numL ~= length(OPTIONS.SourceOrient))
    error('The number of elements in the HeadModel structure should equal the length of the cell array OPTIONS.SourceOrient.')
end
if ~isempty(OPTIONS.loose) && (OPTIONS.loose>=1 || OPTIONS.loose<=0)
    error('loose value should be smaller than 1 and bigger than 0, or empty for no loose orientations.')
end
if OPTIONS.weightexp>1 || OPTIONS.weightexp<0
    error('weightexp should be a scalar between 0 and 1')
end
if OPTIONS.magreg>1 || OPTIONS.magreg<0
    error('magreg should be a scalar between 0 and 1')
end
if OPTIONS.eegreg>1 || OPTIONS.eegreg<0
    error('eegreg should be a scalar between 0 and 1')
end
if OPTIONS.ecogreg>1 || OPTIONS.ecogreg<0
    error('ecogreg should be a scalar between 0 and 1')
end
if OPTIONS.seegreg>1 || OPTIONS.seegreg<0
    error('seegreg should be a scalar between 0 and 1')
end

%% ===== NOISE COVARIANCE RANK =====
% Get indices of MEG and EEG channels
iMeg  = find(strncmpi(OPTIONS.ChannelTypes,'MEG',3));
iEeg  = find(strncmpi(OPTIONS.ChannelTypes,'EEG',3));
iEcog = find(strncmpi(OPTIONS.ChannelTypes,'ECOG',3));
iSeeg = find(strncmpi(OPTIONS.ChannelTypes,'SEEG',3));
% Diagonal noisecov
if OPTIONS.diagnoise
    C_noise = diag(variances);
    rnkC_noise_meg = length(iMeg);
    rnkC_noise_eeg = length(iEeg);
    rnkC_noise_ecog = length(iEcog);
    rnkC_noise_seeg = length(iSeeg);
    disp('wMNE> Setting off diagonal elements of the noise covariance to zero.');
    disp(['wMNE> Rank of noise covariance is ' num2str(size(C_noise,1))]);
% Full noisecov
else
    % Estimate noise covariance matrix rank separately for sensor types
    if ~isempty(iMeg)
        rnkC_noise_meg = rank(single(C_noise(iMeg,iMeg))); % Rey added this. Separate rank of MEG. 3/23/11
        disp(['wMNE> Rank of MEG part of noise covariance is ' num2str(rnkC_noise_meg)]);
    end
    if ~isempty(iEeg)
        rnkC_noise_eeg = rank(single(C_noise(iEeg,iEeg))); % Rey added this. Separate rank of EEG. 3/23/11
        disp(['wMNE> Rank of EEG part of noise covariance is ' num2str(rnkC_noise_eeg)]);
    end
    if ~isempty(iEcog)
        rnkC_noise_ecog = rank(single(C_noise(iEcog,iEcog)));  % FT added 21-Feb-13
        disp(['wMNE> Rank of ECOG part of noise covariance is ' num2str(rnkC_noise_ecog)]);
    end
    if ~isempty(iSeeg)
        rnkC_noise_seeg = rank(single(C_noise(iSeeg,iSeeg)));  % FT added 21-Feb-13
        disp(['wMNE> Rank of SEEG part of noise covariance is ' num2str(rnkC_noise_seeg)]);
    end
    % Sets off-diagonal terms to zero. Rey added this. 3/23/11
    C_noise_new = 0 * C_noise;
    if ~isempty(iMeg)
        C_noise_new(iMeg,iMeg) = C_noise(iMeg,iMeg);
    end
    if ~isempty(iEeg)
        C_noise_new(iEeg,iEeg) = C_noise(iEeg,iEeg);
    end
    if ~isempty(iEcog)
        C_noise_new(iEcog,iEcog) = C_noise(iEcog,iEcog);
    end
    if ~isempty(iSeeg)
        C_noise_new(iSeeg,iSeeg) = C_noise(iSeeg,iSeeg);
    end
    C_noise = C_noise_new;
end

%% ===== REGULARIZE NOISE COVARIANCE MATRIX =====
% Only if option is selected
if OPTIONS.regnoise
    listTypes = unique(OPTIONS.ChannelTypes);
    % Loop on all the required data types (MEG MAG, MEG GRAD, EEG)
    for iType = 1:length(listTypes)
        % Get channel indices
        iChan = find(strcmpi(OPTIONS.ChannelTypes, listTypes{iType}));
        % Regularize noise covariance matrix
        switch listTypes{iType}
            case 'MEG GRAD', reg = OPTIONS.gradreg;
            case 'MEG MAG',  reg = OPTIONS.magreg;
            case 'MEG',      reg = OPTIONS.gradreg;
            case 'EEG',      reg = OPTIONS.eegreg;
            case 'ECOG',     reg = OPTIONS.ecogreg;
            case 'SEEG',     reg = OPTIONS.seegreg;
        end
        % Original Line 4/5/13:
        % C_noise(iChan,iChan) = C_noise(iChan,iChan) + (reg * mean(variances(iChan)) * eye(length(iChan)));
        % JCM 4/5/13, mods just to be clear
        % mean of the diagonal variances
        % Options could be median, maximum, minimum, etc. Note, this is
        % not the eigenspectrum, but the homoskedastic spectrum.
        % TODO try other forms of matrix norms for noise regualarization
        LAMBDA_REGULARIZER = reg * mean(variances(iChan));
        % Now add this Tikhonov regularizer to the noise diagonal.
        C_noise(iChan,iChan) = C_noise(iChan,iChan) + diag(zeros(length(iChan),1) + LAMBDA_REGULARIZER);

    end
end

%% ===== WHITENING OPERATOR =====
% Rey added all of this, 3/23/11
% Modified FT 21-Feb-2013
% Whitening of each modality separately (MEG,EEG,ECOG,SEEG), which assumes
% zero covariance between them (i.e., a block diagonal noise covariance). This
% was recommended by Matti as EEG does not measure all the signals from the same
% environmental noise sources as MEG.
nChan = size(C_noise,2);
W = zeros(0, nChan);
if ~isempty(iMeg)
    W_meg = CalculateWhitener('MEG', C_noise, iMeg, rnkC_noise_meg, OPTIONS.pca);
    W_tmp = zeros(size(W_meg,1), nChan);
    W_tmp(:,iMeg) = W_meg;
    W = [W; W_tmp];
end
if ~isempty(iEeg)
    W_eeg = CalculateWhitener('EEG', C_noise, iEeg, rnkC_noise_eeg, OPTIONS.pca);
    W_tmp = zeros(size(W_eeg,1), nChan);
    W_tmp(:,iEeg) = W_eeg;
    W = [W; W_tmp];
end
if ~isempty(iEcog)
    W_ecog = CalculateWhitener('ECOG', C_noise, iEcog, rnkC_noise_ecog, OPTIONS.pca);
    W_tmp = zeros(size(W_ecog,1), nChan);
    W_tmp(:,iEcog) = W_ecog;
    W = [W; W_tmp];
end
if ~isempty(iSeeg)
    W_seeg = CalculateWhitener('SEEG', C_noise, iSeeg, rnkC_noise_seeg, OPTIONS.pca);
    W_tmp = zeros(size(W_seeg,1), nChan);
    W_tmp(:,iSeeg) = W_seeg;
    W = [W; W_tmp];
end
% Check for whitener integrity
if any(isnan(W(:))) || any(isinf(W(:)))
    error('Invalid noise covariance matrix.')
end
% Display rank of the whitener
rnkC_noise = size(W,1);
display(['wMNE> Total rank is ' num2str(rnkC_noise) '.'])

%% ===== PROCESSING LEAD FIELD MATRICES, WEIGHTS, AREAS & ORIENTATIONS =====
% Initializing.
spl = zeros(numL,1);
numdipcomp = spl;
for k = 1:numL
    sL = size(HeadModel(k).Gain, 2);
    switch OPTIONS.SourceOrient{k}
        case 'fixed',      numdipcomp(k)=1;
        case 'free',       numdipcomp(k)=3;
        case 'loose',      numdipcomp(k)=3;
        case 'truncated',  numdipcomp(k)=2;
    end
    spl(k) = (sL / 3) * numdipcomp(k); % This is a vector with the total number of dipole components per source space.
end
sspl = sum(spl); % This is the total number of dipole components across all source spaces.
L = zeros(rnkC_noise,sspl);
w = ones(sspl,1);
if isfield(HeadModel, 'area')
    areas = w;
else
    areas = [];
end
itangential = [];
start = 0;
Q_Cortex = [];
for k = 1:numL
    start = start + 1;
    endd = start + spl(k) - 1;
    Lk = HeadModel(k).Gain;
    HeadModel(k).Gain = [];
    %% ===== COMPUTE POWER =====
    szL = size(Lk);
    if OPTIONS.depth
        display(['wMNE> Computing power of gain matrices at each source point for source space ' num2str(k) '.'])
        % Computing power
        % JCM 4/5/2013, this is squared Frobenius norm of the source
        % Options would be the matrix norm (largest singular value, or
        % squared).
        wk = squeeze(sum(sum((reshape(Lk,[szL(1) 3 szL(2)/3])) .^2,1),2));
        wk = repmat(wk',[numdipcomp(k) 1]);
        wk = reshape(wk,[spl(k) 1]);
        w(start:endd) = wk;
        clear wk
    end
    switch OPTIONS.SourceOrient{k}
        case 'fixed'
            display('wMNE> Appying fixed dipole orientations.')
            Lk = bst_gain_orient(Lk,HeadModel(k).GridOrient);
        case 'free'
            display('wMNE> Using free dipole orientations. No constraints.')
        case 'loose'
            display('wMNE> Transforming lead field matrix to cortical coordinate system.')
            [Lk, Q_Cortex] = bst_xyz2lf(Lk, HeadModel(k).GridOrient');
            % Getting indices for tangential dipoles.
            itangentialtmp = start:endd;
            itangentialtmp(1:3:end) = [];
            itangential = [itangential itangentialtmp];  %#ok<AGROW>
        case 'truncated'
            display('wMNE> Truncating the dipole component pointing in the direction with least variance (i.e., silent component for single sphere head model.')
            [Lk, Q_Cortex] = bst_remove_silent(Lk);
    end
    %% ===== WHITEN LEAD FIELD MATRIX =====
    % Whiten lead field.
    display(['wMNE> Whitening lead field matrix for source space ' num2str(k) '.'])
    Lk = W * Lk;
    if isfield(HeadModel(k),'area') && ~isempty(HeadModel(k).area)
        areav = HeadModel(k).area;
        areav = repmat(areav', [numdipcomp(k) 1]);
        areav = reshape(areav, [spl(k) 1]);
        areas(start:endd) = areav;
    end
    L(:,start:endd) = Lk;
    start = endd;
end
% Computing reciprocal of power.
w = 1 ./ w;
% Clear memory
clear Lk endd start itangentialtmp sL szL

%% ===== APPLY AREAS =====
if ~isempty(areas)
    display('wMNE> Applying areas to compute current source density.')
    areas = areas.^2;
    w = w .* areas;
end
clear areas

%% ===== APPLY DEPTH WEIGHTHING =====
if OPTIONS.depth
    % ===== APPLY WEIGHT LIMIT =====
    % Applying weight limit.
    display('wMNE> Applying weight limit.')
    weightlimit2 = OPTIONS.weightlimit .^ 2;
    %limit=min(w(w>min(w)*weightlimit2));  % This is the Matti way.
    limit = min(w) * weightlimit2;  % This is the Rey way (robust to possible weight discontinuity).
    w(w>limit) = limit; %JCM note, 4/5/2013, w = min(w,limit);

    % ===== APPLY WEIGHT EXPONENT =====
    % Applying weight exponent.
    display('wMNE> Applying weight exponent.')
    w = w .^ OPTIONS.weightexp;
    clear limit weightlimit2
end

%% ===== APPLY LOOSE ORIENTATIONS =====
if ~isempty(itangential)
    display(['wMNE> Applying loose dipole orientations. Loose value of ' num2str(OPTIONS.loose) '.'])
    w(itangential) = w(itangential) * (OPTIONS.loose);
end

%% ===== APPLY fMRI PRIORS =====
% Apply fMRI Priors
if ~isempty(OPTIONS.fMRI)
    display('wMNE> Applying fMRI priors.')
    ifmri = (OPTIONS.fMRI < OPTIONS.fMRIthresh);
    w(ifmri) = w(ifmri) * OPTIONS.fMRIoff;
end

%% ===== ADJUSTING SOURCE COVARIANCE MATRIX =====
% Adjusting Source Covariance matrix to make trace of L*C_J*L' equal to number of sensors.
% JCM 4/5/2013, i.e. average signal covariance is one, so trace is number
% of sensors.
display('wMNE> Adjusting source covariance matrix.')
C_J = speye(sspl, sspl);
C_J = spdiags(w, 0, C_J);
trclcl = trace(L * C_J * L');
C_J = C_J * (rnkC_noise / trclcl);
Rc = chol(C_J);
LW = L * Rc;
clear C_J trclcl sspl itangential rnkC_noise

%% BEGIN SOURCE MODELING ROUTINES

%% JCM 4/5/2013 GLS and other routines, using above setup

if any(strcmpi(OPTIONS.InverseMethod, {'gls','gls_p','glsr','glsr_p'})),

   start = 0;
   Kernel = zeros(size(LW,2),size(LW,1)); % note transpose of LW
   lambda2 = OPTIONS.SNR^(-2); % for regularizing the GLS

   for k = 1:numL,
      start = start + 1;
      endd = start + spl(k) - 1;
      Lk = LW(:,start:endd); % this whitened gain matrix
      Rck = Rc(start:endd,start:endd); % covariance priors
      % Now process each source in the gain matrix for it's own inversion
      NumSources = size(Lk,2)/numdipcomp(k); % total number of sources

      for i = 1:NumSources,
         ndx = ((1-numdipcomp(k)):0)+i*numdipcomp(k); % next index
         % SVD the source
         [Ua,Sa,Va] = svd(Lk(:,ndx),0); % svd of this source
         Sa = diag(Sa);
         Tolerance_Source = size(Lk,1)*eps(single(Sa(1)));
         Rank_Source = sum(Sa > Tolerance_Source);
         % Trim decomposition
         Ua = Ua(:,1:Rank_Source);
         Sa = Sa(1:Rank_Source);
         Va = Va(:,1:Rank_Source);

         % calculate the weighted subspace
         Reg_Weights_Source = Sa.^2 ./ (Sa.^2 + lambda2*Sa(1)^2); % regularizer
         % now write the pseudoinverse results back into this same matrix

         switch OPTIONS.InverseMethod
            case 'gls'
               % The true pseudo-inverse
               Lk(:,ndx) = Ua * diag(1./Sa) * Va'*Rck(ndx,ndx);
            case 'gls_p'
               % The model performance
               % Model may be reduced rank
               Lk(:,ndx) = 0; % zero out
               Lk(:,ndx(1:size(Ua,2))) = Ua; % model performance
            case 'glsr'
               % Regularized
               Lk(:,ndx) = Ua * diag(Reg_Weights_Source./Sa) * Va'*Rck(ndx,ndx);
            case 'glsr_p'
               % Regularized model performance
               % Model may be reduced rank
               Lk(:,ndx) = 0; % zero out
               Lk(:,ndx(1:size(Ua,2))) = Ua * diag(sqrt(Reg_Weights_Source)); % model performance
            otherwise
               error('Bad Options String %s',OPTIONS.InverseMethod)
         end
      end

      % Now we have a matrix almost ready for use as an imaging kernel
      % Later, below, the whitener will be added

      Kernel(start:endd,:) = Lk';
      start = endd;
   end

end % JCM GLS routine

if any(strcmpi(OPTIONS.InverseMethod, {'mnej','mnej_p'})), %JCM Min Norm
   % mnej should be identical to Rey, but mnej_p is novel

   start = 0;
   Kernel = zeros(size(LW,2),size(LW,1)); % note transpose of LW
   lambda2 = OPTIONS.SNR^(-2); % for regularizing the MN

   for k = 1:numL,
      start = start + 1;
      endd = start + spl(k) - 1;
      Lk = LW(:,start:endd); % this whitened gain matrix
      Rck = Rc(start:endd,start:endd); % covariance priors

      % Setup the Min Norm
      % First, generate the population data covariance
      wCD = LW*LW' + (diag(zeros(size(LW,1),1) + lambda2)); % whitened data covariance

      % Decompose
      [Ud,Sd] = svd(wCD);

      % Data Whitener
      iWd = Ud*diag(1./sqrt(diag(Sd)))*Ud';

      % Now process each source in the gain matrix for it's own inversion
      NumSources = size(Lk,2)/numdipcomp(k); % total number of sources

      for i = 1:NumSources,
         ndx = ((1-numdipcomp(k)):0)+i*numdipcomp(k); % next index
         % SVD the data whitened source
         [Ua,Sa,Va] = svd(iWd*Lk(:,ndx),0); % svd of this source
         Sa = diag(Sa);
         Tolerance_Source = size(Lk,1)*eps(single(Sa(1)));
         Rank_Source = sum(Sa > Tolerance_Source);
         if Rank_Source < length(Sa),
            fprintf('%.0f ',i);
         end
         % Trim decomposition
         Ua = Ua(:,1:Rank_Source);
         Sa = Sa(1:Rank_Source);
         Va = Va(:,1:Rank_Source);

         SNR_Weights_Source = Sa.^2 ./ (1 - Sa.^2); % SNR
         % now write the pseudoinverse results back into this same matrix

         switch OPTIONS.InverseMethod
            case 'mnej'
               % The true solution
               Lk(:,ndx) = iWd * Ua * diag(Sa) * Va' * Rck(ndx,ndx);
            case 'mnej_p'
               % The model performance
               %Lk(:,ndx) = iWd * Ua * diag(sqrt(SNR_Weights_Source)); % model performance
               Lk(:,ndx) = iWd * Ua; % CHEAT model performance

            otherwise
               error('Bad Options String %s',OPTIONS.InverseMethod)
         end
      end

      % Now we have a matrix almost ready for use as an imaging kernel
      % Later, below, the whitener will be added

      Kernel(start:endd,:) = Lk';
      start = endd;
   end

end %JCM Min norms

%% The dSPM and sLORETA functions rely on 'wmne' being calculated

if any(strcmpi(OPTIONS.InverseMethod, {'wmne','dspm','sloreta'}))

   %% ===== SINGULAR VALUE DECOMPOSITION =====
   % Set regularization parameter based on SNR
   lambda2 = OPTIONS.SNR ^ (-2);
   % Compute SVD.
   display('wMNE> Computing SVD of whitened and weighted lead field matrix.')
   [V,S,U] = svd(LW','econ'); % JCM 4/5/2013 transpose for greater speed
   s = diag(S);
   ss = s ./ (s.^2 + lambda2);
   clear LW lambda2 s

   %% ===== WHITENED MNE IMAGING KERNEL =====
   % Compute whitened MNE operator.
   Kernel = Rc * V * diag(ss) * U';
   clear Rc V ss U
end

%% ===== WHITENED dSPM IMAGING KERNEL =====
% Compute dSPM operator.
if strcmpi(OPTIONS.InverseMethod, 'dspm')
   display('wMNE> Computing dSPM inverse operator.')
   start = 0;
   for k = 1:numL
      start = start+1;
      endd = start + spl(k) - 1;
      dspmdiag = sum(Kernel(start:endd,:) .^2, 2);
      if (numdipcomp(k) == 1)
         dspmdiag = sqrt(dspmdiag);
      elseif (numdipcomp(k)==3 || numdipcomp(k)==2)
         dspmdiag = reshape(dspmdiag, [numdipcomp(k), spl(k)/numdipcomp(k)]);
         dspmdiag = sqrt(sum(dspmdiag)); % Taking trace and sqrt.
         dspmdiag = repmat(dspmdiag, [numdipcomp(k), 1]);
         dspmdiag = reshape(dspmdiag, [spl(k), 1]);
      end
      Kernel(start:endd,:) = bst_bsxfun(@rdivide, Kernel(start:endd,:), dspmdiag);
      start = endd;
   end

   %% ===== WHITENED sLORETA IMAGING KERNEL =====
   % Compute sLORETA operator.
elseif strcmpi(OPTIONS.InverseMethod, 'sloreta')
   display('wMNE> Computing sLORETA inverse operator.')
   start=0;
   for k = 1:numL
      start = start + 1;
      endd = start + spl(k) - 1;
      if (numdipcomp(k) == 1)
         sloretadiag = sqrt(sum(Kernel(start:endd,:) .* L(:,start:endd)', 2));
         Kernel(start:endd,:) = bst_bsxfun(@rdivide, Kernel(start:endd,:), sloretadiag);
      elseif (numdipcomp(k)==3 || numdipcomp(k)==2)
         for spoint = start:numdipcomp(k):endd
            R = Kernel(spoint:spoint+numdipcomp(k)-1,:) * L(:,spoint:spoint+numdipcomp(k)-1);
            SIR = sqrtm(pinv(R));
            Kernel(spoint:spoint+numdipcomp(k)-1,:) = SIR * Kernel(spoint:spoint+numdipcomp(k)-1,:);
         end
      end
      start=endd;
   end
end
disp(' ');

%% JCM, 4/5/2013, this loose orientation may be compatible with GLS and others

%% ===== LOOSE ORIENTATION: RE-ORIENT COMPONENTS =====
% WARNING: Changing the orientations of the dipoles from MNE to Brainstorm
% => Make "Unconstrained" and "loose"/"truncated" models directly comparable
if ~isempty(Q_Cortex)
   % Creating a block diagonal matrix
   N = size(Kernel,1);
   Nout = N / numdipcomp * 3;
   iRow = reshape(repmat(reshape(1:Nout,3,[]), numdipcomp, 1), 1, []);
   iCol = reshape(repmat(1:N,3,[]), 1, []);
   Q_Cortex = sparse(iRow, iCol, Q_Cortex(:));
   % Applying orientations
   Kernel = Q_Cortex * Kernel;
end

%% ===== ASSIGN IMAGING KERNEL =====
% Multiply inverse operator and whitening matrix, so no need to whiten data.
Kernel = Kernel * W;
% Return results structure
Results.ImagingKernel = Kernel;
Results.ImageGridAmp  = [];
Results.Whitener      = W;
if (length(numdipcomp) > 1)
    Results.nComponents = 0;
else
    Results.nComponents = numdipcomp;
end

end

%% ==============================================================================
%  ===== HELPER FUNCTIONS =======================================================
%  ==============================================================================
% TODO JCM 4/5/2013: Cleanup the whitener to use SVD, reduced rank. isPCA
% is actually for truncated SVD of the matrix for small values.
function W = CalculateWhitener(Modality, C_noise, iChannel, rnkC_noise, isPca)
    [V,D] = eig(C_noise(iChannel,iChannel));
    D = diag(D);
    [D,I] = sort(D,'descend');
    V = V(:,I);
    % No PCA case.
    if ~isPca
        display(['wMNE> Not doing PCA for ' Modality '.'])
        D = 1./D;
        W = diag(sqrt(D)) * V';
    % Rey's approach. MNE has been changed to implement this.
    else
        display(['wMNE> Setting small ' Modality ' eigenvalues to zero.'])
        D = 1 ./ D;
        D(rnkC_noise+1:end) = 0;
        W = diag(sqrt(D)) * V';
        W = W(1:rnkC_noise,:); % This line will reduce the actual number of variables in data and leadfield to the true rank. This was not done in the original MNE C code.
    end
end

author: brainstorm-tools


text:
% Derive_pwPoly5.m
%
% This script derives the equations that give the bezier coefficients for a
% curve, given the value, slope, and curvature at each endpoint.
%

clc; clear;

% Control Points:
n = 5;  %order of bezier curve
p = sym('p',[n+1,1]); assume(p,'real');

% Independent variable:
syms t 'real'

% Build up the equation for a 5th-order bezier curve:
x = sym(0);
for i=0:n
    tt = (t.^i).*(1-t).^(n-i);
    binom = nchoosek(n,i);
    x = x + binom*p(i+1)*tt;
end

% Now compute the derivaives:
dx = diff(x,t);
ddx = diff(dx,t);

% Evaluate at endpoints:
X0 = simplify(subs(x,t,0));
dX0 = simplify(subs(dx,t,0));
ddX0 = simplify(subs(ddx,t,0));
X1 = simplify(subs(x,t,1));
dX1 = simplify(subs(dx,t,1));
ddX1 = simplify(subs(ddx,t,1));

% User-Inputs:
syms x0 dx0 ddx0 x1 dx1 ddx1 'real'

% Constraint equations:
eqns = [...
    X0-x0;
    dX0-dx0;
    ddX0-ddx0;
    X1-x1;
    dX1-dx1;
    ddX1-ddx1];

% Solve for unknowns
[A,b] = equationsToMatrix(eqns,p);
pSoln = A\b;

% Substitute solution:
xSoln = simplify(subs(x,p,pSoln));
dxSoln = simplify(subs(dx,p,pSoln));
ddxSoln = simplify(subs(ddx,p,pSoln));

% Write solution to a matlab function:
matlabFunction(xSoln, dxSoln, ddxSoln,...
    'file','autoGen_pwPoly5.m',...
    'vars',{t,x0 dx0 ddx0 x1 dx1 ddx1},...
    'outputs',{'x','dx','ddx'});

author: MatthewPeterKelly

repo: dscTutorials


text:
%This Matlab script can be used to reproduce Figure 7.27 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.1 (Last edited: 2018-10-10)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%% Define parameters
S = [1, 16, 256]; %Number of subarrays
M = 256./S; %Number of antennas per subarray
d0 = 10; %Guard distance to the closest subarray
W = 350; %Cell radii
numberOfSetups = 5000; %Number of random subarray drops

%Propagation model parameters as in (2.3)
alpha = 3.76;
Upsilon = -148.1;
sigma_sf = 8;
pathlossdB = @(d) Upsilon-10*alpha*log10(d/1000);
addShadowing = @(pldB) 10.^((pldB +  sigma_sf*randn(size(pldB)))/10);

%% Monte Carlo simulations
betaValues = zeros(numel(S), numberOfSetups);

%Go through the number of subarrays
for s = 1:numel(S)

    %Go through all drops of subarrays
    for it1 = 1:numberOfSetups

        d = sort(sqrt(rand(S(s),1)*(W^2-d0^2) + d0^2)); %The squared distances of uniformly distributed points between two circles are uniformly distributed
        pldB = pathlossdB(d); %Compute pathloss without shadow fading
        pl = addShadowing(pldB); %Add random shadow fading

        %Compute average channel gain
        betaValues(s,it1) = sum(pl)/S(s);

    end

end

%% Plot the simulation results
figure;
hold on; box on;
plot(sort(10*log10(betaValues(1,:))),linspace(0,1,numberOfSetups), 'LineStyle', '-', 'Color', 'k', 'LineWidth', 1);
plot(sort(10*log10(betaValues(2,:))),linspace(0,1,numberOfSetups), 'LineStyle', '--', 'Color', 'r',  'LineWidth', 1);
plot(sort(10*log10(betaValues(3,:))),linspace(0,1,numberOfSetups), 'LineStyle', '-.', 'Color', 'b',  'LineWidth', 1);

xlabel('Average channel gain [dB]')
ylabel('CDF')
legend({['S=', num2str(S(1))],['S=', num2str(S(2))],['S=', num2str(S(3))]}, 'Location', 'NorthWest')
xlim([-160 -60]);

author: emilbjornson

repo: massivemimobook


text:
function status = test_exp_mono(modus)
%TEST_EXP_MONO tests the different expansions of sound fields in frequency-domain
%
%   Usage: status = test_exp_mono(modus)
%
%   Input parameters:
%       modus   - 0: numerical
%                 1: visual
%
%   Output parameters:
%       status  - true or false
%
%   TEST_EXP_MONO(modus) checks, if the circular basis expansions for plane
%   waves and point sources are working. The circular basis expansions are
%   converted to plane wave decompositions. Additionally, modal weighting
%   functions are tested. Optionally, sound field plots of the plane wave
%   decompositions are used for verification.

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

status = false;

%% ===== Checking of input  parameters ===================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Configuration ===================================================
% Parameters
conf = SFS_config;
conf.plot.loudspeakers = false;  % do not plot loudspeakers
conf.modal_window_parameter = 2.0;  % parameter for kaiser window

c = conf.c;  % speed of sound

% For sound field plots
X = [-2, 2];
Y = [-2, 2];
Z = 0;

f = 500;

% Plane waves with equi-angular distribution
Npw = 1024;
phi0 = (0:Npw-1).'*2*pi/Npw;
x0 = [cos(phi0) sin(phi0)];
x0(:,3) = 0;
x0(:,4:6) = x0(:,1:3);
x0(:,7) = 1;

% Test scenarios
scenarios = { ...
  'pw', [ 0.0 -1.0 0.0],  5, 'rect'   , [0.0 0.0 0.0]
  'pw', [ 0.0 -1.0 0.0], 15, 'rect'   , [0.0 0.0 0.0]
  'pw', [ 0.0 -1.0 0.0], 15, 'kaiser' , [0.0 0.0 0.0]
  'pw', [ 0.0 -1.0 0.0],  5, 'rect'   , [0.5 1.0 0.0]
  'pw', [ 0.0 -1.0 0.0], 15, 'rect'   , [0.5 1.0 0.0]
  'pw', [ 0.0 -1.0 0.0], 15, 'kaiser' , [0.5 1.0 0.0]
  'ps', [ 0.0  2.5 0.0],  5, 'rect'   , [0.0 0.0 0.0]
  'ps', [ 0.0  2.5 0.0], 15, 'rect'   , [0.0 0.0 0.0]
  'ps', [ 0.0  2.5 0.0], 15, 'kaiser' , [0.0 0.0 0.0]
  'ps', [ 0.0  2.5 0.0],  5, 'rect'   , [0.5 1.0 0.0]
  'ps', [ 0.0  2.5 0.0], 15, 'rect'   , [0.5 1.0 0.0]
  'ps', [ 0.0  2.5 0.0], 15, 'kaiser' , [0.5 1.0 0.0]
  };

%% ===== Main ============================================================

for ii=1:size(scenarios)

    src = scenarios{ii,1};  % source type
    xs = scenarios{ii,2};  % source position / direction of plane wave
    Nce = scenarios{ii,3};  % modal order
    conf.modal_window = scenarios{ii,4};  % type of modal weighting function
    xq = scenarios{ii,5};  % expansion centre

    % Circular expansion coefficients
    switch src
    case 'pw'
        Pm = circexp_mono_pw(xs,Nce,f,xq,conf);
        g = 1;
    case 'ps'
        Pm = circexp_mono_ps(xs,Nce,f,xq,conf);
        g = 1./(4*pi*norm(xs-xq));
    end

    % Modal weighting of coefficients
    wm = modal_weighting(Nce, conf);
    Pm = bsxfun(@times, Pm, [wm(end:-1:2) wm]);

    % Conversion to plane wave decomposition
    Ppwd = pwd_mono_circexp(Pm, Npw);

    if modus
        % Sound field plot
        P = sound_field_mono(X,Y,Z,x0,'pw',Ppwd,f,conf);
        plot_sound_field(P.*g,X,Y,Z,[],conf);

        % Title string
        str = 'Plane wave decompostion of modally bandlimited';
        switch src
        case 'pw'
            str = sprintf('%s plane wave', str);
        case 'ps'
            str = sprintf('%s point source', str);
        end
        str = sprintf(['%s ([%1.1f %1.1f %1.1f]):\n%s-window (M=%d), ' ...
          'center of modal expansion at [%1.1f %1.1f %1.1f]'], ...
          str, xs, conf.modal_window, Nce, xq);
        title(str);
    end
end

status = true;

author: sfstoolbox

repo: sfs-matlab


text:
function [ x, seed ] = r4_normal_01 ( seed )

%*****************************************************************************80
%
%% R4_NORMAL_01 returns a unit pseudonormal R4.
%
%  Discussion:
%
%    The standard normal probability distribution function (PDF) has
%    mean 0 and standard deviation 1.
%
%    The Box-Muller method is used, which is efficient, but
%    generates two values at a time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the standard normal PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ r1, seed ] = r4_uniform_01 ( seed );
  [ r2, seed ] = r4_uniform_01 ( seed );
  x = sqrt ( -2.0 * log ( r1 ) ) * cos ( 2.0 * pi * r2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef ErrorMultiClass < dagnn.Loss
% author: Hakan Bilen
% computes multi-class accuracy
% inputs{1}->scores
% inputs{2}->gt labels
  properties
    nImgPerClass = []
    nCorPred = []
    accuracy = []
    resetLayer = false
  end

  methods
    function outputs = forward(obj, inputs, params)

      if numel(inputs)~=2
        error('wrong number of inputs');
      end

      nCls = size(inputs{1},3);

      if obj.resetLayer || isempty(obj.nImgPerClass)
        obj.nImgPerClass = zeros(1,size(inputs{1},3));
        obj.nCorPred = zeros(1,size(inputs{1},3));
        obj.accuracy = zeros(1,size(inputs{1},3));

        if obj.resetLayer
          obj.resetLayer = false ;
          obj.average = 0 ;
        end
      end

      [~,predictions] = max(gather(squeeze(inputs{1})),[],1);

      for c=1:nCls
        obj.nImgPerClass(c) = obj.nImgPerClass(c) + sum(inputs{2}==c);
        obj.nCorPred(c)     = obj.nCorPred(c) + sum(predictions==c & inputs{2}==c);
      end

      ni = obj.nImgPerClass;
      ni(ni==0) = 1;

      obj.accuracy = obj.nCorPred ./ ni;
      obj.average = (1-mean(obj.accuracy));
      outputs{1} =  obj.average;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs = cell(1,2);
      derParams = {} ;
    end

    function reset(obj)
      obj.resetLayer = true ;
%       obj.nImgPerClass = [];
%       obj.nCorPred = [];
%       obj.accuracy = [];
%       obj.average = 0;
    end

    function obj = ErrorMultiClass(varargin)
      obj.load(varargin) ;
      obj.loss = 'error_multi_class' ;
    end
  end
end

author: hbilen

repo: dynamic-image-nets


text:
function pass = test_laplacian( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e7*pref.techPrefs.chebfuneps;

% Test with different parity of m,p
% Example 1:
f = ballfun(@(x,y,z)x.^2+y.^2+z.^2);
V = ballfunv(f,2*f,-f);
g = laplacian(V);
exact = ballfunv(ballfun(@(x,y,z)6),ballfun(@(x,y,z)12),ballfun(@(x,y,z)-6));
pass(1) = norm( g - exact ) < tol;

% Example 2:
f1 = ballfun(@(x,y,z)cos(x.*y));
f2 = ballfun(@(x,y,z)sin(y.*z));
f3 = ballfun(@(x,y,z)z.^3);
v = ballfunv(f1,f2,f3);
g = laplacian(v);
exact = grad(div(v))-curl(curl(v));
pass(2) = norm( g - exact ) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end


text:
 function st = newfft_table_init(st, varargin)
%function st = newfft_table_init(st, varargin)
%|
%| Initialize structure for d-dimension NUFFT using table-based interpolator,
%| This should be called only by newfft for its 'table0' or 'table1' mode!
%| Note: default oversample factor is 2^11 or 2^13
%|
%| in
%|	st		structure initialized in newfft.m
%|			uses st.phasing to choose complex or real kernel
%|
%| option
%|	'how'	char	'slow' | 'ratio' | 'fast'
%|			(default depends on kernel type)
%|
%| out
%|	st	strum	(see newfft)
%|
%| Copyright 2008-7-11, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end

arg.how = '';
arg = vararg_pair(arg, varargin);

if isempty(arg.how)
	if streq(st.ktype, 'kb:', 3) % todo: others?
		arg.how = 'ratio';
	else
		arg.how = 'fast';
	end
end

st.phase_shift = []; % compute on-the-fly

Nd = st.Nd;
Jd = st.Jd;
Kd = st.Kd;
dd = length(Nd);
Ld = st.oversample;
if isempty(Ld)
	switch st.mode
	case 'table0'
		st.order = 0; % 0th order
		Ld = 2^11; % default
	case 'table1'
		st.order = 1; % 1st order
		Ld = 2^9; % default
	otherwise
		fail 'bad mode'
	end
end

% dimensionality of input space (usually 2 or 3)
if dd > 1 && length(Ld) == 1
	Ld = repmat(Ld, [1 dd]); % allow scalar L to be used for all dimensions
end
st.oversample = Ld;
st.Ld = Ld;

if dd ~= length(Jd) || dd ~= length(Kd) || dd ~= length(Ld)
	whos
	fail 'inconsistent dim'
end

% desired scaling factors, using fake omega=[0]
ktype_args = {'ktype', st.ktype}; % user's interpolator
if ~isempty(st.kb_m)
	ktype_args = {ktype_args{:}, 'kb_m', st.kb_m, 'kb_alf', st.kb_alf};
elseif ~isempty(st.alpha)
	ktype_args = {ktype_args{:}, 'alpha', st.alpha, 'beta', st.beta};
end

tmp = newfft(zeros(size(Nd)), st.Nd, 'Jd', st.Jd, 'Kd', st.Kd, ...
	'n_shift', 0*st.n_shift, ktype_args{:}, 'mode', 'sparse', ...
	'phasing', st.phasing);
st.sn = tmp.sn;

if streq(st.phasing, 'flipreal')
	iseven = @(n) mod(n,2) == 0;
	st.flips = iseven(st.Nd); % [1 dd]
end

% 1D tables for each dimension
for id = 1:dd

	ktype_args = {'ktype', st.ktype}; ... % user's interpolator
	if ~isempty(st.kb_m)
		ktype_args = {ktype_args{:}, ...
			'kb_m', st.kb_m(id), 'kb_alf', st.kb_alf(id)};
	elseif ~isempty(st.alpha)
		ktype_args = {ktype_args{:}, ...
			'alpha', {st.alpha{id}}, 'beta', {st.beta{id}}};
	end

	make_args = {Nd(id), Jd(id), Kd(id), Ld(id), ktype_args, st.phasing};
	st.h{id} = newfft_table_make(arg.how, make_args{:});

	if 0 % test fast vs slow
		[h0 t0] = newfft_table_make('slow', make_args{:});
		h = st.h{id};
		jf plc 2 2
		jf sub 1, plot(t0, real(h0), 'c-', t0, real(h), 'y.')
		jf sub 2, plot(t0, real(h0-h), 'y.')
		jf sub 3, plot(t0, imag(h0), 'c-', t0, imag(h), 'y.')
		jf sub 4, plot(t0, imag(h0-h), 'y.')
		max_percent_diff(h0, h)
	prompt
	end

	switch st.phasing
	case 'complex'
		if isreal(st.h{id})
			warn 'real kernel?'
%			st.h{id} = complex(st.h{id});
		end

	case {'real', 'flipreal', 'none'}
		if ~isreal(st.h{id})
			fail 'nonreal kernel bug'
		end

	otherwise
		fail bug
	end

end

st = strum(st, { ...
	'fft', @newfft_approx_for, '(x, [om])';
	'adj', @newfft_approx_adj, '(X, [om])';
	'interp_for', @newfft_table_for, '(Xk, [om])';
	'interp_adj', @newfft_table_adj, '(Xk, [om])';
	'plot', @newfft_table_plot, '()';
	});

% newfft_table_make()
function [h, t0] = newfft_table_make(how, N, J, K, L, ktype_args, phasing)

if streq(phasing, 'flipreal')
	phasing = 'none'; % trick: get pure real kernel (with no sign flips)
end
% note: if phasing is already 'real' then keep it that way!

mode_args = {'mode', 'sparse', 'Jd', J, 'n_shift', 0, 'phasing', phasing};

t0 = [-J*L/2:J*L/2]' / L; % [J*L+1]

switch how

% This is a slow and inefficient (but simple) way to get the table
% because it builds a huge sparse matrix but only uses 1 column!
case 'slow'

	om0 = t0 * 2*pi/K; % gam
	s1 = newfft(om0, N, 'Kd', K, ktype_args{:}, mode_args{:});
	h = full(s1.p(:,1));

% This way is "J times faster" than the slow way, but still not ideal.
% It works for any user-specified interpolator.
case 'fast'

	t1 = J/2-1 + [0:(L-1)]' / L;	% [L]
	om1 = t1 * 2*pi/K;		% * gam
	s1 = newfft(om1, N, 'Kd', K, ktype_args{:}, mode_args{:});
	h = col(full(s1.p(:,J:-1:1))); % [J*L+0]
%	h = [h; 0]; % [J*L+1]
	h = [h; h(1)]; % [J*L+1] - assuming symmetric

% This efficient way uses only "J" columns of sparse matrix!
% The trick to this is to use fake small values for N and K,
% which works for interpolators that depend only on the ratio K/N.
case 'ratio' % e.g., 'minmax:kb' | 'kb:*'

	Nfake = J;
	Kfake = Nfake * K/N;
	t1 = J/2-1 + [0:(L-1)]' / L;	% [L]
	om1 = t1 * 2*pi/Kfake;		% "gam"
	s1 = newfft(om1, Nfake, 'Kd', Kfake, ktype_args{:}, mode_args{:});
	h = col(full(s1.p(:,J:-1:1))); % [J*L+0]
%	h = [h; 0]; % [J*L+1]
	h = [h; h(1)]; % [J*L+1] assuming symmetric

	if streq(phasing, 'complex')
		h = h .* exp(1i * pi * t0 * (1/K - 1/Kfake)); % fix phase
	end

otherwise
	fail('bad type %s', how)
end

% newfft_table_plot()
function newfft_table_plot(st)
color = 'cym';
arg = {};
for id = 1:st.dd
	J = st.Jd(id);
	L = st.Ld(id);
	t0 = [-J*L/2:J*L/2]' / L; % [J*L+1]
	arg = {arg{:}, t0, real(st.h{id}), [color(id) '-']};
	arg = {arg{:}, t0, imag(st.h{id}), [color(id) '--']};
end
plot(arg{:})

% newfft_table_for()
function X = newfft_table_for(st, Xk, om)
if ~isvar('om') || isempty(om)
	om = st.om;
end
if isempty(om)
	fail 'need om or st.om'
end
if st.dd ~= size(om,2), fail('omega needs %d columns', dd), end

X = nufft_table_interp(st, Xk, st.order, st.flips, om);

% newfft_exact_adj()
function Xk = newfft_table_adj(st, X, om)
if ~isvar('om') || isempty(om)
	om = st.om;
end
if isempty(om)
	fail 'need om or st.om'
end
if st.dd ~= size(om,2), fail('omega needs %d columns', dd), end

Xk = nufft_table_adj(st, X, st.order, st.flips, om);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function test_tutorial_natmeg2014_timefrequency

% WALLTIME 00:30:00
% MEM 4gb
% DEPENDENCY

% this script executes the MATLAB content from
% http://www.fieldtriptoolbox.org/tutorial/natmeg2014/timefrequency
%
% it corresponds to the wiki version of 7 October 2014

clear all
close all

cd(dccnpath('/home/common/matlab/fieldtrip/data/ftp/workshop/natmeg2014'));

cfg = [];
cfg.dataset = 'oddball1_mc_downsampled.fif';
cfg.channel = 'MEG';

% define trials based on responses
cfg.trialdef.prestim       = 1.5;
cfg.trialdef.poststim      = 2.0;
cfg.trialdef.stim_triggers = [1 2];
cfg.trialdef.rsp_triggers  = [256 4096];
cfg.trialfun               = 'trialfun_oddball_responselocked';
cfg                        = ft_definetrial(cfg);

% preprocess MEG data
cfg.continuous             = 'yes';
cfg.demean                 = 'yes';
cfg.dftfilter              = 'yes';
cfg.dftfreq                = [50 100];

data_MEG_responselocked    = ft_preprocessing(cfg);

cfg              = [];
cfg.output       = 'pow';
cfg.channel      = 'all';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.toi          = [-1 : 0.10 : 1.5];
cfg.foi          = 1:40;
cfg.t_ftimwin    = ones(size(cfg.foi)) * 0.5;

cfg.trials       = find(data_MEG_responselocked.trialinfo(:,1) == 256);
TFR_left_MEG     = ft_freqanalysis(cfg, data_MEG_responselocked);

cfg.trials       = find(data_MEG_responselocked.trialinfo(:,1) == 4096);
TFR_right_MEG    = ft_freqanalysis(cfg, data_MEG_responselocked);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.zlim         = [-2e-26 2e-26];
cfg.showlabels   = 'yes';
cfg.layout       = 'neuromag306mag.lay';
cfg.channel      = 'MEG*1';

figure;
ft_multiplotTFR(cfg, TFR_left_MEG);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.maskstyle    = 'saturation';
cfg.zlim         = [-1e-26 1e-26];
cfg.channel      = 'MEG1041';

figure;
ft_singleplotTFR(cfg, TFR_left_MEG);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.xlim         = [0.4 0.8];
cfg.zlim         = [-4e-27 4e-27];
cfg.ylim         = [15 25];
cfg.marker       = 'on';
cfg.layout       = 'neuromag306mag.lay';
cfg.channel      = 'MEG*1';

figure;
ft_topoplotTFR(cfg, TFR_left_MEG);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.xlim         = [0.4 0.8];
cfg.zlim         = [-4e-27 4e-27];
cfg.ylim         = [15 25];
cfg.marker       = 'on';
cfg.layout       = 'neuromag306mag.lay';
cfg.channel      = 'MEG*1';

figure;
ft_topoplotTFR(cfg, TFR_right_MEG);

cfg = [];
cfg.parameter = 'powspctrm';
cfg.operation = '(x1-x2)/(x1+x2)';

TFR_diff_MEG = ft_math(cfg, TFR_right_MEG, TFR_left_MEG);

cfg = [];
cfg.xlim         = [0.4 0.8];
cfg.zlim         = [-0.4 0.4];
cfg.ylim         = [15 25];
cfg.marker       = 'on';
cfg.layout       = 'neuromag306mag.lay';
cfg.channel      = 'MEG*1';

figure;
ft_topoplotTFR(cfg, TFR_diff_MEG);

cfg = [];
cfg.dataset = 'oddball1_mc_downsampled.fif';

% define trials based on responses
cfg.trialdef.prestim       = 1.5;
cfg.trialdef.poststim      = 2.0;
cfg.trialdef.stim_triggers = [1 2];
cfg.trialdef.rsp_triggers  = [256 4096];
cfg.trialfun               = 'trialfun_oddball_responselocked';
cfg                        = ft_definetrial(cfg);

% preprocess EEG data
cfg.channel                = 'EEG';
cfg.continuous             = 'yes';
cfg.demean                 = 'yes';
cfg.dftfilter              = 'yes';
cfg.dftfreq                = [50 100];

data_EEG_responselocked    = ft_preprocessing(cfg);

if false
  % skip the interactive section
  % select bad channels
  cfg = [];
  cfg.metric  = 'var';
  temp        = ft_rejectvisual(cfg, data_EEG_responselocked);
  % with this little trick we get the names of the selected channels
  badchannels = setdiff(data_EEG_responselocked.label,temp.label);
else
  badchannels = {
    'EEG001'
    'EEG002'
    'EEG003'
    'EEG004'
    'EEG005'
    'EEG006'
    'EEG007'
    'EEG008'
    'EEG015'
    'EEG016'
    'EEG017'
    'EEG097'
    'EEG098'
    'EEG099'
    'EEG100'
    'EEG101'
    'EEG102'
    'EEG111'
    'EEG112'};
end

% determine neighbours structure
cfg            = [];
cfg.method     = 'triangulation';
cfg.senstype   = 'EEG'; % Our data still contains information from the MEG channels, we want to make sure ft_prepare_neighbours does not get confused
neighbours_EEG = ft_prepare_neighbours(cfg, data_EEG_responselocked);

% plotting neighbours
cfg            = [];
cfg.neighbours = neighbours_EEG;
cfg.senstype   = 'EEG';
ft_neighbourplot(cfg, data_EEG_responselocked);

% fix channels
cfg = [];
cfg.method                    = 'spline';
cfg.neighbours                = neighbours_EEG;
cfg.badchannel                = badchannels;
cfg.senstype                  = 'EEG';
data_clean_EEG_responselocked = ft_channelrepair(cfg, data_EEG_responselocked);

cfg = [];
cfg.reref                  = 'yes';
cfg.refchannel             = 'all';
data_clean_EEG_responselocked = ft_preprocessing(cfg, data_clean_EEG_responselocked);

cfg              = [];
cfg.output       = 'pow';
cfg.channel      = 'all';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.toi          = [-1 : 0.10 : 1.5];
cfg.foi          = 1:40;
cfg.t_ftimwin    = ones(size(cfg.foi)) * 0.5;

cfg.trials       = find(data_clean_EEG_responselocked .trialinfo(:,1) == 256);
TFR_left_EEG     = ft_freqanalysis(cfg, data_clean_EEG_responselocked );

cfg.trials       = find(data_clean_EEG_responselocked .trialinfo(:,1) == 4096);
TFR_right_EEG    = ft_freqanalysis(cfg, data_clean_EEG_responselocked );

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.xlim         = [0.5 1.0];
cfg.zlim         = [-4e-12 4e-12];
cfg.ylim         = [15 25];
cfg.marker       = 'on';
cfg.layout       = 'natmeg_customized_eeg1005.lay';

figure;
ft_topoplotTFR(cfg, TFR_left_EEG);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'relchange';
cfg.ylim         = [15 25];
cfg.xlim         = [0.5 1.0];
cfg.zlim         = [-1.2 1.2];
cfg.layout       = 'natmeg_customized_eeg1005.lay';

figure;
ft_topoplotTFR(cfg, TFR_left_EEG);

cfg = [];
cfg.parameter    = 'powspctrm';
cfg.operation    = '(x1-x2)/(x1+x2)';
TFR_diff_EEG = ft_math(cfg, TFR_right_EEG, TFR_left_EEG);

if false
  % skip the interactive section
  % if ft_math didn't work, then just do it by hand - its exactly the same:
  TFR_diff_EEG = TFR_right_EEG;
  TFR_diff_EEG.powspctrm = (TFR_right_EEG.powspctrm - TFR_left_EEG.powspctrm) ./ (TFR_right_EEG.powspctrm + TFR_left_EEG.powspctrm);
end

cfg = [];
cfg.xlim         = [0.4 0.8];
cfg.ylim         = [15 25];
cfg.zlim         = [-0.2 0.2];
cfg.marker       = 'on';
cfg.layout       = 'natmeg_customized_eeg1005.lay';

figure;
ft_topoplotTFR(cfg, TFR_diff_EEG);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.xlim         = [0.4 0.8];
cfg.ylim         = [15 25];
cfg.zlim         = [-1e-24 1e-24];
cfg.marker       = 'on';
cfg.layout       = 'neuromag306planar.lay';

figure;
ft_topoplotTFR(cfg, TFR_left_MEG);

TFR_left_MEG_comb  = ft_combineplanar([],TFR_left_MEG);
TFR_right_MEG_comb = ft_combineplanar([],TFR_right_MEG);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.xlim         = [0.4 0.8];
cfg.ylim         = [15 25];
cfg.zlim         = [-4e-24 4e-24];
cfg.marker       = 'on';
cfg.layout       = 'neuromag306cmb.lay';

figure;
ft_topoplotTFR(cfg, TFR_left_MEG_comb);

cfg = [];
cfg.parameter = 'powspctrm';
cfg.operation = '(x1-x2)/(x1+x2)';

TFR_diff_MEG_comb = ft_math(cfg, TFR_right_MEG_comb, TFR_left_MEG_comb);

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.xlim         = [0.4 0.8];
cfg.ylim         = [15 25];
cfg.zlim         = [-0.3 0.3];
cfg.marker       = 'on';
cfg.layout       = 'neuromag306cmb.lay';

figure;
ft_topoplotTFR(cfg, TFR_diff_MEG_comb);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
report_this_filefun(mfilename('fullpath'));

str = [];
[newmatfile, newpath] = uiputfile([hodi ], 'Save As'); %Syntax change Matlab Version 7, no window positioning on macs

s = [xt2  ; meand ; er   ];
fid = fopen([newpath newmatfile],'w') ;
fprintf(fid,'%6.2f  %6.2f %6.2f\n',s);
fclose(fid)

return

author: CelsoReyes


text:
function jsubs = cfg2jobsubs(item, subs)

% function jsubs = cfg2jobsubs(item, subs)
% Return the subscript into the job tree for a given subscript vector into
% the val part of the cfg tree. In a cfg_branch, this is a struct reference
% to a field with the name of the tag of the corresponding child node.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg2jobsubs.m 1862 2008-06-30 14:12:49Z volkmar $

rev = '$Rev: 1862 $'; %#ok

% Only de-reference subscripts into .val{x} field
if isequal(subs(1), substruct('.','val')) && strcmp(subs(2).type, '{}')
    % return possible indices into .val{1} field
    citem = subsref(item, subs(1:2));
    if numel(subs) > 2
        jsubs1 = cfg2jobsubs(citem, subs(3:end));
    else
        jsubs1 = [];
    end;
    jsubs = [substruct('.', gettag(citem)) jsubs1];
else
    cfg_message('matlabbatch:cfg2jobsubs:wrongsubs', 'Inappropriate subscript reference in item ''%s''.', item.tag);
    jsubs = struct('type',{},'subs',{});
end;

author: spm


text:
function test_failed=test_demos
%TEST_DEMOS  Test if all the demos runs without errors.

  test_failed=0;

  s=dir([ltfatbasepath,filesep,'demos',filesep,'demo_*.m']);

  for ii=1:numel(s)
     filename = s(ii).name;

     disp(filename);

     % The demo is run in separate function to avoid
     % variable name clash
     rundemo(filename(1:end-2));

  end

function rundemo(demoname)
   close all;
   eval(demoname);
   clear all;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function r82col_print_part ( n, a, max_print, title )

%*****************************************************************************80
%
%% R82COL_PRINT_PART prints "part" of an R82COL.
%
%  Discussion:
%
%    An R82COL is an (N,2) array of R8's.
%
%    The user specifies MAX_PRINT, the maximum number of lines to print.
%
%    If N, the size of the vector, is no more than MAX_PRINT, then
%    the entire vector is printed, one entry per line.
%
%    Otherwise, if possible, the first MAX_PRINT-2 entries are printed,
%    followed by a line of periods suggesting an omission,
%    and the last entry.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries of the vector.
%
%    Input, real A(N,2), the vector to be printed.
%
%    Input, integer MAX_PRINT, the maximum number of lines
%    to print.
%
%    Input, string TITLE, a title.
%
  if ( max_print <= 0 )
    return
  end

  if ( n <= 0 )
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  if ( n <= max_print )

    for i = 1 : n
      fprintf ( 1, '  %4d  %14g  %14g\n', i, a(i,1), a(i,2) );
    end

  elseif ( 3 <= max_print )

    for i = 1 : max_print-2
      fprintf ( 1, '  %4d  %14g  %14g\n', i, a(i,1), a(i,2) );
    end
    fprintf ( 1, '  ....  ..............  ..............\n' );
    i = n;
    fprintf ( 1, '  %4d  %14g  %14g\n', i, a(i,1), a(i,2) );

  else

    for i = 1 : max_print - 1
      fprintf ( 1, '  %4d  %14g  %14g\n', i, a(i,1), a(i,2) );
    end
    i = max_print;
    fprintf ( 1, '  %4d  %14g  %14g  ...more entries...\n', i, a(i,1), a(i,2) );

  end

  return
end