text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  battistn@tcnj.edu
% Date Modified: April 2021
% Current Institution: TCNJ
% Date Created: May 27th, 2015
% Institution Created: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (torsional springs or non-invariant beams)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (battistn@tcnj.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the RUBBERBAND-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Rubberband(Nx,Ny)

%------------------------------------------------------------
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%------------------------------------------------------------
%Nx =  128;      % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
%Ny =  128;      % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
Ly = 1.0;        % Length of Eulerian Grid in y-Direction
dx=Lx/Nx;

%---------------------------------------------------------------
% Necessary grid info for ** Concentration ** input file
%  (doesn't need to be changed unless rectangular grid is used)
%---------------------------------------------------------------
x=0:dx:Lx-dx;
y=x;

%------------------------------------------------------------
% Immersed Structure Geometric / Dynamic Parameters %
%------------------------------------------------------------
a = 0.4;         % Length of semi-major axis.
b = 0.2;         % Length of semi-minor axis.
ds_Rest = 0;     % Resting length of springs

%------------------------------------------------------------
% Rubberband overall perimeter length and Lagrangian spacing
%------------------------------------------------------------
L=1.93768964;
ds=Lx/(2*Nx);
s = 0:ds:L;
N=length(s);

%------------------------------------------------------------
% Structure name
%------------------------------------------------------------
struct_name = 'rubberband'; % Name for .vertex, .spring, etc files.

%------------------------------------------------------------
% Call function to construct geometry
%------------------------------------------------------------
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(N,a,b);

%------------------------------------------------------------
% Plot Geometry to test
%------------------------------------------------------------
plot(xLag,yLag,'r-'); hold on;
plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis square;

%------------------------------------------------------------
% Prints .vertex file!
%------------------------------------------------------------
print_Lagrangian_Vertices(xLag,yLag,struct_name);

%------------------------------------------------
% Prints .spring file!
%       ->> find k_Spring for certain resolution
%------------------------------------------------
ds128 = 0.5*Lx/128;   % Lag. Spacing in 128x128 Resolution
F_Spr_128 = 4e5;      % Force used in 128x128 Resolution
%
% Get Spring Stiffness Coefficient void of any particular resolution
k_Spring = F_Spr_128 * ds128^3 / ds128;
%
% Get Spring Force based on particular resolution
F_Spring = k_Spring * ds / ds^3;
%
print_Lagrangian_Springs(xLag,yLag,F_Spring,ds_Rest,struct_name);

%------------------------------------------------------------
% Call function for initial concentration
%------------------------------------------------------------
[Concentration,X,Y]= give_Me_Initial_Concentration(Lx,Ly,Nx,Ny,dx,x,y);

%------------------------------------------------------------
% Prints .concentration file!
%------------------------------------------------------------
print_Concentration_Info(Nx,Nx,Concentration,struct_name);

%------------------------------------------------------------
% Prints .geo_connect file!
%------------------------------------------------------------
print_Geometry_Connections(xLag,struct_name);

% Prints .beam file!
%k_Beam = 0.5; C = 0.0;
%print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .target file!
%k_Target = 1e7;
%print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints GEOMETRY CONNECTIONS to a file called
%                               'struct_name'.geo_connect
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Geometry_Connections(xLag,struct_name)

    N = length(xLag); % Total # of Lag. Pts

    geo_fid = fopen([struct_name '.geo_connect'], 'w');

    % Loops over all Lagrangian Pts.
    for i = 1:N

        if i<N

            s1 = i;
            s2 = i+1;
            fprintf(geo_fid, '%d %d\n', s1, s2);
            fprintf(geo_fid, '%d %d\n', s2, s1);

        elseif i==N

            s1 = N;
            s2 = 1;
            fprintf(geo_fid, '%d %d\n', s1, s2);
            fprintf(geo_fid, '%d %d\n', s2, s1);

        end

    end

    fclose(geo_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called rubberband.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called rubberband.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            else
                %Case s=N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints initial CONCENTRATION to rubberband.concentration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Concentration_Info(Nx,Ny,C,struct_name)

    con_fid = fopen([struct_name '.concentration'], 'w');

    for i=1:Ny
        for j=1:Nx
            fprintf(con_fid, '%1.16e ', C(i,j) );
        end
        fprintf(con_fid,'\n');
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(N,rmin,rmax)

% The immsersed structure is an ellipse %
for i=1:N

    xLag(i) = 0.5 + rmax * cos( 2*pi/N*(i-1) );
    yLag(i) = 0.5 + rmin * sin( 2*pi/N*(i-1) );

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: gives initial concentration in Eulerian domain
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [C,X,Y] = give_Me_Initial_Concentration(Lx,Ly,Nx,Ny,dx,x,y)

[X,Y]=meshgrid(x,y);

C=zeros(Nx,Nx);

author: nickabattista


text:
function v = fibreVolume(odf,h,r,radius,varargin)
% ratio of orientations with a certain orientation
%
% Description
% returns the ratio of mass of the odf that is within a certain
% distance from a given fibre
%
% Syntax
%   v = fibreVolume(odf,h,r,radius)
%
% Input
%  odf    - @SO3Fun
%  h      - @Miller
%  r      - @vector3d
%  radius - double
%
% Options
%  resolution - resolution of discretization
%
% See also
% SO3Fun/volume SO3Fun/entropy SO3Fun/textureindex

% check input
argin_check(h,{'Miller','vector3d'});
if isa(h,'Miller'), h = odf.CS.ensureCS(h);end
argin_check(r,'vector3d');
argin_check(radius,'double');

% get resolution
res = get_option(varargin,'RESOLUTION',min(2.5*degree,radius/50),'double');

% discretisation
sR = odf.CS.fundamentalSector;
S2G = equispacedS2Grid(sR,'resolution',res,varargin{:});
lS2G = length(S2G);
S2G = S2G(angle(h,Miller(S2G,odf.CS))<radius);

% estimate volume portion of odf space
f = length(S2G)/lS2G;

% eval odf
if f==0
  v = 0;
else
  v = min(1,mean(odf.radon(S2G,r)) * f);
end

author: mtex-toolbox

repo: mtex

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


text:
lw=2;
fz=20;

data_fb_real =load('much_longer_fb_real.log');
data_fac_real=load('much_longer_fac_real.log');
data_ola_real=load('much_longer_ola_real.log');

% Columns in data, fb : a M L W gl time
% Columns in data, fac: a M L W time
% Columns in data, ola: a M L W gl bl time

Ls=data_fac_real(:,3);
t_fb_real =data_fb_real(:,6);
t_fac_real=data_fac_real(:,5);
t_ola_real=data_ola_real(:,7);

if 0
  % Color legend
  l1='b';
  l2='b--';
  l3='r';
  l4='r--';
else
  % bw legend
  l1='b';
  l2='b--';
  l3='b-.';
  l4='b:';
end;

figure(1);

M=data_fac_real(1,2);
L=data_fac_real(:,3);

plot(Ls,t_fb_real,l1,...
     Ls,t_fac_real,l2,...
     Ls,t_ola_real,l3,'LineWidth',lw);
set(gca,'Fontsize',fz);

h=legend('Portnoff, real','Fac, real','Fac-OLA, real',...
       'Location','NorthWest');

% Grow the box a little, otherwise the export to .eps is messed up.
q=get(h,'Position');
set(h,'Position',[q(1)*.9 q(2)*.95 q(3)*1.8 q(4)]);

xlabel('Signal length / samples','Fontsize',fz);
ylabel('Running time / seconds','Fontsize',fz);

print -deps plot_much_longer_1.eps

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function Z = plus(X,Y)
% PLUS (overloaded)

dim = [];
if isa(X,'ndsdpvar')
    if isa(Y,'double') && numel(Y)==1
        % Quick return for scalar  + ndsdpvar
        Z = X;Z.basis(:,1) = Z.basis(:,1) + Y;
        Z.conicinfo = [0 0];
        Z.extra.opname='';
        return
    else
        dim = X.dim;
        X = sdpvar(X);
    end
elseif isa(X,'double')
    X = X(:);
end

if isa(Y,'ndsdpvar')
    if isa(X,'double') && numel(X)==1
        % Quick return for scalar  + ndsdpvar
        Z = Y;Z.basis(:,1) = Z.basis(:,1) + X;
        Z.conicinfo = [0 0];
        Z.extra.opname='';
        return
    else
        if isempty(dim)
            dim = Y.dim;
        else
            if isequal(dim,Y.dim)
            else
                error('Dimension mismatch in nD addition')
            end
        end
        Y = sdpvar(Y);
    end
elseif isa(Y,'double')
    Y = Y(:);
end

Z = X + Y;
Z = reshape(Z,dim);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef tdropout < matlab.unittest.TestCase
    % tdropout   Unit tests for transformer.layer.dropout

    % Copyright 2020 The MathWorks, Inc.
    properties(TestParameter)
        Size = struct(...
            'Vector',[10000,1],...
            'Matrix',[100,100],...
            'Tensor',[5,60,70])
        Probability = struct(...
            'Zero',0,...
            'Half',0.5,...
            'OneIsh',0.95)
    end

    methods(TestMethodSetup)
        function resetGlobalSeed(~)
            rng(0);
        end
    end

    methods(Test)
        function doesDropout(test,Size,Probability)
            x = ones(Size);
            y = transformer.layer.dropout(x,Probability);
            % Expectation is that mean(y) is about 1 since we use inverted
            % dropout.
            actMean = mean(y,'all');
            expMean = 1;
            tol = 1e-1;
            test.verifyEqual(actMean,expMean,'AbsTol',tol);
        end

        function isRandom(test)
            % Verify dropout is random on repeated calls
            x = ones([5,4,3]);
            p = 0.5;
            y1 = transformer.layer.dropout(x,p);
            y2 = transformer.layer.dropout(x,p);
            test.verifyNotEqual(y1,y2);
        end

        function supportsDlarrayAndAutodiff(test)
            x = dlarray(ones([5,4,3]));
            p = 0.5;
            f = @(x) transformer.layer.dropout(x,p);
            test.verifyWarningFree(@() f(x));
            function [val,df_val] = df(x)
                val = f(x);
                df_val = dlgradient(sum(val,'all'),x);
            end
            [f_val,df_val] = dlfeval(@df,x);
            % Gradient of dropout is itself
            test.verifyEqual(df_val,f_val);
        end
    end
end

author: matlab-deep-learning

repo: transformer-models


text:
function y = vector_norm(x,dim)
%VECTOR_NORM norm for the rows/columns of a matrix
%
%   Usage: y = vector_norm(x,dim)
%
%   Input parameters:
%       x   - matrix [n x m]
%       dim - dimension along the norm should be calculated
%
%   Output parameter:
%       y   - norm of the vectors [1 x m] or [n x 1]
%
%   VECTOR_NORM(x,dim) calculates the p-norm (with p=2) for the vectors
%   given within the matrix along the dimension dim.
%
%   See also: norm, vector_product

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

%% ===== Checking of input  parameters ==================================
% NOTE: this is disabled due to performance issues in HRTF extrapolation, where
% this function is called multiple times.
%nargmin = 2;
%nargmax = 2;
%narginchk(nargmin,nargmax)
%isargmatrix(x)
%isargpositivescalar(dim)

%% ===== Computation =====================================================
y = sum(abs(x).^2,dim).^(1/2);

author: sfstoolbox

repo: sfs-matlab


text:
function turnDir=turnOrientation(v1,v2,v3)
%%TURNORIENTATION Given 3 two-dimensional vertices in order, determine
%                 whether they form a left turn (go counterclockwise).
%                 Such a determination is necessary for determining
%                 whether a polygon is convex and plays a role in
%                 finding the convex hull of a set of points.
%
%INPUTS: v1, v2, v3 A set of 3 2XN matrices, of N vertex sets in the order
%                   [x;y] These are 2-dimensional vertices in order,
%                   whereby one wishes to determine whether the angle
%                   v1-v2-v3 is going counterclockwise (left) or clockwise
%                   (right) for each set of vertices.
%
%OUTPUTS: turnDir  An NX1 vector where the ith element is is 1 if the ith
%                  set of vertices form a counterclockwise angle, -1 if
%                  they are going clockwise and 0 if they are collinear or
%                  if two of them coincide.
%
%An implementation is given ehre in Matlab. However, the implementation in
%C++ might be more accurate, assuming that one's compilter implements the
%long double datatype as an IEEE-754 souble extended precision datatype.
%A notable exception is Microsoft's Visual Studio (as of 2017), which maps
%long doubles to doubles (64 bits rather than the 80-bit size of Intel/
%AMD floating point registers) as noted in
%https://docs.microsoft.com/en-us/cpp/c-language/type-long-double
%Additionally, the function exactSignOfSumCPP is used to provide the sign
%of an error-free sum. Unfortunately, since the IEEE-754 double extended
%precision datatype only provides a 64-bit mantissa, which is not much of
%an extension over the 54-bit mantissa of a standard double precision
%floating point number, and is not enough to assure no finite precision
%rounding in the products that go into the sum, one cannot guarantee zero
%finite precision errors in this function. On the other hand, if the
%inputs are single-precision (promoted to doubles, since the function does
%not take singles) floating point numbers, the result is exact.
%
%The cross product rule for 3D vectors a and b says that
%norm(cross(a,b))=norm(a)*norm(b)*sin(theta)
%where theta is the positive angle between the vectors. However, If the
%vectors are 2D, (set the z-components to zero), then the cross product
%only have one nonzero component in the z-direction and that component is
%equal to det([a,b]) (for 2D a and b). The interesting thing now, is that
%the sign of the determinant will tell you whether the subsequent vectors
%are going counterclockwise or clockwise. This function returns -1 if the
%vectors are going counterclockwise, 0 if they are exactly collinear and 1
%if they are clockwise.
%
%The Matlab implementation can be used without a compilation. The compiled
%version of the algorithm can be obtained using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%turnDir=turnOrientation(v1,v2,v3);
%
%December 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    S=zeros(6,1);
    S(1)=v3(1)*v1(2);
    S(2)=v1(1)*v2(2);
    S(3)=v2(1)*v3(2);
    S(4)=-v3(1)*v2(2);
    S(5)=-v1(1)*v3(2);
    S(6)=-v2(1)*v1(2);

    turnDir=exactSignOfSum(S);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function out = iszero(f)
%ISZERO    True for zero CHEBTECH objects.
%   ISZERO(F) returns logical TRUE is F.COEFFS has only zero entries and logical
%   FALSE otherwise.

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

out = ~any(f.coeffs, 1);

if ( any(out) )
    % We need this as any(NaN) = 0, which will pass the test above.
    out = out & ~any(isnan(f.coeffs), 1);
    return
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef prtClassMapLog < prtClass
 %prtClassMapLog  Maximum a Posteriori classifier with loglikelihoods
 %  returned

 %    CLASSIFIER = prtClassMap returns a Maximum a Posteriori classifier
 %
 %    CLASSIFIER = prtClassMap(PROPERTY1, VALUE1, ...) constructs a
 %    prtClassMAP object CLASSIFIER with properties as specified by
 %    PROPERTY/VALUE pairs.
 %
 %    A prtClassMap object inherits all properties from the abstract class
 %    prtClass. In addition is has the following property:
 %
 %    rvs    - A prtRv object. This property describes the random variable
 %             model used for Maximum a Posteriori classification.
 %
 %    A prtClassMap object inherits inherits the TRAIN, RUN, CROSSVALIDATE
 %    and KFOLDS methods from prtClass.
 %
 %    Example:
 %
 %    TestDataSet = prtDataGenUnimodal;       % Create some test and
 %    TrainingDataSet = prtDataGenUnimodal;   % training data
 %    classifier = prtClassMap;               % Create a classifier
 %    classifier = classifier.train(TrainingDataSet);    % Train
 %    classified = run(classifier, TestDataSet);         % Test
 %
 %    subplot(2,1,1); classifier.plot;  % Plot results
 %    subplot(2,1,2); prtScoreRoc(classified,TestDataSet);
 %    set(get(gca,'Children'), 'LineWidth',3)
 %
 %    See also prtClass, prtClassLogisticDiscriminant, prtClassBagging,
 %    prtClassMap, prtClassFld, prtClassBinaryToMaryOneVsAll, prtClassDlrt,
 %    prtClassPlsda, prtClassFld, prtClassRvm, prtClassGlrt,  prtClassSvm,
 %    prtClassTreeBaggingCap, prtClassKmsd, prtClassKnn

    properties (SetAccess=private)
        % Required by prtAction
        name = 'Maximum a Posteriori Log'   % Maximum a Posteriori
        nameAbbreviation = 'MAPLog'        % MAP
        isNativeMary = true;            % True
    end

    properties
        rvs = prtRvMvn; % Random variable object containing mean and variance
    end

    methods
        % Constructor
        function Obj = prtClassMapLog(varargin)

            Obj.classTrain = 'prtDataInterfaceCategoricalTargets';
            Obj.classRun = 'prtDataSetBase';
            Obj.classRunRetained = false;

            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end
        % Set function
        function Obj = set.rvs(Obj,val)
            if ~(isa(val, 'prtRv') || isa(val, 'prtBrv'))
                error('prtClassMapLog:rvs','Rvs parameter must be of class prtRv');
            else
                Obj.rvs = val;
            end
        end
    end

    methods (Access = protected, Hidden = true)

        function Obj = trainAction(Obj,DataSet)
            % Repmat the rv objects to get one for each class
            Obj.rvs = repmat(Obj.rvs(:), (DataSet.nClasses - length(Obj.rvs)+1),1);
            Obj.rvs = Obj.rvs(1:DataSet.nClasses);

            % Get the ML estimates of the RV parameters for each class
            for iY = 1:DataSet.nClasses
                Obj.rvs(iY) = train(Obj.rvs(iY), DataSet.retainClassesByInd(iY));
            end
        end

        function DataSet = runAction(Obj,DataSet)

            logLikelihoods = zeros(DataSet.nObservations, length(Obj.rvs));
            for iY = 1:length(Obj.rvs)
                logLikelihoods(:,iY) = getObservations(run(Obj.rvs(iY), DataSet));
            end

            DataSet = prtDataSetClass(logLikelihoods);
        end

        function xOut = runActionFast(Obj,xIn,ds) %#ok<INUSD>

            xOut = zeros(size(xIn,1),length(Obj.rvs));
            for iY = 1:length(Obj.rvs)
                xOut(:,iY) = runFast(Obj.rvs(iY), xIn);
            end

            %xOut = exp(bsxfun(@minus, logLikelihoods, prtUtilSumExp(logLikelihoods.').'));
        end
    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%Finite Element Method 101-2
%National Taiwan University
%2D Elasticity problem

%%clear memory
close all; clc; clear all;

%% Load Mesh for Lab 10
[nodeCoordinates,elementNodes]=Mesh_Lab10('Q12');
% node coordinates are given in mm
NodePerElement=12;
numberNodes=size(nodeCoordinates,1);
numberElements=size(elementNodes,1);

%% Import BCs

% Essential BC's
GDof = 2*numberNodes;
prescribedDof = [1,2,3,4,5,6,7,8];

% Natural BC's
force = zeros(GDof,1);
force(end) = -10; % 10 [kN]

%% Import material and section properties
E = 3E7; % [GPa]
poisson = 0.3; %[-]
thickness = 1; % [mm]

%% Evalute force vector
%force=formForceVectorQ4(GDof,naturalBCs,surfaceOrientation,...
%    elementNodes,nodeCoordinates,thickness);

%% Construct Stiffness matrix for Q12 element
D=E/(1-poisson^2)*[1 poisson 0;poisson 1 0;0 0 (1-poisson)/2];

stiffness=formStiffness2D(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,D,thickness);

%% solution
displacements=solution(GDof,prescribedDof,stiffness,force);

%% output displacements
outputDisplacements(displacements, numberNodes, GDof);
scaleFactor=1.E5;
drawingMesh(nodeCoordinates+scaleFactor*[displacements(1:2:2*numberNodes) ...
    displacements(2:2:2*numberNodes)],elementNodes,'Q4','r--');

%% B matrix & strain
% 4 by 4 quadrature
[gaussWeights,gaussLocations]=gauss2d('4x4');

for e=1:numberElements
  numNodePerElement = length(elementNodes(e,:));
  numEDOF = 2*numNodePerElement;
  elementDof=zeros(1,numEDOF);
  for i = 1:numNodePerElement
      elementDof(2*i-1)=2*elementNodes(e,i)-1;
      elementDof(2*i)=2*elementNodes(e,i);
  end

  % cycle for Gauss point
  for q=1:size(gaussWeights,1)
      GaussPoint=gaussLocations(q,:);
      xi=GaussPoint(1);
      eta=GaussPoint(2);

% shape functions and derivatives
    [shapeFunction,naturalDerivatives]=shapeFunctionQ12(xi,eta);

% Jacobian matrix, inverse of Jacobian,
% derivatives w.r.t. x,y
    [Jacob,invJacobian,XYderivatives]=...
        Jacobian(nodeCoordinates(elementNodes(e,:),:),naturalDerivatives);

%  B matrix
    B=zeros(3,numEDOF);
    B(1,1:2:numEDOF)       = XYderivatives(:,1)';
    B(2,2:2:numEDOF)  = XYderivatives(:,2)';
    B(3,1:2:numEDOF)       = XYderivatives(:,2)';
    B(3,2:2:numEDOF)  = XYderivatives(:,1)';

  elementNodes(e,:);
  dis(1:2:24)=displacements([2*(elementNodes(e,:))-1],1);
  dis(2:2:24)=displacements([2*(elementNodes(e,:))],1);
  strain=B*dis';
  stress(:,q)=D*strain;

  end
  for q=1:size(gaussWeights,1)
      GaussPoint=gaussLocations(q,:);
      xi=1/GaussPoint(1);
      eta=1/GaussPoint(2);
      [shapeFunction,naturalDerivatives]=shapeFunctionQ12(xi,eta);
      stressxx(e,q)=stress(1,1:12)*shapeFunction;
      stressyy(e,q)=stress(2,1:12)*shapeFunction;
      stressxy(e,q)=stress(3,1:12)*shapeFunction;
      vonmises(e,q)=sqrt(0.5*((stressxx(e,q)-stressyy(e,q))^2+(stressyy(e,q))^2+(stressxx(e,q))^2+6*(stressxy(e,q))^2));
      fprintf('\nStress in element %u\n',e)
      fprintf('\nStress in node %u\n',q)
      fprintf('Sigma_xx : %0.6f\n',stressxx(e,q))
      fprintf('Sigma_yy : %0.6f\n',stressyy(e,q))
      fprintf('Sigma_xy : %0.6f\n',stressxy(e,q))
      fprintf('Vonmises : %0.6f\n',vonmises(e,q))
  end

end

repo: Aero-matlab


text:
function [ know, x ] = p09_sol ( n )

%*****************************************************************************80
%
%% P09_SOL returns the solution for problem 9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 0;

  x = zeros ( n, 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function SpeedTuning=SpeedTuningCalcultation(SpeedTrain,SpikeTrain,FrameRate,SpeedRange,MinSpanTime,extend)
% FrameRate=FullNeuronBehaviorDataSet.FrameRate./double(FullNeuronBehaviorDataSet.ImagingPlane);
% Frameselect=find(~isnan(FullNeuronBehaviorDataSet.FullNeuronBehaviorMatrix{1,2}(:,12+4*i)));
% SpeedTrain=FullNeuronBehaviorDataSet.FullNeuronBehaviorMatrix{1,2}(Frameselect,5);
% SpikeTrain=FullNeuronBehaviorDataSet.FullNeuronBehaviorMatrix{1,2}(Frameselect,4*i+12);
% SpeedRange=[0 4 8 12];
% MinSpanTime=10;

SpeedTuning=struct;
SpeedTuning.Rate=zeros(size(SpeedRange));
SpeedTuning.Occ=zeros(size(SpeedRange));
SpeedTuning.SpikeCount=zeros(size(SpeedRange));
SpeedTuning.SpeedSelectivity=0;
SpeedTuning.SpeedScore=0;
SpeedTuning.PeakRate=0;
SpeedTuning.PeakSpeed=1;
SpeedTuning.SpeedRange=SpeedRange;
SpeedRangeExpend=[SpeedRange,extend];
for j=1:1:length(SpeedTrain)
    for k=1:1:size(SpeedRangeExpend,2)-1
        if SpeedTrain(j)>=SpeedRangeExpend(k)&&SpeedTrain(j)<SpeedRangeExpend(k+1)
            SpeedTuning.Occ(k)=SpeedTuning.Occ(k)+1;
            SpeedTuning.SpikeCount(k)=SpeedTuning.SpikeCount(k)+SpikeTrain(j);
        else
        end
    end
end
SpeedTuning.Occ=SpeedTuning.Occ./FrameRate;% frame to second;
SpeedTuning.Occ(SpeedTuning.Occ==0)=0.0000001;
SpeedTuning.Rate=SpeedTuning.SpikeCount./SpeedTuning.Occ;
% SpeedTuning.Rate=smooth(SpeedTuning.Rate,speedsmooth);
SpeedTuning.Rate(SpeedTuning.Occ<MinSpanTime)=nan;

[SpeedTuning.PeakRate,SpeedTuning.PeakSpeed]=max(SpeedTuning.Rate(~isnan(SpeedTuning.Rate)));
% SpeedTuning.PeakSpeed=SpeedRange(SpeedTuning.PeakSpeed);
MinRate=min(SpeedTuning.Rate(~isnan(SpeedTuning.Rate)));
SpeedTuning.SpeedSelectivity=(SpeedTuning.PeakRate-MinRate)./(SpeedTuning.PeakRate+MinRate);

% Con1=find((SpeedRange<=ConditionSeperation).*(~isnan(SpeedTuning.Rate))==1);
% Con2=find((SpeedRange>ConditionSeperation).*(~isnan(SpeedTuning.Rate))==1);
% Condtion1=mean(SpeedTuning.Rate(Con1));
% Condtion2=mean(SpeedTuning.Rate(Con2));
% SpeedTuning.SpeedScore=(Condtion2-Condtion1)./(Condtion1+Condtion2);

end

% end

author: kavli-ntnu


text:
%  udc Undetermined coefficients solution algorithm for DSGE models.
%  The procedure can find all possible solutions for a constant-parameter
%  DSGE model or a regime-switching DSGE model with diagonal transition
%  matrix.
%
%  ::
%
%    [Tz_pb,eigvals,retcode]=udc(A,B,C,Q,T0,TolFun,maxiter)
%
%  Args:
%
%     - **A** [n x n] : Coefficient matrix on lead variables
%
%      - **B** [n x n] : Coefficient matrix on current variables
%
%      - **C** [n x n] : Coefficient matrix on lagged variables
%
%      - **slvOpts.allSols** [true|{false}] : flag for finding all solutions
%
%      - **slvOpts.msvOnly** [{true}|false] : return only MSV solutions
%
%      - **slvOpts.debug** [true|{false}] : slvOpts.debug or not
%
%  Returns:
%     :
%
%      - **Tz_pb** [n x n x h x k array] : Solution set (k solutions)
%
%      - **eigvals** [empty] : Eigenvalues (Not computed)
%
%      - **retcode** [numeric] : 0 if there is no problem
%
%  See also :  groebner
%

author: jmaih


text:
function [ xmin, ixmin ] = r8row_min ( m, n, x )

%*****************************************************************************80
%
%% R8ROW_MIN returns the minimums of rows of an R8ROW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, real X(M,N), the R8ROW.
%
%    Output, real XMIN(M), the minimums of the rows of X.
%
%    Output, integer IXMIN(M); IXMIN(I) is the column of X in which
%    the minimum for row I occurs.
%
  for i = 1 : m

    ixmin(i) = 1;
    xmin(i) = x(i,1);
    for j = 2 : n
      if ( x(i,j) < xmin(i) )
        ixmin(i) = j;
        xmin(i) = x(i,j);
      end
    end

  end

  return
end


text:
function value = circle_segment_contains_point ( r, c, omega1, omega2, xy )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_CONTAINS_POINT reports whether a point is in a circle segment.
%
%  Discussion:
%
%    Begin with a circle of radius R.  Choose two points P1 and P2 on the
%    circle, and draw the chord P1:P2.  This chord divides the circle
%    into two pieces, each of which is called a circle segment.
%    Consider one of the pieces.  The "angle" of this segment is the angle
%    P1:C:P2, where C is the center of the circle.  Let Q be the point on
%    the chord P1:P2 which is closest to C.  The "height" of the segment
%    is the distance from Q to the perimeter of the circle.
%
%    In this function, we allow the circle to have an arbitrary center C,
%    arbitrary radius R, and we describe the points P1 and P2 by specifying
%    their angles OMEGA1 and OMEGA2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%    0 < R.
%
%    Input, real C(2,1), the center of the circle.
%
%    Input, real OMEGA1, OMEGA2, the angles of the two points on the
%    circumference of the circle that define the circle segment.
%    OMEGA1 < OMEGA2 <= OMEGA1 + 2 * PI
%
%    Input, real XY(2,1), a point.
%
%    Output, integer VALUE, is TRUE if the point is inside the circle segment.
%
  if ( r <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CIRCLE_SEGMENT_CONTAINS_POINT - Fatal error!\n' );
    fprintf ( 1, '  R <= 0.0.\n' );
    error ( 'CIRCLE_SEGMENT_CONTAINS_POINT - Fatal error!' )
  end

  while ( omega2 < omega1 )
    omega2 = omega2 + 2.0 * pi;
  end
%
%  Destroy all row vectors.
%
  c = c(:);
  xy = xy(:);
%
%  Compute the vector V = XY - C:
%
  v(1:2,1) = xy(1:2,1) - c(1:2,1);
%
%  a: Point must be inside the circle, so ||V|| <= R.
%
  v_r = norm ( v );
  if ( r < v_r )
    value = 0;
    return
  end
%
%  b: Angle made by the vector V must be between OMEGA1 and OMEGA2.
%
  v_omega = atan2 ( v(2), v(1) );

  while ( omega1 <= v_omega + 2.0 * pi )
    v_omega = v_omega - 2.0 * pi;
  end

  while ( v_omega + 2.0 * pi <= omega1 )
    v_omega = v_omega + 2.0 * pi;
  end

  if ( omega2 < v_omega )
    value = 0;
    return
  end
%
%  c: Projection of V onto unit centerline must be at least R-H.
%
  omegah = 0.5 * ( omega1 + omega2 );
  v_project = v(1) * cos ( omegah ) + v(2) * sin ( omegah );

  theta = omega2 - omega1;
  h = circle_segment_height_from_angle ( r, theta );

  if ( v_project < r - h )
    value = 0;
    return
  end

  value = 1;

  return
end


text:
function test_tutorial_connectivity20130308

% MEM 7gb
% WALLTIME 00:10:00
% DEPENDENCY

% Simulated data with directed connections
% We will first simulate some data with a known connectivity structure built in. This way we know what to expect in terms of connectivity. To simulate data we use ft_connectivitysimulation. We will use an order 2 multivariate autoregressive model. The necessary ingredients are a set of NxN coefficient matrices, one matrix for each time lag. These coefficients need to be stored in the cfg.param field. Next to the coefficients we have to specify the NxN covariance matrix of the innovation noise. This matrix needs to be stored in the cfg.noisecov field. The model we are going to use to simulate the data is as follows:
%
% x(t) = 0.8*x(t-1) - 0.5*x(t-2)
% y(t) = 0.9*y(t-1) + 0.5*z(t-1) - 0.8*y(t-2)
% z(t) = 0.5*z(t-1) + 0.4*x(t-1) - 0.2*z(t-2)

cfg             = [];
cfg.ntrials     = 500;
cfg.triallength = 1;
cfg.fsample     = 200;
cfg.nsignal     = 3;
cfg.method      = 'ar';

cfg.params(:,:,1) = [ 0.8    0    0 ;
                        0  0.9  0.5 ;
                      0.4    0  0.5];

cfg.params(:,:,2) = [-0.5    0    0 ;
                        0 -0.8    0 ;
                        0    0 -0.2];

cfg.noisecov      = [ 0.3    0    0 ;
                        0    1    0 ;
                        0    0  0.2];

data              = ft_connectivitysimulation(cfg);

% The simulated data consists of 3 channels in 500 trials. You can easily visualize the data for example in the first trial using
%
figure
plot(data.time{1}, data.trial{1})
legend(data.label)
xlabel('time (s)')

% or browse through the complete data using
%
cfg = [];
cfg.viewmode = 'vertical';  % you can also specify 'butterfly'
ft_databrowser(cfg, data);

%
%
% Computation of the multivariate autoregressive model
% To be able to compute spectrally resolved Granger causality, or other frequency-domain directional measures of connectivity, we have to fit an autoregressive model to the data. This is done using the ft_mvaranalysis function.
%
% For the actual computation of the autoregressive coefficients FieldTrip makes use of an implementation from third party toolboxes. At present ft_mvaranalysis supports the biosig and bsmart toolboxes for these computations.
%
% In this tutorial we will use the bsmart toolbox. The relevant functions have been included in the FieldTrip release in the fieldtrip/external/bsmart directory.
%
cfg         = [];
cfg.order   = 5;
cfg.toolbox = 'bsmart';
mdata       = ft_mvaranalysis(cfg, data);

% mdata =
%          dimord: 'chan_chan_lag'
%           label: {3x1 cell}
%          coeffs: [3x3x5 double]
%        noisecov: [3x3 double]
%             dof: 500
%     fsampleorig: 200
%             cfg: [1x1 struct]
%
% The resulting variable mdata contains a description of the data in terms of a multivariate autoregressive model. For each time-lag up to the model order (which is 5 in this case), a 33 matrix of coefficients is outputted. The noisecov-field contains covariance matrix of the model's residuals.
%
% Exercise 1
% Compare the parameters specified for the simulation with the estimated coefficients and discuss.
%
% Computation of the spectral transfer function
% From the autoregressive coefficients it is now possible to compute the spectral transfer matrix, for which we use ft_freqanalysis.
%
cfg        = [];
cfg.method = 'mvar';
mfreq      = ft_freqanalysis(cfg, mdata);

% mfreq =
%         label: {3x1 cell}
%          freq: [1x101 double]
%        dimord: 'chan_chan_freq'
%      transfer: [3x3x101 double]
%      noisecov: [3x3 double]
%     crsspctrm: [3x3x101 double]
%           dof: 500
%           cfg: [1x1 struct]
%
% The resulting mfreq data structure contains the pairwise transfer function between the 3 channels for 101 frequencies.
%
% It is also possible to compute the spectral transfer function using non-parametric spectral factorization of the cross-spectral density matrix. For this, we need a Fourier decomposition of the data. This is done in the following section.
%
%
% Non-parametric computation of the cross-spectral density matrix
% Some connectivity metrics can be computed from a non-parametric spectral estimate (i.e. after the application of the FFT-algorithm and conjugate multiplication to get cross-spectral densities), such as coherence, phase-locking value and phase slope index. The following part computes the fourier-representation of the data using ft_freqanalysis. It is not necessary to compute the cross-spectral density at this stage, because the function used in the next step, ft_connectivityanalysis, contains functionality to compute the cross-spectral density from the fourier coefficients.
%
cfg           = [];
cfg.method    = 'mtmfft';
cfg.taper     = 'dpss';
cfg.output    = 'fourier';
cfg.tapsmofrq = 2;
freq          = ft_freqanalysis(cfg, data);

% freq =
%             label: {3x1 cell}
%            dimord: 'rpttap_chan_freq'
%              freq: [1x101 double]
%     fourierspctrm: [1500x3x101 double]
%         cumsumcnt: [500x1 double]
%         cumtapcnt: [500x1 double]
%               cfg: [1x1 struct]
% The resulting freq structure contains the spectral estimate for 3 tapers in each of the 500 trials (hence 1500 estimates), for each of the 3 channels and for 101 frequencies.
%
%
% Computation and inspection of the connectivity measures
% The actual computation of the connectivity metric is done by ft_connectivityanalysis. This function is transparent to the type of input data, i.e. provided the input data allows the requested metric to be computed, the metric will be calculated. Here, we provide an example for the computation and visualization of the coherence coefficient.
%
cfg           = [];
cfg.method    = 'coh';
coh           = ft_connectivityanalysis(cfg, freq);
cohm          = ft_connectivityanalysis(cfg, mfreq);
% Subsequently, the data can be visualized using ft_connectivityplot.
%
cfg           = [];
cfg.parameter = 'cohspctrm';
cfg.zlim      = [0 1];
ft_connectivityplot(cfg, coh, cohm);

%
% The coherence measure is a symmetric measure, which means that it does not provide information regarding the direction of information flow between any pair of signals. In order to analyze directionality in interactions, measures based on the concept of granger causality can be computed. These measures are based on an estimate of the spectral transfer matrix, which can be computed in a straightforward way from the multivariate autoregressive model fitted to the data.
%
cfg           = [];
cfg.method    = 'granger';
granger       = ft_connectivityanalysis(cfg, mfreq);

cfg           = [];
cfg.parameter = 'grangerspctrm';
cfg.zlim      = [0 1];
ft_connectivityplot(cfg, granger);

% Instead of plotting it with ft_connectivityplot, you can use the following low-level MATLAB plotting code which gives a better understanding of the numerical representation of the results.
%
figure
for row=1:3
for col=1:3
  subplot(3,3,(row-1)*3+col);
  plot(granger.freq, squeeze(granger.grangerspctrm(row,col,:)))
  ylim([0 1])
end
end

%
% Exercise 2
% Discuss the differences between the granger causality spectra, and the coherence spectra.
% Exercise 3
% Compute the following connectivity measures from the mfreq data, and visualize and discuss the results: partial directed coherence (pdc), directed transfer function (dtf), phase slope index (psi)
%
% Simulated data with common pick-up and different noise levels
%  this is under progress
%
% When working with electrophysiological data (EEG/MEG/LFP) the signals that are picked up by the individual channels invariably consist of instantaneous mixtures of the underlying source signals. This mixing can severely affect the outcome of connectivity analysis, and thus affects the interpretation. We will demonstrate this by simulating data in 2 channels, where each of the channels consists of a weighted combination of temporally white noise unique to each of the channels, and a common input of a band-limited signal (filtered between 15 and 25 Hz). We will compute connectivity between these channels, and show that the common input can give rise to spurious estimates of connectivity.
%
% % create some instantaneously mixed data
%
% define some variables locally
nTrials  = 100;
nSamples = 1000;
fsample  = 1000;

% mixing matrix
mixing   = [0.8 0.2 0;
              0 0.2 0.8];

data       = [];
data.trial = cell(1,nTrials);
data.time  = cell(1,nTrials);
for k = 1:nTrials
  dat = randn(3, nSamples);
  dat(2,:) = ft_preproc_bandpassfilter(dat(2,:), 1000, [15 25]);
  dat = 0.2.*(dat-repmat(mean(dat,2),[1 nSamples]))./repmat(std(dat,[],2),[1 nSamples]);
  data.trial{k} = mixing * dat;
  data.time{k}  = (0:nSamples-1)./fsample;
end
data.label = {'chan1' 'chan2'}';

figure;plot(dat'+repmat([0 1 2],[nSamples 1]));
title('original ''sources''');

figure;plot((mixing*dat)'+repmat([0 1],[nSamples 1]));
axis([0 1000 -1 2]);
set(findobj(gcf,'color',[0 0.5 0]), 'color', [1 0 0]);
title('mixed ''sources''');

% do spectral analysis
cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'fourier';
cfg.foilim    = [0 200];
cfg.tapsmofrq = 5;
freq          = ft_freqanalysis(cfg, data);
fd            = ft_freqdescriptives(cfg, freq);

figure;plot(fd.freq, fd.powspctrm);
set(findobj(gcf,'color',[0 0.5 0]), 'color', [1 0 0]);
title('powerpectrum');

%
% compute connectivity
cfg = [];
cfg.method = 'granger';
g = ft_connectivityanalysis(cfg, freq);
cfg.method = 'coh';
c = ft_connectivityanalysis(cfg, freq);
% visualize the results
cfg = [];
cfg.parameter = 'grangerspctrm';
ft_connectivityplot(cfg, g);
cfg.parameter = 'cohspctrm';
ft_connectivityplot(cfg, c);

% Exercise 4
% Simulate new data using the following mixing matrix:
% [0.9 0.1 0;0 0.2 0.8]
% and recompute the connectivity measures. Discuss what you see.
%
% Exercise 5
% Play a bit with the parameters in the mixing matrix and see what is the effect on the estimated connectivity.
% Exercise 6
% Simulate new data where the 2 mixed signals are created from 4 underlying sources, and where two of these sources are common input to both signals, and where these two sources are temporally shifted copies of one another.
% Hint: the mixing matrix could look like this:
%
% [a b c 0; 0 d e f];
% and the trials could be created like this:
%
% for k = 1:nTrials
%   dat = randn(4, nSamples+10);
%   dat(2,:) = ft_preproc_bandpassfilter(dat(2,:), 1000, [15 25]);
%   dat(3,1:(nSamples)) = dat(2,11:(nSamples+10));
%   dat = dat(:,1:1000);
%   dat = 0.2.*(dat-repmat(mean(dat,2),[1 nSamples]))./repmat(std(dat,[],2),[1 nSamples]);
%   data.trial{k} = mixing * dat;
%   data.time{k}  = (0:nSamples-1)./fsample;
% end
% Compute connectivity between the signals and discuss what you observe. In particular, also compute measures of directed interaction.
%
%
% Connectivity between MEG virtual channel and EMG
% The previous two examples were using simulated data, either with a clear directed connectivity structure, or with a trivial pick-up of a common source in two channels. We will now continue with connectivity analysis on real MEG data. The dataset is the same as the one used in the Analysis of corticomuscular coherence tutorial.
%
% In short, the dataset consists of combined MEG and EMG recordings while the subject lifted his right hand. The coherence tutorial introduction contains a more elaborate description of the experiment and the dataset and a detailed analysis can be found in the corresponding paper 1). Due to the long distance between the EMG and the MEG, there is no volume conduction and hence no common pick-up. Hence this dataset lends itself well for connectivity analysis. But rather than doing an analysis between the EMG and one of the MEG channels (as in the original study), we will extract the cortical activity using a beamformer virtual channel.
%
%
% Compute the spatial filter for the region of interest
% We start with determining the motor cortex as region of interest. At the end of the coherence tutorial it is demonstrated how to make a 3-D reconstruction of the cortico-muscolar coherence (CMC) using the DICS algorithm. That source reconstruction serves as starting point for this analysis.
%
% You can download the result from the DICS reconstruction from the FieldTrip ftp server (source.mat)
%
% We will first determine the position on which the cortico-muscular coherence is the largest.
%

cd(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/connectivity'));
load source

[maxval, maxindx] = max(source.avg.coh);
maxpos = source.pos(maxindx,:);

% maxpos =
%     4 -3 12
% The cortical position is expressed in individual subject head-coordinates and in centimeter. Relative to the center of the head (in between the ears) the position is 4 cm towards the nose, -3 towards the left side (i.e., 3 cm towards the right!) and 12 cm towards the vertex.
%
% The ft_sourceanalysis methods are usually applied to the whole brain using a regular 3-D grid or using a triangulated cortical sheet. You can also just specify the location of a single or multiple points of interest with cfg.sourcemodel.pos and the LCMV beamformer will simply be performed at the location of interest.
%
% The LCMV beamformer spatial filter for the location of interest will pass the activity at that location with unit-gain, while optimally suppressing all other noise and other source contributions to the MEG data. The LCMV implementation in FieldTrip requires the data covariance matrix to be computed with ft_timelockanalysis.
%
% Rather than doing all the preprocessing again, you can download the preprocessed data from the FieldTrip ftp server (data.mat)
%
load data
%
%% compute the beamformer filter
cfg                   = [];
cfg.covariance        = 'yes';
cfg.channel           = 'MEG';
cfg.vartrllength      = 2;
cfg.covariancewindow  = 'all';
timelock              = ft_timelockanalysis(cfg, data);

cfg             = [];
cfg.method      = 'lcmv';
cfg.hdmfile     = 'SubjectCMC.hdm';
cfg.sourcemodel.pos = maxpos;
cfg.keepfilter  = 'yes';
source          = ft_sourceanalysis(cfg, timelock);
% The source reconstruction contains the estimated power and the source-level time-series of the averaged ERF, but here we are not interested in those. The cfg.keepfilter option results in the spatial filter being kept in the output source structure. That spatial can be used to reconstruct the single-trial time series as a virtual channel by multiplying it with the original MEG data.
%
%
% Extract the virtual channel time-series
%% construct the 3-D virtual channel at the location of interest
beamformer = source.avg.filter{1};

chansel = ft_channelselection('MEG', data.label); % find the names
chansel = match_str(data.label, chansel);         % find the indices

sourcedata = [];
sourcedata.label = {'x', 'y', 'z'};
sourcedata.time = data.time;
for i=1:length(data.trial)
  sourcedata.trial{i} = beamformer * data.trial{i}(chansel,:);
end
% The LCMV spatial filter is computed here without applying any time-domain filters. Consequently, it will have to suppress all noise in the data in all frequency bands. The spatial filter derived from the broadband data allows us to compute a broadband source level time-series.
% If you would know that the subsequent analysis would be limited to a specific frequency range in the data (e.g. everything above 30 Hz), you could first apply a filter using ft_preprocessing (e.g. cfg.hpfilter=yes and cfg.hpfreq=30) prior to computing the covariance and the spatial filter.
%
% The sourcedata structure resembles the raw-data output of ft_preprocessing and consequently can be used in any follow-up function. You can for example visualize the single-trial virtual channel time-series using ft_databrowser:
%
cfg = [];
cfg.viewmode = 'vertical';  % you can also specify 'butterfly'
ft_databrowser(cfg, sourcedata);

%
% Notice that the reconstruction contains three channels, for the x-, the y- and the z-component of the equivalent current dipole source at the location of interest.
%
%
% Project along the strongest dipole direction
% The interpretation of connectivity is facilitated if we can compute it between two plain channels rather than between one channel versus a triplet of channels. Therefore we will project the time-series along the dipole direction that explains most variance. This projection is equivalent to determining the largest (temporal) eigenvector and can be computationally performed using the singular value decomposition (svd).
%
%% construct a single virtual channel in the maximum power orientation
timeseries = cat(2, sourcedata.trial{:});

[u, s, v] = svd(timeseries, 'econ');

% whos u s v
%   Name           Size              Bytes  Class     Attributes
%
%   s              3x3                  72  double
%   u              3x3                  72  double
%   v         196800x3             4723200  double
% Matrix u contains the spatial decomposition, matrix v the temporal and on the diagonal of matrix s you can find the eigenvalues. See "help svd" for more details.
%
% We now recompute the virtual channel time-series, but now only for the dipole direction that has the most power.
%
% this is equal to the first column of matrix V, apart from the scaling with s(1,1)
timeseriesmaxproj = u(:,1)' * timeseries;
virtualchanneldata = [];
virtualchanneldata.label = {'cortex'};
virtualchanneldata.time = data.time;
for i=1:length(data.trial)
  virtualchanneldata.trial{i} = u(:,1)' * beamformer * data.trial{i}(chansel,:);
end
% Rather than using a sourcemodel in the beamformer that consists of all three (x, y, z) directions, you can also have the beamformer compute the filter for only the optimal source orientation. This is implemented using the cfg.lcmv.fixedori='yes' option.
% Recompute the spatial filter for the optimal source orientation and using that spatial filter (a 1151 vector) recompute the time-series.
%
% Investigate and describe the difference between the two time-series. What is the difference between the two dipole orientations?
%
% Note that one orientation is represented in the SVD matrix "u" and the other is in the source.avg.ori field.
%
%
% Combine the virtual channel with the EMG
% The raw data structure containing one (virtual) channel can be combined with the two EMG channels from the original preprocessed data.
%
%% select the two EMG channels
cfg = [];
cfg.channel = 'EMG';
emgdata = ft_selectdata(cfg, data);

%% combine the virtual channel with the two EMG channels
cfg = [];
combineddata = ft_appenddata(cfg, virtualchanneldata, emgdata);

% save combineddata combineddata
%
% Compute the connectivity
% The resulting combined data structure has three channels: the activity from the cortex, the left EMG and the right EMG. We can now continue with regular channel-level connectivity analysis.
%
%% compute the spectral decomposition
cfg            = [];
cfg.output     = 'fourier';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channel    = {'cortex' 'EMGlft' 'EMGrgt'};
freq    = ft_freqanalysis(cfg, combineddata);

cfg = [];
cfg.method = 'coh';
coherence = ft_connectivityanalysis(cfg, freq);
% This computes the spectral decomposition and the coherence spectrum between all channel pairs, which can be plotted with
%
cfg = [];
cfg.zlim = [0 0.2];
figure
ft_connectivityplot(cfg, coherence);
title('coherence')

%
% To look in more detail into the numerical representation of the coherence results, you can use
%
figure
plot(coherence.freq, squeeze(coherence.cohspctrm(1,2,:)))
title(sprintf('connectivity between %s and %s', coherence.label{1}, coherence.label{2}));
xlabel('freq (Hz)')
ylabel('coherence')

%
% The spectrum reveals coherence peaks at 10 and 20 Hz (remember that the initial DICS localizer was done at beta). Furthermore, there is a broader plateau of coherence in the gamma range from 40-50 Hz.
%
% The spectral decomposition was performed with mutitapering and 5 Hz spectral smoothing (i.e. 5Hz in both directions). Recompute the spectral decomposition and the coherence with a hanning taper. Recompute it with mutitapering and 10 Hz smoothing. Plot the three coherence spectra and look at the differences.
% Rather than looking at undirected coherence, the virtual channel level data can now also easily be submitted to directed connectivity measures. Compute the spectrally resolved granger connectivity and try to assess whether the directionality is from cortex to EMG or vice versa.
%
% Summary and further reading
% This tutorial demonstrates how to compute connectivity measures between two time series. If you want to learn how to make a distributed representation of connectivity throughout the whole brain, you may want to continue with the corticomuscular coherence tutorial.
%
% This tutorial was last tested by Robert with revision 6026 of FieldTrip (~20120611) on a 64-bit Mac OS X machine using MATLAB 2011b.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [vol, sens] = ft_prepare_vol_sens(vol, sens, varargin)

% FT_PREPARE_VOL_SENS does some bookkeeping to ensure that the volume
% conductor model and the sensor array are ready for subsequent forward
% leadfield computations. It takes care of some pre-computations that can
% be done efficiently prior to the leadfield calculations.
%
% Use as
%   [vol, sens] = ft_prepare_vol_sens(vol, sens, ...)
% with input arguments
%   sens   structure with gradiometer or electrode definition
%   vol    structure with volume conductor definition
%
% The vol structure represents a volume conductor model, its contents
% depend on the type of model. The sens structure represents a sensor
% array, i.e. EEG electrodes or MEG gradiometers.
%
% Additional options should be specified in key-value pairs and can be
%   'channel'    cell-array with strings (default = 'all')
%   'order'      number, for single shell "Nolte" model (default = 10)
%
% The detailled behaviour of this function depends on whether the input
% consists of EEG or MEG and furthermoree depends on the type of volume
% conductor model:
% - in case of EEG single and concentric sphere models, the electrodes are
%   projected onto the skin surface.
% - in case of EEG boundary element models, the electrodes are projected on
%   the surface and a blilinear interpoaltion matrix from vertices to
%   electrodes is computed.
% - in case of MEG and a multispheres model, a local sphere is determined
%   for each coil in the gradiometer definition.
%  - in case of MEG with a singleshell Nolte model, the volume conduction
%    model is initialized
% In any case channel selection and reordering will be done. The channel
% order returned by this function corresponds to the order in the 'channel'
% option, or if not specified, to the order in the input sensor array.
%
% See also FT_READ_VOL, FT_READ_SENS, FT_TRANSFORM_VOL, FT_TRANSFORM_SENS, FT_COMPUTE_LEADFIELD

% Copyright (C) 2004-2009, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: ft_prepare_vol_sens.m 2833 2011-02-06 20:19:54Z crimic $

% get the options
% fileformat = keyval('fileformat',  varargin);
channel = keyval('channel',  varargin);  % cell-array with channel labels
order   = keyval('order',    varargin);  % order of expansion for Nolte method; 10 should be enough for real applications; in simulations it makes sense to go higher

% set the defaults
if isempty(channel),  channel = sens.label;   end
if isempty(order),    order = 10;             end

% determine whether the input contains EEG or MEG sensors
iseeg = ft_senstype(sens, 'eeg');
ismeg = ft_senstype(sens, 'meg');

% determine the skin compartment
if ~isfield(vol, 'skin_surface')
  if isfield(vol, 'bnd')
    vol.skin_surface   = find_outermost_boundary(vol.bnd);
  elseif isfield(vol, 'r') && length(vol.r)<=4
    [dum, vol.skin_surface] = max(vol.r);
  end
end

% determine the inner_skull_surface compartment
if ~isfield(vol, 'inner_skull_surface')
  if isfield(vol, 'bnd')
    vol.inner_skull_surface  = find_innermost_boundary(vol.bnd);
  elseif isfield(vol, 'r') && length(vol.r)<=4
    [dum, vol.inner_skull_surface] = min(vol.r);
  end
end

% otherwise the voltype assignment to an empty struct below won't work
if isempty(vol)
  vol = [];
end

% this makes them easier to recognise
sens.type = ft_senstype(sens);
vol.type  = ft_voltype(vol);

if isfield(vol, 'unit') && isfield(sens, 'unit') && ~strcmp(vol.unit, sens.unit)
  error('inconsistency in the units of the volume conductor and the sensor array');
end

if ismeg && iseeg
  % this is something that could be implemented relatively easily
  error('simultaneous EEG and MEG not yet supported');

elseif ~ismeg && ~iseeg
  error('the input does not look like EEG, nor like MEG');

elseif ismeg

  % keep a copy of the original sensor array, this is needed for the MEG multisphere model
  sens_orig = sens;

  % always ensure that there is a linear transfer matrix for combining the coils into gradiometers
  if ~isfield(sens, 'tra');
    Nchans = length(sens.label);
    Ncoils = size(sens.pnt,1);
    if Nchans~=Ncoils
      error('inconsistent number of channels and coils');
    end
    sens.tra = sparse(eye(Nchans, Ncoils));
  end

  % select the desired channels from the gradiometer array
  % order them according to the users specification
  [selchan, selsens] = match_str(channel, sens.label);

  % first only modify the linear combination of coils into channels
  sens.label = sens.label(selsens);
  sens.tra   = sens.tra(selsens,:);
  % subsequently remove the coils that do not contribute to any sensor output
  selcoil  = find(sum(sens.tra,1)~=0);
  sens.pnt = sens.pnt(selcoil,:);
  sens.ori = sens.ori(selcoil,:);
  sens.tra = sens.tra(:,selcoil);

  switch ft_voltype(vol)
    case 'infinite'
      % nothing to do

    case 'singlesphere'
      % nothing to do

    case 'concentric'
      % nothing to do

    case 'neuromag'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % if the forward model is computed using the external Neuromag toolbox,
      % we have to add a selection of the channels so that the channels
      % in the forward model correspond with those in the data.
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      [selchan, selsens] = match_str(channel, sens.label);
      vol.chansel = selsens;

    case 'multisphere'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % If the volume conduction model consists of multiple spheres then we
      % have to match the channels in the gradiometer array and the volume
      % conduction model.
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % use the original sensor array instead of the one with a subset of
      % channels, because we need the complete mapping of coils to channels
      sens = sens_orig;

      % remove the coils that do not contribute to any channel output
      % since these do not have a corresponding sphere
      selcoil  = find(sum(sens.tra,1)~=0);
      sens.pnt = sens.pnt(selcoil,:);
      sens.ori = sens.ori(selcoil,:);
      sens.tra = sens.tra(:,selcoil);

      % the initial multisphere volume conductor has a local sphere per
      % channel, whereas it should have a local sphere for each coil
      if size(vol.r,1)==size(sens.pnt,1) && ~isfield(vol, 'label')
        % it appears that each coil already has a sphere, which suggests
        % that the volume conductor already has been prepared to match the
        % sensor array
        return
      elseif size(vol.r,1)==size(sens.pnt,1) && isfield(vol, 'label')
        if ~isequal(vol.label(:), sens.label(:))
          % if only the order is different, it would be possible to reorder them
          error('the coils in the volume conduction model do not correspond to the sensor array');
        else
          % the coil-specific spheres in the volume conductor should not have a label
          % because the label is already specified for the coils in the
          % sensor array
          vol = rmfield(vol, 'label');
        end
        return
      end

      % the CTF way of representing the headmodel is one-sphere-per-channel
      % whereas the FieldTrip way of doing the forward computation is one-sphere-per-coil
      Nchans   = size(sens.tra,1);
      Ncoils   = size(sens.tra,2);
      Nspheres = size(vol.label);

      if isfield(vol, 'orig')
        % these are present in a CTF *.hdm file
        singlesphere.o(1,1) = vol.orig.MEG_Sphere.ORIGIN_X;
        singlesphere.o(1,2) = vol.orig.MEG_Sphere.ORIGIN_Y;
        singlesphere.o(1,3) = vol.orig.MEG_Sphere.ORIGIN_Z;
        singlesphere.r      = vol.orig.MEG_Sphere.RADIUS;
        % ensure consistent units
        singlesphere = ft_convert_units(singlesphere, vol.unit);
        % determine the channels that do not have a corresponding sphere
        % and use the globally fitted single sphere for those
        missing = setdiff(sens.label, vol.label);
        if ~isempty(missing)
          warning('using the global fitted single sphere for %d channels that do not have a local sphere', length(missing));
        end
        for i=1:length(missing)
          vol.label(end+1) = missing(i);
          vol.r(end+1,:)   = singlesphere.r;
          vol.o(end+1,:)   = singlesphere.o;
        end
      end

      multisphere = [];
      % for each coil in the MEG helmet, determine the corresponding channel and from that the corresponding local sphere
      for i=1:Ncoils
        coilindex = find(sens.tra(:,i)~=0); % to which channel does this coil belong
        if length(coilindex)>1
          % this indicates that there are multiple channels to which this coil contributes,
          % which happens if the sensor array represents a synthetic higher-order gradient.
          [dum, coilindex] = max(abs(sens.tra(:,i)));
        end

        coillabel = sens.label{coilindex};                    % what is the label of this channel
        chanindex = strmatch(coillabel, vol.label, 'exact');  % what is the index of this channel in the list of local spheres
        multisphere.r(i,:) = vol.r(chanindex);
        multisphere.o(i,:) = vol.o(chanindex,:);
      end
      vol = multisphere;

      % finally do the selection of channels and coils
      % order them according to the users specification
      [selchan, selsens] = match_str(channel, sens.label);

      % first only modify the linear combination of coils into channels
      sens.label = sens.label(selsens);
      sens.tra   = sens.tra(selsens,:);
      % subsequently remove the coils that do not contribute to any sensor output
      selcoil  = find(sum(sens.tra,1)~=0);
      sens.pnt = sens.pnt(selcoil,:);
      sens.ori = sens.ori(selcoil,:);
      sens.tra = sens.tra(:,selcoil);
      % make the same selection of coils in the multisphere model
      vol.r = vol.r(selcoil);
      vol.o = vol.o(selcoil,:);

    case 'nolte'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % if the forward model is computed using the code from Guido Nolte, we
      % have to initialize the volume model using the gradiometer coil
      % locations
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % compute the surface normals for each vertex point
      if ~isfield(vol.bnd, 'nrm')
        fprintf('computing surface normals\n');
        vol.bnd.nrm = normals(vol.bnd.pnt, vol.bnd.tri);
      end

      % estimate center and radius
      [center,radius] = fitsphere(vol.bnd.pnt);

      % initialize the forward calculation (only if gradiometer coils are available)
      if size(sens.pnt,1)>0
        vol.forwpar = meg_ini([vol.bnd.pnt vol.bnd.nrm], center', order, [sens.pnt sens.ori]);
      end

    case 'openmeeg'
        % nothing ?

    otherwise
      error('unsupported volume conductor model for MEG');
  end

elseif iseeg
  % select the desired channels from the electrode array
  % order them according to the users specification
  if ~isfield(sens, 'tra')
     Nchans    = length(sens.label);
     Ncontacts = length(sens.label);
  else
     Nchans    = size(sens.tra,1);
     Ncontacts = size(sens.tra,2);
  end;

  % In case of Nchans~=Ncontacts it is difficult to determine
  % how to deal with contacts positions (keep the original positions)
  if Nchans == Ncontacts
    [selchan, selsens] = match_str(channel, sens.label);
    sens.label = sens.label(selsens);
    sens.pnt = sens.pnt(selsens,:);
  else
    warning('A sub-selection of channels will not be taken into account')
  end

  % create a 2D projection and triangulation
  try
    sens.prj   = elproj(sens.pnt);
    sens.tri   = delaunay(sens.prj(:,1), sens.prj(:,2));
  catch
    warning('2D projection not done')
  end

  switch ft_voltype(vol)
    case 'infinite'
      % nothing to do

    case 'halfspace'
      pnt    = sens.pnt;
      if ft_voltype(vol,'halfspace')
        d = dist(pnt);
        % scan the electrodes and reposition the ones which are in the
        % wrong halfspace (projected on the plane)
        for i=1:size(pnt,1)
          P = pnt(i,:);
          is_in_empty = acos(dot(vol.ori,(P-vol.pnt)./norm(P-vol.pnt))) < pi/2;
          if is_in_empty
            d = dist(P);
            dPplane = -dot(vol.ori, vol.pnt-P, 2);
            if dPplane>median(d(:))
              error('Some electrodes are too distant from the plane: consider repositioning them')
            else
              % project point on plane
              Ppr  = [0 0 0];
              line = [P vol.ori];
              % get indices of line and plane which are parallel
              par = abs(dot(vol.ori, line(:,4:6), 2))<1e-14;
              % difference between origins of plane and line
              dp = vol.pnt - line(:, 1:3);
              % Divide only for non parallel vectors (DL)
              t = dot(vol.ori(~par,:), dp(~par,:), 2)./dot(vol.ori(~par,:), line(~par,4:6), 2);
              % compute coord of intersection point
              Ppr(~par, :) = line(~par,1:3) + repmat(t,1,3).*line(~par,4:6);
              pnt(i,:) = Ppr;
            end
          end
        end
        sens.pnt = pnt;
      else
        error('Wrong volume type')
      end

    case {'singlesphere', 'concentric'}
      % ensure that the electrodes ly on the skin surface
      radius = max(vol.r);
      pnt    = sens.pnt;
      if isfield(vol, 'o')
        % shift the the centre of the sphere to the origin
        pnt(:,1) = pnt(:,1) - vol.o(1);
        pnt(:,2) = pnt(:,2) - vol.o(2);
        pnt(:,3) = pnt(:,3) - vol.o(3);
      end
      distance = sqrt(sum(pnt.^2,2)); % to the center of the sphere
      if any((abs(distance-radius)/radius)>0.005)
        warning('electrodes do not lie on skin surface -> using radial projection')
      end
      pnt = pnt * radius ./ [distance distance distance];
      if isfield(vol, 'o')
        % shift the center back to the original location
        pnt(:,1) = pnt(:,1) + vol.o(1);
        pnt(:,2) = pnt(:,2) + vol.o(2);
        pnt(:,3) = pnt(:,3) + vol.o(3);
      end
      sens.pnt = pnt;

    case {'bem', 'dipoli', 'asa', 'avo', 'bemcp', 'openmeeg'}
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % do postprocessing of volume and electrodes in case of BEM model
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % project the electrodes on the skin and determine the bilinear interpolation matrix
      if ~isfield(vol, 'tra')
        % determine boundary corresponding with skin and inner_skull_surface
        if ~isfield(vol, 'skin_surface')
          vol.skin_surface = find_outermost_boundary(vol.bnd);
          fprintf('determining skin compartment (%d)\n', vol.skin_surface);
        end
        if ~isfield(vol, 'source')
          vol.source = find_innermost_boundary(vol.bnd);
          fprintf('determining source compartment (%d)\n', vol.source);
        end
        if size(vol.mat,1)~=size(vol.mat,2) && size(vol.mat,1)==length(sens.pnt)
          fprintf('electrode transfer and system matrix were already combined\n');
        else
          fprintf('projecting electrodes on skin surface\n');
          % compute linear interpolation from triangle vertices towards electrodes
          [el, prj] = project_elec(sens.pnt, vol.bnd(vol.skin_surface).pnt, vol.bnd(vol.skin_surface).tri);
          tra       = transfer_elec(vol.bnd(vol.skin_surface).pnt, vol.bnd(vol.skin_surface).tri, el);

          % replace the original electrode positions by the projected positions
          sens.pnt = prj;

          if size(vol.mat,1)==size(vol.bnd(vol.skin_surface).pnt,1)
            % construct the transfer from only the skin vertices towards electrodes
            interp = tra;
          else
            % construct the transfer from all vertices (also inner_skull_surface/outer_skull_surface) towards electrodes
            interp = [];
            for i=1:length(vol.bnd)
              if i==vol.skin_surface
                interp = [interp, tra];
              else
                interp = [interp, zeros(size(el,1), size(vol.bnd(i).pnt,1))];
              end
            end
          end

          % incorporate the linear interpolation matrix and the system matrix into one matrix
          % this speeds up the subsequent repeated leadfield computations
          fprintf('combining electrode transfer and system matrix\n');
          if strcmp(ft_voltype(vol), 'openmeeg')
            nb_points_external_surface = size(vol.bnd(vol.skin_surface).pnt,1);
            vol.mat = vol.mat((end-nb_points_external_surface+1):end,:);
            vol.mat = interp(:,1:nb_points_external_surface) * vol.mat;
          else
            % convert to sparse matrix to speed up the subsequent multiplication
            interp  = sparse(interp);
            vol.mat = interp * vol.mat;
            % ensure that the model potential will be average referenced
            avg = mean(vol.mat, 1);
            vol.mat = vol.mat - repmat(avg, size(vol.mat,1), 1);
          end
        end
      end

    otherwise
      error('unsupported volume conductor model for EEG');
  end

  % FIXME this needs carefull thought to ensure that the average referencing which is now done here and there, and that the linear interpolation in case of BEM are all dealt with consistently
  % % always ensure that there is a linear transfer matrix for
  % % rereferencing the EEG potential
  % if ~isfield(sens, 'tra');
  %   sens.tra = sparse(eye(length(sens.label)));
  % end

end % if iseeg or ismeg

author: PatternRecognition

repo: OpenBMI


text:
function X_zero=assignPerimeterPositions(radius,mesh);
% Define a set of perimeter points on the circle.
%
%   X_zero=assignPerimeterPositions(perimeterDists,mesh);
%
% orderMeshPerimeterPointsAll is the routine that gets the perimeter points
% into the proper order for the mapping here.   The mapping here simply
% puts these points on a circle.  In principle, we could put the points at
% these angles but adjust for distance.  In that case, however, we might
% not have a convex set.  So, we use a circle.
%
% The most anterior point in the sub mesh should map to the top of the
% circle. The most posterior point should be on the left or right.
%
% Stanford.

numPerimPoints=length(mesh.orderedUniquePerimeterPoints);

% Find the maximum y value in the perimeter
[maxHeight,maxHIndex]=max(mesh.uniqueVertices(mesh.orderedUniquePerimeterPoints,2));

angles=linspace(0,2*pi,numPerimPoints)';
angles=shift(angles(:),[-maxHIndex,0]);

X_zero=zeros(numPerimPoints,2);

% Perimeter must be convex for Tutte's algorithm  to work properly -
X_zero(:,1)=radius.*cos(angles);
X_zero(:,2)=radius.*sin(angles);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ p, p2 ] = plot_diffusion_profiles( S, dash )
%UNTITLED5 Summary of this function goes here
%   Detailed explanation goes here

x = 0:1E-5:3E-3;

if nargin == 1
    dash = 0;
end

if dash
    dd = '--';
else
    dd = '-';
end

hold on;

if isfield(S, 'strexp')
    u = S.strexp.D;
    beta = S.strexp.gamma;

    k = 0:5000;

    tww = u^(1/beta);
    l = zeros(size(x));
    for t_counter = 1:length(x)

        t = x(t_counter);

        % see wikipedia page on stretched exponential function
        %tosum = ((-1).^k)./factorial(k).*sin(pi*beta*k).*gamma(beta*k+1).*(t/tww).^(beta*k+1);
        tosum = ((-1).^k).*sin(pi*beta*k).*exp(gammaln(beta*k + 1) - gammaln(k + 1) + (beta*k+1)*log(t/tww));
        tosum(isnan(tosum)) = 0;

        l(t_counter) = -tww/(pi*t)*sum(tosum);
    end
    Y = l./x;
    Y(x < 0) = 0;
    Y(Y <= 0) = 0;
    Y(x > 2.5E-3) = 0; % to avoid numerical errors
    plot(x*1000, Y, dd, 'Color', [0.5 0.5 0.5], 'LineWidth', 2)

end

if isfield(S, 'DKI')
    D = S.DKI.D;
    K = S.DKI.K;

    sigma = sqrt(1/3 * D^2 * K);
    Y = normpdf(x, D, sigma);

    plot(x*1000, Y, [dd 'g'], 'LineWidth', 2);

end

if isfield(S, 'stat')
    mu = S.stat.mu;
    sigma = S.stat.sigma;
    if isfield(S.stat, 'Dmin')
        Dmin = S.stat.Dmin;
    else
        Dmin = 0;
    end
    if isfield(S.stat, 'Dmax')
        Dmax = S.stat.Dmax;
    else
        Dmax = inf;
    end

    Y = normpdf(x, mu, sigma) ...
      / (normcdf(Dmax, mu, sigma) - normcdf(Dmin, mu, sigma));
    Y(x < Dmin) = 0;
    Y(x > Dmax) = 0;

    plot(x*1000, Y, [dd 'k'], 'LineWidth', 2);

end

if isfield(S, 'gamma')
%     D = S.gamma.D;
%     sigma = S.gamma.sigma;
%     kappa = D^2 / sigma^2;
%     theta = sigma^2 / D;

    kappa = S.gamma.kappa;
    theta = S.gamma.theta;

    Ygamma = x.^(kappa - 1) .* exp(- x / theta) ./ (gamma(kappa) * theta^kappa);
    Ygamma(x < 0) = 0;

    plot(x*1000, Ygamma, dd, 'Color', [1 0.69 0.39], 'LineWidth', 2);

end

if 0%isfield(S, 'doubletrunc')

    mu = S.doubletrunc.mu;
    sigma = S.doubletrunc.sigma;
    D_upper = S.doubletrunc.D_upper;

    Y = normpdf(x, mu, sigma) / (normcdf(D_upper, mu, sigma) - normcdf(0, mu, sigma));
    Y(x < 0) = 0;
    Y(x > D_upper) = 0;

    plot(x*1000, Y, 'k--', 'LineWidth', 2);

end

if isfield(S, 'beta')
    a = S.beta.alpha;
    b = S.beta.beta;
    D0 = S.beta.D0;

    P = 1/D0 * gamma(a + b) / (gamma(a) * gamma(b)) *  (x/D0).^(a - 1) .* (1 - x / D0).^(b - 1);
    P(x < 0) = 0;
    P(x > D0) = 0;

    p = plot(x*1000, real(P), [dd 'r'], 'LineWidth', 2);
end

if isfield(S, 'monoexp')
    D = S.monoexp;
    plot([D D]*1000, [0 1500], dd, 'LineWidth', 2);

end

if isfield(S, 'biexp')
    Dfast = S.biexp.Dfast;
    Dslow = S.biexp.Dslow;
    f = S.biexp.f;
    dd = '-';
    p2 = plot([Dfast Dfast NaN Dslow Dslow]*1000, [0 1500*f NaN 0 1500*(1-f)], [dd 'c'], 'LineWidth', 2);

end

%legend('Mono-exponential', 'Bi-exponential', 'DKI (Gaussian)', 'Truncated Gaussian', 'Gamma', 'Beta')

set(gca, 'FontSize', 30)
set(gca, 'LineWidth', 2)
set(gca, 'XLim', [0 3])

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function treepack_test04 ( )

%*****************************************************************************80
%
%% TREEPACK_TEST04 tests TREE_ARC_CENTER.
%
%  Discussion:
%
%    The tree is
%
%    2---3---6---8---1---9
%       /       / \
%      7       5   4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
  nnode = 9;

  inode = [ 2, 3, 3, 6, 8, 8, 8, 1 ];
  jnode = [ 3, 7, 6, 8, 4, 5, 1, 9 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TREEPACK_TEST04\n' );
  fprintf ( 1, '  TREE_ARC_CENTER computes the center of a tree.\n' );

  graph_arc_print ( nnode-1, inode, jnode, '  The edge list of the tree:' );

  [ center, eccent, parity ] = tree_arc_center ( nnode, inode, jnode );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Parity = %d\n', parity );
  fprintf ( 1, '  Eccentricity is %d\n', eccent );

  if ( parity == 0 )
    fprintf ( 1, '  No center node (degenerate case).\n' );
  elseif ( parity == 1 )
    fprintf ( 1, '  Center node: %d\n', center(1) );
  else
   fprintf ( 1, '  Center nodes: %d %d\n', center(1), center(2) );
  end

  return
end


text:
% Genetic Algorithm Toolbox.
% Version 1.2 15-Apr-94
% Department of Automatic Control and Systems Engineering
% University of Sheffield, England
%
% Creating populations
%   crtbase	- create a base vector
%   crtbp	- create a binary population
%   crtrp	- create a real-valued population
%
% Fitness assignment
%   ranking	- rank-based fitness assignment
%   scaling	- proportional fitness-scaling
%
% Selection and reinsertion
%   reins	- uniform random and fitness-based reinsertion
%   rws		- roulette wheel selection
%   select	- high-level selection routine
%   sues		- stochastic universal sampling
%
% Mutation operators
%   mut		- discrete mutation
%   mutate	- high-level mutation function
%   mutbga	- real-value mutation
%
% Crossover operators
%   recdis	- discrete recombination
%   recint	- intermediate recombination
%   reclin	- line recombination
%   recmut	- line recombination with mutation features
%   recombin	- high-level recombination function
%   xovdp	- double-point crossover
%   xovdprs	- double-point reduced surrogate crossover
%   xovmp	- general multi-point crossover
%   xovsh	- shuffle crossover
%   xovshrs	- shuffle reduced surrogate crossover
%   xovsp	- single-point crossover
%   xovsprs	- single-point reduced surrogate crossover
%
% Subpopulation support
%   migrate	- exchange individuals between subpopulations
%
%
% Utility functions
%   bs2rv	- binary string to real-value conversion
%   rep		- matrix replication
%
% Demonstration and other functions
%   mpga	- multi-population genetic algorithm demonstration
%   objfun1	- De Jongs first test function (used by sga)
%   objharv	- harvest function (used in mpga)
%   resplot	- result plotting (used in mpga)
%   sga		- simple genetic algorithm demonstration

author: vonsylvia


text:
function x = p07_root ( i )

%*****************************************************************************80
%
%% P07_ROOT returns a root for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the requested root.
%
%    Output, real X, the value of the root.
%
  if ( i == 1 )
    x = 0.0;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P07_ROOT - Fatal error!\n' );
    fprintf ( 1, '  Illegal root index = %d\n', i );
    error ( 'P07_ROOT - Fatal error!' );
  end

  return
end


text:
function [abbox_scores] = recognize_bboxes_of_all_imgs(...
    model, image_paths, all_bbox_proposals, dst_directory, image_set_name, varargin)
% recognize_bboxes_of_all_imgs: given a bounding box recognition model and
% a set of images with their input candidate bounding boxes, for each
% image it assigns a confidence score to each candidate bounding box that
% represents the likelihood to tightly enclose an object of interest.
%
% INPUTS:
% 1) model: (type struct) the bounding box recognition model
% 2) image_paths: NI x 1 cell array with the image paths of the set in the
% form of strings. NI is the number of images.
% 3) all_bbox_proposals: a cell array with the candidate bounding
% boxes that are going to be processed. It can be of two forms:
%   a) (the flag is_per_class must being set to false) a NI x 1 cell array
%   where NI is the number of images. The all_bbox_proposals{i} element
%   contains the candidate bounding boxes of the i-th image in the form of a
%   NB_i x 4 array where NB_i is the number of candidate bounding boxes of
%   the i-th image. The 4 columns of this array contain the bounding box
%   coordinates in the form of [x0,y0,x1,y1] (where the (x0,y0) and (x1,y1)
%   are the top-left and bottom-right corners).
%   b) (the flag is_per_class must being set to true) a C x 1 cell array
%   where C is the number of categories. The all_bbox_proposals{j} element
%   in this case is NI x 1 cell array, where NI is the number of images,
%   with the candidate bounding boxes of the j-th category for each image.
%   The element all_bbox_proposals{j}{i} is a NB_{i,j} x 4 array, where
%   NB_{i,j} is the number candidate bounding box of the i-th
%   image for the j-th category. The 4 columns of this array are the
%   bounding boxes in the form of [x0,y0,x1,y1] (where the (x0,y0) and (x1,y1)
%   are the top-left and bottom-right corners).
% 4) dst_directory: string with the path of the destination directory where
% the results are going to be cached.
% 5) image_set_name: string with the name of the image set that is being
% processed. It will be used on the name of the file where the results will
% be saved.
% 6) The rest input arguments are given in the form of Name,Value pair
% arguments and are:
% 'is_per_class': boolean value that if set to false, then the 3.a) form
% of the all_bbox_proposals input parameter will be expected; otherwise, if set to
% true then the 3.b) form of the all_bbox_proposals input parameter will be
% expected. Default value: false
% 'suffix': string with a suffix that it will be used on the name of the
% file where the results will be saved. Default value: ''
% 'all_bbox_gt': (optional) a NI cell array with the ground truth bounding
% boxes of each image; it will be used (if given) for printing the mAP on
% regular (see param 'checkpoint_step') step during processing the image set
% 'checkpoint_step': scalar value; during processing the NI images, after each
% checkpoint_step number of images the results are cached and the mAP
% results are printed. Default value: 500
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

ip = inputParser;
ip.addParamValue('is_per_class', false, @islogical);
ip.addParamValue('suffix', '', @ischar);
ip.addParamValue('all_bbox_gt', {},     @iscell);
ip.addParamValue('checkpoint_step', 500, @isnumeric);

ip.parse(varargin{:});
opts = ip.Results;

mkdir_if_missing(dst_directory);
filepath_results     = [dst_directory, filesep, 'scores', '_boxes_', opts.suffix, image_set_name, '.mat'];
in_progress_filepath = [dst_directory, filesep, 'scores', '_boxes_', opts.suffix, image_set_name, '_in_progress.mat'];

timestamp       = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
log_file        = fullfile(dst_directory, 'output', ...
    ['log_file_recognition',opts.suffix, image_set_name, '_', timestamp, '.txt']);
mkdir_if_missing(fileparts(log_file));

t_start = tic();
try
    abbox_scores = load_bboxes_scores(filepath_results);
catch
    diary(log_file);
    abbox_scores = process_images(...
        model, image_paths, all_bbox_proposals, opts.all_bbox_gt, ...
        opts.is_per_class, in_progress_filepath, opts.checkpoint_step);
    save_bboxes_scores(filepath_results, abbox_scores);
    diary off;
    delete(in_progress_filepath);
end
fprintf('Recognize bounding box proposals in %.4f minutes.\n', toc(t_start)/60);
end

function save_bboxes_scores(filename, bboxes_scores)
save(filename, 'bboxes_scores', '-v7.3');
end

function bboxes_scores = load_bboxes_scores(filename)
load(filename, 'bboxes_scores');
end

function all_bboxes_out = process_images(...
    model, image_paths, all_bbox_proposals, all_bbox_gt, ...
    is_per_class, in_progress_filepath, checkpoint_step)

num_imgs      = length(image_paths);
nms_iou_thrs  = 0.3;
max_per_set   = 5 * num_imgs;
max_per_image = 100;
num_classes   = length(model.classes);
in_progress_filepath_prev = [in_progress_filepath, '.prev'];

try
    % load (if any) cached subset of results from images that have being
    % already processed from a previous run of the script that was stopped
    % before finish processing all the images.
    try
        [fist_img_idx, aboxes, all_bboxes_out, thresh] = load_progress(in_progress_filepath);
    catch
        [fist_img_idx, aboxes, all_bboxes_out, thresh] = load_progress(in_progress_filepath_prev);
    end
    mAP   = check_progress_on_mAP(aboxes, all_bbox_gt, fist_img_idx-1, model.classes);
    mAP_i = fist_img_idx-1;
catch exception
    fprintf('Exception message %s\n', getReport(exception));
    aboxes = cell(num_classes, 1);

    for i = 1:num_classes, aboxes{i} = cell(num_imgs, 1); end
    if is_per_class
        all_bboxes_out = cell(num_classes,1);
        for i = 1:num_classes, all_bboxes_out{i} = cell(num_imgs, 1); end
    else
        all_bboxes_out = cell(num_imgs,1);
    end
    thresh = -3.5 * ones(num_classes, 1);
    fist_img_idx = 1; mAP = 0; mAP_i = 0;
end

total_el_time = 0;
model.max_rois_num_in_gpu = 1000; %find_max_rois_num_in_gpu(model);
num_chars = 0;
for i = fist_img_idx:num_imgs
    th = tic;
    % get the bounding box proposals of this image
    bbox_proposals = get_this_img_bbox_proposals(all_bbox_proposals, i, is_per_class);
    image          = get_image(image_paths{i}); % read the image
    % score the bounding box proposals with the recognition model
    bboxes_scores  = recognize_bboxes_of_image(model, image, bbox_proposals(:,1:4));
    recognition_time = toc(th);

    th = tic;
    % prepare the candidate boundng box detections of this image
    bbox_cand_dets = prepare_bbox_cand_dets(bbox_proposals, bboxes_scores, num_classes, is_per_class);
    % prepare the output for this image
    all_bboxes_out = prepare_this_img_output(all_bboxes_out, i, is_per_class, bbox_cand_dets);
    % perform the post-processing step of non-max-suppresion on the
    % candidate boundng box detections of this image
    bbox_detections = post_process_candidate_detections(bbox_cand_dets, ...
        'thresholds',thresh, 'nms_iou_thrs',nms_iou_thrs,'use_gpu',true,...
        'max_per_image',max_per_image);

    for j = 1:num_classes, aboxes{j}{i} = bbox_detections{j}; end
    postprocessing_time = toc(th);

    if mod(i, checkpoint_step) == 0
        th = tic;
        for j = 1:num_classes, [aboxes{j}, thresh(j)] = keep_top_k(aboxes{j}, i, max_per_set, thresh(j)); end
        disp(thresh(:)');

        % save the till now progress
        save_progress(aboxes, all_bboxes_out, thresh, i, in_progress_filepath, in_progress_filepath_prev);

        % evaluate the mAP of the images 1 till i
        mAP = check_progress_on_mAP(aboxes, all_bbox_gt, i, model.classes);
        mAP_i = i;
        diary; diary; % flush diary

        checkpoint_time = toc(th);
        fprintf('check_point time %.2fs\n',checkpoint_time)
        num_chars = 0;
    end

    elapsed_time = recognition_time + postprocessing_time;
    [total_el_time, avg_time, est_rem_time] = timing_process(elapsed_time, total_el_time, fist_img_idx, i, num_imgs);
    fprintf(repmat('\b',[1, num_chars]));

    num_chars = fprintf('%s: bbox rec %d/%d:| ET %.3fs + %.3fs | AT: %.3fs | TT %.4fmin | ERT %.4fmin | mAP[%d/%d] = %.4f\n', ...
        procid(), i, num_imgs, recognition_time, postprocessing_time, avg_time, ...
        total_el_time/60, est_rem_time/60, mAP_i, num_imgs, mAP);

end

delete(in_progress_filepath_prev);
end

function [total_el_time, ave_time, est_rem_time] = timing_process(...
    elapsed_time, total_el_time, fist_img_idx, i, num_imgs)

total_el_time   = total_el_time + elapsed_time;
ave_time        = total_el_time / (i-fist_img_idx+1);
est_rem_time    = ave_time * (num_imgs - i);
end

function bbox_proposals = get_this_img_bbox_proposals(all_bboxes_in, img_id, is_per_class)
if is_per_class
    bbox_proposals_per_class = cellfun(@(x) x{img_id}(:,1:4), all_bboxes_in, 'UniformOutput', false);
    num_bbox_per_class = cellfun(@(x) size(x,1), bbox_proposals_per_class,  'UniformOutput', true);
    bbox_proposals = cell2mat(bbox_proposals_per_class(num_bbox_per_class>0));

    class_indices = single([]);
    for c = 1:length(num_bbox_per_class)
        class_indices = [class_indices; ones(num_bbox_per_class(c),1,'single')*c];
    end

    bbox_proposals = single([bbox_proposals, class_indices]);
    if isempty(bbox_proposals), bbox_proposals = zeros(0,5,'single'); end
else
    bbox_proposals = all_bboxes_in{img_id}(:,1:4);
end
end

function bbox_cand_dets = prepare_bbox_cand_dets(bbox_proposals, bbox_scores, num_classes, is_per_class)
if is_per_class
    class_indices  = bbox_proposals(:,5);
    bbox_cand_dets = cell(num_classes,1);
    for c = 1:num_classes
        this_cls_mask = class_indices==c;
        bbox_cand_dets{c} = single([bbox_proposals(this_cls_mask,1:4),bbox_scores(this_cls_mask,c)]);
        if isempty(bbox_cand_dets{c}), bbox_cand_dets{c} = zeros(0,5,'single'); end
    end
else
    bbox_cand_dets = single([bbox_proposals(:,1:4), bbox_scores]);
end
end

function all_bboxes_out = prepare_this_img_output(all_bboxes_out, img_idx, is_per_class, bbox_this_img)
if is_per_class
    for j = 1:length(bbox_this_img), all_bboxes_out{j}{img_idx} = bbox_this_img{j}; end
else
    all_bboxes_out{img_idx} = bbox_this_img;
end
end

function [fist_img_idx, aboxes, abbox_scores, thresh] = load_progress(in_progress_filepath)
ld = load(in_progress_filepath);
fist_img_idx = ld.progress_state.img_idx + 1;
aboxes       = ld.progress_state.aboxes;
abbox_scores = ld.progress_state.abbox_scores;
thresh       = ld.progress_state.thresh;
end

function save_progress(aboxes, abbox_scores, thresh, img_idx, in_progress_filepath, in_progress_filepath_prev)
progress_state              = struct;
progress_state.img_idx      = img_idx;
progress_state.aboxes       = aboxes;
progress_state.abbox_scores = abbox_scores;
progress_state.thresh       = thresh;

if exist(in_progress_filepath, 'file')
    % in case it crash during updating the in_progress_filepath file
    copyfile(in_progress_filepath, in_progress_filepath_prev);
end
save(in_progress_filepath, 'progress_state', '-v7.3');
end

function mAP = check_progress_on_mAP(aboxes, all_bbox_gt, img_idx, classes)
mAP = 0;
if ~isempty(all_bbox_gt)
    aboxes      = cellfun(@(x) x(1:img_idx), aboxes, 'UniformOutput', false);
    mAP_result  = evaluate_average_precision_pascal( all_bbox_gt(1:img_idx), aboxes, classes );
    printAPResults(classes, mAP_result);
    mAP = mean([mAP_result(:).ap]');
end
end

function [boxes, thresh] = keep_top_k(boxes, end_at, top_k, thresh)
% ------------------------------------------------------------------------
% Keep top K
X = cat(1, boxes{1:end_at});
if isempty(X), return; end

scores = sort(X(:,end), 'descend');
thresh = scores(min(length(scores), top_k));
for image_index = 1:end_at
    bbox = boxes{image_index};
    keep = find(bbox(:,end) >= thresh);
    boxes{image_index} = bbox(keep,:);
end

end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet