text:
function h = entropy(vec1)
%=========================================================
%
%This is a prog in the MutualInfo 0.9 package written by
% Hanchuan Peng.
%
%Disclaimer: The author of program is Hanchuan Peng
%      at <penghanchuan@yahoo.com> and <phc@cbmv.jhu.edu>.
%
%The CopyRight is reserved by the author.
%
%Last modification: April/19/2002
%
%========================================================
%
% h = entropy(vec1)
% calculate the entropy of a variable vec1
%
% demo:
%  a=[1 2 1 2 1]';b=[2 1 2 1 1]';
%  fprintf('entropy(a) = %d\n',entropy(a));
%
% the same as entropycond(vec1)
%
% By Hanchuan Peng, April/2002
%

if nargin<1,

  disp('Usage: h = entropy(vec1).');
  h = -1;

else,

  [p1] = estpa(vec1);
  h = estentropy(p1);

end;

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function [ c, ierror ] = r8mat_cholesky_factor ( n, a )

%*****************************************************************************80
%
%% R8MAT_CHOLESKY_FACTOR computes the Cholesky factor of a symmetric matrix.
%
%  Discussion:
%
%    The matrix must be symmetric and positive semidefinite.
%
%    For a positive semidefinite symmetric matrix A, the Cholesky factorization
%    is a lower triangular matrix L such that:
%
%      A = L * L'
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the matrix.
%
%    Output, real C(N,N), the lower triangular Cholesky factor.
%
%    Output, integer IERROR, error flag.
%    0, no error.
%    1, warning, the matrix is positive semidefinite.  The factorization
%    was carried out, but the matrix is singular.
%    2, error, the matrix has at least one negative eigenvalue.  The
%    factorization could not be completed.
%
  ierror = 0;

  c(1:n,1:n) = a(1:n,1:n);

  for j = 1 : n

    c(1:j-1,j) = 0.0;

    for i = j : n

      ctc = c(j,i) - c(i,1:j-1) c(j,1:j-1)';

      if ( i == j )
        if ( ctc < 0.0 )
          ierror = 2;
          return
        elseif ( ctc == 0.0 )
          ierror = 1;
        else
          c(i,j) = sqrt ( ctc );
        end
      else
        if ( c(j,j) ~= 0.0 )
          c(i,j) = ctc / c(j,j);
        else
          c(i,j) = 0.0;
        end
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out = Dstar(R, n, e, P, fun)
%finds excess distortion for a given rate at a given blocklength for
%transmission of a GMS over an AWGN (noise variance 1)
%R - rate
%n - block length (scalar)
%e - excess probability
%P - channel SNR
%fun - which function to use for calculation

%
%   Created in 2012 by Victoria Kostina (vkostina@caltech.edu)
%

%starting points - 'persistent' to make optimization faster
persistent x0;
persistent y0;
persistent nprev;

if isempty(nprev)
    nprev = 0;
end

tol = 1e-12;

options = optimset('TolX',tol, 'MaxFunEvals', 500, 'Display', 'off');
dbar = 1/(1+P);
V1 = 2*dbar^2; %no coding (amplifiers)
V = 2*dbar^2*(2 - dbar^2); %optimal blocklength-n code

Vc = P*(P+2)/2/(P+1)^2;
Vs = .5;

deltan = n - nprev; %distance between consecutive n's
x0correction = sqrt(V1/n)*Qinv(e)/2/n^(3/2)*deltan;
y0slack = .1;

k = ceil(R*n);

if ~isempty(x0)
    x0 = x0 - x0correction;
    d0 = x0;
else
    d0 = 1.8*dbar;
end

switch lower(fun)
    case 'aseparate'
        out = ASeparate();
    case 'ajoint'
        %Achievability via JSCC coding
        options = optimset('TolX',tol, 'MaxFunEvals', 500, 'Algorithm', 'active-set', 'Display', 'off');
        out = AJoint();
    case 'cht'
        %Converse via hypothesis testing
        out = Cht();
    case 'nocoding'
        %no coding
        out = ANoCoding();
    case 'approxnocoding'
        out = ApproxNoCoding();
    case 'approx'
        %gaussian approximation
        out = Approx();
    otherwise
        disp('Error: unknown type.')
end

%--------------------------------------------------------------------------
    function out = ANoCoding()
        %No coding

        out = fzero(@(d)1 - e - chi2cdf(n*d/dbar, n), d0);
    end

%--------------------------------------------------------------------------
    function out = ApproxNoCoding()
        %No coding
        out = dbar + sqrt(V1/n)*Qinv(e);
    end
%--------------------------------------------------------------------------
    function out = Approx()
        %No coding
        %C = .5*log(1 + P);
       % Rd = .5*log(1/d);
        %out = fzero(@(d)C + .5*log(d) - sqrt(1/n*(Vc + Vs))*Qinv(e), d0);
        out = dbar + sqrt(V/n)*Qinv(e);
    end

%--------------------------------------------------------------------------
    function out = Cht()
        %Converse via meta-converse

        tau = thres(k);

        out = fzero(@(d)Qrob(d, tau) - 1, d0);
        updatex0(out);

        function out = Qrob(d, tau)

            b = sqrt( tau*n/k/d );
            out = quad(@density, 0, b);
            function out = density(r)
                % v - noncentral chi square (n) /n
                % v0 - central chi square (k) /n

                out = ncx2cdf(n*(tau - k/n*d*r.^2)/P, n, n*(1 + 1/P)).*r.^(k-1).*k;

            end
        end

        function out = thres(k)
            %find tau:
            tau0 = 1 + k/n + 2*sqrt( (n*Vc + k*Vs)/n^2 )*Qinv(e);
            %out = fzero(@(tau)Prob(tau) - 1 + e, tau0, options);
            out = fzero(@(tau)Prob(tau) - 1 + e, tau0, options);

            function out = Prob(tau)

                if n < 200
                    b = 10*k/n; %Ev0 = k/n
                else
                    b = 2*k/n;
                end

                out = quad(@density, 0, b);
                function out = density(v0)
                    % v - noncentral chi square (n) /n
                    % v0 - central chi square (k) /n

                    out = ncx2cdf(n*(tau - v0)*(1+P)/P, n, n/P).*chi2pdf(n*v0,k)*n;

                end
            end
        end
    end

%--------------------------------------------------------------------------
    function out = AJoint()
        %dP/dQ is bounded by a constant. Find that constant.
        [~, gamma2] = fminbnd(@(x)-ncx2pdf(n*x, n, n*P)./chi2pdf(n*x/(1+P), n)*(1+P), 0, 10);
        gamma2 = -gamma2;

        out = fzero(@(d)bound(d) - e, d0);
        updatex0(out);
%        fprintf('AJoint: n = %i, R = %f \n', n, out);

        function out = bound(d)
            %optimize with respect to r and g; y0 = [r gg] - TOO SLOW
            if d > 1
                out = NaN;
                return;
            end
            out = goal([d 3]);

            function out = goal(dg)
                d = dg(1);

                gamma1 = log(n)/2*dg(2);

                r0 = sqrt(1 - d);
                a = max(0, r0 - sqrt(d));
                b = r0 + sqrt(d);

                A = n/2*log(1+P) + n/2 - log(gamma1) - log(gamma2);

                Fouter = @(v0)chi2pdf(k*v0, k).*k.*Iinner(v0);
                Ftail = @(v0)(1 - ncx2cdf(thresv(v0)*n,n,n/P)).*chi2pdf(k*v0,k).*k;
                %   Ftail = @(v0) exp(n/P*t./(1-2*t))./(1-2*t).^(n/2).*exp(-t*thresv(v0)).*chi2pdf(k*v0,k).*k;
                out = quad(Fouter, a^2 + tol, b^2- tol)...                      %exponent < 1
                    + quad(Ftail, a^2 + tol, b^2 - tol)...                      %exponent = 1
                    + chi2cdf((a^2 + tol)*k, k) + 1 - chi2cdf((b^2 - tol)*k, k)...      %nontypical source realization
                    + exp(1-gamma1);

                %out = quad(Ftail, a^2 + tol, b^2 - tol); %most significant
                %term

                function out = Iinner(v0)
                    %exponent < 1
                    out = zeros(1, numel(v0));
                    for j0 = 1:numel(v0)
                        Finner = @(v) ncx2pdf(n*v, n, n/P).*n.*exp(P/2/(1+P)*n*(v - thresv(v0(j0))));
                        out(j0) = quad(Finner, 0, thresv(v0(j0)) );
                    end
                    if any(isnan(out))
                        return;
                    end
                end

                function out = thresv(x)
                    %threshold for v as a function of v0
                    out = max( 0, 2*(1+P)/P/n*(A + logPdball(x)) );
                end

                function out = logPdball(x)
                    %x - distance squared from the origin (x*k central chi square k)
                    %lower bound

                    cosa = (r0^2 + x - d)./(2*r0.*sqrt(x));
                    sina = (1 - cosa.^2).^(1/2);
                    out = loggamma(k/2+1, 'LB') - loggamma((k-1)/2 +1, 'UB')- .5*log(pi)-log(k) + (k-1)*log(sina);

                    %handle small distances from the origin
                    out(x <= sqrt(d) - r0) = -Inf;
                end
            end
        end
    end

%--------------------------------------------------------------------------
    function out = ASeparate()
        %achievability via separate coding
        %R = 1

        if isempty(y0)
            %first pass
            startL = 1;
            endL = n;
        else
            startL = max( ceil((y0 - 4*y0slack)*n), 1);
            endL = min( ceil((y0 + 4*y0slack)*n), n );
        end

        %precompute channel error probability
        PeCh = ones(1, endL);
        for L = startL:endL %exp(r) - the number of messages passed from the source enc to the channel enc
            PeCh(L) = PeubChannel(P, n, L);
        end

        Peubseparate(k, d0, PeCh, e);
        out = fzero(@(d)Peubseparate(k, d, PeCh, e) - e, d0);
        updatex0(out);
        [err, curL] = Peubseparate(k, out, PeCh, e);
        updatey0(curL/n);
    end

%--------------------------------------------------------------------------
%Utilities for faster computation
    function updatex0(x)
        if ~isnan(x) && n > 100
            x0 = x;
            nprev = n;
        end
    end

%--------------------------------------------------------------------------
%Utilities for faster computation
    function updatey0(y)
        if ~isnan(y) && n > 100
            y0 = y;
        end
    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function all_bbox_proposals = extract_selective_search_boxes_from_dataset(image_db, ss_boxes_dst_file)
%
% 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]
% ---------------------------------------------------------

try
    ld = load(ss_boxes_dst_file);
    all_bbox_proposals = ld.all_bbox_proposals;
catch
    chunk_size = 1000;
    num_imgs   = numel(image_db.image_paths);
    num_chunks = ceil(num_imgs/chunk_size);

    ss_boxes_dst_file_in_progress1 = regexprep(ss_boxes_dst_file, '.mat', '_in_progress.mat');
    ss_boxes_dst_file_in_progress2 = regexprep(ss_boxes_dst_file, '.mat', '_in_progress_prev.mat');

    try
        try
            ld = load(ss_boxes_dst_file_in_progress1);
            all_bbox_proposals = ld.all_bbox_proposals;
            first_chunk = ld.chunk + 1;
        catch
            ld = load(ss_boxes_dst_file_in_progress2);
            all_bbox_proposals = ld.all_bbox_proposals;
            first_chunk = ld.chunk + 1;
        end
    catch exception
        fprintf('Exception message %s\n', getReport(exception));
        all_bbox_proposals = cell(num_imgs,1);
        first_chunk = 1;
    end

    total_num_elems = 0;
    total_time = 0;
    count = 0;
    for chunk = first_chunk:num_chunks
        start_idx = (chunk-1) * chunk_size + 1;
        stop_idx  = min(chunk * chunk_size, num_imgs);
        th = tic;
        all_bbox_proposals(start_idx:stop_idx) = extract_selective_search_prposlas(image_db.image_paths(start_idx:stop_idx));
        for i = start_idx:stop_idx
            count = count + 1;
            total_num_elems = total_num_elems + numel(all_bbox_proposals{i});
        end
        elapsed_time = toc(th);
        total_time = total_time + elapsed_time;
        est_rem_time = (total_time / count) * (num_imgs - stop_idx);
        est_num_bytes = (total_num_elems / count) * num_imgs * 4 / (1024*1024*1024);
        fprintf('Extract Selective Search boxes %s %d/%d: ET %.2fmin | ETA %.2fmin | EST. NUM BYTES %.2f giga\n', ...
            image_db.image_set_name, stop_idx, num_imgs, ...
            total_time/60, est_rem_time/60, est_num_bytes);

        if (exist(ss_boxes_dst_file_in_progress1,'file')>0)
            copyfile(ss_boxes_dst_file_in_progress1,ss_boxes_dst_file_in_progress2);
        end

        save(ss_boxes_dst_file_in_progress1, 'all_bbox_proposals', 'chunk', '-v7.3');
    end

    save(ss_boxes_dst_file, 'all_bbox_proposals', '-v7.3');
end
end

function all_box_proposals = extract_selective_search_prposlas(image_paths)
fast_mode = true;
num_imgs = length(image_paths);
all_box_proposals = cell(num_imgs,1);
parfor (i = 1:num_imgs)
%     th = tic;
    img = imread(image_paths{i});
    all_box_proposals{i} = selective_search_boxes(img, fast_mode);
    all_box_proposals{i} = single(all_box_proposals{i}(:,[2 1 4 3]));
%     fprintf(' image %d/%d: elapsed time %.2f\n', i, num_imgs, toc(th))
end
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function im_rec = tomo_recon_myart_im(im, angles, n_it)
%TOMO_RECON_ART_IM   Tomographic reconstruction using ART method.
%
%   Phymhan
%   09-Aug-2013 09:15:41

if nargin < 3
    n_it = 100;
end
siz = size(im);
[W, p, ~, ~] = buildWeightMatrix(im,angles);
f = myart_solve(W,p,n_it);
im_rec = reshape(f,siz);
end

author: Sable

repo: mcbench-benchmarks


text:
function f = flops_inv(n)
% FLOPS_INV    Flops for matrix inversion.
% FLOPS_INV(n) returns the number of flops for inv(rand(n,n)).
% For n=1,2 the number of flops is exact.
% For n>2, the number of flops is estimated by assuming the matrix is
% decomposed via LU or Cholesky and then inverted by back-substitution.
% See flops_solve.

% Written by Tom Minka

if n == 2
  % inv([a b; c d]) = [d -b; -c a]/(ad-bc)
  % 3 flops for determinant, 1 divide, 4 multiplies, 1 negation
  % thanks to Zhang Xiaoying for pointing out this special case.
  f = 3+flops_div + 5;
else
  f = flops_solve(n,n,n);
end

author: andrewssobral

repo: lrslibrary


text:
function mae = CalMAE(smap, gtImg)
% Code Author: Wangjiang Zhu
% Email: wangjiang88119@gmail.com
% Date: 3/24/2014
[m,n,~] = size(gtImg);
smap = im2double(smap(:,:,1));
smap = imresize(smap,[m,n]);
smap = (smap-min(smap(:)))/(max(smap(:))-min(smap(:)));

if size(smap, 1) ~= size(gtImg, 1) || size(smap, 2) ~= size(gtImg, 2)
    error('Saliency map and gt Image have different sizes!\n');
end

if ~islogical(gtImg)
    gtImg = gtImg(:,:,1) > 128;
end

fgPixels = smap(gtImg);
fgErrSum = length(fgPixels) - sum(fgPixels);
bgErrSum = sum(smap(~gtImg));
mae = (fgErrSum + bgErrSum) / numel(gtImg);

author: ArcherFMY


text:
function [DI, k] = detectChange(obj, t1, t2)
% A naive implementation of the Multivariate Alteration Detection Algorithm
% Partially refer to http://people.compute.dtu.dk/alan/software.html
% and https://github.com/mortcanty/earthengine
[rows, cols, chns1] = size(t1);
chns2 = size(t2, 3);
chns = min(chns1, chns2);
x1 = reshape(double(t1), cols*rows, chns1);
x1 = x1 - mean(x1, 1);  % Subtract mean
x2 = reshape(double(t2), cols*rows, chns2);
x2 = x2 - mean(x2, 1);
sigma = cov([x1 x2]);
sigma11 = sigma(1:chns1, 1:chns1);
sigma12 = sigma(1:chns1, chns1+1:end);
sigma21 = sigma(chns1+1:end, 1:chns1);
sigma22 = sigma(chns1+1:end, chns1+1:end);

[a, e1] = eig(sigma12 / sigma22 * sigma21, sigma11);
e1 = sqrt(diag(e1));
[b, e2] = eig(sigma21 / sigma11 * sigma12, sigma22);
e2 = sqrt(diag(e2));
[~, idx1] = sort(e1, 'descend');
a = a(:,idx1);
a = a(:,chns:-1:1); % My choice is simply to drop the part of (q-p)
[~, idx2] = sort(e2, 'descend');
b = b(:,idx2);
b = b(:,chns:-1:1);

% Normalize a for unit dispersion
% Ensure that a'*s11*a=I to meet the constraints
vars1 = (a'*sigma11*a);
a = a ./ sqrt(diag(vars1))';
% Similiar operations on b
vars2 = (b'*sigma22*b);
b = b ./ sqrt(diag(vars2))';

% Ensure sum of positive correlations between x1 and x1*a is positive
% I have no idea what this is used for???
invStd1 = diag(1./std(x1));
sgn = diag(sign(sum(invStd1*sigma11*a)));
a = a*sgn;

% Assure positive correlation between pair of canonical variates
b = b .* diag(sign(a'*sigma12*b))';

mad = x1*a - x2*b;

% Normalize mad
% This is no regular operation, yet appears to improve the result
mad = (mad - mean(mad)) ./ (std(mad)+eps);

chi2 = sum(mad.^2, 2);
prob = chi2cdf(chi2, chns);
DI = reshape(abs(mad), rows, cols, chns);

% Select bands
k = min(chns, obj.nMADUsed);
DI = DI(:,:,1:k);
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function a = myiswt2(swc, lo_r, hi_r)
% ISWT2 Inverse discrete stationary wavelet transform 2-D.
% ISWT2 performs a multilevel 2-D stationary wavelet
% reconstruction using either a specific orthogonal wavelet
% ('wname', see WFILTERS for more information) or specific
% reconstruction filters (Lo_r and Hi_r).
%
% Filters lo_r and hi_r must be normalized by dividing by sqrt(2)
% (necessary for implmentation of Stationary Discrete WT)
%
% For X = MYISWT2(swc,Lo_r,Hi_r)
% swc contains H(1:N) V(1:N) D(1:N) A(N)
% lo_r is the reconstruction low-pass filter.
% hi_r is the reconstruction high-pass filter.
%
% See also IDWT2, SWT2, WAVEREC2.

% based on iswt2.m
% Sathish Ramani, October 11, 2009

%% Set DWT_Mode to 'per'.
old_modeDWT = dwtmode('status','nodisp');
modeDWT = 'per';
dwtmode(modeDWT,'nodisp');

%% Load coefficients
n = (size(swc, 3)-1)/3; %% Number of levels
h = swc(:,:,1:n);
v = swc(:,:,n+1:2*n);
d = swc(:,:,2*n+1:3*n);
a = swc(:,:,3*n+1);

[rx, cx, dump] = size(h);
for k = n:-1:1
	step = 2^(k-1);
	last = step;
	for first1 = 1:last
		iRow = first1:step:rx;
		lR   = length(iRow);
		for first2 = 1:last
			iCol = first2:step:cx;
			lC   = length(iCol);
			sR   = iRow;
			sC   = iCol;
			a(iRow,iCol) = idwt2LOC(a(sR,sC), h(sR,sC,k), v(sR,sC,k), d(sR,sC,k), ...
				lo_r, hi_r, [lR lC]);
		end
	end
end

% Restore DWT_Mode.
dwtmode(old_modeDWT, 'nodisp');

%===============================================================%
% INTERNAL FUNCTIONS
%===============================================================%
function y = idwt2LOC(a, h, v, d, lo_r, hi_r, sy)

y = upconvLOC(a,lo_r,lo_r,sy) ... % Approximation
	+ upconvLOC(h,hi_r,lo_r,sy) ... % Horizontal Detail
	+ upconvLOC(v,lo_r,hi_r,sy) ... % Vertical Detail
	+ upconvLOC(d,hi_r,hi_r,sy); % Diagonal Detail

%---------------------------------------------------------------%
function y = upconvLOC(x, f1, f2, s)

lf = length(f1);
% y  = dyadup(x,'mat',0,1);
y = x;
y  = wextend('2D', 'per', y, [lf/2,lf/2]);
y  = wconv2('col', y, f1);
y  = wconv2('row', y, f2);
y  = wkeep2(y, s, [lf lf]);
%===============================================================%

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function varargout = lvd_EditVectorPlaneAngleGUI(varargin)
% LVD_EDITVECTORPLANEANGLEGUI MATLAB code for lvd_EditVectorPlaneAngleGUI.fig
%      LVD_EDITVECTORPLANEANGLEGUI, by itself, creates a new LVD_EDITVECTORPLANEANGLEGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITVECTORPLANEANGLEGUI returns the handle to a new LVD_EDITVECTORPLANEANGLEGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITVECTORPLANEANGLEGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITVECTORPLANEANGLEGUI.M with the given input arguments.
%
%      LVD_EDITVECTORPLANEANGLEGUI('Property','Value',...) creates a new LVD_EDITVECTORPLANEANGLEGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditVectorPlaneAngleGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditVectorPlaneAngleGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_EditVectorPlaneAngleGUI

% Last Modified by GUIDE v2.5 23-Feb-2021 08:38:27

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditVectorPlaneAngleGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditVectorPlaneAngleGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_EditVectorPlaneAngleGUI is made visible.
function lvd_EditVectorPlaneAngleGUI_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to lvd_EditVectorPlaneAngleGUI (see VARARGIN)

    % Choose default command line output for lvd_EditVectorPlaneAngleGUI
    handles.output = hObject;

    centerUIFigure(hObject);

    angle = varargin{1};
    setappdata(hObject, 'angle', angle);

    lvdData = varargin{2};
    setappdata(hObject, 'lvdData', lvdData);

	populateGUI(handles, lvdData, angle);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_EditVectorPlaneAngleGUI wait for user response (see UIRESUME)
    uiwait(handles.lvd_EditVectorPlaneAngleGUI);

function populateGUI(handles, lvdData, angle)
    set(handles.nameText,'String',angle.getName());

    [vectorsListBoxStr, ~] = lvdData.geometry.vectors.getListboxStr();
    [planesListBoxStr, ~] = lvdData.geometry.planes.getListboxStr();

    handles.vectorCombo.String = vectorsListBoxStr;
    handles.vectorCombo.Value = lvdData.geometry.vectors.getIndsForVectors(angle.vector);

    handles.planeCombo.String = planesListBoxStr;
    handles.planeCombo.Value = lvdData.geometry.planes.getIndsForPlanes(angle.plane);

    handles.lineColorCombo.String = ColorSpecEnum.getListboxStr();
    handles.lineColorCombo.Value = ColorSpecEnum.getIndForName(angle.lineColor.name);

    handles.lineSpecCombo.String = LineSpecEnum.getListboxStr();
    handles.lineSpecCombo.Value = LineSpecEnum.getIndForName(angle.lineSpec.name);

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditVectorPlaneAngleGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = false;
    else
        lvdData = getappdata(hObject, 'lvdData');
        angle = getappdata(hObject, 'angle');

        angle.setName(handles.nameText.String);

        angle.vector = lvdData.geometry.vectors.getVectorAtInd(handles.vectorCombo.Value);
        angle.plane = lvdData.geometry.planes.getPlaneAtInd(handles.planeCombo.Value);

        str = handles.lineColorCombo.String{handles.lineColorCombo.Value};
        angle.lineColor = ColorSpecEnum.getEnumForListboxStr(str);

        str = handles.lineSpecCombo.String{handles.lineSpecCombo.Value};
        angle.lineSpec = LineSpecEnum.getEnumForListboxStr(str);

        varargout{1} = true;
        close(handles.lvd_EditVectorPlaneAngleGUI);
    end

% --- Executes on button press in saveAndCloseButton.
function saveAndCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveAndCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    errMsg = validateInputs(handles);

    if(isempty(errMsg))
        uiresume(handles.lvd_EditVectorPlaneAngleGUI);
    else
        msgbox(errMsg,'Invalid Angle Inputs','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    if(isempty(strtrim(handles.nameText.String)))
        errMsg{end+1} = 'Angle name must contain more than white space and must not be empty.';
    end

% --- Executes on button press in cancelButton.
function cancelButton_Callback(hObject, eventdata, handles)
% hObject    handle to cancelButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    close(handles.lvd_EditVectorPlaneAngleGUI);

function nameText_Callback(hObject, eventdata, handles)
% hObject    handle to nameText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of nameText as text
%        str2double(get(hObject,'String')) returns contents of nameText as a double

% --- Executes during object creation, after setting all properties.
function nameText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to nameText see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on lvd_EditVectorPlaneAngleGUI or any of its controls.
function lvd_EditVectorPlaneAngleGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditVectorPlaneAngleGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_EditVectorPlaneAngleGUI);
    end

% --- Executes on selection change in lineSpecCombo.
function lineSpecCombo_Callback(hObject, eventdata, handles)
% hObject    handle to lineSpecCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns lineSpecCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from lineSpecCombo

% --- Executes during object creation, after setting all properties.
function lineSpecCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to lineSpecCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in lineColorCombo.
function lineColorCombo_Callback(hObject, eventdata, handles)
% hObject    handle to lineColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns lineColorCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from lineColorCombo

% --- Executes during object creation, after setting all properties.
function lineColorCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to lineColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in vectorCombo.
function vectorCombo_Callback(hObject, eventdata, handles)
% hObject    handle to vectorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns vectorCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from vectorCombo

% --- Executes during object creation, after setting all properties.
function vectorCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to vectorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in planeCombo.
function planeCombo_Callback(hObject, eventdata, handles)
% hObject    handle to planeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns planeCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from planeCombo

% --- Executes during object creation, after setting all properties.
function planeCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to planeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [node_ele_assignment] = revert_assignment(elements,sparse_flag)

% REVERT_ASSIGNMENT
%
% $Id$

N_nodes = max(max(elements));
N_elements = max(size(elements));
[sort_elements_vec, sort_ind] = sort(elements(:));
[dummy, m, dummy] = unique(sort_elements_vec,'last');
occurence = diff([0;m]);

node_ele_assignment = zeros(max(occurence),N_nodes);
sort_ind = mod(sort_ind,N_elements);
sort_ind(sort_ind == 0) = N_elements;
flip = diff([sort_elements_vec;sort_elements_vec(end)]);
col = 1;
row = 1;

for i=1:length(flip)
    node_ele_assignment(row,col) = sort_ind(i);
    row = row+1;
    if(flip(i))
        col = col+1;
        row = 1;
    end
end

s1 = whos('node_ele_assignment');
sp = sparse(node_ele_assignment);
s2 = whos('sp');
if(s2.bytes/s1.bytes < sparse_flag)
    node_ele_assignment = sp;
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% trigradient
% Below is a demonstration of the features of the |trigradient| function

%%
clear; close all; clc;

%% Syntax
% |[ux,uy]=trigradient(TRI,V,INT);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


function [region_meta_info]=preprocess_mcg_candidates(imnames, mcgdir, gtdir, ovoutdir, sptextdir, regspimgdir, N)
if(~exist(ovoutdir, 'file'))
    mkdir(ovoutdir);
end
if(~exist(sptextdir, 'file'))
    mkdir(sptextdir);
end
if(~exist(regspimgdir, 'file'))
    mkdir(regspimgdir);
end

for i=1:numel(imnames)
    tmp=load(fullfile(mcgdir, [imnames{i} '.mat']));
    candidates=tmp.candidates; clear tmp;
    %get reg2sp
    sp=candidates.superpixels;
    N1=min(N, numel(candidates.labels));
    reg2sp=false(max(sp(:)),N1);
    for j=1:N1
        reg2sp(candidates.labels{j},j)=true;
    end

    %compress to remove irrelevant sps
    [sp, reg2sp]=compressSp2reg(sp, reg2sp);

    if(~exist(fullfile(ovoutdir, [imnames{i} '.mat']), 'file'))
    %compute overlaps between everything
    ov=double(reg2sp')*double(reg2sp);
    ov=ov>0;
    save(fullfile(ovoutdir, [imnames{i} '.mat']), 'ov');
    end
    %load gt
    if(isempty(gtdir))
        %no gt
        inst=zeros(size(sp));
        categories=[];
    else
    [cls, inst, categories]=load_gt(gtdir, imnames{i});
    end
    %compute overlaps with ground truth
    overlap=get_gt_overlaps(reg2sp, sp, inst);
    box_overlap=get_gt_overlaps_box(reg2sp, sp, inst);
    %populate meta info
    region_meta_info.num_regs(i)=size(reg2sp,2)+numel(categories);
    region_meta_info.overlaps{i}=overlap;
    region_meta_info.box_overlaps{i}=box_overlap;
    region_meta_info.gt{i}=[categories(:)' zeros(1,size(reg2sp,2))];

    %add ground truth to reg2sp
    gt_sp=double(inst+1);
    gt_reg2sp=[zeros(1,numel(categories)); eye(numel(categories))];
    gt_reg2sp=logical(gt_reg2sp);
    [sp, reg2sp]=intersect_partitions(gt_sp, sp, gt_reg2sp, reg2sp);
    region_meta_info.boxes{i}=get_region_boxes(sp, reg2sp);
    %write sprep

    textfile=fullfile(sptextdir, [imnames{i} '.txt']);
    if(~exist(textfile, 'file'))
    write_sprep_text(sp, textfile);
    reg2spfile=fullfile(regspimgdir, [imnames{i} '.png']);
    imwrite(uint8(reg2sp),reg2spfile);
    end
    fprintf('Done %d\n', i);
end

function write_sprep_text(sp, filename)
fid=fopen(filename, 'w');
fprintf(fid,'%d %d\n', size(sp,1), size(sp,2));
fprintf(fid, '%d ', sp(:));
fclose(fid);

function [sp, sp2reg]=intersect_partitions(sp1, sp2, sp2reg1, sp2reg2)
%do a unique
[unique_pairs, junk, sp]=unique([sp1(:) sp2(:)], 'rows');
sp=reshape(sp, size(sp1));

%create new sp2reg
sp2reg=[sp2reg1(unique_pairs(:,1),:) sp2reg2(unique_pairs(:,2),:)];

%break disconnected sp
 spNew = zeros(size(sp));
  sp2regNew = false(size(sp2reg));
  cnt = 0;
  %% Check connectivity of superpixels
  for i = 1:max(sp(:))
    tt = bwconncomp(sp == i);
    for j = 1:tt.NumObjects,
      cnt = cnt+1;
      spNew(tt.PixelIdxList{j}) = cnt;
      sp2regNew(cnt, :) = sp2reg(i, :);
    end
  end

sp=spNew;
sp2reg=sp2regNew;

function [spNew sp2regNew] = compressSp2reg(sp, sp2reg)
  assert(islogical(sp2reg), 'sp2reg should be logical');
  [sp2reg bb spSame] = unique(sp2reg, 'rows');
  sp = spSame(sp);

  spNew = zeros(size(sp));
  sp2regNew = false(size(sp2reg));
  cnt = 0;
  %% Check connectivity of superpixels
  for i = 1:max(sp(:))
    tt = bwconncomp(sp == i);
    for j = 1:tt.NumObjects,
      cnt = cnt+1;
      spNew(tt.PixelIdxList{j}) = cnt;
      sp2regNew(cnt, :) = sp2reg(i, :);
    end
  end


text:
function [y,S,k] = spm_csd_fmri_mar(P,M,U)
% Prediction of MAR coefficients for DCM
% FORMAT [y,S,K] = spm_csd_fmri_mar(P,M,U)
%
% P - model parameters
% M - model structure
% U - model inputs (expects U.csd as complex cross spectra)
%
% y - y(nw,nn,nn} - cross-spectral density for nn nodes
%                 - for nw frequencies in M.Hz
% K - Volterra kernels
% S - directed transfer functions (complex)
%
% This routine computes the spectral response of a network of regions
% driven by  endogenous fluctuations and exogenous (experimental) inputs.
% It returns the complex cross spectra of regional responses as a
% three-dimensional array. The endogenous innovations or fluctuations are
% parameterised in terms of a (scale free) power law, in frequency space.
%
% When the observer function M.g is specified, the CSD response is
% supplemented with observation noise in sensor space; otherwise the CSD
% is noiseless.
%
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_csd_fmri_mar.m 5618 2013-08-17 10:36:56Z karl $

% number of nodes and endogenous (neuronal) fluctuations
%--------------------------------------------------------------------------
np   = M.p;                                   % number of MAR lags
nn   = M.l;                                   % number of nodes (regions)

% cross-covaraince functions of neuronal fluctations (Vu) and noise (Vn)
%==========================================================================

% experimental inputs
%--------------------------------------------------------------------------
for i = 1:nn
    for j = 1:nn
        if any(any(P.C))
            for k = 1:M.N
                V(k) = P.C(i,:)*squeeze(U.ccf(k,:,:))*P.C(j,:)';
            end
            Vu{i,j} = toeplitz(V);
        else
            Vu{i,j} = sparse(M.N,M.N);
        end
    end
end

% neuronal inputs
%--------------------------------------------------------------------------
for i = 1:nn
    Vu{i,i} = Vu{i,i} + exp(P.a(1,i))*spm_Q(exp(P.a(2,i))/2,M.N);
end
Vu    = spm_cat(Vu);

% observation noise
%--------------------------------------------------------------------------
for i = 1:nn

    % global component
    %----------------------------------------------------------------------
    for j = 1:nn
        V       = exp(P.b(1,1))*spm_Q(exp(P.b(2,1))/2,np + 1)/64;
        Vn{i,j} = V((1:np),(1:np));
        Rn{i,j} = V((1:np) + 1,1);
    end

    % region specific
    %----------------------------------------------------------------------
    V       = exp(P.c(1,i))*spm_Q(exp(P.c(2,i))/2,np + 1)/8;
    Vn{i,i} = Vn{i,j} + V((1:np),(1:np));
    Rn{i,i} = Rn{i,j} + V((1:np) + 1,1);

end
Vn    = spm_cat(Vn);
Rn    = spm_cat(Rn);

% first-order Volterra kernel
%==========================================================================
P.C   = speye(nn,nn);
[S,k] = spm_dcm_mtf(P,M);

% matix form
%--------------------------------------------------------------------------
for i = 1:nn
    for j = 1:nn
        K{i,j} = k(:,i,j);
    end
end
K     = spm_cat(K);

% lagged matix form
%--------------------------------------------------------------------------
for i = 1:nn
    for j = 1:nn
        KK{i,j}    = zeros(nn,np);
        for p = 1:np
            t = (1 + p):size(k,1);
            KK{i,j}(t,p) = k(t - p,i,j);
        end
    end
end
KK    = spm_cat(KK);

% predicted MAR coefficients
%--------------------------------------------------------------------------
y     = spm_inv(KK'*Vu*KK + Vn)*(KK'*Vu*K + Rn);

author: spm


text:
function layer = genNetworkFeedForward_obj_auto(inputDim, hiddenLayerSize, outputDim, costFn, LastActivation4MSE)
if nargin<5
    LastActivation4MSE = 'linear';
end
inputDim = double(inputDim);

inputStreamIdx = 1;
layer{1} = InputNode(inputStreamIdx,inputDim);

for i=1:length(hiddenLayerSize)
    layer{end+1} = AffineNode(hiddenLayerSize(i));
    layer{end+1} = SigmoidNode(hiddenLayerSize(i));
end
layer{end+1} = AffineNode(outputDim);

if strcmpi(costFn, 'CrossEntropy')
    layer{end+1} = SoftmaxNode(outputDim);
else    % there are a few options for MSE cost function
    switch LastActivation4MSE
        case 'linear'
            % do nothing
        case 'tanh'
            layer{end+1} = TanhNode(outputDim);
        case 'sigmoid'
            layer{end+1} = SigmoidNode(outputDim);
        case 'relu'
            layer{end+1} = ReluNode(outputDim);
        otherwise
            fprintf('Error: unknown activation for last hidden layer, skipped!\n');
    end
end
layer = ConnectNodesLinear(layer);

inputStreamIdx = inputStreamIdx +1;
if strcmpi(costFn, 'CrossEntropy')
    layer{end+1} = InputNode(inputStreamIdx, 1);
    layer{end+1} = CrossEntropyNode();
else
    layer{end+1} = InputNode(inputStreamIdx, outputDim);
    layer{end+1} = MeanSquareErrorNode(1);
end
layer{end}.prev = [-2 -1];
layer{end}.dim = [1 outputDim];
layer = FinishLayer_obj(layer);
end

author: singaxiong

repo: SignalGraph


text:
function failed_old_specest_vs_oldimplementation

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

% read in data from /home/common/testdata (I just my own testdata)
% data = .....

%% COHERENCE specest_mtmfft
cfg = [];
cfg.keeptrials = 'yes';
cfg.foilim     = [2 100];
cfg.output     = 'powandcsd';
cfg.taper      = 'dpss';
cfg.tapsmofrq  = 4;
cfg.pad        = 'maxperlen';
cfg.calcdof    = 'yes';
cfg.keeptrials = 'no';
% old style
cfg.method     = 'mtmfft_old';
tic; freqold = ft_freqanalysis(cfg,data); toc
% new style
cfg.method     = 'mtmfft';
tic; freqnew = ft_freqanalysis(cfg,data); toc

% calculate coherence
cfg = [];
cfg.method = 'coh';
connold = ft_connectivityanalysis(cfg,freqold);
connnew = ft_connectivityanalysis(cfg,freqnew);

% plot 'a' chancomb
figure
subplot(3,1,1)
plot(connold.freq,squeeze(connold.cohspctrm(10,:)))
axis tight
xlabel('frequency (Hz)')
ylabel('coh')
title('old implementation')
subplot(3,1,2)
plot(connnew.freq,squeeze(connnew.cohspctrm(10,:)))
axis tight
xlabel('frequency (Hz)')
ylabel('coh')
title('new implementation')
subplot(3,1,3)
hold on
plot(connold.freq,squeeze(connold.cohspctrm(10,:)))
plot(connnew.freq,squeeze(connnew.cohspctrm(10,:)),'r')
legend
axis tight
xlabel('frequency (Hz)')
ylabel('coh')
title('both overlayed, using freq axis of new implementation')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% POWER specest_mtmfft
cfg = [];
cfg.trials     = 1:3;
cfg.keeptrials = 'yes';
cfg.foilim     = [0 256/2];
cfg.output     = 'pow';
cfg.taper      = 'dpss';
cfg.tapsmofrq  = 4;
cfg.pad        = 'maxperlen';
cfg.calcdof    = 'yes';
% old style
cfg.method     = 'mtmfft_old';
tic; freqold = ft_freqanalysis(cfg,data); toc
% new style
cfg.method     = 'mtmfft';
tic; freqnew = ft_freqanalysis(cfg,data); toc

% plot first trial first channel, single timepoint
figure
subplot(3,1,1)
plot(freqold.freq,squeeze(freqold.powspctrm(1,1,:)))
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('old implementation')
subplot(3,1,2)
plot(freqnew.freq,squeeze(freqnew.powspctrm(1,1,:)))
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('new implementation')
subplot(3,1,3)
hold on
plot(freqnew.freq,squeeze(freqold.powspctrm(1,1,:)))
plot(freqnew.freq,squeeze(freqnew.powspctrm(1,1,:)),'r')
legend
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('both overlayed, using freq axis of new implementation')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% COHERENCE specest_mtmconvol
cfg = [];
cfg.trials     = 1:2;
cfg.keeptrials = 'no';
cfg.keeptapers = 'no';
cfg.channel    = 'all';
cfg.channelcmb = {'all' 'all'};
%cfg.pad        = 4.24;
cfg.foi        = 2:2:20;
cfg.toi        = data.time{1}(1):0.050:data.time{1}(end);
cfg.output     = 'powandcsd';
cfg.taper      = 'hanning';
cfg.t_ftimwin  = 4 ./ cfg.foi; %ones(length(cfg.foi),1) .* 1.5; %
cfg.tapsmofrq  = ones(length(cfg.foi),1) * 8;
cfg.calcdof    = 'yes';
cfg.correctt_ftimwin = 'no';
% old style
cfg.method     = 'mtmconvol_old';
tic; freqold = ft_freqanalysis(cfg,data); toc
% new style
cfg.method     = 'mtmconvol';
tic; freqnew = ft_freqanalysis(cfg,data); toc

% calculate coherence
cfg = [];
cfg.method = 'coh';
connold = ft_connectivityanalysis(cfg,freqold);
connnew = ft_connectivityanalysis(cfg,freqnew);

% plot 'a' chancomb at a single timepoint
timepoint = 20;
figure
subplot(3,1,1)
plot(connold.freq,squeeze(connold.cohspctrm(10,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('coh')
title('old implementation')
subplot(3,1,2)
plot(connnew.freq,squeeze(connnew.cohspctrm(10,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('coh')
title('new implementation')
subplot(3,1,3)
hold on
plot(connold.freq,squeeze(connold.cohspctrm(10,:,timepoint)))
plot(connnew.freq,squeeze(connnew.cohspctrm(10,:,timepoint)),'r')
legend
axis tight
xlabel('frequency (Hz)')
ylabel('coh')
title('both overlayed, using freq axis of new implementation')

% plot 'a' chancomb at a all timepoints
figure
subplot(3,1,1)
imagesc(connold.time,connold.freq,squeeze(connold.cohspctrm(10,:,:)))
axis tight
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old implementation')
colorbar; axis xy; caxis([0 1])
subplot(3,1,2)
imagesc(connnew.time,connnew.freq,squeeze(connnew.cohspctrm(10,:,:)))
axis tight
ylabel('frequency (Hz)')
xlabel('time (s)')
title('new implementation')
colorbar; axis xy; caxis([0 1]);
subplot(3,1,3)
hold on
imagesc(connold.time,connold.freq,squeeze(connold.cohspctrm(10,:,:))-squeeze(connnew.cohspctrm(10,:,:)))
legend
axis tight
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old minus new')
colorbar; caxis([0 1]); axis xy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% POWER specest_mtmconvol
cfg = [];
cfg.implementation = 'old';
cfg.trials     = 1:3;
cfg.keeptrials = 'yes';
cfg.keeptapers = 'no';
%cfg.pad        = 4.24;
cfg.foi        = 0:2:128;
cfg.toi        = data.time{1}(1):0.0050:data.time{1}(end);
cfg.output     = 'pow';
cfg.taper      = 'hanning';
cfg.t_ftimwin  = 4 ./ cfg.foi; %ones(length(cfg.foi),1) .* 1.5; %
cfg.tapsmofrq  = ones(length(cfg.foi),1) * 8;
cfg.calcdof    = 'yes';
cfg.correctt_ftimwin = 'no';
% new style
cfg.method     = 'mtmconvol';
tic; freqnew = ft_freqanalysis(cfg,data); toc
% old style
cfg.method     = 'mtmconvol_old';
cfg.foi = freqnew.freq;
cfg.t_ftimwin  = 4 ./ cfg.foi; %ones(length(cfg.foi),1) .* 1.5; %
tic; freqold = ft_freqanalysis(cfg,data); toc

% plot first trial first channel, single timepoint
timepoint = 20;
figure
subplot(3,1,1)
plot(freqold.freq,squeeze(freqold.powspctrm(1,1,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('old implementation')
subplot(3,1,2)
plot(freqnew.freq,squeeze(freqnew.powspctrm(1,1,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('new implementation')
subplot(3,1,3)
hold on
plot(freqnew.freq,squeeze(freqold.powspctrm(1,1,:,timepoint)))
plot(freqnew.freq,squeeze(freqnew.powspctrm(1,1,:,timepoint)),'r')
legend
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('both overlayed, using freq axis of new implementation')

% plot first trial first channel, all timepoints
figure
subplot(3,1,1)
imagesc(freqold.time,freqold.freq,squeeze(freqold.powspctrm(1,1,:,:)))
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old implementation')
subplot(3,1,2)
imagesc(freqold.time,freqnew.freq,squeeze(freqnew.powspctrm(1,1,:,:)))
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('new implementation')
subplot(3,1,3)
hold on
imagesc(freqold.time,freqold.freq,squeeze(freqold.powspctrm(1,1,:,:)) - squeeze(freqnew.powspctrm(1,1,:,:)))
legend
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old minus new')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% POWER specest_wavelet
cfg = [];
cfg.trials     = 1;
cfg.keeptrials = 'yes';
cfg.keeptapers = 'no';
%cfg.pad        = 4.24;
cfg.foi        = 2:1:40;
cfg.toi        = 0:0.05:4; %data.time{1}(1):(1/data.fsample):data.time{1}(end);
cfg.output     = 'pow';
cfg.calcdof    = 'yes';
cfg.correctt_ftimwin = 'no';
% new style
cfg.method     = 'wavelet';
tic; freqnew = ft_freqanalysis(cfg,data); toc
% old style
cfg.foi = freqnew.freq;
cfg.method     = 'wltconvol_old';
tic; freqold = ft_freqanalysis(cfg,data); toc

% plot first trial first channel, single timepoint
timepoint = 40;
figure
subplot(3,1,1)
plot(freqold.freq,squeeze(freqold.powspctrm(1,1,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('old implementation')
subplot(3,1,2)
plot(freqnew.freq,squeeze(freqnew.powspctrm(1,1,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('new implementation')
subplot(3,1,3)
hold on
plot(freqnew.freq,squeeze(freqold.powspctrm(1,1,:,timepoint)))
plot(freqnew.freq,squeeze(freqnew.powspctrm(1,1,:,timepoint)),'r')
legend
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('both overlayed, using freq axis of new implementation')

% plot first trial first channel, all timepoints
figure
subplot(3,1,1)
imagesc(freqold.time,freqold.freq,squeeze(freqold.powspctrm(1,1,:,:)))
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old implementation')
subplot(3,1,2)
imagesc(freqold.time,freqnew.freq,squeeze(freqnew.powspctrm(1,1,:,:)))
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('new implementation')
subplot(3,1,3)
hold on
imagesc(freqold.time,freqold.freq,squeeze(freqold.powspctrm(1,1,:,:)) - squeeze(freqnew.powspctrm(1,1,:,:)))
legend
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old minus new')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% CSD specest_wavelet
cfg = [];
cfg.trials     = 'all';
cfg.keeptrials = 'yes';
cfg.keeptapers = 'no';
%cfg.pad        = 4.24;
cfg.foi        = 2:2:40;
cfg.toi        = data.time{1}(1):(1/data.fsample):data.time{1}(end);
cfg.output     = 'powandcsd';
cfg.calcdof    = 'yes';
cfg.correctt_ftimwin = 'no';
% new style
cfg.method     = 'wavelet';
tic; freqnew = ft_freqanalysis(cfg,data); toc
% old style
cfg.foi = freqnew.freq;
cfg.method     = 'wltconvol_old';
tic; freqold = ft_freqanalysis(cfg,data); toc

% plot first trial first channel, single timepoint
timepoint = 200;
figure
subplot(3,1,1)
plot(freqold.freq,abs(squeeze(freqold.crsspctrm(1,1,:,timepoint))))
axis tight
xlabel('frequency (Hz)')
ylabel('')
title('old implementation')
subplot(3,1,2)
plot(freqnew.freq,abs(squeeze(freqnew.crsspctrm(1,1,:,timepoint))))
axis tight
xlabel('frequency (Hz)')
ylabel('')
title('new implementation')
subplot(3,1,3)
hold on
plot(freqnew.freq,abs(squeeze(freqold.crsspctrm(1,1,:,timepoint))))
plot(freqnew.freq,abs(squeeze(freqnew.crsspctrm(1,1,:,timepoint))),'r')
legend
axis tight
xlabel('frequency (Hz)')
ylabel('')
title('both overlayed, using freq axis of new implementation')

% plot first trial first channel, all timepoints
figure
subplot(3,1,1)
imagesc(freqold.time,freqold.freq,abs(squeeze(freqold.crsspctrm(1,1,:,:))))
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old implementation')
subplot(3,1,2)
imagesc(freqold.time,freqnew.freq,abs(squeeze(freqnew.crsspctrm(1,1,:,:))))
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('new implementation')
subplot(3,1,3)
hold on
imagesc(freqold.time,freqold.freq,abs(squeeze(freqold.crsspctrm(1,1,:,:)) - squeeze(freqnew.crsspctrm(1,1,:,:))))
legend
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old minus new')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% COHERENCE specest_wavelet
cfg = [];
cfg.trials     = 'all';
cfg.keeptrials = 'yes';
cfg.keeptapers = 'no';
cfg.channel    = 'all';
cfg.channelcmb = {'all' 'all'};
%cfg.pad        = 4.24;
cfg.foi        = 2:2:20;
cfg.toi        = data.time{1}(1):0.050:data.time{1}(end);
cfg.output     = 'powandcsd';
cfg.calcdof    = 'yes';
% new style
cfg.method     = 'wavelet';
tic; freqnew = ft_freqanalysis(cfg,data); toc
% old style
cfg.foi = freqnew.freq;
cfg.method     = 'wltconvol_old';
tic; freqold = ft_freqanalysis(cfg,data); toc

% calculate coherence
cfg = [];
cfg.method = 'coh';
connold = ft_connectivityanalysis(cfg,freqold);
connnew = ft_connectivityanalysis(cfg,freqnew);

% plot 'a' chancomb at a single timepoint
timepoint = 20;
figure
subplot(3,1,1)
plot(connold.freq,squeeze(connold.cohspctrm(10,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('coh')
title('old implementation')
subplot(3,1,2)
plot(connnew.freq,squeeze(connnew.cohspctrm(10,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('coh')
title('new implementation')
subplot(3,1,3)
hold on
plot(connold.freq,squeeze(connold.cohspctrm(10,:,timepoint)))
plot(connnew.freq,squeeze(connnew.cohspctrm(10,:,timepoint)),'r')
legend
axis tight
xlabel('frequency (Hz)')
ylabel('coh')
title('both overlayed, using freq axis of new implementation')

% plot 'a' chancomb at a all timepoints
figure
subplot(3,1,1)
imagesc(connold.time,connold.freq,squeeze(connold.cohspctrm(10,:,:)))
axis tight
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old implementation')
colorbar; axis xy; caxis([0 1])
subplot(3,1,2)
imagesc(connnew.time,connnew.freq,squeeze(connnew.cohspctrm(10,:,:)))
axis tight
ylabel('frequency (Hz)')
xlabel('time (s)')
title('new implementation')
colorbar; axis xy; caxis([0 1]);
subplot(3,1,3)
hold on
imagesc(connold.time,connold.freq,squeeze(connold.cohspctrm(10,:,:))-squeeze(connnew.cohspctrm(10,:,:)))
legend
axis tight
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old minus new')
colorbar; caxis([0 1]); axis xy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% POWER specest_hilbert
cfg = [];
cfg.trials     = 1:3;
cfg.keeptrials = 'yes';
cfg.keeptapers = 'no';
cfg.width      = 1;
cfg.pad        = 'maxperlen';
cfg.foi        = 10:2:40;
cfg.toi        = 0:0.05:4; %data.time{1}(1):(1/data.fsample):data.time{1}(end);
cfg.output     = 'pow';
cfg.calcdof    = 'yes';
cfg.correctt_ftimwin = 'no';
% new style
cfg.method     = 'hilbert';
cfg.filtorder  = round((3 ./ cfg.foi) .* data.fsample);
cfg.filttype   = 'fir';
cfg.width      = 5;

tic; freqnew = ft_freqanalysis(cfg,data); toc
% old style
cfg.foi = freqnew.freq;
cfg.method     = 'mtmconvol';
cfg.taper      = 'hanning';
cfg.t_ftimwin  = 3 ./ cfg.foi; %ones(length(cfg.foi),1) .* 1.5; %
tic; freqold = ft_freqanalysis(cfg,data); toc

% plot first trial first channel, single timepoint
timepoint = 40;
figure
subplot(3,1,1)
plot(freqold.freq,squeeze(freqold.powspctrm(1,1,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('mtmconvol')
subplot(3,1,2)
plot(freqnew.freq,squeeze(freqnew.powspctrm(1,1,:,timepoint)))
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('hilbert')
subplot(3,1,3)
hold on
plot(freqnew.freq,squeeze(freqold.powspctrm(1,1,:,timepoint)))
plot(freqnew.freq,squeeze(freqnew.powspctrm(1,1,:,timepoint)),'r')
legend
axis tight
xlabel('frequency (Hz)')
ylabel('power')
title('both overlayed, using freq axis of hilbert')

% plot first trial first channel, all timepoints
figure
subplot(3,1,1)
imagesc(freqold.time,freqold.freq,squeeze(freqold.powspctrm(1,1,:,:)))
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('mtmconvol')
subplot(3,1,2)
imagesc(freqold.time,freqnew.freq,squeeze(freqnew.powspctrm(1,1,:,:)))
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('hilbert')
subplot(3,1,3)
hold on
imagesc(freqold.time,freqold.freq,squeeze(freqold.powspctrm(1,1,:,:)) - squeeze(freqnew.powspctrm(1,1,:,:)))
legend
axis tight; axis xy; colorbar
ylabel('frequency (Hz)')
xlabel('time (s)')
title('old minus new')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% FOURIER specest_hilbert
cfg = [];
cfg.trials     = 1:3;
cfg.keeptrials = 'yes';
cfg.keeptapers = 'yes';
cfg.width      = 1;
cfg.pad        ='maxperlen';
cfg.foi        = 10:2:40;
cfg.toi        = 0:0.05:4; %data.time{1}(1):(1/data.fsample):data.time{1}(end);
cfg.output     = 'fourier';
cfg.calcdof    = 'yes';
cfg.correctt_ftimwin = 'no';
% new style
cfg.method     = 'hilbert';
cfg.filtorder  = round((3 ./ cfg.foi) .* data.fsample);
cfg.filttype   = 'fir';
cfg.width      = 1;
tic; freqnew = ft_freqanalysis(cfg,data); toc
% old style
cfg.foi = freqnew.freq;
cfg.method     = 'mtmconvol';
cfg.taper      = 'hanning';
cfg.t_ftimwin  = 3 ./ cfg.foi; %ones(length(cfg.foi),1) .* 1.5; %
tic; freqold = ft_freqanalysis(cfg,data); toc

% plot first trial first channel, real part
freq = 3;
figure
subplot(3,1,1)
plot(freqold.time,real(squeeze(freqold.fourierspctrm(1,1,freq,:))))
axis tight
xlabel('frequency (Hz)')
ylabel('amplitude')
title('mtmconvol')
subplot(3,1,2)
plot(freqnew.time,real(squeeze(freqnew.fourierspctrm(1,1,freq,:))))
axis tight
xlabel('frequency (Hz)')
ylabel('amplitude')
title('hilbert')
subplot(3,1,3)
hold on
plot(freqnew.time,real(squeeze(freqold.fourierspctrm(1,1,freq,:))))
plot(freqnew.time,real(squeeze(freqnew.fourierspctrm(1,1,freq,:))),'r')
legend
axis tight
xlabel('frequency (Hz)')
ylabel('amplitude')
title('both overlayed, using freq axis of hilbert')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% FOURIER specest_wavelet VS specest_tfr
cfg = [];
cfg.trials     = 1;
cfg.keeptrials = 'yes';
cfg.keeptapers = 'yes';
cfg.pad        = 6;
cfg.foi        = 2:5:20;
cfg.toi        = 0:0.05:4; %data.time{1}(1):(1/data.fsample):data.time{1}(end);
cfg.output     = 'fourier';
cfg.calcdof    = 'yes';
cfg.correctt_ftimwin = 'no';
% new style
cfg.method     = 'wavelet';
tic; freqwav = ft_freqanalysis(cfg,data); toc
% old style
cfg.foi = freqwav.freq;
cfg.method     = 'tfr';
tic; freqtfr = ft_freqanalysis(cfg,data); toc

% plot first trial first channel, single timepoint
timepoint = 50;
figure
subplot(4,1,1)
plot(freqtfr.freq,real(squeeze(freqtfr.fourierspctrm(1,1,:,timepoint))))
axis tight
xlabel('frequency (Hz)')
ylabel('amplitude')
title('tfr implementation')
subplot(4,1,2)
plot(freqwav.freq,real(squeeze(freqwav.fourierspctrm(1,1,:,timepoint))))
axis tight
xlabel('frequency (Hz)')
ylabel('amplitude')
title('wavelet implementation')
subplot(4,1,3)
plot(freqtfr.freq,angle(squeeze(freqtfr.fourierspctrm(1,1,:,timepoint))))
axis tight
xlabel('frequency (Hz)')
ylabel('phase')
title('tfr implementation')
subplot(4,1,4)
plot(freqwav.freq,angle(squeeze(freqwav.fourierspctrm(1,1,:,timepoint))))
axis tight
xlabel('frequency (Hz)')
ylabel('phase')
title('wavelet implementation')


text:
% This function add a full matrix with a sparse matrix

function out = AddSpMatMat(w1,spMat, w2, Mat, sp_elements_only)
[m,n] = size(spMat);

if w1==-1
    spMat = -spMat;
elseif w1~=1
    spMat = spMat * w1;
end
if w2==-1
    Mat = -Mat;
elseif w2~=1
    Mat = Mat * w2;
end

out = Mat;
if 0
    idx = find(spMat);
    out(idx) = double(out(idx)) + spMat(idx);
elseif 1
    if m<n
        non_zero_col = find(sum(abs(spMat),1)>0);
        out(:,non_zero_col) = out(:,non_zero_col) + full(spMat(:,non_zero_col));
    else
        non_zero_row = find(sum(abs(spMat),2)>0);
        out(non_zero_row,:) = out(non_zero_row,:) + full(spMat(non_zero_row,:));
    end
else
    out = full(spMat) + Mat;
end


text:
function mean = genlogistic_mean ( a, b, c )

%*****************************************************************************80
%
%% GENLOGISTIC_MEAN returns the mean of the Generalized Logistic PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < B,
%    0.0 < C.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = a + b * ( euler_constant ( ) + digamma ( c ) );

  return
end


text:
%% Contents
%  buzcode\analysis\

%% lfp
%    bz_Filter:             Filter - Filter samples.
%    bz_DownsampleLFP:      takes a buzcode lfp structure and returns a buzcode lfp structure downsampled
%    bz_LFPPowerDist:       calculates the power distribution of an LFP signal.
%    bz_LFPSpecToExternalVar: calculates the relationship of the LFP spectrum to an external variable
%    bz_PowerSpectrumSlope:  calculates the slope of the power spectrum
%    bz_Comodulogram        calculates the power-power comodulogram for an lfp file.
%    bz_GradDescCluster     clusters recording sites given a pairwise similarity matrix using gradient descent
%
%%% CrossFrecuencyCoupling
%    bz_ModIndex            calculates LFP phase-amplitude modulation index
%    bz_PhaseAmplitudeDist  calculates mean amplitude of higher freq bands at a the phase for a given lower freq band signal
%    PhaseAmpCouplingByAmp  calculates phase amplitude coupling - NOT BUZCODE FORMAT
%
%%% CurrentSourceDensity
%   bz_CSD                  Calculates the 1D approximation of CSD from a linear array of LFPs
%   bz_eventCSD             Calculates event-triggered CSD map from a linear array of LFPs
%
%%% SharpWaveRipples
%   bz_GetBestRippleChan    Not working ???
%   bz_RippleStats          Compute descriptive stats for ripples
%   bz_PlotRippleStats      Plot descriptive stats for ripples
%
%%% SpectralAnalyses
%   bz_WaveSpec             calculates the wavelet transform of a signal
%   bz_MTCoherogram         Compute LFP coherogram by multi-taper estimation - NOT BUZCODE FORTMAT
%
%
%%% LFPConnectivity
%
%
%% lfp_spikes

%% spikes

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function mono_upto_next_grevlex_test ( )

%*****************************************************************************80
%
%% MONO_UPTO_NEXT_GREVLEX_TEST tests MONO_UPTO_NEXT_GREVLEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONO_UPTO_NEXT_GREVLEX_TEST\n' );
  fprintf ( 1, '  MONO_UPTO_NEXT_GREVLEX can list the monomials\n' );
  fprintf ( 1, '  in M variables, of total degree up to N,\n' );
  fprintf ( 1, '  one at a time, in graded reverse lexicographic order.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We start the process with (0,0,...0,0).\n' );
  fprintf ( 1, '  The process ends with (N,0,...,0,0)\n' );

  test_num = 2;
  n_test = [ 4, 3 ];
  m_test = [ 3, 4 ];

  for test = 1 : test_num

    n = n_test(test);
    m = m_test(test);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Let M = %d\n', m );
    fprintf ( 1, '      N = %d\n', n );
    fprintf ( 1, '\n' );

    x = zeros ( m, 1 );
    i = 1;

    while ( 1 )

      fprintf ( 1, '  %2d:', i );
      for j = 1 : m
        fprintf ( 1, '  %1d', x(j) );
      end
      fprintf ( 1, '\n' );

      if ( x(1) == n )
        break
      end

      x = mono_upto_next_grevlex ( m, n, x );
      i = i + 1;

    end

  end

  return
end


text:
function dist = pntdist(p1, p2)

% PNTDIST returns the euclidian distance between two points
%
%  [dist] = pntdist(pnt1, pnt2)
%
% where pnt1 and pnt2 must be Npnt x 3
% or either one can be Npnt x 1

% Copyright (C) 2002, Robert Oostenveld
%
% $Log: pntdist.m,v $
% Revision 1.2  2003/03/04 21:46:19  roberto
% added CVS log entry and synchronized all copyright labels
%

if size(p1,1)==1
  p1 = repmat(p1, size(p2,1), 1);
elseif size(p2,1)==1
  p2 = repmat(p2, size(p1,1), 1);
end

dist = sqrt(sum((p1-p2).^2, 2));


text:
function func = ft_getuserfun(func, prefix)

% FT_GETUSERFUN will search the MATLAB path for a function with the
% appropriate name, and return a function handle to the function.
% Considered are, in this order:
%  - the name itself, i.e. you get exactly the same func back as you put in;
%  - the name with the specified prefix;
%  - the name with 'ft_' and the specified prefix.
%
% For example, calling FT_GETUSERFUN('general', 'trialfun') might return a
% function named 'general', 'trialfun_general', or 'ft_trialfun_general',
% whichever of those is found first and is not a compatibility wrapper.
%
% func can be a function handle, in which case it is returned as-is.
%
% If no appropriate function is found, the empty array [] will be returned.

% Copyright (C) 2012, Eelke Spaak
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if isa(func, 'function_handle')
  % treat function handle as-is
elseif isfunction(func) && ~iscompatwrapper(func)
  func = str2func(func);
elseif isfunction([prefix '_' func]) && ~iscompatwrapper([prefix '_' func])
  func = str2func([prefix '_' func]);
elseif isfunction(['ft_' prefix '_' func])
  func = str2func(['ft_' prefix '_' func]);
else
  ft_warning(['no function by the name ''' func ''', ''' prefix '_' func...
    ''', or ''ft_' prefix '_' func ''' could be found']);
  func = [];
end


text:
function dy = lotka(t,y,a,b,d,g)
dy = [
    a*y(1) - b*y(1)*y(2);
    d*y(1)*y(2) - g*y(2);
    ];

author: eurika-kaiser

repo: SINDY-MPC


text:
function [model] = tapas_linear_model(model, pars)
%% Define the model.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

model = struct('graph', []);
model.graph = cell(4, 1);

for i = 1:4
    model.graph{i} = struct('llh', [], 'htheta', []);
end

model.graph{1}.llh = @tapas_linear_llh;
model.graph{2}.llh = @tapas_linear_hier_llh;
model.graph{3}.llh = @tapas_linear_hier_llh;
model.graph{4}.llh = [];

model.graph{1}.htheta = struct('pe', 0.5, 'T', pars.T);
model.graph{2}.htheta = struct('pe', 0.5, 'T', ones(size(pars.T)));
model.graph{3}.htheta = struct('pe', 0.5, 'T', ones(size(pars.T)));

% The last level is a dummy used to store the hyperpriors.

model.graph{4}.htheta = struct('y', [], 'u', []);

end

author: translationalneuromodeling

repo: tapas


text:
function [ cx, cy ] = cswap ( n, cx, incx, cy, incy )

%*****************************************************************************80
%
%% CSWAP interchanges two complex vectors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, complex CX(*), one of the vectors to swap.
%
%    Input, integer INCX, the increment between successive entries of CX.
%
%    Input, complex CY(*), one of the vectors to swap.
%
%    Input, integer INCY, the increment between successive elements of CY.
%
%    Output, complex CX(*), the input vector CX, with some elements swapped.
%
%    Input, integer INCX, the increment between successive entries of CX.
%
%    Output, complex CY(*), the input vector CY, with some elements swapped.
%
  temp                    = cx(1:incx:1+(n-1)*incx);
  cx(1:incx:1+(n-1)*incx) = cy(1:incy:1+(n-1)*incy);
  cy(1:incy:1+(n-1)*incy) = temp;

  return
end


text:
function ssmult_demo
%SSMULT_DEMO simple demo for ssmult.
%
% Example:
%   ssmult_demo
%
% See also ssmult, ssmult_unsorted, ssmultsym, sstest, sstest2.

% Copyright 2007-2009, Timothy A. Davis, University of Florida

type ssmult_demo
load west0479
A = west0479 ;
B = sprand (A) ;
C = A*B ;
D = ssmult (A,B) ;
err = norm (C-D,1) / norm (C,1) ;
fprintf ('ssmult west0479 error: %g\n', err) ;


text:
% Demo for aggregate channel features object detector on Inria dataset.
%
% (1) Download data and helper routines from Caltech Peds Website
%  www.vision.caltech.edu/Image_Datasets/CaltechPedestrians/
%  (1a) Download INRIA files: set00.tar, set01.tar, and annotations.zip
%  (1b) Copy above three files to dataDir/ and untar/unzip contents
%  (1c) Download evaluation code (routines necessary for extracting images)
% (2) Set dataDir/ variable below to point to location of INRIA data.
% (3) Launch "matlabpool open" for faster training if available.
% (4) Run demo script and enjoy your newly minted fast ped detector!
%
% Note: pre-trained model files are provided (delete to re-train).
% Re-training may give slightly variable results on different machines.
%
% Piotr's Image&Video Toolbox      Version 3.22
% Copyright 2013 Piotr Dollar & Ron Appel.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

%% extract training and testing images and ground truth
dataDir = 'D:\code\research\detectorAcf\data\Inria\';
for s=1:2
  if(s==1), set='00'; type='train'; else set='01'; type='test'; end
  if(exist([dataDir type '/posGt'],'dir')), continue; end
  seqIo([dataDir 'set' set '/V000'],'toImgs',[dataDir type '/pos']);
  seqIo([dataDir 'set' set '/V001'],'toImgs',[dataDir type '/neg']);
  V=vbb('vbbLoad',[dataDir 'annotations/set' set '/V000']);
  vbb('vbbToFiles',V,[dataDir type '/posGt']);
end

%% set up opts for training detector (see acfTrain)
opts=acfTrain(); opts.modelDs=[100 41]; opts.modelDsPad=[128 64];
opts.posGtDir=[dataDir 'train/posGt']; opts.nWeak=[32 128 512 2048];
opts.posImgDir=[dataDir 'train/pos']; opts.pJitter=struct('flip',1);
opts.negImgDir=[dataDir 'train/neg']; opts.pBoost.pTree.fracFtrs=1/16;
opts.pLoad={'squarify',{3,.41}}; opts.name='models/AcfInria';

%% train detector (see acfTrain)
detector = acfTrain( opts );

%% modify detector (see acfModify)
detector = acfModify(detector,'cascThr',-1,'cascCal',0);

%% run detector on a sample image (see acfDetect)
imgNms=bbGt('getFiles',{[dataDir 'test/pos']});
I=imread(imgNms{1}); tic, bbs=acfDetect(I,detector); toc
figure(1); im(I); bbApply('draw',bbs); pause(.1);

%% test detector and plot roc (see acfTest)
[miss,~,gt,dt]=acfTest('name',opts.name,'imgDir',[dataDir 'test/pos'],...
  'gtDir',[dataDir 'test/posGt'],'pLoad',opts.pLoad,'show',2);

%% optional timing test for detector (should be ~30 fps)
if( 0 )
  detector1=acfModify(detector,'pad',[0 0]); n=60; Is=cell(1,n);
  for i=1:n, Is{i}=imResample(imread(imgNms{i}),[480 640]); end
  tic, for i=1:n, acfDetect(Is{i},detector1); end;
  fprintf('Detector runs at %.2f fps on 640x480 images.\n',n/toc);
end

%% optionally show top false positives ('type' can be 'fp','fn','tp','dt')
if( 0 ), bbGt('cropRes',gt,dt,imgNms,'type','fn','n',50,...
    'show',3,'dims',opts.modelDs([2 1])); end

author: zouchuhang

repo: LayoutNet


text:
function varargout = getsamples(filename, varargin)
% Read samples from acq file
%   filename - acq file
%   varargin - channel names
% Returns column vectors of all channel data samples of requested channels.
% Example:
%   [SourceR, EMG_L, EMG_R] = getsamples('myfile.acq', 'SourceR', 'EMG_L', 'EMG_R');
% Note:
% Developed and tested for version 3.5 AcqKnowledge files, but should
% probably work with other versions as well.
% Multiple sampling rates are not supported.
% No copyright. Feel free to use this code snippet in any way you want.

n = length(varargin);        % Number of channels to read

% Step through each file section until Channel Data Section is reached.
fp = fopen(filename, 'r');
% Graph Header Section
fseek(fp, 6, 'bof');
lExtItemHeaderLen = fread(fp, 1, 'int32');
nChannels = fread(fp, 1, 'int16');
% Per Channel Data Section
sectionstart = lExtItemHeaderLen;
fseek(fp, sectionstart, 'bof');
lChanHeaderLen = fread(fp, 1, 'int32');
fseek(fp, sectionstart+88, 'bof');
lBufLength = fread(fp, 1, 'int32');     % Number of samples, same for all channels
for i=1:nChannels
    sectionstart = lExtItemHeaderLen + (i-1)*lChanHeaderLen;
    fseek(fp, sectionstart+6, 'bof');
    szCommentText(i,1:40) = fscanf(fp, '%40c', 1);  % This is the channel label
    fseek(fp, sectionstart+92, 'bof');
    dAmplScale(i) = fread(fp, 1, 'double');
    dAmplOffset(i) = fread(fp, 1, 'double');
end
% Foreign Data Section
fseek(fp, lExtItemHeaderLen+lChanHeaderLen*nChannels, 'bof');
nLength = fread(fp, 1, 'int16');
% Per Channel Data Types Section
sectionstart = lExtItemHeaderLen + lChanHeaderLen*nChannels + nLength;
fseek(fp, sectionstart, 'bof');
for i=1:nChannels
    nSize(i) = fread(fp, 1, 'int16');
    nType(i) = fread(fp, 1, 'int16');
end
choffset = cumsum([0 nSize(1:end-1)]);  % Interleave offset
% Channel Data Section
% First find indexes of channels to read by searching channel labels, ...
for k=1:n
    for i=1:nChannels
        chlabel = char(varargin{k});
        if strncmp(chlabel, szCommentText(i,:), length(chlabel))
            chindex(k) = i;
            break;
        end
    end
end
% ...then read samples, one channel at a time.
sectionstart = lExtItemHeaderLen + lChanHeaderLen*nChannels + nLength + 4*nChannels;
fseek(fp, sectionstart, 'bof');
blocksize = sum(nSize);
for k=1:n
    fseek(fp, sectionstart + choffset(chindex(k)), 'bof');
    switch nType(chindex(k))
        case 1  % double
            varargout(k) = {fread(fp, lBufLength, 'double', blocksize-8)};
        case 2  % int
            varargout(k)= {dAmplScale(chindex(k)) * fread(fp, lBufLength, 'int16', blocksize-2) + ...
                                    dAmplOffset(chindex(k))};
    end
end

fclose(fp);


text:
function out=declus_wintec(mCatalog, eMethod)
    % DECLUS_WINTEC execute declustering by windowing technique
    %
    % DECLUS_WINTEC(catalog, nMethod)
    %
    % Incoming variables:
    % mCatalog : EQ catalog in ZMAP format
    % eMethod  : a choice from DeclusterWindowingMethods
    %
    %
    % out : struct of declustered details [output from CALC_DECLUSTER]
    % J. Woessner, woessner@seismo.ifg.ethz.ch
    % updated: 14.08.02
    %
    % see DECLUS_INP, CALC_DECLUSTER

    report_this_filefun();

    %%% Decluster catalog using window technique
    out = struct;
    % [out.mCatDecluster, out.mCatAfter, out.vCluster, out.vCl, out.vMain] = calc_decluster(mCatalog,eMethod);
    [out.declusteredCatalog, out.aftershockCatalog, out.aftershockClusterIdx, out.allClusterIdx, out.mainshockClusterIdx] = calc_decluster(mCatalog,eMethod);
    %vSel = (out.vMain(:,1) > 0); % Selects mainshocks of clusters
    %mCluster = mCatalog.subset(vSel);
    out.description = ["Decluster results have been written to gk_decluster_output with the following fields";...
    	"  declusteredCatalog: Declustered earthquake catalog (with mainshocks, no fore/aftershocks)";...
        "  aftershockCatalog: Catalog of aftershocks (and foreshocks)";...
        "  aftershockClusterIdx : Vector indicating only aftershocks/foreshocks in cluster using a cluster number (background seismicity == 0)";...
        "  allClusterIdx: Vector indicating all events in clusters using a cluster number (background seismicity == 0)";...
        "  mainshockClusterIdx: index into original catalog, indicating mainshocks, where any non-zero value is the cluster number"];

    %%% Plot comparison to window length
    %[vMags, vClusTime, vDist]= plot_cluscomp(vMain, vCluster, mCatalog, eMethod);
    % plot_cluscomp(out.mainshockClusterIdx, out.aftershockClusterIdx, mCatalog, eMethod); % output arguments weren't being used
    return % CGR - Sep 1 2020

    %%% Plot seismicity map, clusters and mainshocks
    replaceMainCatalog(out.mCatDecluster);
    zmap_update_displays();
    plot(mCluster.Longitude, mCluster.Latitude,'m+');

    describe_clusters(mCatalog, out)
    plot_mag_histogram(catalog, mCatAfter)

end

function describe_clusters(mCatalog, out)
    %%% Calculate moment release [local variables], only used for description
    [fMomentCluster, vMomentCluster] = calc_moment(out.mCatAfter);

    [fMomentorg, vMomentorg] = calc_moment(out.mCatalog);

    fMomentpercentage = 100*fMomentCluster/fMomentorg;

    fEventpercentage = 100*out.mCatAfter.Count/mCatalog.Count; % Percentage of events in clusters

    %% Setup message box
    sInfost1 = sprintf(...
        [' The declustering found %d clusters of earthquakes, a total of %d (%g\%) events out of %d. ',...
        ' The map window now displays the declustered catalog containing %d events as blue dots.', ....
        ' The individual clusters are displayed as magenta pluses. The seismic moment released by the clusters'...
        ' is %g Nm which is about %g\% of the total seismic moment (%g Nm) of the catalog.'],...
        max(out.vMain) , out.mCatAfter.Count, fEventpercentage, mCatalog.Count,...
        out.mCatDecluster.Count,...
        fMomentCluster, fMomentpercentage, fMomentorg);

    msgbox(sInfost1,'Declustering Information')
end

function plot_mag_histogram(origCatalog, declusteredCatalog)
    %%% Plotting magnitude histogram
    if exist('hd1_declus_wintec','var') && ishandle(hd1_declus_wintec)
        set(0,'Currentfigure',hd1_declus_wintec);
        disp('Figure exists');
    else
        hd1_win_fig=figure('tag','fig_declus_wintec','Name','Histogram',...
            'Units','normalized','Nextplot','add',...
            'Numbertitle','off','Position',[0.4 0.2 .4 .6],'Menubar','none');
    end

    set(gca,'tag','ax_declus_wintec_mag','Nextplot','replace','box','on','Xticklabel', [0 10 100]);
    axs1=findobj('tag','ax_declus_wintec_mag');
    axes(axs1(1));
    maxMagLimit=max(origCatalog.Magnitude);
    magEdges = 0:0.1:maxMagLimit;
    histogram(declusteredCatalog.Magnitude, magEdges);
    xlabel('Magnitude (events of all clusters)');
end

author: CelsoReyes


text:
function pass = test_hermpoly( prefs )
% Test the HERMPOLY command.

if ( nargin < 1 )
    prefs = chebfunpref();
end
tol = 1e6*prefs.techPrefs.chebfuneps;

% Physicists' Hermite polynomials
h0 = hermpoly(0);
err = norm( h0 - chebfun(@(x) 1+0*x, [-inf,inf] ) );
pass(1) = err < tol;

h1 = hermpoly(1);
err = norm( h1 - chebfun(@(x) 2*x, [-inf,inf] ) );
pass(2) = err < tol;

h2 = hermpoly(2);
err = norm( h2 - chebfun(@(x) 4*x.^2-2, [-inf,inf] ) );
pass(3) = err < tol;

h3 = hermpoly(3);
x = linspace(-1,1);
err = norm( feval(h3,x) - 8*x.^3+12*x );
pass(4) = err < tol;

h4 = hermpoly(4);
x = linspace(-1,1);
err = norm( feval(h4,x) - (16*x.^4-48*x.^2+12) );
pass(5) = err < tol;

h5 = hermpoly(5);
x = linspace(-1,1);
err = norm( feval(h5,x) - (32*x.^5-160*x.^3+120*x) );
pass(6) = err < 10*tol;

h6 = hermpoly(6);
x = linspace(-1,1);
err = norm( feval(h6,x) - (64*x.^6-480*x.^4+720*x.^2-120) );
pass(7) = err < 100*tol;

% probabilists' Hermite polynomials:
h0 = hermpoly(0, 'prob');
err = norm( h0 - chebfun(@(x) 1+0*x, [-inf,inf] ) );
pass(8) = err < tol;

h1 = hermpoly(1, 'prob');
err = norm( h1 - chebfun(@(x) x, [-inf,inf] ) );
pass(9) = err < tol;

h2 = hermpoly(2, 'prob');
g2 = chebfun(@(x) x.^2-1, [-inf,inf] );
err = norm( h2 - chebfun(@(x) x.^2-1, [-inf,inf] ) );
pass(10) = err < tol;

h3 = hermpoly(3, 'prob');
x = linspace(-1,1);
err = norm( feval(h3,x) - (x.^3-3*x) );
pass(11) = err < tol;

h4 = hermpoly(4, 'prob');
x = linspace(-1,1);
err = norm( feval(h4,x) - (x.^4-6*x.^2+3) );
pass(12) = err < tol;

h5 = hermpoly(5, 'prob');
x = linspace(-1,1);
err = norm( feval(h5,x) - (x.^5-10*x.^3+15*x) );
pass(13) = err < tol;

h6 = hermpoly(6, 'prob');
x = linspace(-1,1);
norm( feval(h6,x) - (x.^6-15*x.^4+45*x.^2-15) );
pass(14) = err < 10*tol;

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Population = GDE3_EnvironmentalSelection(Population,Offspring,N)
% The environmental selection of CCGDE3

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    % Select the offsprings dominating their corresponding parents firstly,
    % and then select the offspring population non-dominated with the
    % parent population
    PopObj = Population.objs;
    OffObj = Offspring.objs;
    % The offsprings which can replace its parent
    updated = all(PopObj>=OffObj,2);
    % The offsprings which can add to the population
    selected = any(PopObj<OffObj,2) & any(PopObj>OffObj,2);
    % Update the population
    Population(updated) = Offspring(updated);
    Population          = [Population,Offspring(selected)];

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,N);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% (Internal) Plot the area under the ROC curve
%
%       plot_auc( areas, pvalues, strTitle, strLegend, bPrint, fig_hdl)
%
%
% Arguments:
%
%      + areas: string to parse
%
%      + pvalues: string to parse
%
%      + strTitle: string to parse
%
%      + strLegend: string to parse
%
%      + bPrint: string to parse
%
%      + fig_hdl: string to parse
%
% Example:
%
%
% See also ECGwrapper
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate  : 30/7/2014
% Last update: 30/7/2014
% Copyright 2008-2015
%
function plot_auc( areas, pvalues, strTitle, strLegend, bPrint, fig_hdl)

if( nargin < 4  )
    strLegend = [];
end

if( nargin < 3  )
    strTitle = [];
end

if( nargin < 5 || isempty(bPrint) )
    bPrint = false;
end

if( size(areas,1) ~= 3 )
    error('areas should have lower CI - mean Area - upper CI in the rows')
end

if( nargin < 6 || isempty(fig_hdl) )
    fig_hdl = figure();
else
    fig_hdl = figure(fig_hdl);
%     clf(fig_hdl);
end

areas2plot = size(areas,2);
Xvalues = [ ];
Yvalues = [ ];
significance_levels = [0.05 0.01 0.001];
significance_levels_str = {'0.05' '0.01' '0.001'};

box_dist = 1;
box_width = box_dist* 2/3;
box_loc = 1:box_dist:box_dist*areas2plot;
box_height_k = 0.75;
colors = my_colormap(areas2plot);
areas_min = min(min(areas));
areas_max = max(max(areas));
areas_range = areas_max - areas_min;

if(bPrint)
    strFontName = 'LMRoman12';
    FontSizeTitle = 22;
    FontSizeLabel = 20;
    FontSizeTick = 18;
else
    strFontName = 'Helvetica';
    FontSizeTitle = 12;
    FontSizeLabel = 12;
    FontSizeTick = 10;
end

axes_hdl = gca();
% cla(axes_hdl);

xlim([0 (box_dist*areas2plot)+box_dist])
% ylim([areas_min - 0.1*areas_range  areas_max + 0.2*areas_range ])

hold on

for ii = 1:areas2plot

    box_height = box_height_k * (areas(3,ii) - areas(1,ii));
    this_colour = colors( 1+rem(ii-1, size(colors,1)),:);
    plot( box_loc(ii) + [-box_width/4 box_width/4], [areas(3,ii) areas(3,ii)], 'Color', this_colour, 'LineWidth', 2 );
    plot( box_loc(ii) + [-box_width/4 box_width/4], [areas(1,ii) areas(1,ii)], 'Color', this_colour, 'LineWidth', 2 )
    plot( [box_loc(ii) box_loc(ii)], [areas(3,ii) areas(1,ii)], 'Color', this_colour, 'LineWidth', 2 )
    fill( box_loc(ii) + [-box_width/2 -box_width/2 box_width/2 box_width/2], areas(2,ii) + [ -box_height/2 box_height/2 box_height/2 -box_height/2], this_colour, 'EdgeColor',  this_colour );
    text( box_loc(ii), areas(2,ii), num2str(areas(2,ii), '%3.2f'), 'HorizontalAlignment', 'center', 'FontName', strFontName, 'FontSize', FontSizeTick, 'Color', (1-this_colour));

end

for ii = 1:(areas2plot-1)
    for jj = ii+1:areas2plot
        if( pvalues(ii,jj) < 0.05 )
            y_loc = max(areas(3,ii), areas(3,jj)) + (0.05*areas_range);
            if( bPrint )
                plot( [box_loc(ii) box_loc(ii) box_loc(jj) box_loc(jj)], y_loc + (0.02*areas_range) * [0 1 1 0], 'LineWidth', 1 , 'Color', [0 0 0] );
            else
                plot( [box_loc(ii) box_loc(ii) box_loc(jj) box_loc(jj)], y_loc + (0.02*areas_range) * [0 1 1 0], 'LineWidth', 2 , 'Color', [0 0 0] );
            end
        end
    end

    % overprint the p-value
    for jj = ii+1:areas2plot
        aux_idx = find(pvalues(ii,jj) < significance_levels, 1, 'last');
        if( ~isempty(aux_idx) )
            y_loc = max(areas(3,ii), areas(3,jj)) + (0.05*areas_range);

            if( bPrint )
                aux_hdl = text( box_loc(ii)+(box_loc(jj) - box_loc(ii))/2, y_loc + (0.07*areas_range), [ 'p < ' significance_levels_str{aux_idx} ], 'HorizontalAlignment', 'center', 'FontName', strFontName, 'FontSize', FontSizeTick, 'Visible', 'off');
                aux_extent = get(aux_hdl, 'Extent');
                delete(aux_hdl)
                aux_extent(1) = aux_extent(1) + 0.1*aux_extent(3);
                aux_extent(3) = 0.8*aux_extent(3);
                aux_extent(2) = aux_extent(2) + 0.25*aux_extent(4);
                aux_extent(4) = 0.7*aux_extent(4);
                fill( aux_extent(1)+ [ 0 0 aux_extent(3) aux_extent(3) ], aux_extent(2)+ [ 0 aux_extent(4) aux_extent(4) 0 ], [1 1 1], 'EdgeColor', [1 1 1]  );
                text( box_loc(ii)+(box_loc(jj) - box_loc(ii))/2, y_loc + (0.07*areas_range), [ 'p < ' significance_levels_str{aux_idx} ], 'HorizontalAlignment', 'center', 'FontName', strFontName, 'FontSize', FontSizeTick);
            else
                text( box_loc(ii)+(box_loc(jj) - box_loc(ii))/2, y_loc + (0.07*areas_range), [ 'p < ' significance_levels_str{aux_idx} ], 'HorizontalAlignment', 'center', 'FontName', strFontName, 'FontSize', FontSizeTick, 'BackgroundColor', [1 1 1]);
            end
        end
    end

end

hold off

set(axes_hdl, 'Box', 'off' );
set(axes_hdl, 'Xtick', box_loc );
set(axes_hdl, 'XtickLabel', strLegend );
rotateticklabel(axes_hdl, 20);
set(axes_hdl, 'Ytick', linspace(areas_min, areas_max, 5) );
set(axes_hdl, 'YtickLabel', num2str(colvec(linspace(areas_min, areas_max, 5)), '%3.2f') );

set(axes_hdl, 'FontName', strFontName );
set(axes_hdl, 'FontSize', FontSizeTick );

title(strTitle, 'FontName', strFontName, 'FontSize', FontSizeTitle)

ylabel('AUC', 'FontName', strFontName, 'FontSize', FontSizeLabel)

author: marianux

repo: ecg-kit

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


text:
function [dVVectECI] = getNTW2ECIdvVect(dVVectNTW, rVect, vVect)
%getNTWdvVect Summary of this function goes here
%   Detailed explanation goes here
    rVect = reshape(rVect, 3,1);
    vVect = reshape(vVect, 3,1);
    dVVectNTW = reshape(dVVectNTW, 3,1);

    tHat = vVect/norm(vVect);

    wHatC = crossARH(rVect,vVect);
    wHat = wHatC/norm(wHatC);

    nHatC = crossARH(tHat,wHat);
    nHat = nHatC/norm(nHatC);

    ECI2TWNRotMat = [tHat,wHat,nHat];
%     TWN2ECIRotMat = inv(ECI2TWNRotMat);
    TWN2ECIRotMat = ECI2TWNRotMat'; %rotation matrix, inv = transpose

    dVVectECI = TWN2ECIRotMat \ dVVectNTW;
end


text:
function Ft=frametight(F);
%FRAMETIGHT  Construct the canonical tight frame
%   Usage: Ft=frametight(F);
%
%   `Ft=frametight(F)` returns the canonical tight frame of *F*.
%
%   The canonical tight frame can be used to get perfect reconstruction if
%   it is used for both analysis and synthesis. This is demonstrated in the
%   following example:::
%
%     % Create a frame and its canonical tight
%     F=frame('dgt','hamming',32,64);
%     Ft=frametight(F);
%
%     % Compute the frame coefficients and test for perfect
%     % reconstruction
%     f=gspi;
%     c=frana(Ft,f);
%     r=frsyn(Ft,c);
%     norm(r(1:length(f))-f)
%
%   See also: frame, framepair, framedual

complainif_notenoughargs(nargin,1,'FRAMETIGHT');
complainif_notvalidframeobj(F,'FRAMETIGHT');

% Default operation, works for a lot of frames
Ft=F;

% Handle the windowed transforms
switch(F.type)
  case {'dgt','dgtreal','dwilt','wmdct','filterbank','ufilterbank',...
        'nsdgt','unsdgt','nsdgtreal','unsdgtreal'}

    Ft=frame(F.type,{'tight',F.g},F.origargs{2:end});

  case {'filterbankreal','ufilterbankreal'}
    Ft=frame(F.type,{'realtight',F.g},F.origargs{2:end});

  case 'gen'
    [U,sv,V] = svd(F.g,'econ');
    Ft=frame('gen',U*V');

  case 'tensor'
    for ii=1:F.Nframes
        tight_frames{ii}=frametight(F.frames{ii});
    end;
    F=frame('tensor',tight_frames{:});

  case 'fusion'
    tight_w=1./F.w;
    for ii=1:F.Nframes
        tight_frames{ii}=frametight(F.frames{ii});
    end;
    Ft=frame('fusion',tight_w,tight_frames{:});

  case 'ufwt'
    % The canonical tight made from ufwt might not keep the iterated
    % filterbank structure
    [g,a] = wfbt2filterbank({F.g,F.J,'dwt'});
    g = comp_filterbankscale(g,a,F.flags.scaling);

    Ft = frametight(frame('filterbank',g,ones(numel(g),1),numel(g)));

  case 'uwfbt'
    % The canonical tight made from uwfbt might not keep the iterated
    % filterbank structure
    [g,a] = wfbt2filterbank(F.g,F.J);
    g = comp_filterbankscale(g,a,F.flags.scaling);

    Ft = frametight(frame('filterbank',g,ones(numel(g),1),numel(g)));

  case 'uwpfbt'
    % The canonical tight made from uwpfbt might not keep the iterated
    % filterbank structure
    [g, a] = wpfbt2filterbank(F.g,F.flags.interscaling);
    g = comp_filterbankscale(g,a,F.flags.scaling);

    Ft = frametight(frame('filterbank',g,ones(numel(g),1),numel(g)));

  case 'fwt'
    is_basis = abs(sum(1./F.g.a)-1)<1e-6;
    is_tight = F.info.istight;

    if is_basis && is_tight
        Ft = F;
    else
        error(['%s: Cannot create the canonical tight frame with the ',...
               'same structure. Consider casting the system to an ',...
               'uniform filterbank.'],...
               upper(mfilename));
    end

  case 'wfbt'
    is_basis = all(cellfun(@(nEl) abs(sum(1./nEl.a)-1)<1e-6,F.g.nodes));
    is_tight = F.info.istight;

    if is_basis && is_tight
        Ft = F;
    else
        error(['%s: Cannot create the canonical tight frame with the ',...
               'same structure. Consider casting the system to an ',...
               'uniform filterbank.'],...
               upper(mfilename));
    end

   case 'wpfbt'
     % WPFBT is too wierd.
     error(['%s: Canonical tight frame of wpfbt might not keep the ',...
           'same structure. '],upper(mfilename))

end;

switch(F.type)
    case {'ufwt','uwfbt','uwpfbt'}
      warning(sprintf(['%s: The canonical tight system does not preserve ',...
                     'the iterated filterbank structure.'],...
                      upper(mfilename)));
end

% Treat the fixed length frames
if isfield(F,'fixedlength') && F.fixedlength && isfield(F,'L')
   Ft = frameaccel(Ft,F.L);
   Ft.fixedlength = 1;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function combo_test36 ( )

%*****************************************************************************80
%
%% COMBO_TEST36 tests SETPART_ENUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST36\n' );
  fprintf ( 1, '  Set partitions:\n' );
  fprintf ( 1, ' \n' );
  fprintf ( 1, '  SETPART_ENUM enumerates,\n' );
  fprintf ( 1, ' \n' );
%
%  Enumerate.
%
  for n = 1 : 6
    npart = setpart_enum ( n );
    fprintf ( 1, '  %4d  %4d\n', n, npart );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function s = tt_size2str(sz)
%TT_SIZE2STR Convert size to a string that can be printed.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if isempty(sz)
    s = sprintf('[empty tensor]');
    return;
end

if numel(sz) == 1
    s = sprintf('%d',sz);
else
    s = [sprintf('%d x ',sz(1:end-1)) sprintf('%d', sz(end)) ];
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt