function [dm] = directmethods(y,lags,options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'directmethods' generates local projeciton IRF and direct forecasts

% Core Inputs:
% - y, data columns variables
% - lags, lag order of the VAR

% Additonal Inputs collected options
% (...) see below
% See the Hitchhiker's guide for more details.
% https://github.com/naffe15/BVAR_/blob/master/HitchhikerGuide_.pdf

% Filippo Ferroni, 27/02/2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 2
    error('the BVAR toolbox needs at least two inputs: data and number of lags');
end
if lags < 1
    error('lags cannot be zero or negative');
end
% number of observable variables
[T,ny]       = size(y);

%********************************************************
%* DEFAULT SETTINGS
%********************************************************
% Control random number generator
if isOctave == 0
    isMatlab = 1;
    rng('default');
    rng(999);
else
    isMatlab = 0;
    % pkg load optim
     pkg load statistics
    randn('state',999);
    rand('state',999);
end

hor          = 24;
conf_sig     = 0.9;
controls_    = 0;
robust_se_   = 1; % by default robust SE
proxy_       = 0;
noconstant   = 0;
ns           = 0;
Q            = eye(ny);
dummy        = 0;
K            = 5000;
max_prior_tau_ = 0;
lb             = -1e10;
ub             = 1e10;
clags         = 0;

%********************************************************
%* CUSTOMIZED SETTINGS
%********************************************************
if nargin>2
    %======================================================================
    % Various options
    %======================================================================
    if isfield(options,'hor') ==1 % horizon of IRF and Forecasts
        hor = options.hor;
    end
    if isfield(options,'conf_sig') ==1 % CI of OLS estimation of LP and DF
        conf_sig = options.conf_sig;
    end
    if isfield(options,'Q') == 1 % orthonormal matrix
        Q = options.Q;
    end
    if isfield(options,'K') == 1 % # of draws for the BLP
        K = options.K;
    end
    %======================================================================
    % options: adding controls
    %======================================================================
    if isfield(options,'controls') ==1
        controls_  = 1;
        controls   = options.controls;
        if T~=size(controls,1)
            error('Control variables and observables must have the same time length')
        end
        if isfield(options,'clags') ==1
            clags = options.clags;
        end
    end
    %======================================================================
    % options: adding proxy variable for identification
    %======================================================================
    if isfield(options,'proxy') ==1
        proxy_  = 1;
        proxy   = options.proxy;
        ns      = size(proxy,2);
        if T~=size(proxy,1)
            error('Shocks proxies and observables must have the same time length')
        end
    end
    %======================================================================
    % options: Compuing Robust Standard Errors
    %======================================================================
    if isfield(options,'robust_se_') == 1
        robust_se_ = options.robust_se_;
        % robust_se_ = 0    unadjusted SE
        % robust_se_ = 1    NW Robust SE:  Hamilton (1994), Ch 10 pag 282, eq (10.5.20)
        % robust_se_ = 5    Matlab HAC function: need Matlab Econ Toolbox
    end
    %======================================================================
    % Conjugate/Hierachical prior options
    %======================================================================
    if (isfield(options,'priors')==1 && strcmp(options.priors.name,'Conjugate')==1) || (isfield(options,'priors')==1 && strcmp(options.priors.name,'conjugate')==1) || ...
            (isfield(options,'prior')==1 && strcmp(options.prior.name,'Conjugate')==1) || (isfield(options,'prior')==1 && strcmp(options.prior.name,'conjugate')==1)

        dummy = 2;
        prior.name= 'Conjugate';
        % Priors for the AR parameters
        if isfield(options.priors,'Phi') == 1
            % mean
            if isfield(options.priors.Phi,'mean') == 1
                prior.Phi.mean  = options.priors.Phi.mean;
                if size(prior.Phi.mean) ~= [ny*lags+(1-noconstant)    ny]
                    error('Size mismatch')
                end
            else
                warning(['You did not provide a prior mean for the AR coeff. ' ...
                    'Assume zeros everywhere.'])
                prior.Phi.mean  = zeros(ny*lags+(1-noconstant) , ny);
            end
            % variance
            if isfield(options.priors.Phi,'cov') == 1
                prior.Phi.cov   = options.priors.Phi.cov;
                if length(prior.Phi.cov) ~= (ny*lags+(1-noconstant) ) || size(prior.Phi.cov,1 )~=size(prior.Phi.cov,2)
                    error('Size mismatch: Covariance Phi should be square, e.g. size(Phi.mean,1)x size(Phi.mean,1)x')
                end
            else
                warning(['You did not provide a Covariance for the AR coeff. ' ...
                    'Assume 10 times Identity Matrix.'])
                prior.Phi.cov  = 10 * eye((ny*lags+(1-noconstant) ));
            end
        else
            warning(['You did not provide prior mean and covariance for the AR coeff ' ...
                'Assume zeros everywhere with covariance 10 times Identity Matrix.'])
            %             prior.Phi.cov   = 10 * eye((ny*lags+(1-noconstant) +timetrend ) * ny);
            prior.Phi.cov   = 10 * eye((ny*lags+(1-noconstant)  ));
            prior.Phi.mean  = zeros(ny*lags+(1-noconstant) ,  ny);
        end
        % Priors for the Residual Covariance
        if isfield(options.priors,'Sigma') == 1
            % scale
            if isfield(options.priors.Sigma,'scale') == 1
                prior.Sigma.scale = options.priors.Sigma.scale;
                if size(prior.Sigma.scale) ~= [ny ny]
                    error('Size mismatch')
                end
            else
                warning(['You did not provide a prior mean for the Residual Covariance. ' ...
                    'Assume identity matrix.'])
                prior.Sigma.scale = eye(ny);
            end
            % degrees of freedom
            if isfield(options.priors.Sigma,'df') == 1
                prior.Sigma.df = options.priors.Sigma.df;
                if length(prior.Sigma.df) ~= 1
                    error('Size mismatch')
                end
                if prior.Sigma.df < ny
                    error('Too few degrees of freedom - prior on variance residuals')
                end
            else
                warning(['You did not provide the degrees of freedom for the Residual Covariance. ' ...
                    'Assume ny+1.'])
                prior.Sigma.df = ny + 1;
            end
        else
            warning(['You did not provide prior mean and variance for the Residual Covariance. ' ...
                'Assume an identity matrix matrix with N+1 degrees of freedom.'])
            prior.Sigma.scale = eye(ny);
            prior.Sigma.df    = ny + 1;
        end
        if isfield(options.priors,'tau') == 1 && isnumeric(options.priors.tau) == 1
            prior.tau = options.priors.tau;
            if length(options.priors.tau) ~= hor
                error('tau must be a vector of size hor')
            end
        elseif isfield(options.priors,'max_tau') == 1 && options.priors.max_tau ==1
            max_prior_tau_ = 1;
            max_compute    = 3;
            prior.tau = ones(hor,1);
        else
            prior.tau = ones(hor,1);
        end
    end
    % options for the maximization
    if isfield(options,'max_compute') == 1
        max_compute    = options.max_compute;
        max_prior_tau_ = 1;
    end
    if isfield(options,'ub') == 1
        ub    = options.ub;
        if length(ub) ~=  length(prior.tau)
            error('Mismatch between the size of upper bounds and the param vector');
        end
    end
    if isfield(options,'lb') == 1
        lb    = options.lb;
        if length(lb) ~=  length(prior.tau)
            error('Mismatch between the size of lower bounds and the param vector');
        end
    end
end

% Confidence Interval Points and Index for confidence serts
alpha  = 1 - conf_sig;
talpha = abs(tinv(alpha/2,T-ny-1));
sort_idx   = round((0.5 + [-conf_sig, conf_sig, 0]/2) * K);

%**************************************************
% Construct the RHS matrix
%**************************************************

X1 = lagX(y,1:lags);
positions_nylags = 1  : size(X1,2);
nylags           = length(positions_nylags);

% forecast launching point
fdata_initval            = lagX(y(end-lags+1:end, :),0:lags-1);
fdata_initval(1:end-1,:) = [];
% add controls if any
if controls_ == 1
    a  = size(X1,2) + 1;
    X1 = [X1 lagX(controls,0:clags)];
    b  = size(X1,2);
    position_controls = a:b;
    % forecast launching point
    %lastvalc = lagX(controls(end-lags+1:end, :),0:lags-1);
    fdata_initval  = [fdata_initval, controls(end,:)];

else
    position_controls = [];
end
% add proxy shocks if any
if proxy_ == 1
    a     = size(X1,2) + 1;
    X1    = [X1 proxy];
    b     = size(X1,2);
    position_proxy = a:b ;
    % forecast launching point (assume shocks are zero)
    fdata_initval  = [fdata_initval,  zeros(1,ns)];
else
    position_proxy = [];
end
X_ = X1;
if noconstant == 0
    X_ = [X1 ones(T,1)];
    position_constant = size(X_,2);
    % forecast launching point
    fdata_initval  = [fdata_initval, 1];
else
    position_constant = [];
end
nx           = 1 - noconstant;
% other checks
if (lags + hor) >= size(X_,1)
    error('More parameters than observations: consider reducing ''hor'' or ''lags''')
end

%**************************************************
% pre allocation
%**************************************************
Omegaproxy  = eye(ny);
Omega       = eye(ny);
ir_lp       = nan(ny,hor+1,ny,3);  % variable, horizon, shock and mean upper lower
irproxy_lp  = nan(ny,hor+1,1,3);
forecasts   = nan(hor,ny,3);
if dummy == 2
    ir_blp                  = nan(ny,hor+1,ny,K);
    irproxy_blp             = nan(ny,hor+1,1,K);
    bforecasts_no_shocks    = nan(hor,ny,K);
    bforecasts_with_shocks  = nan(hor,ny,K);
    log_dnsty               = nan(hor,1);
end

%**************************************************
%* Computing the LP and DF
%**************************************************
wb = waitbar(0, 'Direct Methods');
for hh = 0 : hor % iteration over horizon
    ytmp = lagX(y, -hh);
    % Reduced Form estimations
    if robust_se_ ~= 0 % Robust SE
        options.robust_se_ = robust_se_;
        options.L     = lags + hh + 1;
        olsreg_(hh+1) = ols_reg(ytmp, X_, options);
    else
        olsreg_(hh+1) = ols_reg(ytmp, X_);
    end
    % Proxy IV identifications
    if proxy_
        irproxy_lp(:, hh+1, :, 2) = olsreg_(hh+1).beta(position_proxy, :)'; % mean
        irproxy_lp(:, hh+1, :, 3) = irproxy_lp(:, hh+1, :, 2) + talpha * olsreg_(hh+1).se(position_proxy, :)'; % upper
        irproxy_lp(:, hh+1, :, 1) = irproxy_lp(:, hh+1, :, 2) - talpha * olsreg_(hh+1).se(position_proxy, :)'; % lower
        if hh == 0
            Omegaproxy(:,1) = olsreg_(hh+1).beta(position_proxy, :)';
            Omegasproxy(:,1,:) =  repmat(Omegaproxy(:,1), 1, 1, K) + ...
                olsreg_(hh+1).se(position_proxy, :)'.* randn(ny,1,K);
        end
    end
    % Choleski/Rotated identification
    if hh == 0
        try
            Omega = chol(olsreg_(hh+1).Serror,'Lower') * Q;
        catch
            warning('Covariance of the error term is almost singular')
            Omega = chol(cov(olsreg_(hh+1).error),'Lower') * Q;
        end
        % generate uncertainty on S (flat prior)
        Sbar = olsreg_(hh+1).error'*olsreg_(hh+1).error;
        df   = olsreg_(hh+1).N - olsreg_(hh+1).K;
        if df > 0
        [~, Omegas] = generateOmegas(Sbar,df,K,Q);
            Omegasort = sort(Omegas,3);
        else
            warning('not enought degrees of freedom - too many regressors')
            Omegasort = repmat(Omega,1,1,K);
        end
        ir_lp(:, hh+1, :, 2) = Omega; % mean
        ir_lp(:, hh+1, :, 3) = Omegasort(:,:,sort_idx(2)); % UPPER
        ir_lp(:, hh+1, :, 1) = Omegasort(:,:,sort_idx(1)); % LOWER
    else
        [ir] =  iresponse(olsreg_(hh).beta(positions_nylags, :), eye(ny) , 2, Omega);
        ir_lp(:, hh+1, :, 2) = ir(:, 2, :);  % mean
        [ir3] =  iresponse(olsreg_(hh).beta(positions_nylags, :) + talpha * olsreg_(hh).se(positions_nylags, :), eye(ny) , 2, Omega);
        ir_lp(:, hh+1, :, 3) = ir3(:, 2, :); % upper
        [ir1] =  iresponse(olsreg_(hh).beta(positions_nylags, :) - talpha * olsreg_(hh).se(positions_nylags, :), eye(ny) , 2, Omega);
        ir_lp(:, hh+1, :, 1) = ir1(:, 2, :); % lower
    end

    % forecast part
    forecasts(hh+1, :, 2) = (fdata_initval * olsreg_(hh+1).beta);
    forecasts(hh+1, :, 3) = forecasts(hh+1, :, 2) + talpha * diag(olsreg_(hh+1).Serror)' ;
    forecasts(hh+1, :, 1) = forecasts(hh+1, :, 2) - talpha * diag(olsreg_(hh+1).Serror)' ;

    %======================================================================
    % Baysian DM
    if dummy == 2 % activating Bayesian Direct methods.
        if hh == 0
            % cholesky IRF on impact
            ir_blp(:, hh+1, :, :)  = Omegas;
            % proxy IRF on impact
            if proxy_
                irproxy_blp(:, hh+1, :, :) = Omegasproxy;
            end
            % one-step ahead forecast
            bforecasts_no_shocks(hh+1, :, :)   = repmat(forecasts(hh+1, :, 2),1,1,K);
            fnoise = Omega * randn(ny,K);
            bforecasts_with_shocks(hh+1, :, :) = bforecasts_no_shocks(hh+1, :, :) + reshape(fnoise,1,ny,K);

            % computing the companion matrix
            F       = [prior.Phi.mean(1 : ny * lags, :)'; eye(ny*(lags-1), ny*lags)];
            % constant
            Fo        = [prior.Phi.mean(end, :)'; zeros(ny * (lags-1), 1)];
            % Shocks Companion
            G       = eye(ny * lags, ny);

        else % hh > 0
            %********************************************************
            % Conjugate Prior: MN-IW
            %********************************************************
            if max_prior_tau_ == 1 % Maximize the shrinkage on the VAR coefficients
                try
                    disp(['***********************************************'])
                    disp(['***********************************************'])
                    disp(['Optimization at horizon ' num2str(hh)])
                    x0   = log(prior.tau(hh));
                    switch max_compute
                        case 1
                            optim_options = optimset('display','iter','MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6);
                            [xh,fh,~,~,~,~] = ...
                                fminunc('blp_opt_hyperpara',x0,optim_options,...
                                hh,prior,olsreg_(hh),F,G,Fo,positions_nylags,position_constant);
                            %=====================================================================
                        case 2 % constraint
                            % Set default optimization options for fmincon.
                            optim_options = optimset('display','iter', 'LargeScale','off', 'MaxFunEvals',100000, 'TolFun',1e-8, 'TolX',1e-6);
                            [xh,fh,~,~,~,~,~] = ...
                                fmincon('blp_opt_hyperpara',x0,[],[],[],[],lb(hh),ub(hh),[],optim_options,y,lags,options);
                            %=====================================================================
                        case 3 % Sims
                            crit = 10e-5;
                            nit  = 10e-4;
                            [fh, xh, ~, ~, ~, ~, ~] = ...
                                csminwel('blp_opt_hyperpara',x0,.1*eye(length(x0)),[],crit,nit,hh,prior,olsreg_(hh),F,G,Fo,positions_nylags,position_constant);
                            %=====================================================================
                        case 7 % Matlab's simplex (Optimization toolbox needed).
                            optim_options = optimset('display','iter','MaxFunEvals',30000,'MaxIter',10000,'TolFun',1e-3,'TolX',1e-3);
                            [xh,fh,~,~] = fminsearch('blp_opt_hyperpara',x0,optim_options,y,lags,options);
                    end
                    fprintf('%s = %0.5g\n','Hyper-parameter Mode ',exp(xh))
                    fprintf('%s = %0.5g\n','Marginal Likelihood ',-fh)
                catch
                    warning('Maximization NOT Successful')
                    disp('Using hyper parameter default values')
                    xh = x0;
                end
                prior.tau(hh) = exp(xh);
                disp(['***********************************************'])
            end
            % constructing the posterior moments given the (optimal) shrinkage
            [posterior_, ~]   = p2p(hh,prior.tau(hh),prior,olsreg_(hh),F,G,Fo,positions_nylags,position_constant);
            % computing the marginal likelihood
            log_dnsty(hh)     = blp_ml(prior.tau(hh),hh,prior,olsreg_(hh),F,G,Fo,positions_nylags,position_constant);
            % Second moments
            S_inv_upper_chol  = chol(inv(posterior_.S));
            XXi_lower_chol    = chol(posterior_.XXi)';
            % number of regressors
            nk = nylags + nx;
            %**************************************************
            %* Generating draws form the Posterior Distribution
            %**************************************************
            for  d =  1 : K %  Gibbs Sampler
                %======================================================================
                % Inferece: Drawing from the posterior distribution
                % Step 1: draw from the Covariance
                Sigma = rand_inverse_wishart(ny, posterior_.df, S_inv_upper_chol);

                % Step 2: given the Covariance Matrix, draw from the AR parameters
                Sigma_lower_chol = chol(Sigma)';
                Phi1 = randn(nk * ny, 1);
                Phi2 = kron(Sigma_lower_chol , XXi_lower_chol) * Phi1;
                Phi3 = reshape(Phi2, nk, ny);
                Phi  = Phi3 + posterior_.PhiHat;

                % Step 3: compute IRF
                blp  =  iresponse(Phi, eye(ny) , 2, Omega);
                ir_blp(:, hh+1, :, d) = blp(:, 2, :);  % mean
                if proxy_
                    blpproxy  =  iresponse(Phi, eye(ny) , 2, Omegaproxy);
                    irproxy_blp(:, hh+1, :, d) = blpproxy(:,2,1);
                end
                % Step 4: compute Forecasts
                bforecasts_no_shocks(hh+1, :, d) = (fdata_initval(1, [positions_nylags position_constant]) * Phi);
                bforecasts_with_shocks(hh+1, :, d) = (fdata_initval(1, [positions_nylags position_constant]) * Phi) + (Sigma_lower_chol * randn(ny,1))';
            end
        end
    end
    waitbar(hh/hor, wb);
end
close(wb);

%********************************************************
%* Storing the resutls
%*******************************************************
% dm.olsreg_     = olsreg_;
dm.forecasts   = forecasts;
dm.ir_lp       = ir_lp;
dm.irproxy_lp  = irproxy_lp;
if dummy == 2
    dm.ir_blp                 = ir_blp;
    dm.irproxy_blp            = irproxy_blp;
    dm.bforecasts.no_shocks   = bforecasts_no_shocks;         % trajectories of forecasts without shocks
    dm.bforecasts.with_shocks = bforecasts_with_shocks;       % trajectories of forecasts with shocks
    dm.logmlike               = log_dnsty;
    dm.prior                  = prior;
else
    dm.ir_blp                 = [];
    dm.irproxy_blp            = [];
    dm.bforecasts.no_shocks   = [];
    dm.bforecasts.with_shocks = [];
    dm.logmlike               = [];
    dm.prior                  = [];
end

%**************************************************************************
%**************************************************************************
%**************************************************************************
%**************************************************************************
function [S,Omegas] = generateOmegas(Sbar,df,K,Q)
% generate uncertainty on S
ny      = size(Sbar,1);
S       = nan(ny,ny,K);
Omegas  = nan(ny,ny,K);

S_inv_upper_chol = chol(inv(Sbar));

for  d = 1:K
    S(:,:,d)      = rand_inverse_wishart(ny, df, S_inv_upper_chol);
    Omegas(:,:,d) = chol(S(:,:,d),'Lower') * Q;
end


text:
function H = cross(F, G)
%CROSS   Cross product of CHEBFUN3V objects.
%   CROSS(F, G) returns a CHEBFUN3V representing the 3D cross product of
%   the CHEBFUN3V objects F and G.

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

% Empty check:
if ( isempty(F) || isempty(G) )
    H = chebfun3v;
    return
end

% Get number of components:
Fc = F.components;
Gc = G.components;

if ( F.nComponents == 3 && G.nComponents == 3 )
    H = [ Fc{2} .* Gc{3} - Fc{3} .* Gc{2};
          Fc{3} .* Gc{1} - Fc{1} .* Gc{3};
          Fc{1} .* Gc{2} - Fc{2} .* Gc{1}];
else
    error('CHEBFUN:CHEBFUN3V:cross:components', ...
        'CHEBFUN3V objects must be both 3-vectors.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% quantile() - computes the quantiles of the data sample from a distribution X
%
% Description:
%      If F is the cumulative distribution function (CDF) of X,
%      the p-th-quantile Qp of distribution X is the value for which holds
%                        F(x) <  p,  for x < Qp, and
%                        F(x) >= p,  for x >= Qp.
%      for example, for p = 0.5, Qp is the median of X. p must be in [0..1].
%
% Usage:
%   >>  q = quantile( data, pc );
%
% Inputs:
%   data   - vector of observations
%   pc     - the quantiles will be estimated at the values in pc [0..1]
%
% Outputs:
%   q      - pc-th-quantiles of the distribution generating the observation
%
% Authors: Scott Makeig & Luca Finelli, CNL/Salk Institute-SCCN, August 21, 2002
%
% Note: this function overload the function from the statistics toolbox. In
%       case the statistic toolbox is present, the function from the
%       statistics toolbox is being called instead.
%
% See also:
%   pop_sample(), eeglab()

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

function q = quantile(data,varargin);

% detect overloaded method in stat toolbox
curPath = fileparts(which('quantile'));
rmpath(curPath);
path2   = fileparts(which('quantile'));
addpath(curPath);
if ~isempty(path2)
    addpath(path2);
    q = quantile(data,varargin{:});
    return;
else
    pc = varargin{1};
end;

if nargin < 2
	help quantile;
	return;
end;

[prows pcols] = size(pc);
if prows ~= 1 & pcols ~= 1
    error('pc must be a scalar or a vector.');
end
if any(pc > 1) | any(pc < 0)
    error('pc must be between 0 and 1');
end
[i,j] = size(data);
sortdata = sort(data);
if i==1 | j==1 % if data is a vector
    i = max(i,j); j = 1;
    if i == 1,
        fprintf('  quantile() note: input data is a single scalar!\n')
        q = data*ones(length(pc),1); % if data is scalar, return it
        return;
    end
    sortdata = sortdata(:);
end
pt = [0 ((1:i)-0.5)./i 1];
sortdata = [min(data); sortdata; max(data)];
q = interp1(pt,sortdata,pc);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function out = scimat_resample(inp,out,kind)
% SCIMAT_RESAMPLE resamples a scimat file given an input and output scimat.
%
% Intended for mapping 2D slices in 3D volume or extract 2D
% slices from 3D volume. Also works nicely between 1D and 2D to extract
% profiles, for example.
%
% OUT = SCIMAT_RESAMPLE(INP, OUT, KIND)
%
%   INP is the input scimat file
%
%   OUT is the target scimat file
%
%   KIND is a string, deciding on the kind of interpolation/averaging
%
% See also: scimat_resize3, scimat_downsample.

% Author: Nicolas Basty <nicolas.basty@eng.ox.ac.uk>
% Copyright © 2016 University of Oxford
% Version: 0.1.2
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
% _________________________________________________________________

[x,y,z] = inp.axis.spacing; %get in and output resolutions
[x_o,y_o,z_o] = out.axis.spacing;

if x_o>x %taking lower resolution
    x = x_o;
    y = y_o;
    z = z_o;
end

sigmaa = [x,y,inf]/2;
ws =  sigmaa*3; %window size
ws(3) = z/2;
rrrr=r;
if strcmp(kind,'gaussian') % get the weighting function
    W = @(xl,yl,zl) exp(-(xl.^2)/(2*sigmaa(1)^2)...
        -(yl.^2)/(2*sigmaa(2)^2)...
        -(zl.^2)/(2*sigmaa(3)^2));
elseif strcmp(kind,'sinc')
    W = @(xl,yl,zl) (sin(sigmaa(1)*xl)./(sigmaa(1)*xl))...
        .* (sin(sigmaa(2)*yl)./(sigmaa(2)*yl))...
        .* (sin(sigmaa(3)*zl)./(sigmaa(3)*zl));
elseif strcmp(kind,'nearest')
    %does something later, l84
else
    fprintf('Weighting method expected to be gaussian, sinc or nearest\n')
    return
end

To = scimat2transform(inp); %gets transform from input image mapping IDX to RW
Ti = scimat2transform(out);  %gets transform from target image mapping IDX to RW
T = pinv(To) * Ti; %combining gives IDXin to IDXout

% create grid of output�s indices
if numel(out.data)>2
    [X, Y, Z] = ndgrid(1:size(out.data,1), 1:size(out.data,2), 1:size(out.data,3));
elseif  numel(out.data)<3
    [X, Y, Z] = ndgrid(1:size(out.data,1), 1:size(out.data,2), 1);
end
% for each of the output index, find its corresponding input index
XYZin = T * [X(:), Y(:), Z(:), ones(size(X(:)))]';

% round them to do nearest neightbour assignment
Xin = round(XYZin(1,:));
Yin = round(XYZin(2,:));
Zin = round(XYZin(3,:));

% find which ones are valid (input space)
valid_indices = (Xin > 0) & (Xin <= size(inp.data,1)) & (Yin > 0) & (Yin <= size(inp.data,2)) & (Zin > 0) & (Zin <= size(inp.data,3));

%%
inpt = double(inp.data); %get data out of the structures
outp = double(out.data);

allind = 1:numel(valid_indices);
allind = allind(valid_indices);

[xi,yi,zi] = scimat_ndgrid(inp); % grid of input rw
[xo,yo,zo] = scimat_ndgrid(out); % grid of output rw

vox_out = [xo(:),yo(:),zo(:)];

[Xi, Yi, Zi] = ndgrid(1:size(inp.data,1), 1:size(inp.data,2), 1:size(inp.data,3));
xyzin = round(pinv(T)*[Xi(:) Yi(:) Zi(:) ones(size(Xi(:)))]');
xin = xyzin(1,:);
yin = xyzin(2,:);
zin = xyzin(3,:);
%using indices to get valid voxels, but the other part works better in
%rw... window size gets messed up when converting to index
valid_indices_out = (xin > 0) & (xin <= size(out.data,1)) & (yin > 0) & (yin <= size(out.data,2)) & (zin > 0) & (zin <= size(out.data,3));

inpt_val = inpt(valid_indices_out);
%keep voxels relevant to the averaging
xi = xi(valid_indices_out);
yi = yi(valid_indices_out);
zi = zi(valid_indices_out);

if strcmp(kind,'nearest')
    source_indices = sub2ind(size(inp.data), Xin(valid_indices), Yin(valid_indices), Zin(valid_indices));
    valid_intensities_input = inp.data(source_indices);
    out.data(valid_indices) = valid_intensities_input;
    return
end
outpdummy = zeros(numel(allind),1);
% xyz = [xi(:) yi(:) zi(:)];
% par parfor is actually slower for the data I am using so I commented it
% out
for r = 1:numel(allind)
    p =  vox_out(allind(r),:)';
    % get a neighbourhood of vox_out around this point (rw)
    mask = (abs(xi-p(1)) <= ws(1)) & (abs(yi - p(2)) <= ws(2)) & (abs(zi - p(3)) <= ws(3));
    %     mask = bsxfun(@le, (abs(bsxfun(@minus, xyz, p'))), ws);
    %     mask = sum(mask,2)==3;
    % get a shortlist of voxels in the window & distances for weighting
    xlist = xi(mask) - p(1);
    ylist = yi(mask) - p(2);
    zlist = zi(mask) - p(3);
    Ilist = inpt_val(mask);

    % get a weighted average
    outpdummy(r) = sum(W(xlist,ylist,zlist).*Ilist) / (sum(W(xlist,ylist,zlist)+eps));
end
outp(allind) = outpdummy;
outp(isnan(outp)) = 0;
out.data = outp;
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [J,residuals] = deconvlucydenoiseNew(varargin) %CHANGED - added meanvec
%DECONVLUCY Deblur image using Lucy-Richardson method.
%   J = DECONVLUCY(I,PSF) deconvolves image I using Lucy-
%   Richardson algorithm, returning deblurred image J. The assumption is
%   that the image I was created by convolving a true image with a
%   point-spread function PSF and possibly by adding noise.
%
%   I can be an N-Dimensional array.
%
%   To improve the restoration, additional parameters can be passed in
%   (use [] as a place holder if an intermediate parameter is unknown):
%   J = DECONVLUCY(I,PSF,NUMIT)
%   J = DECONVLUCY(I,PSF,NUMIT,DAMPAR)
%   J = DECONVLUCY(I,PSF,NUMIT,DAMPAR,WEIGHT)
%   J = DECONVLUCY(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT)
%   J = DECONVLUCY(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT,SUBSMPL), where
%
%   NUMIT   (optional) is the number of iterations (default is 10).
%
%   DAMPAR  (optional) is an array that specifies the threshold deviation
%   of the resulting image from the image I (in terms of the standard
%   deviation of Poisson noise) below which the damping occurs. The
%   iterations are suppressed for the pixels that deviate within the
%   DAMPAR value from their original value. This suppresses the noise
%   generation in such pixels, preserving necessary image details
%   elsewhere. Default is 0 (no damping).
%
%   WEIGHT  (optional) is assigned to each pixel to reflect its recording
%   quality in the camera. A bad pixel is excluded from the solution by
%   assigning it zero weight value. Instead of giving a weight of one for
%   good pixels, you can adjust their weight according to the amount of
%   flat-field correction. Default is a unit array of the same size as
%   input image I.
%
%   READOUT (optional) is an array (or a value) corresponding to the
%   additive noise (e.g., background, foreground noise) and the variance
%   of the read-out camera noise. READOUT has to be in the units of the
%   image. Default is 0.
%
%   SUBSMPL (optional) denotes subsampling and is used when the PSF is
%   given on a grid that is SUBSMPL times finer than the image. Default
%   is 1.
%
%   Note that the output image J could exhibit ringing introduced by the
%   discrete Fourier transform used in the algorithm. To reduce the
%   ringing use I = EDGETAPER(I,PSF) prior to calling DECONVLUCY.
%
%   Note also that DECONVLUCY allows you to resume deconvolution starting
%   from the results of an earlier DECONVLUCY run. To initiate this
%   syntax, the input image I has to be passed in as cell array, {I}.
%   Then the output J becomes a cell array and can be passed as the input
%   array into the next DECONVLUCY call. The input cell array can contain
%   one numeric array (on initial call), or four numeric arrays (when it
%   is the output from a previous run of DECONVLUCY). The output J
%   contains four elements, where J{1}=I, J{2} is the image resulted from
%   the last iteration, J{3} is the image from one before last iteration,
%   J{4} is an array used internally by the iterative algorithm.
%
%   Class Support
%   -------------
%   I and PSF can be uint8, uint16, int16, double, or single. DAMPAR and
%   READOUT must have the same class as the input image. Other inputs have to
%   be double. The output image (or the first array of the output cell) has
%   the same class as the input image.
%
%   Example
%   -------
%
%      I = checkerboard(8);
%      PSF = fspecial('gaussian',7,10);
%      V = .0001;
%      BlurredNoisy = imnoise(imfilter(I,PSF),'gaussian',0,V);
%      WT = zeros(size(I));WT(5:end-4,5:end-4) = 1;
%      J1 = deconvlucy(BlurredNoisy,PSF);
%      J2 = deconvlucy(BlurredNoisy,PSF,20,sqrt(V));
%      J3 = deconvlucy(BlurredNoisy,PSF,20,sqrt(V),WT);
%      subplot(221);imshow(BlurredNoisy);
%                     title('A = Blurred and Noisy');
%      subplot(222);imshow(J1);
%                     title('deconvlucy(A,PSF)');
%      subplot(223);imshow(J2);
%                     title('deconvlucy(A,PSF,NI,DP)');
%      subplot(224);imshow(J3);
%                     title('deconvlucy(A,PSF,NI,DP,WT)');
%
%   See also DECONVWNR, DECONVREG, DECONVBLIND, EDGETAPER, IMNOISE, PADARRAY,
%            PSF2OTF, OTF2PSF.

%   Copyright 1993-2011 The MathWorks, Inc.
%   $Revision: 1.6.4.10 $
%

%   References
%   ----------
%   "Acceleration of iterative image restoration algorithms, by D.S.C. Biggs
%   and M. Andrews, Applied Optics, Vol. 36, No. 8, 1997.
%   "Deconvolutions of Hubble Space Telescope Images and Spectra",
%   R.J. Hanisch, R.L. White, and R.L. Gilliland. in "Deconvolution of Images
%   and Spectra", Ed. P.A. Jansson, 2nd ed., Academic Press, CA, 1997.

% Parse inputs to verify valid function calling syntaxes and arguments
[J,PSF,NUMIT,DAMPAR,READOUT,WEIGHT,SUBSMPL,sizeI,classI,numNSdim]=...
  parse_inputs(varargin{:});

% 1. Prepare PSF. If PSF is known at a higher sampling rate, it has to be
% padded with zeros up to sizeI(numNSdim)*SUBSMPL in all non-singleton
% dimensions. Or its OTF could take care of it:
sizeOTF = sizeI;
sizeOTF(numNSdim) = SUBSMPL*sizeI(numNSdim);
H = psf2otf(PSF,sizeOTF);

% 2. Prepare parameters for iterations
%
% Create indexes for image according to the sampling rate
idx = repmat({':'},[1 length(sizeI)]);
for k = numNSdim,% index replicates for non-singleton PSF sizes only
  idx{k} = reshape(repmat(1:sizeI(k),[SUBSMPL 1]),[SUBSMPL*sizeI(k) 1]);
end

wI = max(WEIGHT.*(READOUT + J{1}),0);% at this point  - positivity constraint
J{2} = J{2}(idx{:});
scale = real(ifftn(conj(H).*fftn(WEIGHT(idx{:})))) + sqrt(eps);
clear WEIGHT;
DAMPAR22 = (DAMPAR.^2)/2;

if SUBSMPL~=1,% prepare vector of dimensions to facilitate the reshaping
  % when the matrix is binned within the iterations.
  vec(2:2:2*length(sizeI)) = sizeI;
  vec(2*numNSdim-1) = -1;
  vec(vec==0) = [];
  num = fliplr(find(vec==-1));
  vec(num) = SUBSMPL;
else
  vec = [];
  num = [];
end

% 3. L_R Iterations
%
residuals=zeros(1,NUMIT); %CHANGED - initializing residualvector

lambda = 2*any(J{4}(:)~=0);
% h = waitbar(0,'Applying PVE correction to input volume ... please wait'); % waitbar
tic
for k = lambda + (1:NUMIT)

  % 3.a Make an image predictions for the next iteration
  if k > 2,
    lambda = (J{4}(:,1).'*J{4}(:,2))/(J{4}(:,2).'*J{4}(:,2) +eps);
    lambda = max(min(lambda,1),0);% stability enforcement
  end
  Y = max(J{2} + lambda*(J{2} - J{3}),0);% plus positivity constraint

  % 3.b  Make core for the LR estimation
  fprintf('Iteration %u ... ',k)
  [CC,residuals(k)] = corelucydenoiseNew(Y,H,DAMPAR22,wI,READOUT,SUBSMPL,idx,vec,num); %CHANGED to "denoise" and to call meanvec

  % 3.c Determine next iteration image & apply positivity constraint
  J{3} = J{2};
  J{2} = max(Y.*real(ifftn(conj(H).*CC))./scale,0);
  clear CC;
  J{4} = [J{2}(:)-Y(:) J{4}(:,1)];
  % waitbar(k/NUMIT,h,'Applying PVE correction to input volume ... please wait') % waitbar
end
% waitbar(1,h,'Applying PVE correction to input volume ... DONE'), pause(1) % waitbar
toc
% close(h)
clear wI H scale Y;

% 4. Convert the right array (for cell it is first array, for notcell it is
% second array) to the original image class & output whole thing
num = 1 + strcmp(classI{1},'notcell');
if ~strcmp(classI{2},'double'),
  J{num} = changeclass(classI{2},J{num});
end

if num==2,% the input & output is NOT a cell
  J = J{2};
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Function: parse_inputs
function [J,PSF,NUMIT,DAMPAR,READOUT,WEIGHT,SUBSMPL,sizeI,classI,numNSdim] = ...
    parse_inputs(varargin)
%
% Outputs:
% I=J{1}   the input array (could be any numeric class, 2D, 3D)
% PSF      operator that distorts the ideal image
% numNSdim non-singleton dimensions of PSF
%
% Defaults:
%
NUMIT = [];NUMIT_d = 10;% Number of  iterations, usually produces good
                        % result by 10.
DAMPAR =[];DAMPAR_d = 0;% No damping is default
WEIGHT =[];             % All pixels are of equal quality, flat-field is one
READOUT=[];READOUT_d= 0;% Zero readout noise or any other
           % back/fore/ground noise associated with CCD camera.
           % Or the Image is corrected already for this noise by user.
SUBSMPL= [];SUBSMPL_d= 1;% Image and PSF are given at equal resolution,
           % no over/under sampling at all.

narginchk(2,7);

% First, assign the inputs starting with the image
%
if iscell(varargin{1}),% input cell is used to resume interrupted iterations
  classI{1} = 'cell';% or interrupt the iteration to resume it later
  J = varargin{1};
else % no-cell array is used to do a single set of iterations
  classI{1} = 'notcell';
  J{1} = varargin{1};% create a cell array in order to do the iterations
end;

% check the Image, which is the first array of the cell
classI{2} = class(J{1});

validateattributes(J{1},{'uint8' 'uint16' 'double' 'int16','single'},...
              {'real' 'nonempty' 'finite'},mfilename,'I',1);

if length(J{1})<2,
    error(message('images:deconvlucy:inputImagesMustHaveAtLeast2Elements'))
elseif ~isa(J{1},'double'),
    J{1} = im2double(J{1});
end

% now since the image is OK&double, we assign the rest of the J cell
len = length(J);
if len == 1,% J = {I} will be reassigned to J = {I,I,0,0}
  J{2} = J{1};
  J{3} = 0;
elseif len ~= 4,% J = {I,J,Jm1,gk} has to have 4 or 1 arrays
    error(message('images:deconvlucy:inputCellMustHave1or4Elements'));
else % check if J,Jm1,gk are double in the input cell
  if ~all([isa(J{2},'double'),isa(J{3},'double'),isa(J{4},'double')]),
    error(message('images:deconvlucy:inputImageCellElementsMustBeDouble'))
  end
end;

% Second, Assign the rest of the inputs:
%
PSF = varargin{2};%      deconvlucy(I,PSF)
switch nargin
case 3,%                 deconvlucy(I,PSF,NUMIT)
  NUMIT = varargin{3};
case 4,%                 deconvlucy(I,PSF,NUMIT,DAMPAR) CHANGED
  NUMIT = varargin{3};
  DAMPAR = varargin{4};
case 5,%                 deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT)
  NUMIT = varargin{3};
  DAMPAR = varargin{4};
  WEIGHT = varargin{5};
case 6,%                 deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT)
  NUMIT = varargin{3};
  DAMPAR = varargin{4};
  WEIGHT = varargin{5};
  READOUT = varargin{6};
case 7,%                 deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT,SUBSMPL)
  NUMIT = varargin{3};
  DAMPAR = varargin{4};
  WEIGHT = varargin{5};
  READOUT = varargin{6};
  SUBSMPL = varargin{7};
end

% Third, Check validity of the input parameters:
%
% NUMIT check number of iterations
if isempty(NUMIT),
  NUMIT = NUMIT_d;
else
  validateattributes(NUMIT,{'double'},{'scalar' 'positive' 'finite'},...
                mfilename,'NUMIT',3);
end

% SUBSMPL check sub-sampling rate
if isempty(SUBSMPL),
  SUBSMPL = SUBSMPL_d;
else
  validateattributes(SUBSMPL,{'double'},{'scalar' 'positive' 'finite'},...
                mfilename,'SUBSMPL',7);
end

% PSF array
[sizeI, sizePSF] = padlength(size(J{1}), size(PSF));
numNSdim = find(sizePSF~=1);
if prod(sizePSF)<2,
  error(message('images:deconvlucy:psfMustHaveAtLeast2Elements'))
elseif all(PSF(:)==0),
  error(message('images:deconvlucy:psfMustNotBeZeroEverywhere'))
elseif any(sizePSF(numNSdim)/SUBSMPL > sizeI(numNSdim)),
  error(message('images:deconvlucy:psfMustBeSmallerThanImage'))
end
if length(J)==3,% assign the 4-th element of input cell now
  J{4}(prod(sizeI)*SUBSMPL^length(numNSdim),2) = 0;
end;

% DAMPAR check damping parameter
if isempty(DAMPAR),
  DAMPAR = DAMPAR_d;
elseif (numel(DAMPAR)~=1) && ~isequal(size(DAMPAR),sizeI),
  error(message('images:deconvlucy:damparMustBeSameSizeAsImage'))
elseif ~isa(DAMPAR,classI{2}),
  error(message('images:deconvlucy:damparMustBeSameClassAsInputImage'))
elseif ~strcmp(classI{2},'double'),
  DAMPAR = im2double(DAMPAR);
end

validateattributes(DAMPAR,{'double'},{'finite'},mfilename,'DAMPAR',4);

% READOUT check read-out noise
if isempty(READOUT),
  READOUT = READOUT_d;
elseif (numel(READOUT)~=1) && ~isequal(size(READOUT),sizeI),
  error(message('images:deconvlucy:readoutMustBeSameSizeAsImage'))
elseif ~isa(READOUT,classI{2}),
  error(message('images:deconvlucy:readoutMustBeSameClassAsInputImage'))
elseif ~strcmp(classI{2},'double'),
  READOUT = im2double(READOUT);
end

validateattributes(READOUT,{'double'},{'finite'},mfilename,'READOUT',6);

% WEIGHT check weighting
if isempty(WEIGHT),
  WEIGHT = ones(sizeI);
else
    validateattributes(WEIGHT,{'double'},{'finite'},mfilename,'WEIGHT',5);
    if (numel(WEIGHT)~=1) && ~isequal(size(WEIGHT),sizeI),
      error(message('images:deconvlucy:weightMustBeSameSizeAsImage'))
    elseif numel(WEIGHT)== 1,
      WEIGHT = repmat(WEIGHT,sizeI);
    end
end

author: mvallieres

repo: radiomics


text:
% Copyright Mohammad SAFEEA, 17th-Aug-2017
function factor=jointLimitAvoidance(jointIndex,jPos,dq)
% Return the displacemnt according to joint limits
clearance=pi/24;
maxLimit=[170,120,170,120,170,120,175]-clearance*ones(1,7);
minLimit=-maxLimit;
i=jointIndex;
    if(abs(jPos{i}+dq)>maxLimit(i))
        factor=0;
    else
        factor=1;
    end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function y = maskrepelem(x, mask)
%MASKREPELEM Replicate elements according to a mask array.
%
%   MASKREPELEM(X, MASK) replicates each element of X so it becomes a block
%   with the size of MASK and the element inserted at the positions where
%   MASK has a non-zero value.  X must be a matrix, but MASK may be any ND
%   array.  X and MASK may be of any class.
%
%   When X and MASK are double arrays, and MASK is a matrix of zeros and
%   ones, MASKREPELEM(X, MASK) returns the same as KRON(X, MASK).
%
%   For example, maskrepelem(magic(2), [0 1;1 0]) returns
%
%      [ 0  1  0  3
%        1  0  3  0
%        0  4  0  2
%        4  0  2  0 ]

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:20:49 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

% some error checking
error(nargchk(2, 2, nargin));
if ndims(x) > 2
   error('X must be a 2D matrix.');
end

% get sizes, number of elements and number of dimensions
sx = size(x);                       % size of x
nx = prod(sx);                      % number of elements in x
sm = size(mask);                    % size of mask
nm = prod(sm);                      % number of elements in mask
dm = ndims(mask);                   % number of dimensions in mask

% initialize output; replicate and insert the elements of x
i = find(mask);                     % index of non-zero elements
x = reshape(x, 1, nx);              % make x a row vector
y(nm, nx) = feval(class(x), 0);     % let y have same class as x
y(i,:) = x(ones(length(i),1),:);    % replicate and insert elements

% manipulate the blocks
y = reshape(y, [ sm sx ]);
y = permute(y, [ 1 1+dm 2 2+dm 3:dm ]);
y = reshape(y, [ sx(1)*sm(1) sx(2)*sm(2) sm(3:end) ]);

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function obj=coord_flip(obj)
% coord_flip Flips plot axes
%
% Using coord_flip flips the plot axes, i.e. x data will be on the vertical axis and y data on
% the horizontal axis. This is mainly useful to draw stat_ visualisations horizontally,
% such as boxplots, violin plots, etc. Coord flip sets the following axes properties
% 'CameraUpVector',[1 0 0],'XDir','reverse','YAxisLocation','Right' and
% sets things up for correct resizing of the figure. You can further modify
% these axe properties with axe_property() in order to get the desired axes
% orientation and location

    %CameraUpVector takes in accound 'XDir': when using 'reverse' the
    %CameraUpVector doesn't need to be negative X [-1 0 0]
    obj.axe_property('CameraUpVector',[1 0 0],'XDir','reverse','YAxisLocation','Right');
    for obj_ind=1:numel(obj)
        obj(obj_ind).is_flipped=true;
    end
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
function [pg, data, imsegs] = ijcvTestImage(im, imsegs, classifiers, smaps, spdata, adjlist, edata)
%  [pg, data, imsegs] = ijcvTestImage(im, imsegs, classifiers, smaps,
%  spdata, adjlist, edata)
%
% Computes the marginals of the geometry for the given input image
% spdata, adjlist, edata are optional inputs
% Note: only first three arguments are required

nsegments = [10 20 30 40 50 60 80 100];

% if imsegs is a cell, it contains the system command to make the
% segmentation image and the filename of the segmentation image

if iscell(imsegs)
%    syscall = imsegs{1};
%    outfn = imsegs{2};
%    system(syscall);
   segs = segment_fh(im, imsegs{1}, imsegs{2}, imsegs{3});
    imsegs = processSuperpixelImage2(segs);
end

vclassifier = classifiers.vclassifier;
hclassifier = classifiers.hclassifier;
sclassifier = classifiers.sclassifier;

if ~exist('spdata', 'var') || isempty(spdata)
    spdata = mcmcGetSuperpixelData(im, imsegs);
end

if ~exist('adjlist', 'var') || ~exist('edata', 'var') || isempty(adjlist) || isempty(edata)
    [edata, adjlist] = mcmcGetEdgeData(imsegs, spdata);
end

imdata = mcmcComputeImageData(im, imsegs);

if ~exist('smaps', 'var') || isempty(smaps)
    eclassifier = classifiers.eclassifier;
    ecal = classifiers.ecal;
    if isfield(classifiers, 'vclassifierSP')
        vclassifierSP = classifiers.vclassifierSP;
        hclassifierSP = classifiers.hclassifierSP;
        [pvSP, phSP, pE] = mcmcInitialize(spdata, edata, ...
            adjlist, imsegs, vclassifierSP, hclassifierSP, eclassifier, ecal, 'none');
    else
        pE = test_boosted_dt_mc(eclassifier, edata);
        pE = 1 ./ (1+exp(ecal(1)*pE+ecal(2)));
    end
    smaps = generateMultipleSegmentations2(pE, adjlist, imsegs.nseg, nsegments);
end

nsp = imsegs.nseg;
segs = cell(nsp, 1);

if exist('pvSP', 'var')
    for k = 1:numel(segs)
        segs{k}{end+1} = k;
    end
    pg = [pvSP(:, 1) repmat(pvSP(:, 2), 1, 5).*phSP pvSP(:, 3)];
else
    pg = zeros(nsp, 7);
end

for k = 1:size(smaps, 2)

    for s = 1:max(smaps(:, k))

        [segs, ind] = checksegs(segs, smaps(:, k), s);

        if ~isempty(ind)

            labdata = mcmcGetSegmentFeatures(imsegs, spdata, imdata, smaps(:, k), s);

            vconf = test_boosted_dt_mc(vclassifier, labdata);
            vconf = 1 ./ (1+exp(-vconf));
            vconf = vconf / sum(vconf);

            hconf = test_boosted_dt_mc(hclassifier, labdata);
            hconf = 1 ./ (1+exp(-hconf));
            hconf = hconf / sum(hconf);

            sconf = test_boosted_dt_mc(sclassifier, labdata);
            sconf = 1 ./ (1+exp(-sconf));

            pgs = [vconf(1) vconf(2)*hconf vconf(3)]*sconf;

            pg(ind, :) = pg(ind, :) + repmat(pgs, numel(ind), 1);
        end

    end

end

pg = pg ./ max(repmat(sum(pg, 2), 1, size(pg, 2)), 0.00001);

data.smaps = smaps;
data.edata = edata;
data.adjlist = adjlist;
data.spdata = spdata;
data.imdata = imdata;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [segs, ind] = checksegs(segs, map, s)
% Checks whether this segment has been seen before

ind = find(map==s);

if isempty(ind)
    return;
end

% if numel(ind)==1 % already accounted for by superpixels
%     ind = [];
%     return;
% end

oldsegs = segs{ind(1)};

for k = 1:numel(oldsegs)
    if (numel(oldsegs{k})==numel(ind)) && all(oldsegs{k}==ind)
        ind = [];
        return;
    end
end

segs{ind(1)}{end+1} = ind;

author: Cloud-CV

repo: object-proposals


text:
function [Obj, log] = SOFAexpand(Obj,VarName)
%SOFAexpand - Expand variables in a SOFA object
%   Usage: Obj = SOFAexpand(Obj,VarName)
%
%   SOFAexpand(Obj,VarName) expands the the singleton dimensions of
%   the variable VarName, e.g., a variable SourcePosition containing a single entry
%   thus size of [I C], will be expanded to the next possible,thus size of [I M].
%   The potential for the expansion is defined in API.Dimensions.VarName.
%   Obj.API.Dimensions will be updated to the new dimensions.
%
%   SOFAexpand(Obj) expands the singleton dimensions of all variables of Obj.
%   Note that Data variables are not expanded.
%
%	Obj = SOFAexpand(Obj,VarName) expands the singleton dimensions of the variable VarName.
%
%   [Obj,log] = SOFAexpand(..) returns a log of expanded variables.

% #Author: Piotr Majdak
% #Author: Michael Mihocic: header documentation updated (28.10.2021)
%
% SOFA Toolbox - function SOFAexpand
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or - as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

%% Initial check
OC = SOFAgetConventions(Obj.GLOBAL_SOFAConventions,'a');
log={''};

%% Update dimensions
Obj=SOFAupdateDimensions(Obj,'nodata');

%% If VarName given, expand a single variable only
if ~exist('VarName','var'),

  %% Expand all variables
    % create field names which should have dimensions
  X=rmfield(Obj,{'Data','API'});
  if isfield(X,'PRIVATE'), X=rmfield(X,'PRIVATE'); end
  Xf=fieldnames(X);

  % Update the dimensions structure w/o data
  for ii=1:length(Xf)
    if ~isempty(strfind(Xf{ii},'_')),	continue; end; % is an attribute --> not expandable
    if ~isfield(OC.API.Dimensions, Xf{ii}), continue; end; % is a used-defined variable --> not expandable
    dim=OC.API.Dimensions.(Xf{ii}); % all possible dimensions
    if ~iscell(dim), continue; end;	% is a variable with a single dimension definition --> not expandable
    if numel(dim)==1, continue; end;	% is a variable with a single dimension definition --> not expandable
    [varNew,dimNew]=expand(Obj,Xf{ii},dim);
    if ~isempty(varNew),
      Obj.(Xf{ii})=varNew;
      Obj.API.Dimensions.(Xf{ii})=dimNew;
      log{end+1}=[Xf{ii} ' expanded to ' dimNew];
    end
  end

  % Expand the dimensions of Data
  Xf=fieldnames(Obj.Data);
  for ii=1:length(Xf)
    if ~isempty(strfind(Xf{ii},'_')),	continue; end; % is an attribute --> not expandable
    if ~isfield(OC.API.Dimensions.Data, Xf{ii}), continue; end; % is a used-defined variable --> not expandable
    dim=OC.API.Dimensions.Data.(Xf{ii}); % all possible dimensions
    if ~iscell(dim), continue; end;	% is a variable with a single dimension definition --> not expandable
    if numel(dim)==1, continue; end;	% is a variable with a single dimension definition --> not expandable
    [varNew,dimNew]=expandData(Obj,Xf{ii},dim);
    if ~isempty(varNew),
      Obj.Data.(Xf{ii})=varNew;
      Obj.API.Dimensions.Data.(Xf{ii})=dimNew;
      log{end+1}=['Data.' Xf{ii} ' expanded to ' dimNew];
    end
  end

else	% Expand a single variable only
	if isempty(strfind(VarName,'_')),	% is an attribute --> not expandable
    if strncmp(VarName,'Data.',length('Data.'))
      % variable within the Data. structure
      VarName=VarName(length('Data.')+1:end);
      if isfield(OC.API.Dimensions.Data, VarName), % is a used-defined variable --> not expandable
        dim=OC.API.Dimensions.Data.(VarName); % all possible dimensions
        if iscell(dim), % is a variable with a single dimension definition --> not expandable
            [varNew,dimNew]=expandData(Obj,VarName,dim);
            if ~isempty(varNew),
                Obj.Data.(VarName)=varNew;
                Obj.API.Dimensions.Data.(VarName)=dimNew;
                log{end+1}=['Data.' VarName ' expanded to ' dimNew];
            end
        end
      end
    else
      if isfield(OC.API.Dimensions, VarName), % is a used-defined variable --> not expandable
        dim=OC.API.Dimensions.(VarName); % all possible dimensions
        if iscell(dim), % is a variable with a single dimension definition --> not expandable
          [varNew,dimNew]=expand(Obj,VarName,dim);
          if ~isempty(varNew),
              Obj.(VarName)=varNew;
              Obj.API.Dimensions.(VarName)=dimNew;
              log{end+1}=[VarName ' expanded to ' dimNew];
          end
        end
      end
		end
	end
end

%% log variable
if length(log)>1, log=log(2:end); else log={}; end;

%% expand a single variable
% Obj: the full SOFA object
% f: name of the variable
% dims: allowed dimensions of that variable (cell array)
% var: expanded variable, or empty if nothing happened
% dN: new dimension, or empty if nothing happened
function [var,dN]=expand(Obj,f,dims)
	d=cell2mat(strfind(dims,'I'));	% all choices for a singleton dimensions
	for jj=1:length(d)	% loop through all expandable dimensions
		len=size(Obj.(f),d(jj)); % size of the considered dimension
		if len>1, continue; end;	% the expandable dimension is already expanded
		dN=dims{cellfun('isempty',strfind(dims,'I'))==1};
		var=bsxfun(@times,Obj.(f),ones([getdim(Obj,dN) 1]));
	end
	if ~exist('var','var'), var=[]; dN=[]; end;
%% Get the sizes of the dimension variables according the dimension variables in str
function vec=getdim(Obj,str)
	vec=arrayfun(@(f)(Obj.API.(f)),upper(str));

%% expand a single Data variable
% Obj: the full SOFA object
% f: name of the Data variable
% dims: allowed dimensions of that variable (cell array)
% var: expanded variable, or empty if nothing happened
% dN: new dimension, or empty if nothing happened
function [var,dN]=expandData(Obj,f,dims)
	d=cell2mat(strfind(dims,'I'));	% all choices for a singleton dimensions
	for jj=1:length(d)	% loop through all expandable dimensions
		len=size(Obj.Data.(f),d(jj)); % size of the considered dimension
		if len>1, continue; end;	% the expandable dimension is already expanded
		dN=dims{cellfun('isempty',strfind(dims,'I'))==1};
		var=bsxfun(@times,Obj.Data.(f),ones([getdim(Obj,dN) 1]));
	end
	if ~exist('var','var'), var=[]; dN=[]; end;

author: sofacoustics

repo: SOFAtoolbox


text:
function pyramid_num_test ( )

%*****************************************************************************80
%
%% PYRAMID_NUM_TEST tests PYRAMID_NUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PYRAMID_NUM_TEST\n' );
  fprintf ( 1, '  PYRAMID_NUM computes the pyramidal numbers.\n' );
  fprintf ( 1, '\n' );

  for n = 1 : 10
    fprintf ( 1, '  %2d  %6d\n', n, pyramid_num ( n ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef featuresObject
    % Version 2: Individual object for train and test
    % Working from M50 - removing unneeded code
    % - Anything unique to first test/train sets
    % - Grads
    % Moving in external functions - compile features and get filelists

    properties
        tt = 'Test' % Train or Test
        type = string('Uncompiled');
        params
        gradsParams = []
        use
        useGrads % Now vector, see set.useGrads
        fileLists
        feaNames
        labels % Training labels
        keepIdx % New keepIdx (no NaN rows)
        applyNewSafeIdx = 1;
        dataSet
        subSegLists % Not used
        divS
        files % Cell: train, test x divS
        SSL % Cell: train, test
        hybrid = 0 % Add subject numbers to general model?
        path % Path to main data directory
    end

    methods
        function obj = featuresObject(params, use)

            if isfield(params, 'tt')
                obj.tt = params.tt;
            end

            if isfield(params, 'divS')
                obj.params = params;
                % Get list of divS to use
                % set method sorts
                obj.divS = params.divS;
            else
                obj.params = 'Auto';
                obj = obj.findDivS;
            end

            if isfield(params, 'applyNewSafeIdx')
                obj.applyNewSafeIdx = params.applyNewSafeIdx;
            end

            % Check hybrid mode
            if isfield(params, 'hybrid')
                obj.hybrid = params.hybrid;
            end

            % Set paths
            obj.path = params.paths;

            % Save file lists and use
            obj = setFileLists(obj, {});
            obj.use = use;

            % Run compile now?
            % obj = compileFeatures(obj);
        end

        function obj = set.divS(obj, divS)
            obj.divS = sort(divS, 'descend');
        end

        function obj = compileFeatures(obj)

            % Check the requested files are available
            obj.files = ...
                obj.checkFiles();

            % Load files
            nD = numel(obj.divS);
            for d = 1:nD
                fn = ['divS', num2str(obj.divS(d))];

                disp(['Loading ', fn])

                data = load(obj.files{1,d});

                allData.(fn).SSL = data.subSegList;

                % Compile this data
                % Moved external function to static
                [allData.(fn).dataSet, obj.feaNames] = ...
                    obj.comboFeatures(data.feas, ...
                    data.subSegList, obj.use);
            end
            clear data

            % Join available data
            if nD>1
                % Create new table
                nt = [];
                for d = 1:nD-1

                    % Join next (d+1) to this (d)
                    d1 = ['divS', num2str(obj.divS(d))];
                    d2 = ['divS', num2str(obj.divS(d+1))];
                    disp(['Joining ', d2 ' to ', d1])

                    % Create join keys
                    % Use pedantic key generation (from M41)
                    % First get list of unique files without using unique
                    % (it sorts)
                    disp('Keygen:')
                    fns = string(allData.(d1).SSL.Files);
                    toDel = zeros(numel(fns),1);
                    for f = 2:numel(fns)
                        % If this row is same as above row, mark to delete
                        toDel(f) = fns(f).contains(fns(f-1));
                    end
                    unqFiles = fns(~toDel);
                    % Now, for this list of files, generate keys
                    % This is slow
                    k1 = [];
                    k2 = [];
                    nFiles = numel(unqFiles);
                    for f = 1:nFiles

                        k1Idx = ...
                            string(allData.(d1).SSL.Files)...
                            .contains(unqFiles(f));
                        nK1 = sum(k1Idx);

                        k2Idx = ...
                            string(allData.(d2).SSL.Files)...
                            .contains(unqFiles(f));
                        nK2 = sum(k2Idx);

                        nextIdx = length(k1)+1;

                        k1Add = nextIdx : (nextIdx-1)+nK1;
                        k2Add = round(linspace(k1Add(1),k1Add(end), nK2));

                        k1 = [k1; k1Add']; %#ok<AGROW> % Can't preallocate
                        k2 = [k2; k2Add']; %#ok<AGROW>
                        if ~mod(f,1000)
                            disp([num2str(f), '/' num2str(nFiles)])
                        end
                    end

                    % Join
                    % Simplify syntax for now (add test/train loop?)
                    % Need to drop subSegID col, and append divs to col names
                    if d==1
                        jd1 = allData.(d1).dataSet(:,2:end);
                    else
                        jd1 = nt;
                    end
                    jd1.Properties.VariableNames = ...
                        cellstr(string(jd1.Properties.VariableNames) + d1);
                    jd1.jk = k1;
                    jd2 = allData.(d2).dataSet(:,2:end);
                    jd2.Properties.VariableNames = ...
                        cellstr(string(jd2.Properties.VariableNames) + d2);
                    jd2.jk = k2;

                    nt = join(jd2, jd1);
                    % Remove key
                    nt.jk = [];

                    clear jd1 jd2

                end
            else % Just 1 divS, no joining
                % But remember to drop subSegID
                nt = ...
                    allData.(['divS', num2str(obj.divS)]).dataSet(:,2:end);
            end

            % Save to object
            obj.dataSet = nt;
            obj.feaNames = nt.Properties.VariableNames;

            % Labels should be last divS
            labs = allData.(['divS', num2str(obj.divS(end))]).SSL.Class;
            obj.labels = labs;

            % Also save final subSegList
            obj.SSL = allData.(['divS', num2str(obj.divS(end))]).SSL;

            % Set feature type
            obj = setType(obj);

            % Set safeIdx (if training)
            switch obj.tt
                case {'train', 'Train'}
                    obj = setSafeIdx(obj);
                case {'test', 'Test'}
                    % Set keepIdx
                    obj = setkeepIdx(obj);
            end

            % Add subject numbers to sets if hybrid
            if obj.hybrid
                obj = addSubs(obj);
            end

        end

        function obj = addSubs(obj)
            switch obj.dataSet.Properties.VariableNames{1}
                case 'Subject'
                    % Subject already added, do nothing
                otherwise
                    obj.dataSet = ...
                        [obj.SSL{1}(:, 'Subject'), obj.dataSet];
            end
        end

        function obj = setType(obj)
            % General/Single/Hybrid | NoGrads/Grads/PartialGrads
            % eg. Hybrid.PartialGrads
            % Model will inherit type

            if numel(unique(obj.SSL.Subject))>1 && obj.hybrid == 1
                s1 = 'Hybrid';
            else
                % Features will contain all 3 subs
                % Set model object to single later if needed
                obj.hybrid = 0;
                s1 = 'General';
            end

            if all(obj.useGrads)
                s2 = 'Grads';
            elseif any(obj.useGrads) && numel(obj.useGrads)>1
                s2 = 'PartialGrads';
            elseif ~all(obj.useGrads)
                s2 = 'NoGrads';
            end

            obj.type = string({s1, s2});

        end

        % Set file lists
        function obj = setFileLists(obj, fileLists)
            if isempty(fileLists)
                obj.fileLists = ...
                    getFileLists(obj);
            else
                obj.fileLists = fileLists;
            end
        end

        % Get file lists - static functions (were external)
        function [fileList, labels] = ...
                getFileLists(obj)
            % NB: Labels not being returned - readd for training

            switch obj.tt
                case {'train', 'Train'}
                    [fileList] = ...
                        obj.genFileListSomeSingles({'1','2','3'}, ...
                        obj.tt, obj.path);
                case {'test', 'Test'}
                    [fileList] = ...
                        obj.genFileListNoConcat({'1','2','3'}, ...
                        obj.tt, obj.path);
            end
        end

        % Not a set method
        function obj = setkeepIdx(obj)
            % Sets keepIdx

            % Set newKeepIdx
            keepIdx1 = obj.newKeepIdx(obj.dataSet);

            obj.keepIdx = keepIdx1;
        end

        % Load and apply safe idx from second release of data
        function obj = setSafeIdx(obj)
            if obj.applyNewSafeIdx

                % Load csv
                % Contains all training data and affected test data
                safe = readtable('train_and_test_data_labels_safe.csv');

                % Drop refernces to test data
                flStr = string(obj.fileLists.File);
                safeStr = string(safe.image);
                trIdx = safeStr.contains(flStr);
                safe = safe(trIdx, :);
                % Remove class column for safety
                safe.class = [];
                % Rename key variable
                safe.Properties.VariableNames{1} = 'File';

                % Readd the new files in the training set (with modified
                % names)
                % Files in fileList
                str = string(obj.fileLists.File);
                % But not in safeStr
                missingIdx = ~str.contains(safeStr);
                % Add thses to safe, along with safe tag
                nAdd = sum(missingIdx);
                add = table(cellstr(str(missingIdx)), ones(nAdd,1));
                add.Properties.VariableNames{1} = 'File';
                add.Properties.VariableNames{2} = 'safe';

                % NB: Union drops dupes, but there shouldn't be any
                safe = union(safe, add);

                % Join
                % To fileList
                nt = join(obj.fileLists, safe);

                % Order as fileList?
                % plot(nt.SubSegID); hold on; plot(obj.fileLists.SubSegID)
                if ~all(nt.SubSegID == obj.fileLists.SubSegID)
                    keyboard
                else
                    obj.fileLists = nt;
                end

                % Join
                % To biggest subSegList
                % (.File is called .Files here)
                safe.Properties.VariableNames{1} = 'Files';
                nSSL = join(obj.SSL, safe);

                % Order as fileList?
                % plot(nSSL.SubSegID); hold on; plot(obj.SSL.SubSegID)
                if ~all(nSSL.SubSegID == obj.SSL.SubSegID)
                    keyboard
                else
                    obj.SSL = nSSL;
                end

                % Set newKeepIdx
                keepIdx2 = nSSL.safe;

            else
                % Don't apply new safeIdx
                keepIdx2 = true(height(obj.SSL,1));
            end

            % Get oklist from training data and combine
            keepIdx1 = obj.newKeepIdx(obj.dataSet);

            obj.keepIdx = keepIdx1 & keepIdx2;
        end

        function divS = findDivS(obj)
            % Find available features files
            files = dir(['*', obj.tt, 'divS*_*.mat']); %#ok<*PROP>
            nFiles = numel(files);

            divS = NaN(1, nFiles);
            for f = 1:nFiles
                s = string(files(f).name);
                divS(f) = s.extractAfter('divS').extractBefore('_').double;
            end
        end

        function files = checkFiles(obj)
            nD = numel(obj.divS);
            files = cell(1, nD);

            for d = 1:nD
                % Check
                fn = [obj.tt, '_divS', num2str(obj.divS(d)), ...
                    '.mat'];
                fn = dir(fn);

                if exist(fn.name, 'file')
                    % Exists
                    files{1,d} = fn.name;
                else
                    % Run generation
                    files{1,d} = [obj.tt, '_divS', num2str(obj.divS(d)), ...
                        '.mat'];

                    obj.params.divS = obj.divS(d);
                    obj.params.fn = files{1,d};

                    [~, ~, ~] = ... Outputs unused here
                        obj.redivideSegmentsConcatSome(obj.fileLists, ...
                        obj.params, obj.tt);

                end
            end
        end
    end

    methods (Static)
        function [allData, feaNames] = ...
                comboFeatures(feas, ...
                subSegList, use)
            % Compile data tables from feature structures

            fld = fieldnames(feas.data);
            allData = NaN(size(feas.data.(fld{1}),1), 0, 'single');
            feaNames = cell(1,0);

            fns = fieldnames(use);
            nF = numel(fns);
            for f = 1:nF
                if use.(fns{f})

                    % Collect selected data
                    allData = [allData, feas.data.(fns{f})]; %#ok<AGROW>

                    % Collect names
                    feaNames = [feaNames, feas.names.(fns{f})]; %#ok<AGROW>
                end
            end

            % Make sure feaNames doesn't contain spaces
            feaNames = strrep(feaNames, ' ', '_');
            % Add SubSegID to use as key later
            feaNames = ['SubSegID', feaNames];

            % Temp bodge to allow for non-multiples (M40)
            % There may be additional row of allTrain (rounding error?)
            aUse = 1:height(subSegList);

            allData = ...
                array2table(...
                [subSegList.SubSegID, double(allData(aUse,:))]);
            allData.Properties.VariableNames = cellstr(feaNames);

        end

        function keepIdx = newKeepIdx(data)
            % All data is prcessed
            % All zero windows will return NaNs for some features
            keepIdx = all(~isnan(data{:,:}),2);
        end

        function [fullFileList] = genFileListSomeSingles(subs, str, paths)
            warning('off', 'MATLAB:table:RowsAddedExistingVars')

            vars = {...
                NaN, 'Subject', 'Subject'; ...
                NaN, 'ID', ''; ...
                NaN, 'Class', ''; ...
                NaN, 'Segment', 'The 10 min segment this is from'; ...
                NaN, 'nInSegment', 'n within segment'; ...
                NaN, 'SegID', 'n within segment'; ...
                NaN, 'SubSegID', 'n within segment'; ...
                cell(1), 'Filename', ''; ...
                cell(1), 'File', ''; ...
                NaN, 'n', 'total n'; ...
                };

            fullFileList = table(vars{:,1});
            fullFileList.Properties.VariableNames = vars(:,2);
            fullFileList.Properties.VariableDescriptions = vars(:,3);

            nSubs = numel(subs);
            % Generate file list
            tn = 0;
            for s = 1:nSubs
                % Create sub table for this subject
                fileList = table(vars{:,1});
                fileList.Properties.VariableNames = vars(:,2);
                fileList.Properties.VariableDescriptions = vars(:,3);

                sDir = [paths.dataDir, str, '_', subs{s}, '\'];

                % files = dir([sDir, subs{s}, '*']);
                files = dir([sDir, '*.mat']);
                nFilesSub = length(files);
                r = 0;
                for n = 1:nFilesSub
                    tn=tn+1;
                    fn = [sDir, files(n).name];
                    disp(['Importing file: ', fn, '(', num2str(n), ...
                        '/', num2str(nFilesSub), ')']);

                    switch str
                        case {'train', 'Train'}
                            Y = str2double(fn(end-4));
                            IDIdx = strfind(files(n).name, '_');
                            ID = files(n).name(IDIdx(1)+1:IDIdx(2)-1);
                            ID = str2double(ID);
                            % ID duplicated between 0 and 1s
                            % Add 6000 to 1 segement IDs
                            % (More than max of each subject, may still be duplicate IDs
                            % between subjects)
                            if Y==1
                                ID = ID+6000;
                            end

                        otherwise % Test
                            Y = NaN;
                            IDIdx1 = strfind(files(n).name, '_');
                            IDIdx2 = strfind(files(n).name, '.');
                            ID = files(n).name(IDIdx1+1:IDIdx2-1);
                            ID = str2double(ID);
                    end

                    r = r+1;
                    fileList.Subject(r,1) = s;
                    fileList.Filename{r,1} = fn;
                    fileList.File{r,1} = files(n).name;
                    fileList.Class(r,1) = Y;
                    fileList.ID(r,1) = ID;
                    % fileList.n(r,1) = tn;
                end

                % Subject has been added, sort by ID and label segements
                fileList = sortrows(fileList, 'ID', 'ascend');

                % In some files,
                % The are 6 files for each segment
                % In others, they are singles
                % Load list of singles for this subject
                singles = load(['singles', num2str(s), '.mat']);

                exSinglesIdx = ~string(fileList.File).contains(singles.safe4);
                nFilesSubEx = sum(exSinglesIdx);
                % Label the grouped data
                fileList.Segment(find(exSinglesIdx)) = ...
                    reshape(repmat((1:nFilesSubEx/6)', 1, 6)',nFilesSubEx,1); %#ok<FNDSB>
                fileList.nInSegment(find(exSinglesIdx)) = repmat((1:6)', ...
                    nFilesSubEx/6, 1); %#ok<FNDSB>

                % Now label the remaing as single segments
                fileList.Segment(find(~exSinglesIdx)) = ...
                    max(fileList.Segment)+1:...
                    (numel(singles.safe4)+ max(fileList.Segment)); %#ok<FNDSB>
                fileList.nInSegment(find(~exSinglesIdx)) = 1; %#ok<FNDSB>

                % Append to full table
                fullFileList = [fullFileList; fileList]; %#ok<AGROW>
                clear fileList
            end

            % Remove empty first row
            fullFileList(1,:) = [];
            fullFileList.n(:,1) = 1:tn;
            segID = cellstr([num2str(fullFileList.Subject), ...
                repmat('000', height(fullFileList),1), ...
                num2str(fullFileList.Segment)]);
            subSegID = cellstr([num2str(fullFileList.Subject), ...
                repmat('000', height(fullFileList),1), ...
                num2str(fullFileList.Segment), ...
                num2str(fullFileList.nInSegment)]);

            cellstr([segID{:}, ...
                ]);

            % FFS
            for s = 1:height(fullFileList)
                segID{s} = str2double(strrep(segID{s}, ' ', '0'));
                subSegID{s} = str2double(strrep(subSegID{s}, ' ', '0'));
            end
            fullFileList.SegID = cell2mat(segID);
            fullFileList.SubSegID = cell2mat(subSegID);
            % plot(fullFileListTrain.ID); hold on; plot(fullFileListTrain.Segment)
        end

        function [fullFileList] = genFileListNoConcat(subs, str, paths)
            warning('off', 'MATLAB:table:RowsAddedExistingVars')

            vars = {...
                NaN, 'Subject', 'Subject'; ...
                NaN, 'ID', ''; ...
                NaN, 'Class', ''; ...
                NaN, 'Segment', 'The 10 min segment this is from'; ...
                NaN, 'nInSegment', 'n within segment'; ...
                NaN, 'SegID', 'n within segment'; ...
                NaN, 'SubSegID', 'n within segment'; ...
                cell(1), 'Filename', ''; ...
                cell(1), 'File', ''; ...
                NaN, 'n', 'total n'; ...
                };

            fullFileList = table(vars{:,1});
            fullFileList.Properties.VariableNames = vars(:,2);
            fullFileList.Properties.VariableDescriptions = vars(:,3);

            nSubs = numel(subs);
            % Generate file list
            tn = 0;
            for s = 1:nSubs
                % Create sub table for this subject
                fileList = table(vars{:,1});
                fileList.Properties.VariableNames = vars(:,2);
                fileList.Properties.VariableDescriptions = vars(:,3);

                sDir = [paths.dataDir, str, '_', subs{s}, '_New\'];

                % files = dir([sDir, subs{s}, '*']);
                files = dir([sDir, '*.mat']);
                nFilesSub = length(files);
                r = 0;
                for n = 1:nFilesSub
                    tn=tn+1;
                    fn = [sDir, files(n).name];
                    disp(['Importing file: ', fn, ...
                        '(', num2str(n), '/', num2str(nFilesSub), ')']);

                    switch str
                        case {'train', 'Train'}
                            Y = str2double(fn(end-4));
                            IDIdx = strfind(files(n).name, '_');
                            ID = files(n).name(IDIdx(1)+1:IDIdx(2)-1);
                            ID = str2double(ID);
                            % ID duplicated between 0 and 1s
                            % Add 6000 to 1
                            % segement IDs (More than max of each subject,
                            % may still be duplicate IDs between subjects)
                            if Y==1
                                ID = ID+6000;
                            end

                        otherwise % Test
                            Y = NaN;
                            IDIdx1 = strfind(files(n).name, '_');
                            IDIdx2 = strfind(files(n).name, '.');
                            ID = files(n).name(IDIdx1+1:IDIdx2-1);
                            ID = str2double(ID);
                    end

                    r = r+1;
                    fileList.Subject(r,1) = s;
                    fileList.Filename{r,1} = fn;
                    fileList.File{r,1} = files(n).name;
                    fileList.Class(r,1) = Y;
                    fileList.ID(r,1) = ID;
                    % fileList.n(r,1) = tn;
                end

                % Subject has been added, sort by ID and label segements
                fileList = sortrows(fileList, 'ID', 'ascend');
                % The are 6 files for each segment
                % fileList.Segment = reshape(repmat((1:nFilesSub/6)', 1, 6)',nFilesSub,1);
                % New for master 33: Each file is individual segment
                fileList.Segment = (1:height(fileList))';
                fileList.nInSegment = repmat(1, nFilesSub, 1);

                % Append to full table
                fullFileList = [fullFileList; fileList]; %#ok<AGROW>
                clear fileList
            end

            % Remove empty first row
            fullFileList(1,:) = [];
            fullFileList.n(:,1) = 1:tn;
            segID = cellstr([num2str(fullFileList.Subject), ...
                repmat('000', height(fullFileList),1), ...
                num2str(fullFileList.Segment)]);
            subSegID = cellstr([num2str(fullFileList.Subject), ...
                repmat('000', height(fullFileList),1), ...
                num2str(fullFileList.Segment), ...
                num2str(fullFileList.nInSegment)]);

            cellstr([segID{:}, ...
                ]);

            % FFS
            for s = 1:height(fullFileList)
                segID{s} = str2double(strrep(segID{s}, ' ', '0'));
                subSegID{s} = str2double(strrep(subSegID{s}, ' ', '0'));
            end
            fullFileList.SegID = cell2mat(segID);
            fullFileList.SubSegID = cell2mat(subSegID);
            % plot(fullFileListTrain.ID); hold on; plot(fullFileListTrain.Segment)
        end

        function [subSegList, fileList, feas] = ...
                redivideSegmentsConcatSome(fileList, params, str)
            % Divide and OK check
            warning('off', 'MATLAB:table:RowsAddedExistingVars')

            % Unpack params
            divS = params.divS;
            OKThresh = 0.5;
            params.plotOn = 0;

            % First check it's necessary to run
            fn = params.fn;
            vars = {'feas', 'subSegList'};
            if exist(fn, 'file') && ~params.force
                disp('Reloading')
                load(fn, vars{:});
            else

                % Prepare table
                vars2 = {...
                    NaN, 'Subject', 'Subject'; ...
                    NaN, 'ID', ''; ...
                    NaN, 'Class', ''; ...
                    NaN, 'Segment', 'The 10 min segment this is from'; ...
                    NaN, 'nInSegment', 'n within segment'; ...
                    NaN, 'SegID', 'n within segment'; ...
                    NaN, 'SubSegID', 'n within segment'; ...
                    cell(1), 'Filenames', ''; ...
                    cell(1), 'Files', ''; ...
                    NaN, 'n', 'total n'; ...
                    NaN, 'SubN', 'n in whole segment'; ...
                    NaN, 'Of', 'total epochs in 6-file concatonated segment. Should always be the same!'; ...
                    NaN, 'Lengh', 'Length of subSeg in S'; ...
                    NaN, 'Prop0', 'Proportion of zeros'; ...
                    true, 'OK', 'Passed OK Check?'; ...
                    };

                subSegList = table(vars2{:,1});
                subSegList.Properties.VariableNames = vars2(:,2);
                subSegList.Properties.VariableDescriptions = vars2(:,3);

                % Set lengths
                % unqSub = unique(fileList.Subject);
                % nSub = numel(unqSub);
                unqSeg = unique(fileList.SegID);
                nSeg = numel(unqSeg);
                % unqSubSeg = unique(fileList.SubSegID);
                % nSubSeg = numel(unqSubSeg);

                % Known parameters - fixed, not discovered or varied.
                fs = 400;
                dl = 240000;
                nChans = 16;
                % nInSeg = 1; % Master33: Now dynamic

                % Consistent paraemters
                divPts = divS*fs;
                % nDivs = floor((dl*nInSeg)/divPts);
                % cut = dl*nInSeg - nDivs*divPts;
                % Changed master33
                %     if nDivs>6
                %         from = reshape(repmat(1:6, nDivs/6, 1), nDivs, 1);
                %     else
                %         from = (1:6)';
                %     end
                % from = ones(nDivs,1)

                % totalToDo = height(fileList);

                % Features processed here
                % tRows = height(fileList)/nInSeg*nDivs;
                % tRows = height(fileList)*ceil((10*60)/divS);

                % Count singles and blocks of 6 here
                % Concatonate .Subject and .ID as chars
                ac = string([num2str(fileList.Subject), num2str(fileList.Segment)]);
                % Convert to string and .replace any spaces, return as doubles
                ac = ac.replace(' ' ,'').double();
                % Use this as accarray index on ones
                counts = accumarray(ac, ones(height(fileList),1));

                n1s = sum(counts == 1);
                n6s = sum(counts == 6);

                row1s = n1s*ceil((10*60*1)/divS);
                row6s = n6s*ceil((10*60*6)/divS);

                tRows = row1s + row6s;

                feas.data.hillsBandsLogAv = NaN(tRows, ...
                    numel(params.HillsBands.Range), 'single');
                feas.data.hillsBandsLog2D = NaN(tRows, ...
                    numel(params.HillsBands.Range)*16, 'single');
                feas.data.maxHills2D = NaN(tRows, ...
                    16, 'single');
                feas.data.maxHillsAv = NaN(tRows, ...
                    1, 'single');
                feas.data.summ32D = NaN(tRows, ...
                    8*nChans, 'single');
                feas.data.summ3Av = NaN(tRows, ...
                    8, 'single');
                feas.data.bandsLin2D = NaN(tRows, ...
                    11*16, 'single');
                feas.data.bandsLinAv = NaN(tRows, ...
                    11, 'single');
                feas.data.maxBands2D = NaN(tRows, ...
                    16, 'single');
                feas.data.maxBandsAv = NaN(tRows, ...
                    1, 'single');
                feas.data.mCorrsT = NaN(tRows, ...
                    4*16, 'single');
                feas.data.mCorrsF = NaN(tRows, ...
                    3*16, 'single');

                sslRow = 0;
                finalFeaRow = 0;
                % for s = 1:nSeg % For each unique segment
                for s = 1:nSeg % For each file individually
                    % Get ID
                    sn = unqSeg(s);

                    % Get index of files
                    segIdx = fileList.SegID == sn;

                    % Get subFileList structure just for this segment
                    % Might be one file, might be 6
                    subFileList = fileList(segIdx,:);

                    % Moved here to make dynamic
                    nInSeg = height(subFileList);
                    nDivs = floor((dl*nInSeg)/divPts);
                    cut = dl*nInSeg - nDivs*divPts;

                    if cut>0
                        % If cut is >0 there are unused values at end
                        % Add one more window achnored to end to fill these values
                        fillWindow = true;

                        % Also add 1 to nDivs, but not yet
                    else
                        fillWindow = false;
                    end

                    from = round(linspace(1,height(subFileList),nDivs))'; % Estimated for non-multiples of length
                    % Previous - didn't work witn non-multiples eg divS = 400
                    % from = reshape(repmat(subFileList.nInSegment', ...
                    %    nDivs/nInSeg, 1), nDivs, 1);

                    % Load and concatonate the data in the segement
                    % If there's more than one file
                    data1 = NaN(dl*nInSeg, nChans, 'single');
                    files = subFileList.File';
                    filenames = subFileList.Filename';
                    nFiles = numel(files);
                    sIdx = 1;
                    eIdx = 240000;
                    for ss = 1:nFiles
                        a = load(filenames{ss});
                        if nInSeg==1
                            disp(['Not-Concatonating ', filenames{ss}])
                        else
                            disp(['Concatonating ', filenames{ss}])
                        end
                        data1(sIdx:eIdx,:) = a.dataStruct.data;

                        sIdx = sIdx + 240000;
                        eIdx = eIdx + 240000;
                    end
                    clear a

                    % Do preprocessing here

                    % Epoch
                    % Cut off any less than full bins, from start
                    % data = data(1+cut:end,:);
                    % Can increase data here by using sliding window instead

                    % Epoch - M40
                    % Now cutting end off, and if any is cut, adding this to process as
                    % start cut off (below)
                    data = data1(1:end-cut,:);

                    % Idiot check
                    % dataEp = reshape(data', 16, divPts, nDivs);
                    % dataEp = permute(dataEp, [2, 1, 3]);
                    % figure
                    % subplot(2,1,1), plot(data(1:4000, 1:6))
                    % subplot(2,1,2), plot(dataEp(:, 1:6, 1))
                    % figure
                    % subplot(2,1,1), plot(data(4001:8000, 1:6))
                    % subplot(2,1,2), plot(dataEp(:, 1:6, 2))
                    data = reshape(data', 16, divPts, nDivs);
                    data = permute(data, [2, 1, 3]);

                    if fillWindow
                        % Add the end cut off window to data
                        % fw = data1(1+cut:end,:);
                        fw = data1(end-(divPts-1):end,:);
                        fw = reshape(fw', 16, divPts, 1);
                        fw = permute(fw, [2, 1, 3]);

                        % Set new nDivs here and add window to process to 3rd dim of
                        % data
                        nDivs = nDivs+1;
                        data(:,:,nDivs) = fw;
                        % Also remake from
                        from = round(linspace(1,height(subFileList),nDivs))';
                    end

                    % For each of the new epochs, do OK check and add to table
                    nInSegment = 0;
                    for n = 1:nDivs
                        sslRow = sslRow + 1;
                        nInSegment = nInSegment+1;
                        % Get the index (of the 6 rows from fileList) that this data came
                        % from
                        frIdx = from(n);

                        % Copy in columns from subFilefileList
                        nc = width(subFileList);
                        for c = 1:nc
                            subSegList(sslRow, c) = subFileList(frIdx, c);
                        end

                        % OK Check
                        prop0 = sum(sum(data(:,:,n)==0)) / (divPts*nChans);
                        chk = prop0 < OKThresh;

                        % Add additional information
                        subSegList.SubN(sslRow,1) = n;
                        subSegList.Of(sslRow,1) = nDivs;
                        subSegList.Length(sslRow,1) = divS;
                        subSegList.Prop0(sslRow,1) = prop0;
                        subSegList.OK(sslRow,1) = chk;
                        subSegList.nInSegment(sslRow,1) = nInSegment;

                    end

                    disp('Computing Hills bands feature')
                    feaRow = finalFeaRow;
                    for n = 1:nDivs
                        feaRow = feaRow+1;
                        [feas.data.hillsBandsLog2D(feaRow, :), ...
                            feas.data.hillsBandsLogAv(feaRow, :), ...
                            feas.data.maxHills2D(feaRow, :), ...
                            feas.data.maxHillsAv(feaRow, :), ...
                            feas.names.hillsBandsLog2D, ...
                            feas.names.hillsBandsLogAv, ...
                            feas.names.maxHills2D, ...
                            feas.names.maxHillsAv] = ...
                            featuresObject.extractHillBands(data(:,:,n), ...
                            params);
                    end

                    disp('Computing summary stats features')
                    feaRow = finalFeaRow;
                    for n = 1:nDivs
                        feaRow = feaRow+1;
                        [feas.data.summ32D(feaRow, :), ...
                            feas.data.summ3Av(feaRow, :), ...
                            feas.names.summ32D, ...
                            feas.names.summ3Av] = ...
                            featuresObject.extractSumm3(data(:,:,n), ...
                            params);
                    end

                    disp('Computing bandsLin feature')
                    feaRow = finalFeaRow;
                    for n = 1:nDivs
                        feaRow = feaRow+1;
                        [feas.data.bandsLin2D(feaRow, :), ...
                            feas.data.bandsLinAv(feaRow, :), ...
                            feas.data.maxBands2D(feaRow, :), ...
                            feas.data.maxBandsAv(feaRow, :), ...
                            feas.names.bandsLin2D, ...
                            feas.names.bandsLinAv, ...
                            feas.names.maxBands2D, ...
                            feas.names.maxBandsAv] = ...
                            featuresObject.extractBandsLin(data(:,:,n), ...
                            params);
                    end

                    disp('Computing channel correlations T feature')
                    feaRow = finalFeaRow;
                    for n = 1:nDivs
                        feaRow = feaRow+1;
                        [feas.data.mCorrsT(feaRow, :), ...
                            feas.names.mCorrsT] = ...
                            featuresObject.extractChannelCorrelationT(data(:,:,n), ...
                            params);
                    end

                    disp('Computing channel correlations F feature')
                    feaRow = finalFeaRow;
                    for n = 1:nDivs
                        feaRow = feaRow+1;
                        [feas.data.mCorrsF(feaRow, :), ...
                            feas.names.mCorrsF] = ...
                            featuresObject.extractChannelCorrelationF(data(:,:,n), ...
                            params);
                    end

                    finalFeaRow = feaRow;

                    disp([str, ' Done: ', num2str(s), '/', num2str(nSeg), ...
                        '(', num2str(s/nSeg*100), '%)'])

                end

                % Reapply subSegID
                subSegList.SubSegID = ...
                    str2double(strrep(cellstr([num2str(subSegList.SegID), ...
                    num2str(subSegList.nInSegment)]), ' ', '0'));

                save(fn, vars{:}, '-v7.3')
            end
        end

        [bandsLin2D, bandsLinAv, ...
            mB2D, mBAv, ...
            names2D, namesAv, ...
            namesmB2D, namesmBAv] = extractBandsLin(data, params)

        [mCorrsF2D, names2D] = ...
            extractChannelCorrelationF(data, params)

        [mCorrsT2D, names2D] = ...
            extractChannelCorrelationT(data, params)

        [HBLog2D, HBLogAv, ...
            mB2D, mBAv, ...
            names2D, namesAv, ...
            namesmB2D, namesmBAv] = ...
            extractHillBands(data, params)

        [summ32D, summ3Av, names32D, names3Av] = extractSumm3(data, params)

    end
end

author: garethjns

repo: Kaggle-EEG


text:
function [trl] = ft_trialfun_emgdetect(cfg)

% Note that there are some parameters, like the EMG channel name and the
% processing that is done on the EMG channel data, which are hardcoded in
% this trial function. You should change these parameters if necessary.

% read the header and determine the channel number corresponding with the EMG
hdr         = ft_read_header(cfg.headerfile);
chanindx    = strmatch('EMGlft', hdr.label);

if length(chanindx)>1
  error('only one EMG channel supported');
end

% read all data of the EMG channel, assume continuous file format
emg = ft_read_data(cfg.datafile, 'header', hdr, ...
              'begsample', 1, 'endsample', hdr.nSamples*hdr.nTrials, ...
              'chanindx', chanindx, 'checkboundary', false);

% apply filtering, hilbert transformation and boxcar convolution (for smoothing)
emgflt      = ft_preproc_highpassfilter(emg, hdr.Fs, 10); % highpassfilter
emghlb      = abs(hilbert(emgflt')');                     % hilbert transform
emgcnv      = conv2([1], ones(1,hdr.Fs), emghlb, 'same'); % smooth using convolution
emgstd      = ft_preproc_standardize(emgcnv, 2);          % z-transform, i.e. mean=0 and stdev=1
emgtrl      = emgstd>0;                                   % detect the muscle activity
emgtrl      = diff(emgtrl, [], 2);

emgon       = find(emgtrl(:)== 1);
emgoff      = find(emgtrl(:)==-1);

trl(:,1) = emgon (:) + hdr.Fs*0.5;  % as a consequence of the convolution with a one-second boxcar
trl(:,2) = emgoff(:) - hdr.Fs*0.5;  % as a consequence of the convolution with a one-second boxcar
trl(:,3) = 0;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% read the NIST format wav files
% Author: Xiao Xiong
% Created: 1 Feb 2005
% Last modified: 1 Feb 2005

function [data] = readNIST(file_name,big_endian);

FILE = fopen( file_name );
if FILE <1
    error('File open failed: %s', file_name);
end
fread(FILE, 1024, 'int8');      %read the header first
if big_endian
    data = fread(FILE, 'int16=>short','b');
else
    data = fread(FILE, 'int16=>short','l');
end
% plot(data);
data = double(data(1:length(data)));        % store the vecors as column vector of wav
% plot(data);
fclose(FILE);

author: singaxiong

repo: SignalGraph


text:
function [X_mean]=wmean(X,W)

% function [X_mean]=wmean(X,W)
% ------------------------------------------------------------------------
% This function calculates the weighted mean of X using the weights defined
% in the vector W.
%
% Uses the following formula: X_mean=sum(X.*W)./sum(W);
%
% 12/09/2008
% ------------------------------------------------------------------------

%%

X=X(:);
W=W(:);
X_mean=sum(X.*W)./sum(W);

%% END

%%
% _*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