text:
%implementation of the TVE SLM model by Martha Banbura and Andries van Vlodrop (2017)
%by Ben Schumann (2020)

function [Ys, Yt, Y, data_endo, const, priorValues, dataValues, sizetraining]=...
    TVESLM_prior(data_endo, data_exo, names, endo, lags, lambda1, lambda2, lambda3, lambda5, ar, bex, dataSLM, namesSLM, datesSLM, const, priorexo, gamma)
%input n, data_endo, data_exo lags, q, lambda1, lambda2, lambda3, lambda5,ar bex,dataSLM, namesSLM, datesSLM, endo, names,n
%output ys, data_endo, Yt, const, priorValues,dataValues. data.values

%% Preliminaries
p=lags;
T = size(data_endo,1);
n = size(data_endo,2);
% compute m, the number of exogenous variables in the model
% if data_exo is empty, set m=0
if isempty(data_exo)==1
m=0;
% if data_exo is not empty, count the number of exogenous variables that will be included in the model
else
m=size(data_exo,2);
% Also, trim a number initial rows equal to the number of lags, as they will be suppressed from the endogenous as well to create initial conditions
data_exo=data_exo(p+1:end,:);
end

if const ==1
    m=m-1;
    const=0; %no constant as the VAR is estimated on the local mean adjusted variables
end

% determine k, the number of parameters to estimate in each equation; it is equal to np+m
k=n*p+m;
% determine q, the total number of VAR parameters to estimate; it is equal to n*k
q=n*k;
%Divide the sample into a training sample and an estimation sample
sizetraining = floor(size(data_endo,1)/5) + lags;
Yt = data_endo(1:sizetraining,:);           %training sample
dataValues.Yt = Yt;                         %training sample

Y = data_endo((sizetraining-lags)+1:end,:); %estimation sample
dataValues.Y = Y;                           %estimation sample

dataValues.data_endo_full = data_endo;
data_endo = data_endo((sizetraining-lags)+1:end,:);

%Determine for which variables there is a survey local mean and build selection matrix

HaveSLM = ismember(endo,namesSLM)'; %first establish for which variables we have survey local mean data
ns = sum(HaveSLM);                  %check how many endogenous variables have a survey local mean
Ys = nan(size(Y,1),ns);
ni = 0;
Ppsi = zeros(ns,n);                 %create selection matrix for measurement equation
% for kk = 1:n
%     if HaveSLM(1,kk) ==1
%         ni = ni+1;
%         Ys(:,kk) = dataSLM(sizetraining-lags+1:end,kk);
%         Ppsi(ni,kk) = 1;            %selection matrix for local mean measurement equation
%     end
% end
for kk = 1:n
    if ismember(endo{kk,1},namesSLM)
        %find entry in ismember
        IndexC= find(strcmp(namesSLM, endo{kk,1}));
        ni = ni+1;
        Ys(:,IndexC) = dataSLM(sizetraining-lags+1:end,IndexC);
        Ppsi(IndexC,kk) = 1;            %selection matrix for local mean measurement equation
    end
end

dataValues.Ppsi = Ppsi;

dataValues.Ys = Ys;                 %Survey local mean after training sample

%%%%%%%%%%%%%%%%%% set the prior values%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% prior for stochastic volatility in VAR residuals priors (H = A^-1'*Lambda_t*A^-1')
%estimate an ar(p) process for setting the prior variances
[arvar]=bear.arloop(Yt,1,lags,n);
arvar = arvar*((size(Yt,1)-lags-lags-1)/(size(Yt,1)-lags-1)); %translate bear into Banbura et al (2017) estimate of ar(4) VCV

priorValues.priorVarAscaling_H =10;  %prior variance for the below diagonal elements of A (constant part of VCV
priorValues.priorMeanAscaling_H=0;   %prior mean for the below diagonal elements of A

priorValues.mean_ln_h0=log(arvar);    %prior mean of the initial condition for stochastic volatilty VAR residuals
priorValues.var_ln_h0 = 10*ones(n,1); %prior variance for the inital condition  for stochastic volatilty VAR residuals

priorValues.phi_h=0.01;               %centering parameter for the variance of the stochastic volatility process in the VAR residuals
priorValues.d_h=10;                   %scaling parameter for the variance of the stochastic volatility process in the VAR residuals

priorValues.offset_c = 0.001;         %constant for log transformation in order to avoid numerical problems of log close to 0

%% Priors for Reduced form VAR coefficients (independent normal wishart)
priorValues.vars  =arvar;
%for bear sampling
priorValues.lambda1=lambda1;        % overall tightness
priorValues.lambda2=lambda2;        % cross variable shrinkage
priorValues.lambda3=lambda3;        % lag shrinkage
priorValues.lambda5=lambda5;        % block exogeneity shrinkage
priorValues.ar = ar;                % prior values for ar(1) coefficients of the variable

%for banbura et al sampling
priorValues.lambda=lambda1^2;       % only used for banbura large model sampling
priorValues.theta=lambda3;          % only used for banbura large model sampling

if size(ar,1) < n              % if no variable specific priors were given
    artemp=zeros(n,1);
    for kk=1:n
        artemp(kk,1)=ar;
    end
    ar = artemp;
end

% start with beta0, defined in (1.3.4)
% it is a q*1 vector of zeros, save for the n coefficients of each variable on their own first lag
beta0=zeros(q,1);

for ii=1:n
beta0((ii-1)*k+ii,1)=ar(ii,1);
end

% if a prior for the exogenous variables is selected put it in here:
for ii=1:n
    beta0(k*ii)=priorexo(ii,1);
end

% next compute omega0, the variance-covariance matrix of beta, defined in (1.3.8)
% set it first as a q*q matrix of zeros
omega0=zeros(q,q);

% set the variance on coefficients related to own lags, using (1.3.5)
for ii=1:n
   for jj=1:p
   omega0((ii-1)*k+(jj-1)*n+ii,(ii-1)*k+(jj-1)*n+ii)=(lambda1/jj^lambda3)^2;
   end
end

%  set variance for coefficients on cross lags, using (1.3.6)
for ii=1:n
   for jj=1:p
      for kk=1:n
      if kk==ii
      else
      omega0((ii-1)*k+(jj-1)*n+kk,(ii-1)*k+(jj-1)*n+kk)=(arvar(ii,1)/arvar(kk,1))*(((lambda1*lambda2)/(jj^lambda3))^2);
      end
      end
   end
end

% finally set the variance for exogenous variables, using (1.3.7)
for ii=1:n
   for jj=1:m
   omega0(ii*k-m+jj,ii*k-m+jj)=arvar(ii,1)*((lambda1*lambda4)^2);
   end
end

% if block exogeneity has been selected, implement it, according to (1.7.4)
if bex==1
   for ii=1:n
      for jj=1:n
         if blockexo(ii,jj)==1
            for kk=1:p
            omega0((jj-1)*k+(kk-1)*n+ii,(jj-1)*k+(kk-1)*n+ii)=omega0((jj-1)*k+(kk-1)*n+ii,(jj-1)*k+(kk-1)*n+ii)*lambda5^2;
            end
         else
         end
      end
   end
% if block exogeneity has not been selected, don't do anything
else
end

priorValues.beta0=beta0;        %prior for the reduced form VAR coefficients in vectorized form
priorValues.omega0=omega0;       %prior variance for reduced form VAR coefficients

%% Priors for the measurement equation (G)

priorValues.mean_ln_g0=zeros(ns,1); %mean of the initial condition in the variance of measurement equation residuals
priorValues.var_ln_g0 = 10*ones(ns,1); %variance of the initial condition in the variance of measurement equation residuals

priorValues.phi_g=0.01;    %centering parameter for the variance of measurement equation residuals
priorValues.d_g=10;        %scaling parameter for the variance of the stochastic volatility process in the VAR residuals

%% Priors for the local mean equation
priorValues.meanTS=mean(Yt)';          %start values for mean corrected Kalman smoother simulation in order to draw from the local mean process
priorValues.kappa = 1000;              %variance of the initial state for the local mean

priorValues.gamma = gamma;             %value for ar coefficient in

priorValues.mean_ln_v0=log(arvar);     %prior mean of the initial condition for stochastic volatilty in the state transition equation
priorValues.var_ln_v0 = 10*ones(n,1);  %prior variance of the initial condition for stochastic volatilty in the state transition equation

priorValues.d_v=10;                    %scaling parameter for the variance of the stochastic volatility process in the VAR residuals
priorValues.phi_v=0.01;                %centering parameter for the variance of transition equation residuals

author: european-central-bank

repo: BEAR-toolbox


text:
function y = subsref( x, varargin )
y = cvxtuple( subsref( x.value_, varargin{:} ) );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


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

    properties

    end

    methods
        normvect = getPlaneNormVectAtTime(obj, time, vehElemSet, inFrame)

        originPt = getPlaneOriginPtAtTime(obj, time, vehElemSet, inFrame)

        name = getName(obj)

        setName(obj, name)

        listboxStr = getListboxStr(obj)

        useTf = openEditDialog(obj)

        tf = isVehDependent(obj)

        tf = usesGroundObj(obj, groundObj)

        tf = usesGeometricPoint(obj, point)

        tf = usesGeometricVector(obj, vector)

        tf = usesGeometricCoordSys(obj, coordSys)

        tf = usesGeometricRefFrame(obj, refFrame)

        tf = usesGeometricAngle(obj, angle)

        tf = usesGeometricPlane(~, ~)

        tf = isInUse(obj, lvdData)
    end

    methods(Sealed)
        function tf = eq(a,b)
            tf = eq@handle(a,b);
        end

        function tf = ne(a,b)
            tf = ne@handle(a,b);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function set_theory_test ( )

%*****************************************************************************80
%
%% SET_THEORY_TEST tests the SET_THEORY library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SET_THEORY_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SET_THEORY library.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Demonstrate some set theory operations that\n' );
  fprintf ( 1, '  can be implemented in MATLAB.\n' );

  set_theory_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SET_THEORY_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [cm] = in2cm(in)
% Convert length from inches to centimeters.
% Chad A. Greene
cm = in*2.54;

author: Sable

repo: mcbench-benchmarks


text:
%% GrabCut segmentation demo
% Interactive foreground extraction using the GrabCut algorithm.
%
% This program demonstrates GrabCut segmentation: select an object in a
% region and then grabcut will attempt to segment it out.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/grabcut.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/grabcut.py>
% * <https://docs.opencv.org/3.2.0/d8/d83/tutorial_py_grabcut.html>
%

%% Theory
%
% GrabCut algorithm was designed by Carsten Rother, Vladimir Kolmogorov and
% Andrew Blake from Microsoft Research Cambridge, UK. in their paper:
%
% * "GrabCut": interactive foreground extraction using iterated graph cuts
%   <http://dl.acm.org/citation.cfm?id=1015720>
%
% An algorithm was needed for foreground extraction with minimal user
% interaction, and the result was *GrabCut*.
%
% How it works from user point of view? Initially user draws a rectangle
% around the foreground region (foreground region should be completely inside
% the rectangle). Then algorithm segments it iteratively to get the best
% result. Done. But in some cases, the segmentation won't be fine, like, it
% may have marked some foreground region as background and vice versa. In that
% case, user need to do fine touch-ups. Just give some strokes on the images
% where some faulty results are there. Strokes basically says: "Hey, this
% region should be foreground, you marked it background, correct it in next
% iteration", or its opposite for background. Then in the next iteration, you
% get better results.
%
% See the image below. First player and football is enclosed in a blue
% rectangle. Then some final touchups with white strokes (denoting foreground)
% and black strokes (denoting background) is made. And we get a nice result.
%
% <<https://docs.opencv.org/3.2.0/grabcut_output1.jpg>>
%
% So what happens in background ?
%
% * User inputs the rectangle. Everything outside this rectangle will be taken
%   as sure background (That is the reason it is mentioned before that your
%   rectangle should include all the objects). Everything inside rectangle is
%   unknown. Similarly any user input specifying foreground and background are
%   considered as hard-labelling which means they won't change in the process.
% * Computer does an initial labelling depeding on the data we gave. It labels
%   the foreground and background pixels (or it hard-labels)
% * Now a Gaussian Mixture Model (GMM) is used to model the foreground and
%   background.
% * Depending on the data we gave, GMM learns and create new pixel
%   distribution. That is, the unknown pixels are labelled either probable
%   foreground or probable background depending on its relation with the other
%   hard-labelled pixels in terms of color statistics (It is just like
%   clustering).
% * A graph is built from this pixel distribution. Nodes in the graphs are
%   pixels. Additional two nodes are added, *Source node* and *Sink node*.
%   Every foreground pixel is connected to source node and every background
%   pixel is connected to sink node.
% * The weights of edges connecting pixels to source node/end node are defined
%   by the probability of a pixel being foreground/background. The weights
%   between the pixels are defined by the edge information or pixel
%   similarity. If there is a large difference in pixel color, the edge
%   between them will get a low weight.
% * Then a mincut algorithm is used to segment the graph. It cuts the graph
%   into two separating source node and sink node with minimum cost function.
%   The cost function is the sum of all weights of the edges that are cut.
%   After the cut, all the pixels connected to source node become foreground
%   and those connected to sink node become background.
% * The process is continued until the classification converges.
%
% It is illustrated in below image
% (Image Courtesy: <http://www.cs.ru.ac.za/research/g02m1682/>)
%
% <<https://docs.opencv.org/3.2.0/grabcut_scheme.jpg>>
%

%% Code
% This is an interactive tool using grabcut. You can also watch this
% <http://www.youtube.com/watch?v=kAwxLTDDAwU youtube video> on how to use it.
%

function varargout = grabcut_demo_gui(im)
    % load an image
    if nargin < 1
        src = imread(fullfile(mexopencv.root(),'test','fruits.jpg'));
    elseif isempty(im)
        fmts = imformats();
        filtspec = strjoin(strcat('*.', [fmts.ext]), ';');
        [fn,fp] = uigetfile(filtspec, 'Select an image');
        if fp==0, error('No file selected'); end
        src = imread(fullfile(fp,fn));
    elseif ischar(im)
        src = imread(im);
    else
        src = im;
    end

    % we expect an 8-bit RGB image
    validateattributes(src, {'uint8'}, ...
        {'ndims',3, 'size',[nan nan 3], 'nonempty'});

    % initialize app state, and create the UI
    app = initApp(src);
    h = buildGUI(src, app);

    % hook event handlers
    opts = {'Interruptible','off', 'BusyAction','cancel'};
    set(h.pop, 'Callback',@onChange);
    set(h.btn(1), 'Callback',@onHelp);
    set(h.btn(2), 'Callback',@onReset);
    set(h.btn(3), 'Callback',@onNext, opts{:});
    set(h.fig, 'WindowKeyPressFcn',@onType, ...
        'WindowButtonDownFcn',@onMouseDown, opts{:});

    % return graphics handles
    if nargout > 0, varargout{1} = h; end

    % ========== Event Handlers ==========

    function onHelp(~,~)
        %ONHELP  Display usage help dialog

        helpdlg({
            'This program demonstrates GrabCut segmentation:'
            'select an object in a region and then grabcut'
            'will attempt to segment it out.'
            ''
            'Select a rectangular area around the object you'
            'want to segment.'
            ''
            'Then press "next" to segment the object (once or a few times).'
            ''
            'For finer touch-ups, set the mode and draw lines on the areas you'
            'want, to mark them as foreground/background (sure or probable),'
            'then press "next" again.'
            ''
            'Hot keys:'
            'ESC - quit the program'
            'r - restore the original image'
            'n - next iteration'
            'left mouse button - First set rectangle, then set pixels'
            '    as either BGD/FGD/PR_BGD/PR_FGD depending on selected'
            '    mode in dropdown menu.'
        });
    end

    function onReset(~,~)
        %ONRESET  Event handler for reset button

        app.mask(:) = 0;
        app.bgdModel(:) = 0;
        app.fgdModel(:) = 0;
        app.rect = zeros(0,4);
        app.rectxy = zeros(0,2);
        app.pts = repmat({zeros(0,2)}, [1 4]);
        app.iterCount = 0;
        app.isInitialized = false;

        set(h.txt, 'String','Iter =  0');
        set(h.img, 'CData',app.img0);
        set(h.rect, 'XData',NaN, 'YData',NaN);
        set(h.line(:), 'XData',NaN, 'YData',NaN);
        drawnow;
    end

    function onNext(~,~)
        %ONNEXT  Event handler for next button

        if app.isInitialized
            % set pixels in GC mask using drawing points
            if any(~cellfun(@isempty, app.pts))
                setLblsInMask();
            end
            % continue using current mask
            tic
            [app.mask, app.bgdModel, app.fgdModel] = cv.grabCut(...
                app.img0, app.mask, 'Mode','Eval', 'IterCount',1, ...
                'BgdModel',app.bgdModel, 'FgdModel',app.fgdModel);
            toc
        elseif any(~cellfun(@isempty, app.pts))
            % set foreground pixels in GC mask using rectangle
            setRectInMask();
            % set pixels in GC mask using drawing points
            setLblsInMask();
            % init using mask
            tic
            [app.mask, app.bgdModel, app.fgdModel] = cv.grabCut(...
                app.img0, app.mask, 'Mode','InitWithMask', 'IterCount',1);
            toc
        elseif ~isempty(app.rect)
            % init using rectangle
            rect = app.rect - [1 1 0 0];
            tic
            [app.mask, app.bgdModel, app.fgdModel] = cv.grabCut(...
                app.img0, rect, 'Mode','InitWithRect', 'IterCount',1);
            toc
        else
            disp('First select object to segment by drawing a rectangle');
            return;
        end

        % mark mask as initialized, and increment counter
        app.isInitialized = true;
        app.iterCount = app.iterCount + 1;

        % show result
        showImage();
    end

    function onType(~,e)
        %ONTYPE  Event handler for key press on figure

        % handle keys
        switch e.Key
            case {'q', 'escape'}
                close(h.fig);

            case 'h'
                onHelp([],[]);

            case 'r'
                onReset([],[]);

            case 'n'
                onNext([],[]);

            case {'add', 'subtract'}
                % adjust brush thickness
                if strcmp(e.Character, '+')
                    app.thick = min(app.thick + 2, 40);
                elseif strcmp(e.Character, '-')
                    app.thick = max(app.thick - 2, 1);
                end
                set(h.line(:), 'MarkerSize',app.thick*5.4);

            case {'1', '2', '3', '4'}
                % set brush value
                app.currIdx = str2double(e.Key);
                set(h.pop, 'Value',app.currIdx);
        end
    end

    function onChange(~,~)
        %ONCHANGE  Event handler for UI controls

        % change current GC mask drawing value: BGD/FGD/PR_BGD/PR_FGD
        app.currIdx = get(h.pop, 'Value');
    end

    function onMouseDown(~,~)
        %ONMOUSEDOWN  Event handler for mouse down on figure

        % ignore anything but left mouse clicks
        if ~strcmp(get(h.fig,'SelectionType'), 'normal')
            return;
        end

        % one of two phases: drawing rectangle, or free-drawing of points
        if isempty(app.rect)
            % select and draw rectangle
            select_rectangle();
            if isempty(app.rect), return; end
            set(h.rect, 'XData',app.rectxy(:,1), 'YData',app.rectxy(:,2));
        else
            % attach event handlers, and change mouse pointer
            set(h.fig, 'Pointer','circle', ...
                'WindowButtonMotionFcn',@onMouseMove, ...
                'WindowButtonUpFcn',@onMouseUp);
        end
    end

    function onMouseMove(~,~)
        %ONMOUSEMOVE  Event handler for mouse move on figure

        % get current point and append it
        app.pts{app.currIdx}(end+1,:) = getCurrentPoint();

        % update corresponding graphic line
        set(h.line(app.currIdx), ...
            'XData',app.pts{app.currIdx}(:,1), ...
            'YData',app.pts{app.currIdx}(:,2));
    end

    function onMouseUp(~,~)
        %ONMOUSEUP  Event handler for mouse up on figure

        % detach event handlers, and restore mouse pointer
        set(h.fig, 'Pointer','arrow', ...
            'WindowButtonMotionFcn','', ...
            'WindowButtonUpFcn','');
    end

    % ========== Helper Functions ==========

    function showImage()
        out = app.img0;
        if app.isInitialized
            % zero-out background pixels
            if true
                binMask = repmat(app.mask == 0 | app.mask == 2, [1 1 3]);
                out(binMask) = 0;
            else
                binMask = (app.mask == 1 | app.mask == 3);
                out = cv.bitwise_and(out, out, 'Mask',binMask);
            end
        end
        set(h.img, 'CData',out);
        set(h.txt, 'String',sprintf('Iter = %2d',app.iterCount));
        drawnow;
    end

    function setRectInMask()
        % convert rectangle to binary mask
        rect_mask = poly2mask(app.rectxy(:,1), app.rectxy(:,2), ...
            app.sz(1), app.sz(2));

        % set foreground pixels in GC mask using rectangle
        app.mask(:) = 0;          % BGD
        app.mask(rect_mask) = 3;  % PR_FGD
    end

    function setLblsInMask()
        % set pixels in GC mask from drawing points: BGD, FGD, PR_BGD, PR_FGD
        for i=1:4
            app.mask = cv.circle(app.mask, app.pts{i}-1, app.thick, ...
                'Color',uint8(i-1), 'Thickness','Filled');
        end

        % clear drawing points after being processed
        app.pts = repmat({zeros(0,2)}, [1 4]);
        set(h.line(:), 'XData',NaN, 'YData',NaN);
    end

    function select_rectangle()
        %TODO: consider IMRECT from image_toolbox
        % create rubberband box to prompt user for a rectangle
        p1 = getCurrentPoint(); % retrieve mouse location before dragging
        rbbox;                  % ignore its output (figure coordinates)
        pause(0.005);           % CP might not get updated if selection was too fast
        p2 = getCurrentPoint(); % retrieve mouse location after dragging

        % form rectangle from two points: [x y w h]
        tl = min([p1;p2]);   % top-left corner
        br = max([p1;p2]);   % bottom-right corner
        if all((br-tl) > 1)  % ignore small rectangles
            app.rect = [tl br-tl];
            app.rectxy = [tl; tl+[app.rect(3) 0]; br; tl+[0 app.rect(4)]; tl];
        end
    end

    function p = getCurrentPoint()
        % retrieve current mouse location
        p = get(h.ax, 'CurrentPoint');
        p = p(1,1:2);

        % clamp to within image coordinates
        p = max(p, [1 1]);
        p = min(p, [app.sz(2) app.sz(1)]);
    end
end

% ========== Initializer functions ==========

function app = initApp(img)
    %INITAPP  Initialize app state

    app = struct();
    app.img0 = img;              % original image
    app.sz = size(img);          % image size
    app.mask = zeros(size(img,1), size(img,2), 'uint8'); % GC mask
    app.bgdModel = zeros(1,64);  % GC background model
    app.fgdModel = zeros(1,64);  % GC foreground model
    app.currIdx = 1;             % drawing value (BGD/FGD/PR_BGD/PR_FGD)
    app.thick = 5;               % drawing thickness
    app.pts = repmat({zeros(0,2)}, [1 4]); % drawing points of each brush
    app.rect = zeros(0, 4);      % rectangle [x,y,w,h]
    app.rectxy = zeros(0,2);     % rectangle points [TL;TR;BR;BL;TL]
    app.iterCount = 0;           % iterations counter
    app.isInitialized = false;   % whethet GC mask is initialized
end

function h = buildGUI(img, app)
    %BUILDGUI  Creates the UI

    % parameters
    sz = size(img);
    sz(2) = max(sz(2), 350);  % minimum figure width

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.fig = figure('Name','GrabCut Demo', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2) sz(1)+29]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax = axes('Parent',h.fig, ...
        'Units','pixels', 'Position',[1 30 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img = imshow(img, 'Parent',h.ax);
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax);
        h.img = imshow(img);
    end
    %axis(h.ax, 'on');
    h.btn(1) = uicontrol('Parent',h.fig, 'Style','pushbutton', ...
        'Position',[5 5 60 20], 'String','Help');
    h.btn(2) = uicontrol('Parent',h.fig, 'Style','pushbutton', ...
        'Position',[70 5 60 20], 'String','Reset');
    h.btn(3) = uicontrol('Parent',h.fig, 'Style','pushbutton', ...
        'Position',[135 5 60 20], 'String','Next');
    h.txt = uicontrol('Parent',h.fig, 'Style','text', 'FontSize',11, ...
        'Position',[200 5 60 20], 'String','Iter =  0');
    h.pop = uicontrol('Parent',h.fig, 'Style','popupmenu', ...
        'Position',[260 5 80 20], 'String',{'BGD','FGD','PR_BGD','PR_FGD'});

    % initialize lines (drawing and rectangle selection)
    clr = 'kwgrb';  % 'rbcmg'
    for i=1:4
        h.line(i) = line(NaN, NaN, 'Color',clr(i), 'Parent',h.ax, ...
            'LineStyle','none', 'Marker','.', 'MarkerSize',app.thick*5.4);
    end
    h.rect = line(NaN, NaN, 'Color',clr(5), 'Parent',h.ax, 'LineWidth',2);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%% Dummy test
lp = [1 1 1 1
      1 2 2 1
      1 3 3 1
      1 1 1 1];
ms = [2 3 4
      1 4 5];

% Get all pairs of neighboring leave regions
[~, idx_neighbors] = seg2gridbmap(lp);
K = max(idx_neighbors.matrix_max(:)) + 1;
neigh_pairs = unique(idx_neighbors.matrix_min+K*idx_neighbors.matrix_max);
neigh_pairs(neigh_pairs==0) = [];
neigh_pairs_min = mod(neigh_pairs,K);
neigh_pairs_max = (neigh_pairs-neigh_pairs_min)/K;

if isrow(neigh_pairs_min)
    neigh_pairs_min = neigh_pairs_min';
end
if isrow(neigh_pairs_max)
    neigh_pairs_max = neigh_pairs_max';
end

n_pairs    = 100;
n_triplets = 100;
curr_cands = mex_get_tree_cands(double(lp)-1, double(ms)-1,...
                                 neigh_pairs_min-1, neigh_pairs_max-1,...
                                 [n_pairs, n_triplets]);

assert(isequal(curr_cands{1},[1 2
                              1 3]))
assert(isempty(curr_cands{2}))

%% Simple test 1
lp = [1 2 3 4 5];
ms = [2 3 6
      4 5 7
      6 1 8
      7 8 9];

% Get all pairs of neighboring leave regions
[~, idx_neighbors] = seg2gridbmap(lp);
K = max(idx_neighbors.matrix_max(:)) + 1;
neigh_pairs = unique(idx_neighbors.matrix_min+K*idx_neighbors.matrix_max);
neigh_pairs(neigh_pairs==0) = [];
neigh_pairs_min = mod(neigh_pairs,K);
neigh_pairs_max = (neigh_pairs-neigh_pairs_min)/K;

if isrow(neigh_pairs_min)
    neigh_pairs_min = neigh_pairs_min';
end
if isrow(neigh_pairs_max)
    neigh_pairs_max = neigh_pairs_max';
end

n_pairs    = 100;
n_triplets = 100;
curr_cands = mex_get_tree_cands(double(lp)-1, double(ms)-1,...
                                 neigh_pairs_min-1, neigh_pairs_max-1,...
                                 [n_pairs, n_triplets]);
assert(isequal(curr_cands{1},[6 7
                              4 8
                              4 6
                              1 2
                              3 7
                              3 4]))
assert(isempty(curr_cands{2}))

%%  Simple test 2
lp = [1 2 3 4];
ms = [1 2 5
      3 4 6
      5 6 7];

% Get all pairs of neighboring leave regions
[~, idx_neighbors] = seg2gridbmap(lp);
K = max(idx_neighbors.matrix_max(:)) + 1;
neigh_pairs = unique(idx_neighbors.matrix_min+K*idx_neighbors.matrix_max);
neigh_pairs(neigh_pairs==0) = [];
neigh_pairs_min = mod(neigh_pairs,K);
neigh_pairs_max = (neigh_pairs-neigh_pairs_min)/K;

if isrow(neigh_pairs_min)
    neigh_pairs_min = neigh_pairs_min';
end
if isrow(neigh_pairs_max)
    neigh_pairs_max = neigh_pairs_max';
end

n_pairs    = 100;
n_triplets = 100;
curr_cands = mex_get_tree_cands(double(lp)-1, double(ms)-1,...
                                neigh_pairs_min-1, neigh_pairs_max-1,...
                                [n_pairs, n_triplets]);
assert(isequal(curr_cands{1},[3 5
                              2 6
                              2 3]))
assert(isempty(curr_cands{2}))

%% Simple with 4-tuples
lp = [1 2 3 4 5
      1 3 3 4 5
      1 4 4 4 5
      1 5 5 5 5];
ms = [1 2 6
      3 6 7
      4 7 8
      5 8 9];

% Get all pairs of neighboring leave regions
[~, idx_neighbors] = seg2gridbmap(lp);
K = max(idx_neighbors.matrix_max(:)) + 1;
neigh_pairs = unique(idx_neighbors.matrix_min+K*idx_neighbors.matrix_max);
neigh_pairs(neigh_pairs==0) = [];
neigh_pairs_min = mod(neigh_pairs,K);
neigh_pairs_max = (neigh_pairs-neigh_pairs_min)/K;

if isrow(neigh_pairs_min)
    neigh_pairs_min = neigh_pairs_min';
end
if isrow(neigh_pairs_max)
    neigh_pairs_max = neigh_pairs_max';
end

n_cands    = [100, 100, 100, 100];
curr_cands = mex_get_tree_cands(double(lp)-1, double(ms)-1,...
                                neigh_pairs_min-1, neigh_pairs_max-1,...
                                n_cands);
assert(isequal(curr_cands{1},[4 5
                              5 7
                              3 4
                              5 6
                              4 6
                              1 5
                              1 4
                              1 3
                              2 3]))
assert(isequal(curr_cands{2},[3 4 5
                              4 5 6
                              1 4 5
                              1 3 5
                              1 3 4
                              2 3 4]))
assert(isequal(curr_cands{3}, [1 3 4 5
                               2 3 4 5]))
assert(isempty(curr_cands{4}))

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function dataSize = getSize(MRSIStruct)
    dataSize = MRSIStruct.sz;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function a=rlogit(p)
% written by Issam El Naqa Spring 2003
% Extracted for generalized use 2005, AJH
%
% Copyright 2010, Joseph O. Deasy, on behalf of the DREES development team.
%
% This file is part of the Dose Response Explorer System (DREES).
%
% DREES development has been led by:  Issam El Naqa, Aditya Apte, Gita Suneja, and Joseph O. Deasy.
%
% DREES has been financially supported by the US National Institutes of Health under multiple grants.
%
% DREES is distributed under the terms of the Lesser GNU Public License.
%
%     This version of DREES is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% DREES is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with DREES.  If not, see <http://www.gnu.org/licenses/>.

    a = log(p ./ (1-p));

return

author: mvallieres

repo: radiomics


text:
function [dtdx,dtrp,vart]=trop_model_prec(time,pos,azel,x)

zhd=tropmodel(pi/2,pos,0); %zenith hydrostatic delay

[mh,mw]=tropmapf(time,pos,azel); %hydrostatic and wet projection coefficient

if azel(2)>0
    %equation: m_w=m_0+m_0*cot(el)*(Gn*cos(az)+Ge*sin(az))
    cotz=1.0/tan(azel(2));
    grad_n=mw*cotz*cos(azel(1)); %north wet projection coefficient
    grad_e=mw*cotz*sin(azel(1)); %east wet projection coefficient
    mw=mw+grad_n*x(2)+grad_e*x(3); %total wet projection coefficient
    dtdx(2)=grad_n*(x(1)-zhd); %north wet delay
    dtdx(3)=grad_e*(x(1)-zhd); %east wet delay
end

dtdx(1)=mw;
dtrp=mh*zhd+mw*(x(1)-zhd);
vart=0.01^2;

return;

repo: GINav


text:
clear; clc;

%% dataset related arguments
train_path = 'data/ml1m/validation/train.txt';
test_path = 'data/ml1m/validation/test.txt';
rate_once = 1;           % whether each item will only be rated once by each user
L = 5;                   % length of sequence
T = 3;                   % number of targets

%% training arguments
seed = 1234;             % random seed
n_iter = 50;             % number of iteration
learning_rate = 1e-3;    % learning rate for ADAM
l2 = 1e-6;               % l2 norm weight
neg_samples = 3;         % negative samples per target
early_stop = true;       % whether to perform early stop during training

%% Caser related arguments
d = 50;                 % number of latent dimensions
nv = 4;                 % number of vertical filters
nh = 16;                % number of horizontal filters
ac_conv = 'relu';       % activation function for convolution layer (i.e., phi_c in paper)
ac_fc = 'relu';         % activation function for fully-connected layer (i.e., phi_a in paper)
drop_rate = 0.5;        % drop ratio when performing dropout

%% Perform training
args = struct(...
        'trainpath',train_path,...
        'testpath',test_path,...
        'rateonce',rate_once,...
        'L',L,...
        'T', T,...
        'seed',seed,...
        'niter',n_iter,...
        'lrate',learning_rate,...
        'l2',l2,...
        'negsample',neg_samples,...
        'earlystop',early_stop,...
        'd',d,...
        'nv',nv,...
        'nh',nh,...
        'acconv',ac_conv,...
        'acfc',ac_fc,...
        'droprate',drop_rate...
        );

caser_train(args);

author: graytowne

repo: caser

save_path: github-repos/MATLAB/graytowne-caser


text:
function cmap3D = MakeDiagonal3Dcolormap(huex,satmin,pw,res,plotdemo)
% draw a square color swatch
% color for upper left corner assigned according to input 'huex'
% lower right: opposite hue of input color, based on the hsv color wheel
% bottem left - upper right diagonal: grayscale
% (bottom left: black; upper right: white)
%%
if ~exist('huex','var'),
    huex = 0;%120/360; % 120/360 for green/magenta; e.g. 0/360 for red/cyan
end
if ~exist('satmin','var'),
    satmin = 0;%0.1;
end
if ~exist('pw','var')
    pw = 1; % scaling of the diagnoal, reminiscent of contrast; 1 is linear
end
if ~exist('res','var'),
    res = 32;
end
%%
cmap3D = zeros(res,res,res,3);
hue = zeros(res,res,res);
sat = zeros(res,res,res);
val = zeros(res,res,res);

for i = 1:res
    for j = 1:res
        for k = 1:res
            % set value
            val(i,j,k) = min((sqrt(i^2+j^2+k^2)/res)^pw,1);

            % set hue
            [~,ix] = max([i,j,k]);
            switch ix
                case 1
                    hue(i,j,k) = huex;
                case 2
                    hue(i,j,k) = huex+120/360;
                case 3
                    hue(i,j,k) = huex+240/360;
            end

            % set saturation
             switch ix
                case 1
                    d = max(0,i-sqrt(j^2+k^2));
                case 2
                    d = max(0,j-sqrt(k^2+i^2));
                case 3
                   d = max(0,k-sqrt(i^2+j^2));
            end
            sat(i,j,k) = max(satmin,d/res);

            cmap3D(i,j,k,:) = hsv2rgb([hue(i,j,k), sat(i,j,k), val(i,j,k)]);
        end
    end
end
%%
if nargin>4
    if plotdemo
        %%
        [X,Y,Z] = ind2sub([res,res,res],1:res^3);
        cmap3D_flat = reshape(cmap3D,res*res*res,3); % for efficient indexing
        figure;scatter3(X,Y,Z,5,cmap3D_flat(1:res^3,:));%,'filled')
%         axis xy
        axis off
        axis equal
    end
end
end

author: xiuyechen

repo: FishExplorer


text:
function mdf = calcMDF(mori,varargin)
% computes an MDF from individuel orientations or misorientations
%
% The function *calcMDF* applies one of the following algorithms to compute
% an MDF from a list of orientations.
%
% # direct kernel density estimation
% # kernel density estimation via Fourier series
% # Bingham estimation
%
% Syntax
%
%   % use kernel density estimation with a 10 degree kernel
%   mori = grains.boundary.misorientation
%   mdf = calcMDF(mori,'halfwidth',10*degree)
%
%   % compute an uncorrelated MDF
%   mdf = calcMDF(grains('phase1').meanorientation)
%
%   % use grain area as weights for the orientations
%   mdf = calcMDF(grains('phase1').meanOrientation,'weights',grains('phase1').diameter)
%
%   % use a specific kernel
%   psi = SO3AbelPoissonKernel('halfwidth',10*degree)
%   mdf = calcMDF(mori,'kernel',psi)
%
%   % compute the MDF as a Fourier series of order 16
%   mdf = calcMDF(mori,'order',16)
%
% Input
%  ori  - @orientation
%  mori - misorientation
%
% Output
%  mdf - @SO3Fun
%
% Options
%  weights    - list of weights for the orientations
%  halfwidth  - halfwidth of the kernel function
%  resolution - resolution of the grid where the MDF is approximated
%  kernel     - @SO3Kernel function (default -- SO3 de la Valee Poussin kernel)
%  order      - order up to which Fourier coefficients are calculated
%
% Flags
%  silent           - no output
%  exact            - no approximation to a corser grid
%  Fourier          - force Fourier method
%  Bingham          - model bingham mdf
%  noFourier        - no Fourier method
%
% See also
% orientation/calcFourierMDF orientation/calcKernelMDF orientation/calcBinghamMDF ebsd_demo EBSD2mdf EBSDSimulation_demo

warning('The command calcODF is depreciated! Please use calcDensity instead.')

% if orientations have been specified
if isa(mori.SS,'specimenSymmetry')

  % compute an ODF first
  odf1 = calcFourierODF(mori,varargin{:});

  if nargin > 1 && isa(varargin{1},'orientation')

    % maybe a second ODF is needed
    odf2 = calcFourierODF(varargin{1},varargin{:});

    % compute the MDF
    mdf = calcMDF(odf1,odf2);

  else

    % compute the MDF
    mdf = calcMDF(odf1);

  end

  return
end

if check_option(varargin,'antipoal') && mori.CS == mori.SS
  mori = [mori(:);inv(mori(:))];
end

mdf = calcODF(mori,'halfwidth',7.5*degree,varargin{:});

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
% test for the tensor structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

path(path, 'toolbox/');

n = 128;
rep = 'faces/';
rep = '';
name = 'polygons_blurred';
name = 'rakan';
name = 'barb';
name = 'lena';
M0 = load_image([rep name]);
M = sum(M0,3);
M = M(end/2-n/2+1:end/2+n/2,end/2-n/2+1:end/2+n/2);
M0 = M0(end/2-n/2+1:end/2+n/2,end/2-n/2+1:end/2+n/2,:);
M = M + 0.1*randn(n);

sigma1 = 1.5;
sigma2 = 4;

options.m_theta = 16;

options.use_anisotropic=0;
options.use_renormalization = 0;
disp('Computing isotropic tensors.');
H = compute_structure_tensor(M,sigma1,sigma2, options);

options.use_anisotropic=0;
options.use_renormalization = 1;
options.sigmat = 0.4;
disp('Computing anisotropic tensors.');
H1 = compute_structure_tensor(M,sigma1,sigma2, options);

clf;
subplot(1,2,1);
plot_tensor_field(H, M0, 5);
title('Isotropic');

subplot(1,2,2);
plot_tensor_field(H1, M0, 5);
title('Anisotropic');

rep = 'results/';
if ~exist(rep)
    mkdir(rep);
end
saveas( gcf, [rep name '_tensor.png'], 'png' );

return;

[e1,e2,l1,l2] = perform_tensor_decomp(H1);

an = l1./l2; ae = l1.*l2;
ae_target = 4; ae = ae_target * ae./mean(ae(:));
an_median = 4;
an = 2.^( log2(an_median)/log2(median(an(:))) * log2(an)   );
l1 = sqrt(ae.*an);
l2 = sqrt(ae./an);
H2 = perform_tensor_recomp(e1,e2,l1,l2);

close all;
plot_tensor_field(H, M0, 5);
figure;
plot_tensor_field(H2, M0, 5);

author: gpeyre

repo: matlab-toolboxes


text:
% StackExchange Signal Processing Q81493
% https://dsp.stackexchange.com/questions/81493
% Applying 2D Sinc Interpolation for Upsampling in the Fourier Domain (DFT / FFT)
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     29/12/2021
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

numRowsI = 5000;
numColsI = 5200;

numRowsO = 10000;
numColsO = 10400;

sincRadius = 5;

%% Generate / Load Data

mX      = GenTest([numRowsI, numColsI], sincRadius);
mYRef   = GenTest([numRowsO, numColsO], sincRadius);

mY = DftReSample2D(mX, [numRowsO, numColsO]);

%% Analysis

disp(['The interpolation error is given by: ', num2str(max(abs(mYRef - mY), [], 'all'))]);

%% Auxilizary Function

function [ mX ] = GenTest( vSize, sincRadius )

vX = linspace(-sincRadius, sincRadius, vSize(2) + 1);
vX(end) = [];
vY = linspace(-sincRadius, sincRadius, vSize(1) + 1);
vY = vY(:);
vY(end) = [];

% mX = abs(vX) + abs(vY) + sinc(sqrt(vX .^2 + vY .^2));
mX = sinc(sqrt(vX .^2 + vY .^2));

end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function ok = write_hdr_raw(fname,hdr,be)
% Write a NIFTI-1 .hdr file.
% FORMAT ok = write_hdr_raw(fname,hdr,be)
% fname - filename of image
% hdr   - a structure containing hdr info
% be    - whether big-endian or not
% ok    - status (1=good, 0=bad)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

[pth,nam,ext] = fileparts(fname);
if isempty(pth), pth = pwd; end

if isfield(hdr,'magic')
    org = niftistruc;
    switch deblank(hdr.magic)
    case {'ni1'}
        hname = fullfile(pth,[nam '.hdr']);
    case {'n+1'}
        hname = fullfile(pth,[nam '.nii']);
    otherwise
        error('Bad header.');
    end;
else
    org   = mayostruc;
    hname = fullfile(pth,[nam '.hdr']);
end;

if nargin >=3
    if be, mach = 'ieee-be';
    else   mach = 'ieee-le';
    end;
else       mach = 'native';
end;

ok = true;
if spm_existfile(hname),
    fp = fopen(hname,'r+',mach);
else
    fp = fopen(hname,'w+',mach);
end
if fp == -1,
    ok = false;
    return;
end

for i=1:length(org)
    if isfield(hdr,org(i).label),
        dat = hdr.(org(i).label);
        if length(dat) ~= org(i).len,
            if length(dat)< org(i).len,
                dat = [dat(:) ; zeros(org(i).len-length(dat),1)];
            else
                dat = dat(1:org(i).len);
            end;
        end;
    else
        dat = org(i).def;
    end;
    % fprintf('%s=\n',org(i).label)
    % disp(dat)
    len = fwrite(fp,dat,org(i).dtype.prec);
    if len ~= org(i).len,
        ok = false;
    end;
end;
fclose(fp);
if ~ok,
     fprintf('There was a problem writing to the header of\n');
     fprintf('"%s"\n', fname);
end;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function transfo = createBasisTransform3d(source, target)
%CREATEBASISTRANSFORM3D Compute matrix for transforming a basis into another basis
%
%   TRANSFO = createBasisTransform3d(SOURCE, TARGET) will create a 4-by-4
%   transformation matrix representing the transformation from SOURCE basis
%   to TARGET basis.
%    SOURCE and TARGET are either standard 1-by-9 geom3d PLANE
%    representations of the form: [x0 y0 z0  ex1 ey1 ez1  ex2 ey2 ez2]
%     OR
%    SOURCE and TARGET may be any string such as 'global' or 'g' in which
%    case they represent the global plane [0 0 0 1 0 0 0 1 0].
%
%   The resulting TRANSFO matrix is such that a point expressed with
%   coordinates of the first basis will be represented by new coordinates
%   P2 = transformPoint3d(P1, TRANSFO) in the target basis.
%
%   Either (or both) SOURCE or TARGET may be an N-by-9 set of N planes. In
%   that case, TRANSFO will be a 4-by-4-by-N array of N transformation
%   matrices.
%
%   Example:
%     % Calculate local plane coords. of a point given in global coords
%     Plane = [10 10 10  1 0 0  0 1 0];
%     Tform = createBasisTransform3d('global', Plane);
%     PT_IN_PLANE = transformPoint3d([3 8 2], Tform)
%     PT_IN_PLANE =
%         13  18  12
%
%   See also
%     transforms3d, transformPoint3d, planePosition, createBasisTransform

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-12-03,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% HISTORY
% 2013-07-03 added support for multiple inputs (Sven Holcombe)
% 2017-10-16 rewrite

% size of input arguments
srcSz = size(source, 1);
tgtSz = size(target, 1);
maxSz = max(srcSz, tgtSz);

% check case of multiple inputs
if maxSz > 1
    [t1, t2] = deal( bsxfun(@times, eye(4), ones(1,1,maxSz)) );
    if srcSz > 1
        source = permute(source, [3 2 1]);
    end
    if tgtSz > 1
        target = permute(target, [3 2 1]);
    end
else
    [t1, t2] = deal(eye(4));
end

% Place source and target planes into t1 and t2 t-form matrices. If either
% input is non-numeric it is assumed to mean 'global', or identity t-form.
if isnumeric(source)
    if maxSz > 1 && srcSz == 1
        source = bsxfun(@times, source, ones(1,1,maxSz));
    end
    t1(1:3, 1, :) = source(1, 4:6, :);
    t1(1:3, 2, :) = source(1, 7:9, :);
    t1(1:3, 3, :) = crossProduct3d(source(1,4:6,:), source(1,7:9,:));
    t1(1:3, 4, :) = source(1, 1:3, :);
end
if isnumeric(target)
    if maxSz > 1 && tgtSz == 1
        target = bsxfun(@times, target, ones(1,1,maxSz));
    end
    t2(1:3, 1, :) = target(1, 4:6, :);
    t2(1:3, 2, :) = target(1, 7:9, :);
    t2(1:3, 3, :) = crossProduct3d(target(1,4:6,:), target(1,7:9,:));
    t2(1:3, 4, :) = target(1, 1:3, :);
end

% compute transform matrix
transfo = zeros(4, 4, maxSz);
for i = 1:maxSz
    % coordinate of four reference points in source basis
    po = t1(1:3, 4, i)';
    px = po + t1(1:3, 1, i)';
    py = po + t1(1:3, 2, i)';
    pz = po + t1(1:3, 3, i)';

    % express coordinates of reference points in the new basis
    t2i = inv(t2(:,:,i));
    pot = transformPoint3d(po, t2i);
    pxt = transformPoint3d(px, t2i);
    pyt = transformPoint3d(py, t2i);
    pzt = transformPoint3d(pz, t2i);

    % compute direction vectors in new basis
    vx = pxt - pot;
    vy = pyt - pot;
    vz = pzt - pot;

    % concatenate result in a 4-by-4 affine transform matrix
    transfo(:,:,i) = [vx' vy' vz' pot' ; 0 0 0 1];
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function cvx_optval = lambda_sum_largest( x, k )

%LAMBDA_SUM_LARGEST   Internal cvx version.

narginchk(2,2);
n = size( x, 1 );
if ndims( x ) > 2 || n ~= size( x, 2 ), %#ok

    error( 'First input must be a square matrix.' );

elseif ~isnumeric( k ) || numel( k ) ~= 1 || ~isreal( k ),

    error( 'Second input must be a constant real scalar.' );

elseif cvx_isconstant( x ),

    cvx_optval = cvx( lambda_max( cvx_constant( x ), k ) );

elseif ~cvx_isaffine( x ),

    error( 'Discipliend convex programming error:\n    LAMBDA_SUM_LARGEST is convex and nonmonotonic, so its input must be affine.' );

elseif k <= 0,

    cvx_optval = 0;

elseif k >= size( x, 1 ),

    cvx_optval = trace( x );

else

    S = [];
    cvx_begin
        variable S(n,n) symmetric
        S == semidefinite(n); %#ok
        minimize( k * lambda_max( x - S ) + trace( S ) );
    cvx_end

end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


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

%% Figure 26.1

% load sample EEG dataset
load sampleEEGdata

% names of the channels you want to synchronize
channel1 = 'p1';
channel2 = 'pz';

% create complex Morlet wavelet
center_freq = 5; % in Hz
time        = -1:1/EEG.srate:1; % time for wavelet
wavelet     = exp(2*1i*pi*center_freq.*time) .* exp(-time.^2./(2*(4/(2*pi*center_freq))^2))/center_freq;
half_of_wavelet_size = (length(time)-1)/2;

% FFT parameters
n_wavelet     = length(time);
n_data        = EEG.pnts;
n_convolution = n_wavelet+n_data-1;

% FFT of wavelet
fft_wavelet = fft(wavelet,n_convolution);

% initialize output time-frequency data
phase_data = zeros(2,EEG.pnts);
real_data  = zeros(2,EEG.pnts);

% find channel indices
chanidx = zeros(1,2); % always initialize!
chanidx(1) = find(strcmpi(channel1,{EEG.chanlocs.labels}));
chanidx(2) = find(strcmpi(channel2,{EEG.chanlocs.labels}));

% run convolution and extract filtered signal (real part) and phase
for chani=1:2
    fft_data = fft(squeeze(EEG.data(chanidx(chani),:,1)),n_convolution);
    convolution_result_fft = ifft(fft_wavelet.*fft_data,n_convolution) * sqrt(4/(2*pi*center_freq));
    convolution_result_fft = convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size);

    % collect real and phase data
    phase_data(chani,:) = angle(convolution_result_fft);
    real_data(chani,:)  = real(convolution_result_fft);
end

% open and name figure
figure, set(gcf,'Name','Movie magic minimizes the mystery.','Number','off');

% draw the filtered signals
subplot(321)
filterplotH1 = plot(EEG.times(1),real_data(1,1),'b');
hold on
filterplotH2 = plot(EEG.times(1),real_data(2,1),'m');
set(gca,'xlim',[EEG.times(1) EEG.times(end)],'ylim',[min(real_data(:)) max(real_data(:))])
xlabel('Time (ms)')
ylabel('Voltage (\muV)')
title([ 'Filtered signal at ' num2str(center_freq) ' Hz' ])

% draw the phase angle time series
subplot(322)
phaseanglesH1 = plot(EEG.times(1),phase_data(1,1),'b');
hold on
phaseanglesH2 = plot(EEG.times(1),phase_data(2,1),'m');
set(gca,'xlim',[EEG.times(1) EEG.times(end)],'ylim',[-pi pi]*1.1,'ytick',-pi:pi/2:pi)
xlabel('Time (ms)')
ylabel('Phase angle (radian)')
title('Phase angle time series')

% draw phase angle differences in cartesian space
subplot(323)
filterplotDiffH1 = plot(EEG.times(1),real_data(1,1)-real_data(2,1),'b');
set(gca,'xlim',[EEG.times(1) EEG.times(end)],'ylim',[-10 10])
xlabel('Time (ms)')
ylabel('Voltage (\muV)')
title([ 'Filtered signal at ' num2str(center_freq) ' Hz' ])

% draw the phase angle time series
subplot(324)
phaseanglesDiffH1 = plot(EEG.times(1),phase_data(1,1)-phase_data(2,1),'b');
set(gca,'xlim',[EEG.times(1) EEG.times(end)],'ylim',[-pi pi]*2.2,'ytick',-2*pi:pi/2:pi*2)
xlabel('Time (ms)')
ylabel('Phase angle (radian)')
title('Phase angle time series')

% draw phase angles in polar space
subplot(325)
polar2chanH1 = polar([phase_data(1,1) phase_data(1,1)]',repmat([0 1],1,1)','b');
hold on
polar2chanH2 = polar([phase_data(1,1) phase_data(2,1)]',repmat([0 1],1,1)','m');
title('Phase angles from two channels')

% draw phase angle differences in polar space
subplot(326)
polarAngleDiffH = polar([zeros(1,1) phase_data(2,1)-phase_data(1,1)]',repmat([0 1],1,1)','k');
title('Phase angle differences from two channels')

% now update plots at each timestep
% Note: in/decrease skipping by 10 to speed up/down the movie
for ti=1:10:EEG.pnts

    % update filtered signals
    set(filterplotH1,'XData',EEG.times(1:ti),'YData',real_data(1,1:ti))
    set(filterplotH2,'XData',EEG.times(1:ti),'YData',real_data(2,1:ti))

    % update cartesian plot of phase angles
    set(phaseanglesH1,'XData',EEG.times(1:ti),'YData',phase_data(1,1:ti))
    set(phaseanglesH2,'XData',EEG.times(1:ti),'YData',phase_data(2,1:ti))

    % update cartesian plot of phase angles differences
    set(phaseanglesDiffH1,'XData',EEG.times(1:ti),'YData',phase_data(1,1:ti)-phase_data(2,1:ti))
    set(filterplotDiffH1,'XData',EEG.times(1:ti),'YData',real_data(1,1:ti)-real_data(2,1:ti))

    subplot(325)
    cla
    polar(repmat(phase_data(1,1:ti),1,2)',repmat([0 1],1,ti)','b');
    hold on
    polar(repmat(phase_data(2,1:ti),1,2)',repmat([0 1],1,ti)','m');

    subplot(326)
    cla
    polar(repmat(phase_data(2,1:ti)-phase_data(1,1:ti),1,2)',repmat([0 1],1,ti)','k');

    drawnow
end

%% Figure 26.2

figure
subplot(221)
polar(repmat(phase_data(2,:)-phase_data(1,:),1,2)',repmat([0 1],1,EEG.pnts)','k');
title([ 'Phase synchronization: ' num2str(abs(mean(exp(1i*(diff(phase_data,1)))))) ])

new_phase_data = phase_data;
for i=2:4
    subplot(2,2,i)

    % add random phase offset
    new_phase_data(1,:) = new_phase_data(1,:)+rand*pi;

    % plot again
    polar(repmat(new_phase_data(2,:)-new_phase_data(1,:)+pi/2,1,2)',repmat([0 1],1,EEG.pnts)','k');
    title([ 'Phase synchronization: ' num2str(abs(mean(exp(1i*(diff(new_phase_data,1)))))) ])
end

%% Figure 26.3

% note: see commented line "time_window_idx..." below for panels C and D

channel1 = 'fz';
channel2 = 'o1';

freqs2use  = logspace(log10(4),log10(30),15); % 4-30 Hz in 15 steps
times2save = -400:20:800;
timewindow = linspace(1.5,3,length(freqs2use)); % number of cycles on either end of the center point (1.5 means a total of 3 cycles))
baselinetm = [-400 -200];

% wavelet and FFT parameters
time          = -1:1/EEG.srate:1;
half_wavelet  = (length(time)-1)/2;
num_cycles    = logspace(log10(4),log10(8),length(freqs2use));
n_wavelet     = length(time);
n_data        = EEG.pnts*EEG.trials;
n_convolution = n_wavelet+n_data-1;

% time in indices
times2saveidx = dsearchn(EEG.times',times2save');
baselineidx   = dsearchn(times2save',baselinetm');

chanidx    = zeros(1,2); % always initialize!
chanidx(1) = find(strcmpi(channel1,{EEG.chanlocs.labels}));
chanidx(2) = find(strcmpi(channel2,{EEG.chanlocs.labels}));

% initialize
ispc = zeros(length(freqs2use),length(times2save));
ps   = zeros(length(freqs2use),length(times2save));

% data FFTs
data_fft1 = fft(reshape(EEG.data(chanidx(1),:,:),1,n_data),n_convolution);
data_fft2 = fft(reshape(EEG.data(chanidx(2),:,:),1,n_data),n_convolution);

for fi=1:length(freqs2use)

    % create wavelet and take FFT
    s = num_cycles(fi)/(2*pi*freqs2use(fi));
    wavelet_fft = fft( exp(2*1i*pi*freqs2use(fi).*time) .* exp(-time.^2./(2*(s^2))) ,n_convolution);

    % phase angles from channel 1 via convolution
    convolution_result_fft = ifft(wavelet_fft.*data_fft1,n_convolution);
    convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
    phase_sig1 = angle(reshape(convolution_result_fft,EEG.pnts,EEG.trials));

    % phase angles from channel 2 via convolution
    convolution_result_fft = ifft(wavelet_fft.*data_fft2,n_convolution);
    convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
    phase_sig2 = angle(reshape(convolution_result_fft,EEG.pnts,EEG.trials));

    % phase angle differences
    phase_diffs = phase_sig1-phase_sig2;

    % compute ICPS over trials
    ps(fi,:) = abs(mean(exp(1i*phase_diffs(times2saveidx,:)),2));

    % compute time window in indices for this frequency
    time_window_idx = round((1000/freqs2use(fi))*timewindow(fi)/(1000/EEG.srate));
%     time_window_idx = round(300/(1000/EEG.srate)); % set 300 to 100 for figure 3c/d

    for ti=1:length(times2save)

        % compute phase synchronization
        phasesynch = abs(mean(exp(1i*phase_diffs(times2saveidx(ti)-time_window_idx:times2saveidx(ti)+time_window_idx,:)),1));

        % average over trials
        ispc(fi,ti) = mean(phasesynch);
    end
end % end frequency loop

figure
contourf(times2save,freqs2use,ispc-repmat(mean(ispc(:,baselineidx(1):baselineidx(2)),2),1,size(ispc,2)),20,'linecolor','none')
set(gca,'clim',[-.08 .08],'yscale','log','ytick',round(logspace(log10(freqs2use(1)),log10(freqs2use(end)),8)))
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

figure
plot(freqs2use,(1000./freqs2use).*timewindow*2,'o-','markerface','k')
hold on
plot(freqs2use,(1000./freqs2use).*timewindow(1)*2,'ro-','markerface','m')
ylabel('Window width (ms)'), xlabel('Frequency (Hz)')
legend({'variable windows';'fixed 3*f window'})

%% Figure 26.4

figure
for i=1:8
    subplot(8,1,i)
    plot(phase_sig1(1:200,i)-phase_sig2(1:200,i))
end

figure
subplot(121)
polar(repmat(phase_sig1(1:200,1)-phase_sig2(1:200,1),2,1),repmat([0 1]',200,1),'k');
title('Phase angle differences over time')

subplot(122)
polar(repmat(phase_sig1(100,1:i)-phase_sig2(100,1:i),2,1),repmat([0 1]',1,i),'k');
title('Phase angle differences over trials')

%% Figure 26.5

figure
contourf(times2save,freqs2use,bsxfun(@minus,ps,mean(ps(:,baselineidx(1):baselineidx(2)),2)),20,'linecolor','none')
set(gca,'clim',[-.2 .2],'yscale','log','ytick',round(logspace(log10(freqs2use(1)),log10(freqs2use(end)),8)),'xlim',[-300 800])
xlabel('Time (ms)'), ylabel('Frequency (Hz)')

%% figure 26.6

time2use = 300; % ms
niterations = 50; % you can decrease this to make the code a bit faster

% initialize
ispcByNandF = zeros(length(freqs2use),EEG.trials);
time2useidx = dsearchn(times2save',time2use);

% data FFTs
data_fft1 = fft(reshape(EEG.data(chanidx(1),:,:),1,n_data),n_convolution);
data_fft2 = fft(reshape(EEG.data(chanidx(2),:,:),1,n_data),n_convolution);

for fi=1:length(freqs2use)

    % create wavelet and take FFT
    s = num_cycles(fi)/(2*pi*freqs2use(fi));
    wavelet_fft = fft( exp(2*1i*pi*freqs2use(fi).*time) .* exp(-time.^2./(2*(s^2))) ,n_convolution);

    % phase angles from channel 1 via convolution
    convolution_result_fft = ifft(wavelet_fft.*data_fft1,n_convolution);
    convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
    phase_sig1 = angle(reshape(convolution_result_fft,EEG.pnts,EEG.trials));

    % phase angles from channel 2 via convolution
    convolution_result_fft = ifft(wavelet_fft.*data_fft2,n_convolution);
    convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
    phase_sig2 = angle(reshape(convolution_result_fft,EEG.pnts,EEG.trials));

    % phase angle differences
    phase_diffs = phase_sig1-phase_sig2;

    for n=1:EEG.trials
        % multiple iterations to select different random sets of trials
        for iteri=1:niterations
            trials2use = randsample(EEG.trials,n);
            ispcByNandF(fi,n) = ispcByNandF(fi,n) + mean(abs(mean(exp(1i*phase_diffs(times2saveidx(time2useidx)-time_window_idx:times2saveidx(time2useidx)+time_window_idx,trials2use)),2)),1);
        end
    end
end

figure
plot(1:EEG.trials,ispcByNandF/iteri)
xlabel('Trials')
ylabel('ICPS-trials')

%% Figure 26.7

% initialize
data4test  = zeros(2,EEG.pnts,EEG.trials);
data4power = zeros(2,EEG.pnts,EEG.trials);

amp_mod = 0.00001;

for triali=1:EEG.trials
    % each trial is a random channel and trial
    trialdata1 = EEG.data(chanidx(1),:,triali);
    trialdata2 = EEG.data(chanidx(2),:,triali);

    % Un/comment the next line of code for band-pass filtered data.
    % This uses the eegfilt function, which is part of the eeglab toolbox.
    % You can also replace this function with your preferred filter method (chapter 14).
    trialdata1 = eegfilt(double(trialdata1),EEG.srate,10,20);
    trialdata2 = eegfilt(double(trialdata2),EEG.srate,10,20);

    % phase angle differences, with and without amplitude dampening
    data4test(1,:,triali) = angle(hilbert(trialdata1)) - angle(hilbert(trialdata2));
    data4test(2,:,triali) = angle(hilbert(trialdata1)) - angle(hilbert(trialdata2*amp_mod));

    data4power(1,:,triali) = abs(hilbert(trialdata2)).^2;
    data4power(2,:,triali) = abs(hilbert(trialdata2*amp_mod)).^2;
end

% compute ITPC
ispc_nomod = abs(mean(exp(1i*data4test(1,:,:)),3));
ispc_mod   = abs(mean(exp(1i*data4test(2,:,:)),3));

% compute power
power = squeeze(mean(data4power,3));

% plot!
figure
subplot(311)
plot(EEG.times,trialdata2)
hold on
plot(EEG.times,trialdata2*amp_mod,'r')
title('Amplitude modulator')

subplot(312)
plot(EEG.times,data4test(1,:,10))
hold on
plot(EEG.times,data4test(2,:,10),'r')
axis tight
set(gca,'ytick',-2*pi:pi:2*pi)
title('Example trials')

subplot(313)
plot(EEG.times,ispc_mod,'ro')
hold on
h=plotyy(EEG.times,ispc_nomod,EEG.times,squeeze(mean(data4power(1,:,:),3)));
legend({'amplitude modulation';'no amp mod'})
set(h(2),'ylim',[12 36])
set(h(1),'ylim',[0 .4])
xlabel('Time (ms)'), ylabel('ICPS')
title('ICPS')

figure
subplot(121)
plot(power(1,:),ispc_nomod,'.')
axis square
xlabel('Power'), ylabel('ICPS')
title('non-modulated power')

subplot(122)
plot(power(2,:),ispc_mod,'.')
axis square
xlabel('Power'), ylabel('ICPS')
title('modulated power')

%% figure 26.8

% select channels
channel1 = 'fz';
channel2 = 'o1';

% wavelet and FFT parameters
time          = -1:1/EEG.srate:1;
half_wavelet  = (length(time)-1)/2;
n_wavelet     = length(time);
n_data        = EEG.pnts*EEG.trials;
n_convolution = n_wavelet+n_data-1;

chanidx    = zeros(1,2); % always initialize!
chanidx(1) = find(strcmpi(channel1,{EEG.chanlocs.labels}));
chanidx(2) = find(strcmpi(channel2,{EEG.chanlocs.labels}));

% data FFTs
data_fft1 = fft(reshape(EEG.data(chanidx(1),:,:),1,n_data),n_convolution);
data_fft2 = fft(reshape(EEG.data(chanidx(2),:,:),1,n_data),n_convolution);

% initialize
spectcoher = zeros(length(freqs2use),length(times2save));

for fi=1:length(freqs2use)

    % create wavelet and take FFT
    s = num_cycles(fi)/(2*pi*freqs2use(fi));
    wavelet_fft = fft( exp(2*1i*pi*freqs2use(fi).*time) .* exp(-time.^2./(2*(s^2))) ,n_convolution);

    % phase angles from channel 1 via convolution
    convolution_result_fft = ifft(wavelet_fft.*data_fft1,n_convolution);
    convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
    sig1 = reshape(convolution_result_fft,EEG.pnts,EEG.trials);

    % phase angles from channel 2 via convolution
    convolution_result_fft = ifft(wavelet_fft.*data_fft2,n_convolution);
    convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
    sig2 = reshape(convolution_result_fft,EEG.pnts,EEG.trials);

    % compute power and cross-spectral power
    spec1 = mean(sig1.*conj(sig1),2);
    spec2 = mean(sig2.*conj(sig2),2);
    specX = abs(mean(sig1.*conj(sig2),2)).^2;

    % alternative notation for the same procedure, using the Euler-like expression: Me^ik
    %spec1 = mean(abs(sig1).^2,2);
    %spec2 = mean(abs(sig2).^2,2);
    %specX = abs(mean( abs(sig1).*abs(sig2) .* exp(1i*(angle(sig1)-angle(sig2))) ,2)).^2;

    % compute spectral coherence, using only requested time points
    spectcoher(fi,:) = specX(times2saveidx)./(spec1(times2saveidx).*spec2(times2saveidx));

    % yet another equivalent notation, just FYI
    %spec1 = sum(sig1.*conj(sig1),2);
    %spec2 = sum(sig2.*conj(sig2),2);
    %specX = sum(sig1.*conj(sig2),2);
    %spectcoher(fi,:) = abs(specX(times2saveidx)./sqrt(spec1(times2saveidx).*spec2(times2saveidx))).^2;

    % imaginary coherence
    %spec1 = sum(sig1.*conj(sig1),2);
    %spec2 = sum(sig2.*conj(sig2),2);
    %specX = sum(sig1.*conj(sig2),2);
    % spectcoher(fi,:) = abs(imag(specX(times2saveidx)./sqrt(spec1(times2saveidx).*spec2(times2saveidx))));
end

figure
subplot(121)
contourf(times2save,freqs2use,spectcoher,20,'linecolor','none') %
set(gca,'clim',[0 .2],'yscale','log','ytick',round(logspace(log10(freqs2use(1)),log10(freqs2use(end)),8)),'xlim',[times2save(1) times2save(end)])
title('"Raw" spectral coherence')

subplot(122)
contourf(times2save,freqs2use,spectcoher-repmat(mean(spectcoher(:,baselineidx(1):baselineidx(2)),2),1,size(spectcoher,2)),20,'linecolor','none') %
set(gca,'clim',[-.1 .1],'yscale','log','ytick',round(logspace(log10(freqs2use(1)),log10(freqs2use(end)),8)),'xlim',[times2save(1) times2save(end)])
xlabel('Time (ms)'), ylabel('Frequency (Hz)')
title('Baseline-subtracted spectral coherence')

%% Figure 26.9

% number of "trials"
n = 100;

figure

subplot(221)
phases = rand(n,1)*pi;
polar([phases; phases],repmat([0 1]',n,1),'k');
pli  = abs(mean(sign(imag(exp(1i*phases)))));
ispc = abs(mean(exp(1i*phases)));
title([ 'PLI=' num2str(pli) ', ISPC=' num2str(ispc) ])

subplot(222)
phases = phases-pi/2;
polar([phases; phases],repmat([0 1]',n,1),'k');
pli  = abs(mean(sign(imag(exp(1i*phases)))));
ispc = abs(mean(exp(1i*phases)));
title([ 'PLI=' num2str(pli) ', ISPC=' num2str(ispc) ])

subplot(223)
phases = rand(n,1)/2+pi/3+.25;
polar([phases; phases],repmat([0 1]',n,1),'k');
pli  = abs(mean(sign(imag(exp(1i*phases)))));
ispc = abs(mean(exp(1i*phases)));
title([ 'PLI=' num2str(pli) ', ISPC=' num2str(ispc) ])

subplot(224)
phases = phases-pi/2;
polar([phases; phases],repmat([0 1]',n,1),'k');
pli  = abs(mean(sign(imag(exp(1i*phases)))));
ispc = abs(mean(exp(1i*phases)));
title([ 'PLI=' num2str(pli) ', ISPC=' num2str(ispc) ])

%% Figure 26.10

% select channels
channel1 = 'fz';
channel2 = 'o1';

% specify some time-frequency parameters
freqs2use  = logspace(log10(4),log10(30),15); % 4-30 Hz in 15 steps
times2save = -400:10:800;
timewindow = linspace(1.5,3,length(freqs2use)); % number of cycles on either end of the center point (1.5 means a total of 3 cycles))
baselinetm = [-400 -200];

% wavelet and FFT parameters
time          = -1:1/EEG.srate:1;
half_wavelet  = (length(time)-1)/2;
num_cycles    = logspace(log10(4),log10(8),length(freqs2use));
n_wavelet     = length(time);
n_data        = EEG.pnts*EEG.trials;
n_convolution = n_wavelet+n_data-1;

% time in indices
times2saveidx = dsearchn(EEG.times',times2save');
baselineidxF  = dsearchn(EEG.times',baselinetm');  % for the full temporal resolution data (thanks to Daniel Roberts for finding/reporting this bug here!)
baselineidx   = dsearchn(times2save',baselinetm'); % for the temporally downsampled data

chanidx = zeros(1,2); % always initialize!
chanidx(1) = find(strcmpi(channel1,{EEG.chanlocs.labels}));
chanidx(2) = find(strcmpi(channel2,{EEG.chanlocs.labels}));

% data FFTs
data_fft1 = fft(reshape(EEG.data(chanidx(1),:,:),1,n_data),n_convolution);
data_fft2 = fft(reshape(EEG.data(chanidx(2),:,:),1,n_data),n_convolution);

% initialize
ispc    = zeros(length(freqs2use),EEG.pnts);
pli     = zeros(length(freqs2use),EEG.pnts);
wpli    = zeros(length(freqs2use),EEG.pnts);
dwpli   = zeros(length(freqs2use),EEG.pnts);
dwpli_t = zeros(length(freqs2use),length(times2save));
ispc_t  = zeros(length(freqs2use),length(times2save));

for fi=1:length(freqs2use)

    % create wavelet and take FFT
    s = num_cycles(fi)/(2*pi*freqs2use(fi));
    wavelet_fft = fft( exp(2*1i*pi*freqs2use(fi).*time) .* exp(-time.^2./(2*(s^2))) ,n_convolution);

    % phase angles from channel 1 via convolution
    convolution_result_fft = ifft(wavelet_fft.*data_fft1,n_convolution);
    convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
    sig1 = reshape(convolution_result_fft,EEG.pnts,EEG.trials);

    % phase angles from channel 2 via convolution
    convolution_result_fft = ifft(wavelet_fft.*data_fft2,n_convolution);
    convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
    sig2 = reshape(convolution_result_fft,EEG.pnts,EEG.trials);

    % cross-spectral density
    cdd = sig1 .* conj(sig2);

    % ISPC
    ispc(fi,:) = abs(mean(exp(1i*angle(cdd)),2)); % note: equivalent to ispc(fi,:) = abs(mean(exp(1i*(angle(sig1)-angle(sig2))),2));

    % take imaginary part of signal only
    cdi = imag(cdd);

    % phase-lag index
    pli(fi,:)  = abs(mean(sign(imag(cdd)),2));

    % weighted phase-lag index (eq. 8 in Vink et al. NeuroImage 2011)
    wpli(fi,:) = abs( mean( abs(cdi).*sign(cdi) ,2) )./mean(abs(cdi),2);

    % debiased weighted phase-lag index (shortcut, as implemented in fieldtrip)
    imagsum      = sum(cdi,2);
    imagsumW     = sum(abs(cdi),2);
    debiasfactor = sum(cdi.^2,2);
    dwpli(fi,:)  = (imagsum.^2 - debiasfactor)./(imagsumW.^2 - debiasfactor);

    % compute time window in indices for this frequency
    time_window_idx = round((1000/freqs2use(fi))*timewindow(fi)/(1000/EEG.srate));

    for ti=1:length(times2save)
        imagsum        = sum(cdi(times2saveidx(ti)-time_window_idx:times2saveidx(ti)+time_window_idx,:),1);
        imagsumW       = sum(abs(cdi(times2saveidx(ti)-time_window_idx:times2saveidx(ti)+time_window_idx,:)),1);
        debiasfactor   = sum(cdi(times2saveidx(ti)-time_window_idx:times2saveidx(ti)+time_window_idx,:).^2,1);
        dwpli_t(fi,ti) = mean((imagsum.^2 - debiasfactor)./(imagsumW.^2 - debiasfactor));

        % compute phase synchronization
        phasesynch     = abs(mean(exp(1i*angle(cdd(times2saveidx(ti)-time_window_idx:times2saveidx(ti)+time_window_idx,:))),1));
        ispc_t(fi,ti)  = mean(phasesynch);
    end
end

% baseline subtraction from all measures
ispc    = bsxfun(@minus,ispc,mean(ispc(:,baselineidxF(1):baselineidxF(2)),2)); % not plotted in the book, but you can plot it for comparison with PLI
ispc_t  = bsxfun(@minus,ispc_t,mean(ispc_t(:,baselineidx(1):baselineidx(2)),2));
pli     = bsxfun(@minus,pli,mean(pli(:,baselineidxF(1):baselineidxF(2)),2));
dwpli   = bsxfun(@minus,dwpli,mean(dwpli(:,baselineidxF(1):baselineidxF(2)),2));
dwpli_t = bsxfun(@minus,dwpli_t,mean(dwpli_t(:,baselineidx(1):baselineidx(2)),2));

figure
subplot(221)
contourf(times2save,freqs2use,pli(:,times2saveidx),40,'linecolor','none')
set(gca,'clim',[-.3 .3],'yscale','log','ytick',round(logspace(log10(freqs2use(1)),log10(freqs2use(end)),8)))
title('PLI over trials')

subplot(222)
contourf(times2save,freqs2use,dwpli(:,times2saveidx),40,'linecolor','none')
set(gca,'clim',[-.2 .2],'yscale','log','ytick',round(logspace(log10(freqs2use(1)),log10(freqs2use(end)),8)))
title('dWPLI over trials')

subplot(223)
contourf(times2save,freqs2use,ispc_t,40,'linecolor','none')
set(gca,'clim',[-.1 .1],'yscale','log','ytick',round(logspace(log10(freqs2use(1)),log10(freqs2use(end)),8)))
title('ICPS over time')

subplot(224)
contourf(times2save,freqs2use,dwpli_t,40,'linecolor','none')
set(gca,'clim',[-.1 .1],'yscale','log','ytick',round(logspace(log10(freqs2use(1)),log10(freqs2use(end)),8)))
title('dWPLI over time')

%% Figure 26.11

trial2plot = 10; % any trial between 1 and 99 (book uses trial 10)
center_freq = 4.6; % Hz (book uses 4.6)

% create wavelet and take FFT
s = 4.5/(2*pi*center_freq);
wavelet_fft = fft( exp(2*1i*pi*center_freq.*time) .* exp(-time.^2./(2*(s^2))) ,n_convolution);
% phase angles from channel 1 via convolution
convolution_result_fft = ifft(wavelet_fft.*data_fft1,n_convolution);
convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
sig1 = reshape(convolution_result_fft,EEG.pnts,EEG.trials);
% phase angles from channel 2 via convolution
convolution_result_fft = ifft(wavelet_fft.*data_fft2,n_convolution);
convolution_result_fft = convolution_result_fft(half_wavelet+1:end-half_wavelet);
sig2 = reshape(convolution_result_fft,EEG.pnts,EEG.trials);
% cross-spectral density
xsd  = sig1 .* conj(sig2);
xsdi = imag(xsd);

dwpli = zeros(size(EEG.times));
ispc  = zeros(size(EEG.times));

[junk,animate_start] = min(abs(EEG.times-0));
[junk,animate_stop]  = min(abs(EEG.times-1000));

time_window_idx = round(100*timewindow(1)/(1000/EEG.srate));

figure
subplot(121)
hpol = polar(repmat(angle(xsd(animate_start:animate_start+time_window_idx-1,trial2plot)),1,2)',[zeros(time_window_idx,1) ones(time_window_idx,1)]','k-o');
subplot(122)
hplo2 = plot(EEG.times,0,'r');
hold on
hplo1 = plot(EEG.times,0,'b');

for idx=animate_start:animate_stop

    % update angles
    for i=1:length(hpol)
        set(hpol(i),'XData',[0; cos(angle(xsd(idx+i,trial2plot)))],'YData',[0; sin(angle(xsd(idx+i,trial2plot)))]);
    end
    title([ num2str(round(EEG.times(idx))) '-' num2str(round(EEG.times(idx+i))) ' ms' ])

    % compute ICPS and dwPLI
    ispc(idx) = abs(mean(exp(1i*angle(xsd(idx:idx+i,trial2plot))),1));

    imagsum        = sum(xsdi(idx:idx+i,trial2plot),1);
    imagsumW       = sum(abs(xsdi(idx:idx+i,trial2plot)),1);
    debiasfactor   = sum(xsdi(idx:idx+i,trial2plot).^2,1);
    dwpli(idx) = mean((imagsum.^2 - debiasfactor)./(imagsumW.^2 - debiasfactor));

    set(hplo1,'XData',EEG.times(1:idx),'YData',ispc(1:idx));
    set(hplo2,'XData',EEG.times(1:idx),'YData',dwpli(1:idx));
    set(gca,'xlim',EEG.times([animate_start animate_stop]),'ylim',[-.1 1.1])
    axis square

    pause(0.01)
end

%% Figure 26.12

% This figure is generated in the code below.

%% Figure 26.13

% generate inline functions (small functions that you define without being
% saved as general Matlab functions)
vtest  = inline('n.*icpcmag*cos(val).*sqrt(2./n)');
gvtest = inline('n.*(icpcmag*exp((-(val).^2)./(4.*pi./n)).*(sqrt(2./n)))');

% Since initially writing this code, Matlab decided to make the inline
% function obsolete in future versions. The following two lines produce
% the identical 'anonymous functions' as the previous two lines.
%vtest  = @(icpcmag,n,val) n.*icpcmag*cos(val).*sqrt(2./n);
%gvtest = @(icpcmag,n,val) n.*(icpcmag*exp((-(val).^2)./(4.*pi./n)).*(sqrt(2./n)));

% figure
clf
subplot(221)
n = 2:100;
plot(n,1-normcdf(vtest(.3,n,pi/10)))
hold on
plot(n,1-normcdf(gvtest(.3,n,pi/10)),'m')
legend({'v-test';'gv-test'})
xlabel('Number of points'), ylabel('P-value')
set(gca,'ylim',[0 .6])
title('angle = pi/10')

subplot(222)
plot(n,1-normcdf(vtest(.3,n,pi/3)))
hold on
plot(n,1-normcdf(gvtest(.3,n,pi/3)),'m')
legend({'v-test';'gv-test'})
xlabel('Number of points'), ylabel('P-value')
set(gca,'ylim',[0 .6])
title('angle = pi/3')

subplot(223)
x=linspace(-pi,pi,50);
n=15;
polar(x,1-normcdf(vtest(.3,n,x-0)))
hold on
polar(x,1-normcdf(gvtest(.3,n,x-0)),'m')
title([ 'N=' num2str(n) ])

subplot(224)
n=600;
polar(x,1-normcdf(vtest(.3,n,x-0)))
hold on
polar(x,1-normcdf(gvtest(.3,n,x-0)),'m')
set(gca,'xtick',round((-pi:pi/4:pi)*100)/100)
title([ 'N=' num2str(n) ])

% number of simulated data points
numUsims = 10000;

u = zeros(2,numUsims);

for i=1:numUsims

    % make some noise
    fake_phase_data = rand(2,EEG.pnts)*pi*2-pi;

    % compute ispc
    ispc_mag = abs  (mean(exp(1i*(diff(fake_phase_data,1)))));
    ispc_phs = angle(mean(exp(1i*(diff(fake_phase_data,1)))));

    % compute statistics
    u(1,i) = vtest (ispc_mag,EEG.pnts,ispc_phs-0);
    u(2,i) = gvtest(ispc_mag,EEG.pnts,ispc_phs-0);
end

% This figure is also figure 26.12 but with no log-scaling
figure
for i=1:2
    subplot(1,2,i)

    [y,x]=hist(u(i,:),100);
    h=bar(x,log10(y),'histc');
    set(h,'linestyle','none')

    title([ num2str(100*sum((1-normcdf(u(i,:)))<.05)/length(u)) '% false positive' ])
end

nrange    = 10:300; % here n is number of datapoints
ispcrange = .05:.01:.7;
pvalmat   = zeros(2,length(nrange),length(ispcrange));

for ni=1:length(nrange)
    for mi=1:length(ispcrange)

        n = nrange(ni);

        pvalmat(1,ni,mi) = 1-normcdf( vtest(ispcrange(mi),nrange(ni),pi/5));
        pvalmat(2,ni,mi) = 1-normcdf(gvtest(ispcrange(mi),nrange(ni),pi/5));
    end
end

figure
subplot(121)
imagesc(ispcrange,nrange,squeeze(pvalmat(1,:,:))), axis xy,
set(gca,'clim',[0 .5])
xlabel('ICPS strength'), ylabel('N')
title('v-test')

subplot(122)
imagesc(ispcrange,nrange,squeeze(pvalmat(2,:,:))), axis xy,
title('gv-test')
xlabel('ICPS strength'), ylabel('N')
set(gca,'clim',[0 .5])

%% end.

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


% some variables needed
u = P.x_trim(4);
v = P.x_trim(5);
w = P.x_trim(6);
Va_trim = sqrt(u^2 + v^2 + w^2);
theta_trim = P.x_trim(8);
alpha_trim = atan(w/u);
delta_e_trim = P.u_trim(1);
delta_t_trim = P.u_trim(4);

% Roll attitude hold
P.delta_a_max = 45 * pi / 180;
P.phi_max = 15 * pi / 180;
P.roll_max = 45 * pi / 180;
a_phi1 = -0.5 * P.rho * Va_trim^2 * P.S_wing * P.b * P.C_p_p * P.b / (2 * Va_trim);
a_phi2 = 0.5 * P.rho * Va_trim^2 * P.S_wing * P.b * P.C_p_delta_a;
P.kp_phi = P.delta_a_max / P.phi_max * sign(a_phi2);
omega_phi = sqrt(abs(a_phi2) * P.delta_a_max / P.phi_max);
zeta_phi = 1.2;   % design parameter
P.kd_phi = (2 * zeta_phi * omega_phi - a_phi1) / a_phi2;
P.ki_phi = 0.2;    % design parameter

% Course hold
W_chi = 8;  % design parameter
omega_chi = 1 / W_chi * omega_phi;
zeta_chi = 0.707;   % design parameter
Vg = P.Va0;
P.kp_chi = 2 * zeta_chi * omega_chi * Vg / P.gravity;
P.ki_chi = omega_chi^2 * Vg / P.gravity;

% Pitch attitude hold
P.delta_e_max = 45 * pi / 180;
P.e_theta_max = 10 * pi / 180;
a_theta1 = -P.rho * Va_trim^2 * P.c * P.S_wing * P.C_m_q * P.c / (2 * P.Jy * 2 * Va_trim);
a_theta2 = -P.rho * Va_trim^2 * P.c * P.S_wing * P.C_m_alpha / (2 * P.Jy);
a_theta3 = P.rho * Va_trim^2 * P.c * P.S_wing * P.C_m_delta_e / (2 * P.Jy);
P.kp_theta = P.delta_e_max / P.e_theta_max * sign(a_theta3);
omega_theta = sqrt(a_theta2 + P.delta_e_max / P.e_theta_max * abs(a_theta3));
zeta_theta = 0.707;     % design parameter
P.kd_theta = (2 * zeta_theta * omega_theta - a_theta1) / a_theta3;
P.theta_max = 20 * pi / 180;
P.pitch_max = 40 * pi / 180;

% Airspeed hold using Throttle
a_V1 = P.rho * Va_trim * P.S_wing / P.mass ...
       * (P.C_D_0 + P.C_D_alpha * alpha_trim + P.C_L_delta_e * delta_e_trim) ...
       + P.rho * P.S_prop / P.mass * P.C_prop * Va_trim;
a_V2 = P.rho * P.S_prop / P.mass * P.C_prop * P.k_motor^2 * delta_t_trim;
a_V3 = P.gravity * cos(theta_trim - alpha_trim);
omega_v = 5;   % design parameter
zeta_v = 0.707;     % design parameter
P.delta_t_max = 1;
P.delta_t_min = 0;
P.ki_v = omega_v^2 / a_V2;
P.kp_v = (2 * zeta_v * omega_v - a_V1) / a_V2;

% Airspeed hold using Pitch
W_v2 = 7;   % design parameter
zeta_v2 = 0.707;    % design parameter
omega_v2 = 1 / W_v2 * omega_theta;
K_theta_dc = P.kp_theta * a_theta3 / (a_theta2 + P.kp_theta * a_theta3);
P.ki_v2 = -omega_v2^2/(K_theta_dc * P.gravity);
P.kp_v2 = (a_V1 - 2 * zeta_v2 * omega_v2) / (K_theta_dc * P.gravity);

% Altitude hold using Pitch
W_h = 10;   % design parameter
omega_h = 1 / W_h * omega_theta;
Va = P.Va0;
zeta_h = 1.2;     % design parameter
P.h_max = 1000;
P.h_min = 0;
P.ki_h = omega_h^2 / (K_theta_dc * Va);
P.kp_h = 2 * zeta_h * omega_h / (K_theta_dc * Va);


text:
%% The SCNLOBJECT cookbook
% This cookbook is designed to familiarize yourself with the scnlobjects used within waveform. There's not really all that much to them, so there'll not be much to this, either.

%% Properties
% Scnlobjects are comprised of only four fields:
%%
%
% * Station
% * Channel
% * Network
% * Location

%%
% STATION and CHANNEL are not optional, while NETWORK and LOCATION are.

%%
% Create a series of scnls for seismic stations on Redoubt volcano
    red = scnlobject('RED','EHZ','AV','--') % full declaration
    rso = scnlobject('RSO','EHZ','AV')      % leaving out location

%%
% 3-component station, this creates rdbw as a 1x3 scnlobject
    rdbw = scnlobject('RDBW',{'BHZ','BHE','BHN'},'AV')

%% Using SET to modify a SCNLOBJECT
% use SET to create another 3 component broadband using rdbw as a template...
    rdjh = set(rdbw,'station','RDJH');

%% Grouping SCNLOBJECTS through concatenation% and finally, create something that will parse through all stations

    redstas = [red rso rdbw rdjh]
%%
% Create a couple scnls with wildcard values
% The value '*' can be entered in place of a field, so that searches can be made.

    anyehz = scnlobject('*','EHZ','*','*') % this one prints as an example...
    anybroadband = scnlobject('*',{'BHZ','BHE','BHN'},'*','*');
    anyNorthComponent = scnlobject('*',{'EHN','SHN','BHN'},'*','*');

%%
% Let's display all the variables we've got, thus far

    whos

%% Methods
% Find out what manipulations can be done with scnlobjects

    methods(redstas)

%%
% I could have typed "methods(scnlobject)", instead

%% using ISMEMBER to see if a particular locale (SCNLOBJECT) is represented
% Ismember will not only tell you if a scnlobject exists, but is also capable of telling WHERE

    [IsInList, whereInList] = ismember(rso,redstas)

%%
% ismember is sensitive to the order in which you place the scnlobjects. The size of the returned values match the size of the first argument. The previous example asks '' where can RSO be found in REDSTAS? '', while the following example essentially asks '' which REDSTAS can be found in RSO (and where)?''

    [IsInList, whereInList] = ismember(redstas,rso)

%% using ISMEMBER with wildcards
% The ismember function looks for a scnlobject in an array of scnlobjects. Lets grab the EHZ component stations using a search, then display each one

    stationsOfInterest = redstas(ismember(redstas,anyehz)) %should get 2

    for n = 1 : numel(stationsOfInterest);
        display(stationsOfInterest(n));
    end

%%
% repeat the exercise for broadband and stations... there should be 6

    mybroadbands = redstas(ismember(redstas,anybroadband));
    size(mybroadbands)

%% Use GET to see which stations we have...
% Get can also retrieve CHANNEL, NETWORK, and LOCATION information.

    get(mybroadbands,'station')

%%
% There are 3 of each because there are three components to each.

    strcat(get(mybroadbands,'station'),'|',get(mybroadbands,'channel'))

%%
% to get a vertical list, I could transpose mybroadbands using (')

    strcat( get(mybroadbands', 'station') , '|' , get(mybroadbands', 'channel') )

%% Using UNIQUE to whittle down a bunch of scnlobjects
% Assume that for some reason, you have a list of scnlobjects with some members being repeated. This may happen when you retrieve scnlobjects from a group of waveforms. You'd like to know which station/channels are actually being represented within your data set. Here's how:
%
% set up the situation as listed:
    manyscnls = [redstas, red, rso, red, redstas]
%%
% Find out which ones we have:
    unique(manyscnls)

author: geoscience-community-codes

repo: GISMO


text:
function [mHitRate , mFalseAlarm, AUC] = performBenchmark(datasetStruct,algStructArray)

evaluateSal = @(sMap,thresholds,gtMap) thresholdBased_HR_FR(sMap,thresholds,gtMap);
switch lower(datasetStruct.datasetName)
    case 'msra'
        evaluateSal = @(sMap,thresholds,gtMap) MSRAthresholdBased_HR_FR(sMap,thresholds,gtMap);
    case 'asd'
    case 'sed1'
    case 'sed2'
    case 'sod'
        %         Do nothing
    otherwise
        fprintf('\nPerforming analysis on unknown dataset (%s)\nMake sure that the ground-truth maps are handled correctly\n',datasetStruct.datasetName);
end

fprintf(['\nEvaluating dataset: ' datasetStruct.datasetName '\n']);
thresholds = 1:-.05:0;
GTfiles=dir([datasetStruct.GTdir '/*.png']);
GTfiles = [GTfiles; dir([datasetStruct.GTdir '/*.jpg'])];
GTfiles = [GTfiles; dir([datasetStruct.GTdir '/*.bmp'])];

numOfFiles = size(GTfiles,1);
numOfAlgos = length(algStructArray);

[hitRate, falseAlarm] = deal(zeros(numOfFiles,length(thresholds),numOfAlgos));
%Iterate over images
for imIndx=1:numOfFiles

%     if imIndx == 3
%         break;
%     end

    fprintf('Processing image %i out of %i\n',imIndx,numOfFiles);
    [~,base_name,ext] = fileparts(GTfiles(imIndx).name);

    gtMap = im2double(imread([datasetStruct.GTdir base_name ext]));
    gtSize = size(gtMap);
    if (length(gtSize) == 3)
        gtMap = rgb2gray(gtMap);
        gtSize(3)= [];
    end
    gtMap = double(gtMap>=0.1);
    for algIdx = 1:numOfAlgos
        sMap = readSaliencyMap(algStructArray(algIdx),base_name,gtSize);
        [hitRate(imIndx,:,algIdx), falseAlarm(imIndx,:,algIdx)] = evaluateSal(sMap,thresholds,gtMap);
    end

end %End of image loop

%Average across images -
mHitRate = permute(mean(hitRate,1),[2 3 1]);
mFalseAlarm = permute(mean(falseAlarm,1),[2 3 1]);

AUC = nan(1,size(mFalseAlarm,2));
for algIdx=1:numOfAlgos
        AUC(algIdx) = trapz(mFalseAlarm(:,algIdx),mHitRate(:,algIdx));
end

end

% Read and resize saliency map
function sMap = readSaliencyMap(algStruct,base_name,gtSize)
file_name = fullfile(algStruct.dir,[algStruct.prefix base_name algStruct.postfix '.' algStruct.ext]);
sMap = imresize(im2double(imread(file_name)),gtSize(1:2));
if (size(sMap,3)==3)
    sMap = rgb2gray(sMap);
end
sMap = sMap./max(sMap(:));
sMap(sMap<0)=0;
end

function [hitRate, falseAlarm] = thresholdBased_HR_FR(sMap,thresholds,gtMap)
numOfThreshs = length(thresholds);
[hitRate, falseAlarm] = deal(zeros(1,numOfThreshs));
for threshIdx=1:numOfThreshs
    cThrsh=thresholds(threshIdx);
    [hitRate(threshIdx) , falseAlarm(threshIdx)] = hitRates((sMap>=cThrsh),gtMap);
end
end

function [hitRate, falseAlarm] = MSRAthresholdBased_HR_FR(sMap,thresholds,gtMap)
numOfThreshs = length(thresholds);
[hitRate, falseAlarm] = deal(zeros(1,numOfThreshs));

for threshIdx=1:numOfThreshs
    cThrsh=thresholds(threshIdx);
    STATS = regionprops(uint8(sMap>=cThrsh),'BoundingBox');
    bMap = zeros(size(sMap));
    bbox = round([STATS.BoundingBox]);
    if (numel(bbox)>0)
        bMap(bbox(2):(bbox(2)+bbox(4)-1),bbox(1):(bbox(1)+bbox(3)-1))=1;
    end
    [hitRate(threshIdx) , falseAlarm(threshIdx)] = hitRates(bMap,gtMap);
end
end

author: ArcherFMY


text:
function parameters = initialise_parameters(ns_ps,Srate,method)
len_val = length(ns_ps);
switch lower(method)
    case 'martin'
        L_val=len_val;
        R_val=len_val/2;
        D_val=150; V_val=15; Um_val=10; Av_val=2.12;
        alpha_max_val=0.96;
        alpha_min_val=0.3;
        beta_max_val=0.8;
        x_val=[1 2 5 8 10 15 20 30 40 60 80 120 140 160];
        Y_M_val=[0 .26 .48 .58 .61 .668 .705 .762 .8 .841 .865 .89 .9 .91];
        Y_H_val=[0 .15 .48 .78 .98 1.55 2.0 2.3 2.52 2.9 3.25 4.0 4.1 4.1];
        xi_val=D_val;
        M_D_val=interp1(x_val,Y_M_val,xi_val);
        H_D_val=interp1(x_val,Y_H_val,xi_val);
        xi_val=V_val;
        M_V_val=interp1(x_val,Y_M_val,xi_val);
        H_V_val=interp1(x_val,Y_H_val,xi_val);
        minact_val(1:L_val,1:Um_val)=max(ns_ps);
        parameters = struct('n',2,'len',len_val,'alpha_corr',0.96,'alpha',0.96*ones(len_val,1),'P',ns_ps,'noise_ps',ns_ps,'Pbar',ns_ps,...
            'Psqbar',ns_ps,'actmin',ns_ps,'actmin_sub',ns_ps,'Pmin_u',ns_ps,'subwc',2,'u',1,'minact',minact_val,'lmin_flag',zeros(len_val,1),...
            'L',L_val,'R',R_val,'D',D_val,'V',V_val,'Um',Um_val,'Av',Av_val,'alpha_max',alpha_max_val,'alpha_min',alpha_min_val,...
            'beta_max',beta_max_val,'Y_M',Y_M_val,'Y_H',Y_H_val,'M_D',M_D_val,'H_D',H_D_val,'M_V',M_V_val,'H_V',H_V_val);
    case 'mcra'
        parameters = struct('n',2,'len',len_val,'P',ns_ps,'Pmin',ns_ps,'Ptmp',ns_ps,'pk',zeros(len_val,1),'noise_ps',ns_ps,...
            'ad',0.95,'as',0.8,'L',round(1000*2/20),'delta',5,'ap',0.2);
    case 'imcra'
        alpha_d_val=0.85;
        alpha_s_val=0.9;
        U_val=8;V_val=15;
        Bmin_val=1.66;gamma0_val=4.6;gamma1_val=3;
        psi0_val=1.67;alpha_val=0.92;beta_val=1.47;
        j_val=0;
        b_val=hanning(3);
        B_val=sum(b_val);
        b_val=b_val/B_val;
        Sf_val=zeros(len_val,1);Sf_tild_val=zeros(len_val,1);
        Sf_val(1) = ns_ps(1);
        for f=2:len_val-1
            Sf_val(f)=sum(b_val.*[ns_ps(f-1);ns_ps(f);ns_ps(f+1)]);
        end
        Sf_val(len_val)=ns_ps(len_val);
        Sf_tild_val = zeros(len_val,1);
        parameters = struct('n',2,'len',len_val,'noise_ps',ns_ps,'noise_tild',ns_ps,'gamma',ones(len_val,1),'Sf',Sf_val,...
            'Smin',Sf_val,'S',Sf_val,'S_tild',Sf_val,'GH1',ones(len_val,1),'Smin_tild',Sf_val,'Smin_sw',Sf_val,'Smin_sw_tild',Sf_val,...
            'stored_min',max(ns_ps)*ones(len_val,U_val),'stored_min_tild',max(ns_ps)*ones(len_val,U_val),'u1',1,'u2',1,'j',2,...
            'alpha_d',0.85,'alpha_s',0.9,'U',8,'V',15,'Bmin',1.66,'gamma0',4.6,'gamma1',3,'psi0',1.67,'alpha',0.92,'beta',1.47,...
            'b',b_val,'Sf_tild',Sf_tild_val);
    case 'doblinger'
        parameters = struct('n',2,'len',len_val,'alpha',0.7,'beta',0.96,'gamma',0.998,'noise_ps',ns_ps,'pxk_old',ns_ps,...
            'pxk',ns_ps,'pnk_old',ns_ps,'pnk',ns_ps);
    case 'hirsch'
        parameters = struct('n',2,'len',len_val,'as',0.85,'beta',1.5,'omin',1.5,'noise_ps',ns_ps,'P',ns_ps);
    case 'mcra2'
        freq_res=Srate/len_val;
        k_1khz=floor(1000/freq_res);
        k_3khz=floor(3000/freq_res);
        %delta_val=[2*ones(k_1khz,1);2*ones(k_3khz-k_1khz,1);5*ones(len_val/2-k_3khz,1);...
        %    5*ones(len_val/2-k_3khz,1);2*ones(k_3khz-k_1khz,1);2*ones(k_1khz,1)];
         delta_val=[2*ones(k_1khz,1);2*ones(k_3khz-k_1khz,1);5*ones(len_val/2-k_3khz,1)];
			delta_val=[delta_val;5;flipud(delta_val(2:end))];

        parameters = struct('n',2,'len',len_val,'ad',0.95,'as',0.8,'ap',0.2,'beta',0.8,'beta1',0.98,'gamma',0.998,'alpha',0.7,...
            'delta',delta_val,'pk',zeros(len_val,1),'noise_ps',ns_ps,'pxk_old',ns_ps,'pxk',ns_ps,'pnk_old',ns_ps,'pnk',ns_ps);

      case 'conn_freq'
        D = 7;
        b = triang(2*D+1)/sum(triang(2*D+1));
        b = b';
        beta_min = 0.7; % for R's recursion
        U = 5;
        V = 8;
        gamma1 = 6;
        gamma2 = 0.5;
        K_tild = 2*sum(b.^2)^2/sum(b.^4);
        alpha_max_val=0.96;
        alpha_min_val=0.3;
        stored_min = max(ns_ps)*ones(len_val,U);

        alpha_c = 0.7;
        noise_ps = ns_ps;
        Rmin_old = 1;
        Pmin_sw = ns_ps;
        Pmin = ns_ps;
        P = ns_ps;
        Decision = zeros(size(P));
        u1 = 1;
        j = 0;
        parameters = struct('len',len_val,'D',D,'b',b,'U',U,'V',V,'gamma1',gamma1,'gamma2',gamma2,'K_tild',K_tild,'alpha_c',alpha_c,...
            'noise_ps',noise_ps,'Rmin_old',Rmin_old,'Pmin_sw',Pmin_sw,'Pmin',Pmin,'SmthdP',P,'u1',u1,'j',j,'alpha',0,'alpha_max',alpha_max_val,...
            'stored_min',stored_min,'beta_min',beta_min,'Decision',Decision);

    otherwise
            error('Method not implemented. Check spelling.');
end

author: jtkim-kaist

repo: Speech-enhancement


text:
function [TrainErr, TestErr, feature_training,feature_testing] = map_FNPAQR(feature_training,class_training,feature_testing,class_testing)
class_training=grp2idx(class_training);
class_testing=grp2idx(class_testing);
[feature_training,ps] = mapminmax(feature_training',0,1);feature_training=feature_training';
feature_testing = mapminmax('apply',feature_testing',ps)';
%%
options.ReducedDim = max(class_training)-1;
[eigvector]=FNPAQR(feature_training, class_training,options);
feature_training = feature_training * eigvector;
feature_testing = feature_testing * eigvector;
%%
TrainPredict = classify(feature_training,feature_training,class_training,'quadratic');
TestPredict = classify(feature_testing,feature_training,class_training,'quadratic');
TrainErr = sum(TrainPredict ~= class_training)/length(class_training)*100;
TestErr = sum(TestPredict ~= class_testing)/length(class_testing)*100;

author: Sable

repo: mcbench-benchmarks


text:
function variance = hypergeometric_variance ( n, m, l )

%*****************************************************************************80
%
%% HYPERGEOMETRIC_VARIANCE returns the variance of the Hypergeometric PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of balls selected.
%    0 <= N <= L.
%
%    Input, integer M, the number of white balls in the population.
%    0 <= M <= L.
%
%    Input, integer L, the number of balls to select from.
%    0 <= L.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = n * m * ( l - m ) * ( l - n ) / ( l * l * ( l - 1 ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% correlate frames with each other
% sort frames by correlations
% take mean image as mean over most correlated frames

function mimg = pick_reg_init(data)
%%
dd = bsxfun(@minus, data, mean(mean(data,1),2));

% WHITENING???
% for i = 1:size(data,3)
%    d0 = dd(:,:,i);
%    fd0 = fft2(d0);
%    dd(:,:,i) = real(ifft2(fd0./abs(fd0)));
% end

dd = reshape(dd, [], size(dd,ndims(dd)));

CC = dd'*dd;
CC = CC./(diag(CC) * diag(CC)').^.5;

% CC = corrcoef(dd);

[CCsort, isort] = sort(CC, 2, 'descend');
nc = size(CCsort,2);
bestCC = mean(CCsort(:, 1:min(nc,20)), 2);
[~, imax] = max(bestCC);

mimg = mean(data(:,:,isort(imax, 1:min(nc,20))), 3);
end

author: cortex-lab


text:
function a = i4col_sort2_a ( m, n, a )

%*****************************************************************************80
%
%% I4COL_SORT2_A ascending sorts the elements of each column of an I4COL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of A.
%
%    Input, integer N, the number of columns of A, and the length
%    of a vector of data.
%
%    Input, integer A(M,N), the array of N columns of M vectors.
%
%    Output, integer A(M,N), the elements of each column of A have been
%    sorted in ascending order.
%
  if ( m <= 1 )
    return
  end

  if ( n <= 0 )
    return
  end
%
%  Initialize.
%
  for col = 1 : n

    i = 0;
    indx = 0;
    isgn = 0;
    j = 0;
%
%  Call the external heap sorter.
%
    while ( 1 )

      [ indx, i, j ] = sort_heap_external ( m, indx, isgn );
%
%  Interchange the I and J objects.
%
      if ( 0 < indx )

        t        = a(i,col);
        a(i,col) = a(j,col);
        a(j,col) = t;
%
%  Compare the I and J objects.
%
      elseif ( indx < 0 )

        if ( a(j,col) < a(i,col) )
          isgn = +1;
        else
          isgn = -1;
        end

      elseif ( indx == 0 )

        break

      end

    end

  end

  return
end


text:
function st_print_some ( row1, row2, col1, col2, nst, ist, jst, ast, title )

%*****************************************************************************80
%
%% ST_PRINT_SOME prints some of an ST file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ROW1, ROW2, the first and last rows to print.
%
%    Input, integer COL1, COL2, the first and last columns to print.
%
%    Input, integer NST, the number of nonzeros.
%
%    Input, integer IST(NST), COL(NST), the row and column indices.
%
%    Input, real A(NST), the nonzero values.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  for k = 1 : nst
    if ( row1 <= ist(k) & ist(k) <= row2 & col1 <= jst(k) & jst(k) <= col2 )
      fprintf ( 1, '  %8d  %8d  %8d  %16.8f\n', k, ist(k), jst(k), ast(k) );
    end
  end

  return
end


text:
%Nice grid by Kris Donhowe, SC Solutions, Inc.
%
h_axes = findobj(get(gcf,'Children'),'Type','axes');
grey = [0.7,0.7,0.7];
set(h_axes,'xcolor',grey,'ycolor',grey, ...
    'GridLineStyle','-','MinorGridLineStyle','-', ...
    'Units','pixels');
grid on
c11=copyobj(h_axes,gcf);
set(c11,'color','none','xcolor','k', ...
    'xgrid','off','ycolor','k', ...
    'ygrid','off');


text:
function [output_canvas] = image_blending_linear(warped_img1,warped_img2)

    c1omask{1} = imfill(im2bw(uint8(warped_img1), 0),'holes');
    c1omask{2} = imfill(im2bw(uint8(warped_img2), 0),'holes');

    c1out{1} = double(warped_img1);
    c1out{2} = double(warped_img2);

    for i = 1 : 2
        if i == 1
            out = c1out{1};
            out_mask = c1omask{1};
            % center of out
            [r, c] = find(c1omask{1}(:, :, 1));
            out_center = [mean(r) mean(c)];
        else % blend out and c1out{i}
            % center of c1out{i}
            [r, c] = find(c1omask{i}(:, :, 1));
            out_i_center = [mean(r) mean(c)];
            % compute weighted mask
            vec = out_i_center - out_center; % vector from out_center to out_i_center
            intsct_mask = c1omask{i}(:, :, 1) & out_mask(:, :, 1); % 1 channel

            [r, c] = find(intsct_mask(:, :, 1));
            idx = sub2ind(size(c1omask{i}(:, :, 1)), r, c);
            out_wmask = zeros(size(c1omask{i}(:, :, 1)));
            proj_val = (r - out_center(1))*vec(1) + (c- out_center(2))*vec(2); % inner product
            out_wmask(idx) = (proj_val - (min(proj_val)+(1e-3))) / ...
                             ((max(proj_val)-(1e-3)) - (min(proj_val)+(1e-3))); % weight map (of overlapped area) for c1out{i}, 1 channel
            % blending
            mask1 = out_mask(:, :, 1)&(out_wmask==0);
            mask2 = out_wmask;
            mask3 = c1omask{i}(:, :, 1)&(out_wmask==0);
            mask1 = cat(3, mask1, mask1, mask1); mask2 = cat(3, mask2, mask2, mask2); mask3 = cat(3, mask3, mask3, mask3);
            out = out.*(mask1+(1-mask2).*(mask2~=0)) + c1out{i}.*(mask2+mask3);
            % update
            out_mask = out_mask | c1omask{i};
            out_center = out_i_center; % update out_center by assign center of c1out{i}
        end
    end

    output_canvas = uint8(out);
end

author: YaqiLYU

repo: AANAP

save_path: github-repos/MATLAB/YaqiLYU-AANAP


text:
%  INTERNAL FUNCTION: Kalman filter for non-regime-switching models
%
%  ::
%
%    f=missing_observations_kalman_filter(y,T,R,Z,H,Q,init,dy,ca,level)
%
%  Args:
%
%       y (p x n matrix): data
%       T (m x m x nT): state matrix, see below
%       R (m x r x nR): state matrix, see below
%       Z (p x m | logical): state matrix or selector, see below
%       H (p x p x nH): Covariances of measurement errors
%       Q (r x r x nQ): Covariances of shocks
%       init (struct):
%       dy (p x ndy):
%       ca (m x nca):
%       level (0 | 1 | 2 | {3}): controls the level output produced.
%
%           - if level==0, only the likelihood is returned
%           - if level==1, filters/forecasts are returned in addition to level 0
%           - if level==2, updates are returned in addition to level 1
%           - if level==3, smoothed values are returned in addition to level 2
%
%       nstep (integer | {1}): Number of forecast steps
%
%  Returns:
%     :
%
%       - **f** [struct]: output structure
%
%           - **retcode** [0|305]: 305 if computation of the likelihood breaks
%             downs
%           - **incr** [vector]: likelihood in each period
%           - **log_lik** [scalar]: log likelihood
%           - **a** [m x n+1 x nstep]: filtered stated (available if level >0)
%           - **P** [m x m x n+1 x nstep]: Covariance matrix of filtered state
%             (available if level >0)
%           - **att** [m x n]: Updated state vector  (available if level >1)
%           - **Ptt** [m x m x n]: Covariance matrix for updates (available if
%             level >1)
%           - **iF** [p x p x n]: Inverses of covariance matrix of forecast
%             errors
%           - **v** [p x n]: Forecast errors
%           - **K** [m x p x n]: Kalman gains
%           - **alpha** [m x n]: smoothed state vector
%           - **V** [m x m x n]: Covariance matrix of smoothed state vector
%           - **r** [m x n+1]: quantity useful for the efficient computation of
%             the state vector
%           - **eta** [r x n]: smoothed shocks
%           - **epsilon** [p x n]: smoothed measurement errors
%
%  Note:
%
%     The system is of the form
%     .. math::
%
%        alpha(t+1) = T*alpha(t) + ca(t) + R(t)*eta(t),  eta ~ N(0, Q)
%        y(t) =   Z*alpha(t) + dy(t) + epsilon(t),  epsilon ~ N(0,H)
%
%

author: jmaih