text:
function imsegs = msCreateSuperpixels(im, imname)

prefix = num2str(floor(rand(1)*1E6));
fn1 = ['./tmpim' prefix '.ppm'];
outfn = ['./tmpimsp' prefix '.ppm'];

imwrite(im, fn1);
segcmd = '/home/dhoiem/cmu/tools/superpixelsFH/segment 0.8 100 100';
syscall = [segcmd ' ' fn1 ' ' outfn];

system(syscall);
imsegs = processSuperpixelImage(outfn);
if exist('imname', 'var')
    imsegs.imname = imname;
end

delete(fn1);
delete(outfn);

author: zouchuhang

repo: LayoutNet


text:
function [P]   = spm_dcm_fit(P,use_parfor)
% Bayesian inversion of DCMs using Variational Laplace
% FORMAT [DCM] = spm_dcm_fit(P)
%
% P          - {N x M} DCM structure array (or filenames) from N subjects
% use_parfor - if true, will attempt to run in parallel (default: false)
%              NB: all DCMs are loaded into memory
%
% DCM  - Inverted (1st level) DCM structures with posterior densities
%__________________________________________________________________________
%
% This routine is just a wrapper that calls the appropriate dcm inversion
% routine for a set a pre-specifed DCMs.
%
% If called with a cell array, each column is assumed to contain 1st level
% DCMs inverted under the same model. Each row contains a different data
% set (or subject).
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_fit.m 7755 2019-12-16 13:19:28Z spm $

if nargin < 2, use_parfor = false; end

% get filenames and set up
%--------------------------------------------------------------------------
if ~nargin
    [P, sts] = spm_select([1 Inf],'^DCM.*\.mat$','Select DCM*.mat files');
    if ~sts, return; end
end
if ischar(P),   P = cellstr(P);  end
if isstruct(P), P = {P};         end

% Number of subjects (data) and models (of those data)
%--------------------------------------------------------------------------
[Ns,Nm] = size(P);

% Find model class and modality
%--------------------------------------------------------------------------
try, load(P{1}); catch, DCM = P{1}; end

model = spm_dcm_identify(DCM);

if isempty(model)
    warning('unknown inversion scheme');
    return
end

% Get data structure for each subject (column)
%--------------------------------------------------------------------------
for i = 1:Ns
    for j = 2:Nm
        switch model
            case{'DEM'}
                P{i, j}.xY = P{i, 1}.Y;
            otherwise
                P{i, j}.xY = P{i, 1}.xY;
        end
    end
end

% Estimate
%--------------------------------------------------------------------------
if use_parfor
    % Estimate DCMs in parallel using parfor
    P = spm_dcm_load(P);

    parfor i = 1:numel(P)
        P{i} = fit_dcm(P{i}, model);
    end
else
    % Estimate DCMs without parfor
    for i = 1:numel(P)
        try
            DCM = load(P{i});
            DCM = DCM.DCM;
        catch
            DCM = P{i};
        end
        P{i} = fit_dcm(DCM, model);
    end
end

% -------------------------------------------------------------------------
function DCM = fit_dcm(DCM, model)
% Inverts a DCM.
% DCM   - the DCM structure
% model - a string identifying the model type (see spm_dcm_identify)

switch model

    % fMRI model
    %----------------------------------------------------------------------
    case{'fMRI'}
        DCM = spm_dcm_estimate(DCM);

    % conventional neural-mass and mean-field models
    %----------------------------------------------------------------------
    case{'fMRI_CSD'}
        DCM = spm_dcm_fmri_csd(DCM);

    % conventional neural-mass and mean-field models
    %----------------------------------------------------------------------
    case{'ERP'}
        DCM = spm_dcm_erp(DCM);

    % cross-spectral density model (complex)
    %----------------------------------------------------------------------
    case{'CSD'}
        DCM = spm_dcm_csd(DCM);

    % cross-spectral density model (steady-state responses)
    %----------------------------------------------------------------------
    case{'TFM'}
        DCM = spm_dcm_tfm(DCM);

    % induced responses
    %----------------------------------------------------------------------
    case{'IND'}
        DCM = spm_dcm_ind(DCM);

    % phase coupling
    %----------------------------------------------------------------------
    case{'PHA'}
        DCM = spm_dcm_phase(DCM);

    % cross-spectral density model (steady-state responses)
    %----------------------------------------------------------------------
    case{'NFM'}
        DCM = spm_dcm_nfm(DCM);

    % behavioural Markov decision process model
    %----------------------------------------------------------------------
    case{'MDP'}
        DCM = spm_dcm_mdp(DCM);

    % generic nonlinear system identification
    %----------------------------------------------------------------------
    case{'NLSI'}
        [Ep,Cp,Eh,F] = spm_nlsi_GN(DCM.M,DCM.xU,DCM.xY);
        DCM.Ep       = Ep;
        DCM.Eh       = Eh;
        DCM.Cp       = Cp;
        DCM.F        = F;

    % hierarchical dynamic mmodel
    %----------------------------------------------------------------------
    case{'DEM'}
        DCM = spm_DEM(DCM);

    % default
    %----------------------------------------------------------------------
    otherwise
        try
            DCM = feval(model, DCM);
        catch
            error('unknown DCM');
        end
end

author: spm


text:
function Psi0 = biasVardistPsi0Compute(biaskern, vardist)

% BIASVARDISTPSI0COMPUTE one line description
% FORMAT
% DESC description
% RETURN Psi0 : description
% ARG biasKern : the kernel structure associated with the white kernel.
% ARG vardist : description
%
%
% SEEALSO : others
%
%
% COPYRIGHT : Michalis K. Titsias, 2009
%

% VARGPLVM

Psi0 = vardist.numData*biaskern.variance;

author: SheffieldML

repo: vargplvm


text:
function ll = gpsimMapLogLikelihood(model, df)

% GPSIMMAPLOGLIKELIHOOD Compute the log likelihood of a GPSIMMAP model.
% FORMAT
% DESC computes the log likelihood of the given Gaussian process
% for use in a single input motif protein network.
% ARG model : the model for which the log likelihood is computed.
% RETURN ll : the log likelihood of the data set.
%
% SEEALSO : gpsimMapCreate, gpsimMapLogLikeGradient, gpsimMapObjective
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% MODIFIED : Pei Gao, 2008

% SHEFFIELDML

numData = length(model.t);

% ll = model.f'*model.invK*model.f ...
%       + log(det(eye(size(model.W))+model.K*model.W));
ll = model.f'*model.invK*model.f ...
    - model.logDetCovf + model.logDetK;

if isfield(model, 'includeNoise') && model.includeNoise
  noiseMat = ones(length(model.t), 1)*model.noiseVar;
  yvar = model.yvar + noiseMat;
else
  yvar = model.yvar;
end

for i = 1:numData
  for j = 1:model.numGenes
    ind = i + (j-1)*numData;
    beta_ij = 1/yvar(ind);
    factor = (model.ypred(model.times_index(i), j)...
              - model.y(ind));
    ll = ll + factor*factor*beta_ij - log(beta_ij);
  end
end

ll = ll + numData*model.numGenes*log(2*pi);

ll = -.5*ll;

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QPSK demonstration packet-based transceiver for Chilipepper
% Toyon Research Corp.
% http://www.toyon.com/chilipepper.php
% Created 10/17/2012
% embedded@toyon.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is the top level entry function for the receiver portion of the
% example. The entire receiver is designed to run at Rate=1 (one clock
% cycle per iteration of the core.
% We follow standard receive practice with frequency offset estimation,
% pulse-shape filtering, time estimateion, and correlation to determine
% tart of packet.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%#codegen

function [r_out, s_f_out, s_c_out, s_t_out, t_est_out, f_est_out] = ...
    qpsk_rx(i_in, q_in, mu_foc_in, mu_toc_in)

    % scale input data coming from the Chilipepper ADC to be purely
    % fractional to avoid scaling issues
    r_in = complex(i_in, q_in);

    % frequency offset estimation. Note that time constant is input as
    % integer
    [s_f_i, s_f_q, fe] = qpsk_rx_foc(i_in, q_in, mu_foc_in);

    % Square-root raised-cosine band-limited filtering
    [s_c_i, s_c_q] = qpsk_rx_srrc(s_f_i, s_f_q);

    % Time offset estimation. Output data changes at the symbol rate.
    [s_t_i, s_t_q, tauh] =  qpsk_rx_toc(s_c_i, s_c_q, mu_toc_in);

    % estimation and correlation values
    t_est_out = tauh;
    f_est_out = fe;

    % original signal out (real version)
    r_out = real(r_in);

    % incremental signal outputs after frequency estimation, filtering, and
    % timining estimation
    s_f_out = complex(s_f_i,s_f_q);
    s_c_out = complex(s_c_i,s_c_q);
    s_t_out = complex(s_t_i,s_t_q);
end

author: Sable

repo: mcbench-benchmarks


text:
function cuboidhyp_omap = generate_cuboid_from_omap(omapmore, vp, OMAP_FACTOR)

[depthorder reglabel regori rpo] = omap_depthorder(omapmore, vp, OMAP_FACTOR);
cuboidhyp_omap = omapobj2cuboidhyp(rpo, reglabel, regori, vp, OMAP_FACTOR);


text:
function D = prtDistanceLNorm(dataSet1,dataSet2,Lnorm)
% prtDistanceLNorm   L Norm distance function.
%
%   DIST = prtDistanceCityBlock(DS1,DS2) calculates the LNorm distance
%   from all the observations in datasets DS1 to DS2, and ouputs a distance
%   matrix of size DS1.nObservations x DS2.nObservations. DS1 and DS2
%   should have the same number of features. DS1 and DS2 should be
%   prtDataSet objects.
%
%   For more information, see:
%
%   http://en.wikipedia.org/wiki/Norm_(mathematics)#p-norm
%
% Example:
%
%   % Create 2 data sets
%   dsx = prtDataSetStandard('Observations', [0 0; 1 1]);
%   dsy = prtDataSetStandard('Observations', [1 0;2 2; 3 3]);
%   % Compute distance
%   distance = prtDistanceLnorm(dsx,dsy)
%
% See also: prtDistanceCityBlock, prtDistanceChebychev
% prtDistanceMahalanobis, prtDistanceSquare, prtDistanceEuclidean

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A portion of IPDM from MATLAB Central is used in this function see
% prtExternal.IPDM.ipdm(). The license information from that file is below.
%
% Copyright (c) 2009, John D'Errico
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

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

[data1,data2] = prtUtilDistanceParseInputs(dataSet1,dataSet2);

% Used to handle memory efficiency paths see IPDM
chunkSize = 2^25;

[nSamples1, nDim1] = size(data1);
[nSamples2, nDim2] = size(data2);

if nDim1 ~= nDim2
    error('Dimensionality of data1 and data2 must be equal')
end

if (nDim1>1) && ((nSamples1*nSamples2*nDim1)<=chunkSize)
    switch Lnorm
        case 1
            D = sum(abs(bsxfun(@minus,reshape(data1,[nSamples1,1,nDim1]),reshape(data2,[1,nSamples2,nDim1]))),3);
        case inf
            D = max(abs(bsxfun(@minus,reshape(data1,[nSamples1,1,nDim1]),reshape(data2,[1,nSamples2,nDim1]))),[],3);
        case 0
            D = min(abs(bsxfun(@minus,reshape(data1,[nSamples1,1,nDim1]),reshape(data2,[1,nSamples2,nDim1]))),[],3);

            % This code has overflow problems for large data1 and data2
            %         case 2
            %             %un-rolled((x-y)^2)) - sqrt below; this takes less time than
            %             %the generic code below for the most common L-norm (2)
            %
            %             %D = repmat(sum((data1.^2), 2), [1 nSamples2]) + repmat(sum((data2.^2),2), [1 nSamples1]).' - 2*data1*(data2.');
            %
            %             %             %Handle overflow issues for large data2
            %             %             muData2 = prtUtilNanMean(data2);
            %             %             data2 = bsxfun(@minus,data2,muData2);
            %             %             data1 = bsxfun(@minus,data1,muData2);
            %
            %             D = bsxfun(@minus,bsxfun(@plus,sum((data1.^2), 2),sum((data2.^2),2).'),2*data1*(data2.'));
        otherwise
            D = sum(bsxfun(@minus,reshape(data1,[nSamples1,1,nDim1]),reshape(data2,[1,nSamples2,nDim1])).^Lnorm,3);
    end
else
    % too big, so that the ChunkSize will have been exceeded, or just 1-d
    if isfinite(Lnorm) && Lnorm ~= 1
        D = bsxfun(@minus,data1(:,1),data2(:,1)').^Lnorm;
    else
        D = abs(bsxfun(@minus,data1(:,1),data2(:,1)'));
    end
    for i=2:nDim1
        switch Lnorm
            case 1
                D = D + abs(bsxfun(@minus,data1(:,i),data2(:,i)'));
            case inf
                D = max(D,abs(bsxfun(@minus,data1(:,i),data2(:,i)')));
            case 0
                D = min(D,abs(bsxfun(@minus,data1(:,i),data2(:,i)')));
            otherwise
                D = D + bsxfun(@minus,data1(:,i),data2(:,i)').^Lnorm;
        end
    end
end

if isfinite(Lnorm) && Lnorm ~= 1
    if Lnorm == 2
        D = sqrt(D);
        if isreal(data1) && isreal(data2)
            D = real(D);
        end
    else
        D = D.^(1./Lnorm);
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function []=TVESLMdisp(beta_median,beta_std,beta_lbound,beta_ubound,sigma_median,sigma_t_median,sigma_t_lbound,sigma_t_ubound,X,Y,n,m,p,k,T,bex,ar,lambda1,lambda2,lambda3,lambda4,lambda5,gamma,IRFt,const,endo,exo,startdate,enddate,stringdates1,decimaldates1,pref, Ylevel, Psi_median, Psi_lbound, Psi_ubound, sizetraining,PriorExcel)

% before displaying and saving the results, start estimating the evaluation measures for the model

% obtain first a point estimate betatilde of the VAR coefficients
% this is simply the median
betatilde=beta_median;
Btilde=reshape(betatilde,k,n);
% use this estimate to produce predicted values for the model, following (a.8.2)
Ytilde=X*Btilde;
% then produce the corresponding residuals, using (a.8.3)
EPStilde=Y-Ytilde;

% check first whether the model is stationary, using (a.7.2)
[stationary eigmodulus]=bear.checkstable(betatilde,n,p,k);

% Compute then the sum of squared residuals
% compute first the RSS matrix, defined in (a.8.4)
RSS=EPStilde'*EPStilde;
% retain only the diagonal elements to get the vector of RSSi values
rss=diag(RSS);

% Go on calculating R2
% generate Mbar
Mbar=eye(T)-ones(T,T)/T;
% then compute the TSS matrix, defined in (a.8.7)
TSS=Y'*Mbar*Y;
% generate the R2 matrix in (a.8.8)
R2=eye(n)-RSS./TSS;
% retain only the diagonal elements to get the vector of R2 values
r2=diag(R2);

% then calculate the adjusted R2, using (a.8.9)
R2bar=eye(n)-((T-1)/(T-k))*(eye(n)-R2);
% retain only the diagonal elements to get the vector of R2bar values
r2bar=diag(R2bar);

% now start displaying and saving the results

% preliminary task: create and open the txt file used to save the results

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'wt');

% print toolbox header

fprintf('%s\n','');
fprintf(fid,'%s\n','');

fprintf('%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
fprintf(fid,'%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
fprintf('%s\n','%                                                                                %%');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    BAYESIAN ESTIMATION, ANALYSIS AND REGRESSION (BEAR) TOOLBOX                                         %');
fprintf(fid,'%s\n','%    BAYESIAN ESTIMATION, ANALYSIS AND REGRESSION (BEAR) TOOLBOX                                     %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    This statistical package has been developed by the external developments division of the ECB.       %');
fprintf(fid,'%s\n','%    This statistical package has been developed by the external developments division of the ECB.   %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    Authors:                                                                                %');
fprintf(fid,'%s\n','%    Authors:                                                                                %');
fprintf('%s\n','%    Romain Legrand  (Romain Legrand <b00148883@essec.edu>)                                              %');
fprintf(fid,'%s\n','%    Romain Legrand  (Romain Legrand <b00148883@essec.edu>)                                          %');
fprintf('%s\n','%    Alistair Dieppe (adieppe@worldbank.org)                                                             %');
fprintf(fid,'%s\n','%    Alistair Dieppe (adieppe@worldbank.org)                                                         %');
fprintf('%s\n','%    Björn van Roye  (Bjorn.van_Roye@ecb.europa.eu)                                                      %');
fprintf(fid,'%s\n','%    Björn van Roye  (Bjorn.van_Roye@ecb.europa.eu)                                                  %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    Version 5,0                                                                                %');
fprintf(fid,'%s\n','%    Version 5.0                                                                                %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    The authors are grateful to Paolo Bonomolo, Marta Banbura, Martin Bruns, Fabio Canova,              %');
fprintf(fid,'%s\n','%    The authors are grateful to Paolo Bonomolo, Marta Banbura, Martin Bruns, Fabio Canova,          %');
fprintf('%s\n','%    Matteo Ciccarelli, Marek Jarocinski, Niccolo Battistini, Gabriel Bobeica                            %');
fprintf(fid,'%s\n','%    Matteo Ciccarelli, Marek Jarocinski, Niccolo Battistini, Gabriel Bobeica                        %');
fprintf('%s\n','%    Michele Lenza, Chiara Osbat, Mirela Miescu, Gary Koop, Giorgio Primiceri                            %');
fprintf(fid,'%s\n','%    Michele Lenza, Chiara Osbat, Mirela Miescu, Gary Koop, Giorgio Primiceri,                       %');
fprintf('%s\n','%    Michal Rubaszek, Barbara Rossi, Ben Schumann, Peter Welz, Hugo Vega de la Cruz and Francesca Loria. %');
fprintf(fid,'%s\n','%  Michal Rubaszek, Barbara Rossi, Ben Schumann, Peter Welz, Hugo Vega de la Cruz and Francesca Loria%');
fprintf('%s\n','%    valuable input and advice which contributed to improve the quality of this work.                    %');
fprintf(fid,'%s\n','%  valuable input and advice which contributed to improve the quality of this work.                  %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%   These programmes are the responsibilities of the authors and not of the ECB and the Worldbank.       %');
fprintf(fid,'%s\n','%   These programmes are the responsibilities of the authors and not of the ECB and the Worldbank.   %');
fprintf('%s\n','%   Errors and ommissions remain those of the authors.                                                   %');
fprintf(fid,'%s\n','%   Errors and ommissions remain those of the authors.                                               %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    Please do not use or quote this work without permission.                                            %');
fprintf(fid,'%s\n','%    Please do not use or quote this work without permission.                                        %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
fprintf(fid,'%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
% print then estimation results

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

toolboxinfo='BEAR toolbox estimates';
fprintf('%s\n',toolboxinfo);
fprintf(fid,'%s\n',toolboxinfo);

time=clock;
datestring=datestr(time);
dateinfo=['Date: ' datestring(1,1:11) '   Time: ' datestring(1,13:17)];
fprintf('%s\n',dateinfo);
fprintf(fid,'%s\n',dateinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

VARtypeinfo='Stochastic volatility BVAR';
fprintf('%s\n',VARtypeinfo);
fprintf(fid,'%s\n',VARtypeinfo);

modelinfo='Survey Local Mean Stochastic Volatility VAR';
fprintf('%s\n',modelinfo);
fprintf(fid,'%s\n',modelinfo);

if IRFt==1
SVARinfo='structural decomposition: none';
elseif IRFt==2
SVARinfo='structural decomposition: choleski factorisation';
elseif IRFt==3
SVARinfo='structural decomposition: triangular factorisation';
elseif IRFt==4
SVARinfo='structural decomposition: sign restrictions';
end
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);

temp='endogenous variables: ';
for ii=1:n
temp=[temp ' ' endo{ii,1} ' '];
end
endoinfo=temp;
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

temp='exogenous variables: ';
if const==0 && m==0
temp=[temp ' none'];
elseif const==1 && m==1
temp=[temp ' constant '];
elseif const==0 && m>0
   for ii=1:m-1
   temp=[temp ' ' exo{ii,1} ' '];
   end
elseif const==1 && m>1
temp=[temp ' constant '];
   for ii=1:m-1
   temp=[temp ' ' exo{ii,1} ' '];
   end
end
exoinfo=temp;
fprintf('%s\n',exoinfo);
fprintf(fid,'%s\n',exoinfo);

sampledateinfo=['estimation sample: ' startdate '-' enddate];
fprintf('%s\n',sampledateinfo);
fprintf(fid,'%s\n',sampledateinfo);

samplelengthinfo=['Size of the Trainingsample: ' num2str(sizetraining)];
fprintf('%s\n',samplelengthinfo);
fprintf(fid,'%s\n',samplelengthinfo);

samplelengthinfo=['sample size (omitting initial conditions and Trainingsample): ' num2str(T)];
fprintf('%s\n',samplelengthinfo);
fprintf(fid,'%s\n',samplelengthinfo);

laginfo=['number of lags included in regression: ' num2str(p)];
fprintf('%s\n',laginfo);
fprintf(fid,'%s\n',laginfo);

if PriorExcel==1
    arprint=[];
    for ii=1:n
        arprint=[arprint num2str(ar(ii,1)) '  '];
    end
 hyperparam2=['autoregressive coefficients (ar):                ' arprint];
else
 hyperparam2=['autoregressive coefficients (ar):                ' num2str(ar(1,1))];
end
fprintf('%s\n',hyperparam2);
fprintf(fid,'%s\n',hyperparam2);

hyperparam3=['overall tightness (lambda1):                    ' num2str(lambda1)];
fprintf('%s\n',hyperparam3);
fprintf(fid,'%s\n',hyperparam3);

hyperparam4=['cross-variable weighting (lambda2):             ' num2str(lambda2)];
fprintf('%s\n',hyperparam4);
fprintf(fid,'%s\n',hyperparam4);

hyperparam5=['lag decay (lambda3):                            ' num2str(lambda3)];
fprintf('%s\n',hyperparam5);
fprintf(fid,'%s\n',hyperparam5);

%hyperparam6=['exogenous variable tightness (lambda4):         ' num2str(lambda4)];
%fprintf('%s\n',hyperparam6);
%fprintf(fid,'%s\n',hyperparam6);

if bex==1
hyperparam7=['block exogeneity shrinkage (lambda5):           ' num2str(lambda5)];
fprintf('%s\n',hyperparam7);
fprintf(fid,'%s\n',hyperparam7);
end

hyperparam8=['AR coefficient of stochastic volatiltiy of VAR residuals(gamma):    ' num2str(gamma)];
fprintf('%s\n',hyperparam8);
fprintf(fid,'%s\n',hyperparam8);

% display coefficient estimates

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

coeffinfo=['VAR coefficients (beta): posterior estimates'];
fprintf('%s\n',coeffinfo);
fprintf(fid,'%s\n',coeffinfo);

for ii=1:n

fprintf('%s\n','');
fprintf(fid,'%s\n','');
if ii~=1
fprintf('%s\n','');
fprintf(fid,'%s\n','');
end

endoinfo=['Endogenous: ' endo{ii,1}];
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

coeffheader=fprintf('%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');
coeffheader=fprintf(fid,'%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');

% handle the endogenous
   for jj=1:n
      for kk=1:p
      values=[beta_median((ii-1)*k+n*(kk-1)+jj,1) beta_std((ii-1)*k+n*(kk-1)+jj,1) beta_lbound((ii-1)*k+n*(kk-1)+jj,1) beta_ubound((ii-1)*k+n*(kk-1)+jj,1)];
      fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
      fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
      end
   end

% handle the exogenous
   % if there is no constant:
   if const==0
      % if there is no exogenous at all, obvioulsy, don't display anything
      if m==0
      % if there is no constant but some other exogenous, display them
      else
         for jj=1:m
         values=[beta_median(ii*k-m+jj,1) beta_std(ii*k-m+jj,1) beta_lbound(ii*k-m+jj,1) beta_ubound(ii*k-m+jj,1)];
         fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
         fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
         end
      end
   % if there is a constant
   else
   % display the results related to the constant
         values=[beta_median(ii*k-m+1,1) beta_std(ii*k-m+1,1) beta_lbound(ii*k-m+1,1) beta_ubound(ii*k-m+1,1)];
         fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
         fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
      % if there is no other exogenous, stop here
      if m==1
      % if there are other exogenous, display their results
      else
         for jj=1:m-1
         values=[beta_median(ii*k-m+jj+1,1) beta_std(ii*k-m+jj+1,1) beta_lbound(ii*k-m+jj+1,1) beta_ubound(ii*k-m+jj+1,1)];
         fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
         fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
         end
      end
   end

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display evaluation measures
rssinfo=['Sum of squared residuals: ' num2str(rss(ii,1),'%.2f')];
fprintf('%s\n',rssinfo);
fprintf(fid,'%s\n',rssinfo);

r2info=['R-squared: ' num2str(r2(ii,1),'%.3f')];
fprintf('%s\n',r2info);
fprintf(fid,'%s\n',r2info);

adjr2info=['adj. R-squared: ' num2str(r2bar(ii,1),'%.3f')];
fprintf('%s\n',adjr2info);
fprintf(fid,'%s\n',adjr2info);

end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display VAR stability results
eigmodulus=reshape(eigmodulus,p,n);
stabilityinfo1=['Roots of the characteristic polynomial (modulus):'];
fprintf('%s\n',stabilityinfo1);
fprintf(fid,'%s\n',stabilityinfo1);
for ii=1:p
temp=num2str(eigmodulus(ii,1),'%.3f');
   for jj=2:n
   temp=[temp,'  ',num2str(eigmodulus(ii,jj),'%.3f')];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end
if stationary==1;
stabilityinfo2=['No root lies outside the unit circle.'];
stabilityinfo3=['The estimated VAR model satisfies the stability condition'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
else
stabilityinfo2=['Warning: at leat one root lies on or outside the unit circle.'];
stabilityinfo3=['The estimated VAR model will not be stable'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display posterior for sigma
sigmainfo1=['sigma (residual covariance matrix): posterior estimates'];
fprintf('%s\n',sigmainfo1);
fprintf(fid,'%s\n',sigmainfo1);
% calculate the (integer) length of the largest number in sigma, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(sigma_median))))));
% add a separator, a potential minus sign, and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:n
temp=[];
   for jj=1:n
   % convert matrix entry into string
   number=num2str(sigma_median(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end
% notice about use of long-run values
sigmainfo2=['Note: Last period observation variance covariance matrix'];
sigmainfo3=['Note: The last period observation variance covariance matrix will also be used for identification.'];
%sigmainfo4=['Note: The historical decomposition and structural shocks are dependend on this choice.'];

fprintf('%s\n',sigmainfo2);
fprintf(fid,'%s\n',sigmainfo2);
fprintf('%s\n',sigmainfo3);
fprintf(fid,'%s\n',sigmainfo3);
%fprintf('%s\n',sigmainfo4);
%fprintf(fid,'%s\n',sigmainfo4);
fclose(fid);

% Finally, display the results in terms of graph
if pref.plot
% plot actual vs. fitted
actualfitted=figure('Tag','BEARresults');
set(actualfitted,'Color',[0.9 0.9 0.9]);
set(actualfitted,'name','model estimation: actual vs fitted')
ncolumns=ceil(n^0.5);
nrows=ceil(n/ncolumns);
for ii=1:n
subplot(nrows,ncolumns,ii)
hold on
plot(decimaldates1,Y(:,ii),'Color',[0 0 0],'LineWidth',2);
plot(decimaldates1,Ytilde(:,ii),'Color',[1 0 0],'LineWidth',2);
hold off
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
title(endo{ii,1},'FontName','Times New Roman','FontSize',10,'FontWeight','normal');
   if ii==1
   plotlegend=legend('actual','fitted');
   set(plotlegend,'FontName','Times New Roman');
   end
end

% plot the residuals
residuals=figure('Tag','BEARresults');
set(residuals,'Color',[0.9 0.9 0.9]);
set(residuals,'name','model estimation: residuals')
for ii=1:n
subplot(nrows,ncolumns,ii)
plot(decimaldates1,EPStilde(:,ii),'Color',[0 0 0],'LineWidth',2)
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
title(endo{ii,1},'FontName','Times New Roman','FontSize',10,'FontWeight','normal');
end

% plot first the time-varying variance and covariance estimates
varcov=figure('Tag','BEARresults');
set(varcov,'Color',[0.9 0.9 0.9]);
set(varcov,'name','model estimation: residual variance and covariance')
for ii=1:n
   for jj=1:ii
   subplot(n,n,n*(ii-1)+jj)
   hold on
   Xpatch=[decimaldates1' fliplr(decimaldates1')];
   Ypatch=[sigma_t_lbound{ii,jj}(p+1:end,1)' fliplr(sigma_t_ubound{ii,jj}(p+1:end,1)')];
   HDpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
   set(HDpatch,'facealpha',0.6);
   set(HDpatch,'edgecolor','none');
   plot(decimaldates1,sigma_t_median{ii,jj}(p+1:end,1),'Color',[0.4 0.4 1],'LineWidth',2);
   plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
   hold off
   set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
      % top labels
      if jj==ii
      title(['var(' endo{ii,1} ')'],'FontWeight','normal');
      else
      title(['cov(' endo{jj,1} ',' endo{ii,1} ')'],'FontWeight','normal');
      end
   end
end

end % pref.plot

% finally, save the results on excel

% compute the cell for actual and fitted
% create the cell that will be saved on excel
afcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
% loop over variables (horizontal dimension)
for ii=1:n
% create cell of actual/fitted for variable ii
temp=['actual and fitted: ' endo{ii,1}];
af_i=[temp {''} {''};{''} {''} {''};{''} {'sample'} {'fitted'};stringdates1 num2cell(Y(:,ii)) num2cell(Ytilde(:,ii))];
afcell=[afcell af_i vertspace];
end
% trim
afcell=afcell(:,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),afcell,'actual fitted','B2');
end

% then compute the cell for the residuals
% create the cell that will be saved on excel
horzspace=repmat({''},1,n);
rescell=[{'residuals'} horzspace;{''} horzspace;{''} endo';stringdates1 num2cell(EPStilde)];
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),rescell,'resids','B2');
end

% finally compute the cell for the time varying variance and covariance
% create the cell that will be saved on excel
varcovcell={};
% build preliminary elements: space between the tables
horzspace=repmat({''},2,5*n);
vertspace=repmat({''},T+3,1);
% loop over variables (vertical dimension)
for ii=1:n
tempcell={};
   % loop over shocks (horizontal dimension)
   for jj=1:ii
   % create cell of hd record for the contribution of shock jj in variable ii fluctuation
      % if a sign restriction identification scheme has been used, use the structural shock labels
      if jj==ii
      temp=['variance of ' endo{ii,1} ' residuals'];
      % otherwise, the shocks are just orthogonalised shocks from the variables: use variable names
      else
      temp=['covariance between ' endo{jj,1} ' and ' endo{ii,1} 'residuals'];
      end
   vc_ij=[temp {''} {''} {''};{''} {''} {''} {''};{''} {'lw. bound'} {'median'} {'up. bound'};stringdates1 num2cell(sigma_t_lbound{ii,jj}(p+1:end,:)) num2cell(sigma_t_median{ii,jj}(p+1:end,:)) num2cell(sigma_t_ubound{ii,jj}(p+1:end,:))];
   tempcell=[tempcell vc_ij vertspace];
   end
   % complete with blanks (for repeated covariance values)
   for jj=1:n-ii
   tempcell=[tempcell cell(T+3,5)];
   end
varcovcell=[varcovcell;horzspace;tempcell];
end
% trim
varcovcell=varcovcell(3:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),varcovcell,'time variation','B2');
end
%trend and cycle decomposition
%plot
% Finally, display the results in terms of graph
if pref.plot
% plot actual vs. fitted
trendcycle=figure('Tag','BEARresults');
set(trendcycle,'Color',[0.9 0.9 0.9]);
set(trendcycle,'name','Local mean (trend) and actual data')
ncolumns=ceil(n^0.5);
nrows=ceil(n/ncolumns);
for ii=1:n
subplot(nrows,ncolumns,ii)
hold on
Xpatch=[decimaldates1' fliplr(decimaldates1')];
Ypatch=[Psi_lbound(p+1:end,ii)' fliplr(Psi_ubound(p+1:end,ii)')];
HDpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
set(HDpatch,'facealpha',0.6);
set(HDpatch,'edgecolor','none');
plot(decimaldates1,Psi_median(p+1:end,ii),'Color',[1 0 0],'LineWidth',2);
plot(decimaldates1,Ylevel(p+1:end,ii),'Color',[0 0 0],'LineWidth',2);
hold off
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
title(endo{ii,1},'FontName','Times New Roman','FontSize',10,'FontWeight','normal');
   if ii==1
   plotlegend=legend('Credible Bands Local Mean', 'Median Estimate Local Mean (Trend component)','Actual Data');
   set(plotlegend,'FontName','Times New Roman');
   % user defined position (ex : right bottom)
   set(plotlegend,'Position',[0.7443    0.1782    0.1518    0.1286])
   % automatic
   %set(plotlegend,'Location','best')
   end
end

%finally record the local mean estimates in excel

% compute the cell for actual and fitted
% create the cell that will be saved on excel
trendcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
% loop over variables (horizontal dimension)
for ii=1:n
% create cell of actual/fitted for variable ii
temp=['Local mean estimates: ' endo{ii,1}];
trendcell_i=[temp {''} {''} {''};{''} {''} {''} {''};{'dates'} {'upper bound'} {'median'} {'lower bound'};stringdates1 num2cell(Psi_ubound(p+1:end,ii)) num2cell(Psi_median(p+1:end,ii)) num2cell(Psi_lbound(p+1:end,ii))];
trendcell=[trendcell trendcell_i vertspace];
end
% trim
trendcell=trendcell(:,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),trendcell,'Local Mean Estimates','B2');
end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function l=lowerbound(dvec, w, alpha, beta, a, b)
% l=lowerbound(dvec, w, alpha, beta, a, b) (1 x I)
% Computes lower bound for variational approximation (Buntine & Jakulin
% DCA 2006)
% dvec: data - each column is one image (JxI)
% w: basis - each column is one basis vector (component, psf); (JxK)
% alpha, beta: parameters of the Gamma distribution for latent h (1xK each)
% a, b: parameters of the approximative Gamma distribution (KxI each)
% J=#pixels
% K=#components
% I=#images.

% this is the same as in update_n ...
% expectation of <log(l_k)>:
e_l = psi(a) - log(b); % (K x I)

% n-update:
[ntmp, z]=update_ntmp(w, a, b); % z (Jx1xI) normalization constant from n update

t1=sum(e_l.*bsxfun(@minus, alpha', a),1);   %(1xI)
t2=sum(dvec.*log(squeeze(z)),1);            %(1xI)

% approximation of the log-gamma function
t31=gammalogapprox(a)-a.*log(b);            % (KxI)
t32=gammalogapprox(alpha)-alpha.*log(beta); % (1xK)
t3=sum(bsxfun(@minus, t31, t32'),1);        % (1xI)
% This t4 is not really necessary as it is only data dependent:
t4=sum(factorialapprox(dvec),1);            % (1xI)
l=t1+t2+t3-t4; %(1xI)

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function f = local_xy ( x , y )

%*****************************************************************************80
%
%% LOCAL_XY computes the local minimum function.
%
%  Discussion:
%
%    This function has a local minimum:
%
%      X* = (  0.28581..., 0.27936...), F(X*) = 5.9225...
%
%    and a global minimum:
%
%      X* = (  -21.026653..., -36.760090...), F(X*) = 0.
%
%    Suggested starting point:
%
%      X = ( 1, 1 ), F(X) = 3.33 * 10^6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Himmelblau,
%    Applied Nonlinear Programming,
%    McGraw Hill, 1972,
%    ISBN13: 978-0070289215,
%    LC: T57.8.H55.
%
%  Parameters:
%
%    Input, real X, Y, the argument of the function.
%
%    Output, real F, the value of the function at (X,Y).
%
  f = ( x^2 + 12 * y - 1 )^2 ...
    + ( 49 * x^2 + 49 * y^2 + 84 * x + 2324 * y - 681 )^2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function val = h5read_safe(fname, name, val)

try
    val = h5read(fname, name);
catch
    switch(class(val))
        case 'char'
            val = '';
        case 'cell'
            val = {};
        otherwise
            val = [];
    end
end

author: BUNPC


text:
function [document, scores] = report_ar_rank(context, experiment, trackers, sequences, varargin)
% report_ar_rank Generate a report based on A-R ranking
%
% Performs A-R ranking analysis and generates a report based on the results.
%
% Input:
% - context (structure): Report context structure.
% - experiment (struct): An experiment structure.
% - trackers (cell): An array of tracker structures.
% - sequences (cell): An array of sequence structures.
% - varargin[UsePractical] (boolean): Use practical difference.
% - varargin[UseTags] (boolean): Rank according to tags (otherwise rank according to sequences).
% - varargin[Average] (string): How to compute average rank.
%     - weighted_mean: Average ranks, average values by taking into account length
%     - mean: Average ranks, average values
%     - pooled: gather all frames and compute ranking on a single combined sequence
% - varargin[Alpha] (double): Statistical significance parameter.
% - varargin[Adaptation] (string): Statistical significance parameter.
% - varargin[HideLegend] (boolean): Hide legend in plots.
%
% Output:
% - document (structure): Resulting document structure.
% - scores (struct): Averaged ranks for entire set.
%

usetags = get_global_variable('report_tags', true);
usepractical = get_global_variable('report_ar_practical', true);
orderingplot = get_global_variable('report_ar_ordering', true);
hidelegend = get_global_variable('report_legend_hide', false);
arplot = get_global_variable('report_ar_arplot', true);
average = get_global_variable('report_ar_average', 'weighted_mean');
adaptation = get_global_variable('report_ar_adaptation', 'mean');
sensitivity = get_global_variable('report_ar_sensitivity', 30);
alpha = get_global_variable('report_ar_alpha', 0.05);
table_format = get_global_variable('report_ar_table_format', 'accrob'); % joined, rankscores, accrob, fragmented
table_orientation = get_global_variable('report_ar_table_orientation', 'trackers'); % trackers, selectors, trackerscores, selectorscores

for i = 1:2:length(varargin)
    switch lower(varargin{i})
        case 'usepractical'
            usepractical = varargin{i+1};
        case 'usetags'
            usetags = varargin{i+1};
        case 'average'
            average = varargin{i+1};
        case 'alpha'
            alpha = varargin{i+1};
        case 'adaptation'
            adaptation = varargin{i+1};
        case 'hidelegend'
            hidelegend = varargin{i+1};
        otherwise
            error(['Unknown switch ', varargin{i}, '!']) ;
    end
end

if ~strcmp(experiment.type, 'supervised')
   error('A-R Ranking analysis only suitable for supervised experiments!');
end

if numel(trackers) < 2
    error('Ranking analysis requires two or more trackers.');
end;

if ~any(strcmp(average, {'mean', 'weighted_mean', 'pooled'}))
   error('Unknown averaging technique "%s"!', average);
end
document = document_create(context, 'ar_rank', 'title', 'AR ranking');

trackers_hash = md5hash(strjoin((cellfun(@(x) x.identifier, trackers, 'UniformOutput', false)), '-'), 'Char', 'hex');
parameters_hash = md5hash(sprintf('%f-%d-%d-%s', alpha, usetags, usepractical, adaptation));

tags = {};

if usetags && isfield(experiment, 'tags')
    tags = union(experiment.tags, {'all'});
    sequences_hash = md5hash(strjoin(tags, '-'), 'Char', 'hex');
else
    sequences_hash = md5hash(strjoin((cellfun(@(x) x.name, sequences, 'UniformOutput', false)), '-'), 'Char', 'hex');
end;

cache_identifier = sprintf('ranking_%s_%s_%s_%s', experiment.name, trackers_hash, sequences_hash, parameters_hash);

result = document_cache(context, cache_identifier, @analyze_accuracy_robustness, experiment, trackers, ...
    sequences, 'tags', tags, 'usepractical', usepractical, ...
    'alpha', alpha, 'adaptation', adaptation);

if usetags
    % When using tags we have inserted a separate one for this
    mask = strcmp('tag_all', result.tags);

    result.accuracy.pooled_values = result.accuracy.values(mask, :);
    result.robustness.pooled_values = result.robustness.values(mask, :);
    result.robustness.pooled_normalized = result.robustness.normalized(mask, :);
    result.accuracy.pooled_ranks = result.accuracy.ranks(mask, :);
    result.robustness.pooled_ranks = result.robustness.ranks(mask, :);

    % Now remove the 'all' tag from results
    result.accuracy.values = result.accuracy.values(~mask, :);
    result.accuracy.ranks = result.accuracy.ranks(~mask, :);
    result.robustness.values = result.robustness.values(~mask, :);
    result.robustness.normalized = result.robustness.normalized(~mask, :);
    result.robustness.ranks = result.robustness.ranks(~mask, :);
    result.lengths = result.lengths(~mask);
    result.tags = result.tags(~mask);

end

result.accuracy.weighted_mean_ranks = nanmean(result.accuracy.ranks, 1);
result.robustness.weighted_mean_ranks = nanmean(result.robustness.ranks, 1);

useable = result.lengths > 0;
result.accuracy.weighted_mean_values = sum(result.accuracy.values(useable, :) ...
    .* repmat(result.lengths(useable), 1, length(trackers)), 1) ./ sum(result.lengths(useable));
result.robustness.weighted_mean_values = sum(result.robustness.values(useable, :) ...
    .* repmat(result.lengths(useable), 1, length(trackers)), 1) ./ sum(result.lengths(useable));
result.robustness.weighted_mean_normalized = sum(result.robustness.normalized(useable, :) ...
    .* repmat(result.lengths(useable), 1, length(trackers)), 1) ./ sum(result.lengths(useable));

result.accuracy.mean_ranks = nanmean(result.accuracy.ranks, 1);
result.robustness.mean_ranks = nanmean(result.robustness.ranks, 1);

result.accuracy.mean_values = nanmean(result.accuracy.values, 1);
result.robustness.mean_values = nanmean(result.robustness.values, 1);
result.robustness.mean_normalized = nanmean(result.robustness.normalized, 1);

switch average

    case 'weighted_mean'

        values = cat(2, result.accuracy.weighted_mean_ranks', ...
            result.robustness.weighted_mean_ranks', ...
            result.accuracy.weighted_mean_values', ...
            result.robustness.weighted_mean_values');

    case 'mean'

        values = cat(2, result.accuracy.mean_ranks', ...
            result.robustness.mean_ranks', ...
            result.accuracy.mean_values', ...
            result.robustness.mean_values');

    case 'pool'

        if usetags
            values = cat(2, result.accuracy.pooled_ranks', ...
                result.robustness.pooled_ranks', ...
                result.accuracy.pooled_values', ...
                result.robustness.pooled_values');
        else
            values = cat(2, result.accuracy.weighted_mean_ranks', ...
                result.robustness.weighted_mean_ranks', ...
                result.accuracy.weighted_mean_values', ...
                result.robustness.weighted_mean_values');
        end
end

scores.name = 'A-R rank';
scores.values = values;
scores.ids = {'arank', 'rrank', 'overlap', 'failures'};
scores.names = {'A-Rank', 'R-Rank', 'Overlap', 'Failures'};
scores.order = {'ascending', 'ascending', 'descending', 'ascending'};

tracker_labels = cellfun(@(x) iff(isfield(x.metadata, 'verified') && x.metadata.verified, [x.label, '*'], x.label), trackers, 'UniformOutput', 0);

print_text('Writing ranking table ...');

document.section('Experiment %s', experiment.name);

if arplot

    ar_and_rank_plot(document, sprintf('%s_mean', experiment.name), ...
        sprintf('experiment %s (mean)', experiment.name), ...
        trackers, result.accuracy.mean_ranks, ...
        result.robustness.mean_ranks, ...
        result.accuracy.mean_values, ...
        result.robustness.mean_normalized, sensitivity, hidelegend);

    ar_and_rank_plot(document, sprintf('%s_weighted_mean', experiment.name), ...
        sprintf('experiment %s (weighted_mean)', experiment.name), ...
        trackers, result.accuracy.weighted_mean_ranks, ...
        result.robustness.weighted_mean_ranks, ...
        result.accuracy.weighted_mean_values, ...
        result.robustness.weighted_mean_normalized, sensitivity, hidelegend);

    if usetags

        ar_and_rank_plot(document, sprintf('%s_pooled', experiment.name), ...
        sprintf('experiment %s (pooled)', experiment.name), ...
        trackers, result.accuracy.pooled_ranks, ...
        result.robustness.pooled_ranks, ...
        result.accuracy.pooled_values, ...
        result.robustness.pooled_normalized, sensitivity, hidelegend);

    end

end;

selector_tags = result.tags;

score_tags = {'A-Rank', 'R-Rank', 'Overlap', 'Failures'};
score_sorting = {'ascending', 'ascending', 'descending', 'ascending'};
raw_scores = cat(3, result.accuracy.ranks', ...
    result.robustness.ranks', ...
    result.accuracy.values', ...
    result.robustness.values');
raw_scores = cat(2, raw_scores, cat(3, result.accuracy.mean_ranks', ...
    result.robustness.mean_ranks', ...
    result.accuracy.mean_values', ...
    result.robustness.mean_values'));
raw_scores = cat(2, raw_scores, cat(3, result.accuracy.weighted_mean_ranks', ...
    result.robustness.weighted_mean_ranks', ...
    result.accuracy.weighted_mean_values', ...
    result.robustness.weighted_mean_values'));

table_selector_tags = selector_tags;
table_selector_tags{end+1} = create_table_cell('Mean', 'class', 'average');
table_selector_tags{end+1} = create_table_cell('Weighted mean', 'class', 'average');

if usetags
    raw_scores = cat(2, raw_scores, cat(3, result.accuracy.pooled_ranks', ...
        result.robustness.pooled_ranks', ...
        result.accuracy.pooled_values', ...
        result.robustness.pooled_values'));
    table_selector_tags{end+1} = create_table_cell('Pooled', 'class', 'average');
end

switch table_format
    case 'joined'
        print_scores_table(document, raw_scores, score_sorting, score_tags, tracker_labels, table_selector_tags, table_orientation, 'Ranks and raw scores');
    case 'rankscores'
        print_scores_table(document, raw_scores(:, :, 1:2), score_sorting(1:2), score_tags(1:2), tracker_labels, table_selector_tags, table_orientation, 'Ranks');
        print_scores_table(document, raw_scores(:, :, 3:4), score_sorting(3:4), score_tags(3:4), tracker_labels, table_selector_tags, table_orientation, 'Raw scores');
    case 'accrob'
        print_scores_table(document, raw_scores(:, :, [1,3]), score_sorting([1,3]), score_tags([1,3]), tracker_labels, table_selector_tags, table_orientation, 'Accuracy');
        print_scores_table(document, raw_scores(:, :, [2,4]), score_sorting([2,4]), score_tags([2,4]), tracker_labels, table_selector_tags, table_orientation, 'Robustness');
    case 'fragmented'
        for t = 1:numel(score_tags)
            print_scores_table(document, raw_scores(:, :, t), score_sorting(t), score_tags(t), tracker_labels, table_selector_tags, table_orientation, score_tags{t});
        end;
end

document.subsection('Detailed plots');

if orderingplot

    h = plot_ordering(trackers, result.accuracy.ranks, selector_tags, ...
        'flip', 1, 'legend', ~hidelegend);
    document.figure(h, sprintf('ordering_accuracy_%s', experiment.name), ...
        'Ranking orderings for accuracy rank');

    close(h);

    h = plot_ordering(trackers, result.accuracy.values, selector_tags, ...
        'scope', [0, 1], 'type', 'Overall overlap', 'legend', ~hidelegend);
    document.figure(h, sprintf('ordering_overlap_%s', experiment.name), ...
        'Orderings for overall overlap');

    close(h);

    h = plot_ordering(trackers, result.robustness.ranks, selector_tags, ...
        'flip', 1, 'legend', ~hidelegend);
    document.figure(h, sprintf('ordering_robustness_%s', experiment.name), ...
        'Ranking orderings for robustness rank');

    close(h);

    robustness = result.robustness.normalized .* sensitivity;

    h = plot_ordering(trackers, robustness, selector_tags, ...
        'scope', [0, max(robustness(:))+eps], 'type', ...
        'Failures', 'legend', ~hidelegend);

    document.figure(h, sprintf('ordering_failures_%s', experiment.name), ...
        'Orderings for failures');

    close(h);

end;

if arplot

    for l = 1:length(selector_tags)

        plot_title = sprintf('Ranking plot for tag %s in experiment %s', ...
            selector_tags{l}, experiment.name);
        plot_id = sprintf('rankingplot_%s_%s', ...
            experiment.name, selector_tags{l});

        hf = plot_ranking(trackers, result.accuracy.ranks(l, :)', ...
            result.robustness.ranks(l, :)', ...
            'title', plot_title, 'limit', numel(trackers), 'legend', ~hidelegend);

        document.figure(hf, plot_id, plot_title);

        close(hf);

        plot_title = sprintf('AR plot for tag %s in experiment %s', ...
            selector_tags{l}, experiment.name);
        plot_id = sprintf('arplot_%s_%s', experiment.name, selector_tags{l});

        hf = plot_ar(trackers, result.accuracy.values(l, :), ...
            result.robustness.normalized(l, :), ...
            'title', plot_title, 'sensitivity', sensitivity, 'legend', ~hidelegend);

        document.figure(hf, plot_id, plot_title);

        close(hf);

    end;

end;

document.write();

end

% --------------------------------------------------------------------- %

function print_scores_table(document, scores, score_sorting, score_tags, tracker_labels, selector_tags, orientation, title)

    % Scores - selectors x trackers x scores

    score_count = numel(score_tags);
    selector_count = numel(selector_tags);
    tracker_count = numel(tracker_labels);

    switch orientation
        case 'trackers'
            row_tags = tracker_labels;
            column_labels = selector_tags;
            row_scores = false;
            sort_columns = false;
        case 'selectors'
            row_tags = selector_tags;
            column_labels = tracker_labels;
            row_scores = false;
            sort_columns = true;
        case 'trackerscores'
            row_tags = tracker_labels;
            column_labels = selector_tags;
            row_scores = true;
            sort_columns = true;
        case 'selectorscores'
            row_tags = selector_tags;
            column_labels = tracker_labels;
            row_scores = true;
            sort_columns = false;
        otherwise
            error('Unknown format %s', orientation);
    end

    column_labels = column_labels(:)';
    row_tags = row_tags(:);

    row_count = numel(row_tags);
    column_count = numel(column_labels);

    if sort_columns
        table_data = cell(tracker_count * score_count, selector_count);

        for s = 1:score_count
            score_table_data = highlight_best_rows(num2cell(scores(:, :, s)), ...
                repmat(score_sorting(s), 1, numel(selector_tags)));
            table_data(s:score_count:end, :) = score_table_data;
        end

        if ~row_scores
            table_data = table_data';
        end
    else
        table_data = cell(tracker_count, selector_count * score_count);

        for s = 1:score_count
            score_table_data = highlight_best_rows(num2cell(scores(:, :, s)), ...
                repmat(score_sorting(s), 1, numel(selector_tags)));
            table_data(:, s:score_count:end) = score_table_data;
        end

        if row_scores
            table_data = table_data';
        end
    end

    if row_scores

        if score_count > 1
            row_tags_exp = cell(score_count * row_count, 2);
            row_tags_exp(:, 1) = repmat({struct()}, 1, row_count * score_count);
            row_tags_exp(1:score_count:end, 1) = cellfun(@(x) create_table_cell(x, 'rows', score_count), row_tags, 'UniformOutput', false);
            row_tags_exp(:, 2) = score_tags(repmat(1:score_count, 1, row_count));
            row_tags = row_tags_exp;
        end

    else

        if score_count > 1
            column_labels_exp = cell(2, score_count * column_count);
            column_labels_exp(1, :) = repmat({struct()}, 1, column_count * score_count);
            column_labels_exp(1, 1:score_count:end) = cellfun(@(x) create_table_cell(x, 'columns', score_count), column_labels, 'UniformOutput', false);
            column_labels_exp(2, :) = score_tags(repmat(1:score_count, 1, column_count));
            column_labels = column_labels_exp;
        end

    end;

    document.table(table_data, 'columnLabels', column_labels, 'rowLabels', row_tags, 'title', title);

end

function ar_and_rank_plot(document, identifier, title, trackers, ...
    accuracy_ranks, robustness_ranks, accuracy_values, ...
    robustness_values, sensitivity, hidelegend)

    plot_title = sprintf('Ranking plot for %s', title);
    plot_id = sprintf('rankingplot_%s', identifier);

    hf = plot_ranking(trackers, accuracy_ranks, ...
        robustness_ranks, ...
        'title', plot_title, 'limit', numel(trackers), 'legend', ~hidelegend);

    document.figure(hf, plot_id, plot_title);

    close(hf);

    plot_title = sprintf('AR plot for %s', title);
    plot_id = sprintf('arplot_%s', identifier);

    hf = plot_ar(trackers, accuracy_values, ...
        robustness_values, ...
        'title', plot_title, 'sensitivity', sensitivity, 'legend', ~hidelegend);

    document.figure(hf, plot_id, plot_title);

    close(hf);

end

author: votchallenge

repo: toolkit-legacy


text:
function dt = GetLastCheckForUpdates()
if ~ispathvalid([getAppDir, 'LastCheckForUpdates.dat'])
    try
        dt = datetime - duration(200,0,0);
    catch
        dt = -1;
    end
else
    fd = fopen([getAppDir, 'LastCheckForUpdates.dat'],'rt');
    dt = fgetl(fd);
    fclose(fd);
end


text:
function [model, B, elapse] = BRE_learn(A, maxbits)
%   This is a wrapper function of Binary Reconstructive Embedding learning.
%
%	Usage:
%	[model, B,elapse] = BRE_learn(A, maxbits)
%
%	      A: Rows of vectors of data points. Each row is sample point
%   maxbits: Code length
%
%     model: Used for encoding a test sample point.
%	      B: The binary code of the input data A. Each row is sample point
%    elapse: The coding time (training time).
%
%
%
%   version 2.0 --Nov/2016
%   version 1.0 --Jan/2013
%
%   Written by  Yue Lin (linyue29@gmail.com)
%               Deng Cai (dengcai AT gmail DOT com)
%

tmp_T = tic;

trainIdx=randperm(size(A,1));
Xtrain=A(trainIdx(1:1000),:);
K = Xtrain*Xtrain';

params=[];
params.disp = 0;
params.n = size(K,1);
params.numbits = maxbits;
params.K = K;
params.hash_size = 20;
hash_inds = zeros(params.hash_size,params.numbits);
for b = 1:params.numbits
	rp = randperm(params.n);
	hash_inds(:,b) = rp(1:params.hash_size)';
end
params.hash_inds = hash_inds;
W0 = .001*randn(params.hash_size,params.numbits);
[W,H] = BRE(W0,params);

model.W=W;
model.X=Xtrain;
model.hash_inds=hash_inds;

B=BRE_compress(A,model);

elapse = toc(tmp_T);
end

author: ZJULearning

repo: MatlabFunc


text:
function output_info=move_output_info_cpu(output_info)

if output_info.is_group_data

    data_child_groups=output_info.data_child_groups;
    output_valid_flags=output_info.child_valid_flags;
    data_child_groups=data_child_groups(output_valid_flags);
    child_num=length(data_child_groups);

    for g_idx=1:child_num
        data_child_groups{g_idx}=move_output_info_cpu(data_child_groups{g_idx});
    end

    output_info.data_child_groups(output_valid_flags)=data_child_groups;

else

    if isfield(output_info, 'x')
        output_info.x=move_to_cpu_values(output_info.x);
    end
    if isfield(output_info, 'aux')
        output_info.aux=move_to_cpu_values(output_info.aux);
    end
    if isfield(output_info, 'dzdx')
        output_info.dzdx=move_to_cpu_values(output_info.dzdx);
    end
    if isfield(output_info, 'dzdw')
        output_info.dzdw=move_to_cpu_values(output_info.dzdw);
    end

end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function out = spm_shoot_divergence(job)
% Compute divergences from velocity fields
% FORMAT spm_shoot_divergence(job)
% job.velocities - Filenames of initial velocity fields
%
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_shoot_divergence.m 7387 2018-08-03 15:13:57Z john $

P = strvcat(job.velocities);
out = cell(size(P,1),1);
for i=1:size(P,1)
    Nii   = nifti(deblank(P(i,:)));
    d     = size(Nii.dat);
    krn   = {[-1;0;1]/2,[-1,0,1]/2,reshape([-1,0,1],[1 1 3])/2};
    dv    = zeros(d(1:3));
    for dm=1:3
        dv = dv - convn(Nii.dat(:,:,:,1,dm),krn{dm},'same');
    end
    dv(:,:,1)=0; dv(:,:,end)=0;
    dv(:,1,:)=0; dv(:,end,:)=0;
    dv(1,:,:)=0; dv(end,:,:)=0;
    Nio            = Nii;
    [pth,nam,ext]  = fileparts(deblank(P(i,:)));
    Nio.dat.fname  = fullfile(pth,['d' nam(1:end) ext]);
    Nio.dat.dim    = d(1:3);
    Nio.descrip    = 'Divergence Field';
    create(Nio);
    Nio.dat(:,:,:) = dv;
    out{i}         = Nio.dat.fname;
end


%% Experiment with the cnn_mnist_fc_bnorm

[net_bn, info_bn] = cnn_mnist(...
  'expDir', 'data/mnist-bnorm', 'batchNormalization', true);

[net_fc, info_fc] = cnn_mnist(...
  'expDir', 'data/mnist-baseline', 'batchNormalization', false);

figure(1) ; clf ;
subplot(1,2,1) ;
semilogy(info_fc.val.objective', 'o-') ; hold all ;
semilogy(info_bn.val.objective', '+--') ;
xlabel('Training samples [x 10^3]'); ylabel('energy') ;
grid on ;
h=legend('BSLN', 'BNORM') ;
set(h,'color','none');
title('objective') ;
subplot(1,2,2) ;
plot(info_fc.val.error', 'o-') ; hold all ;
plot(info_bn.val.error', '+--') ;
h=legend('BSLN-val','BSLN-val-5','BNORM-val','BNORM-val-5') ;
grid on ;
xlabel('Training samples [x 10^3]'); ylabel('error') ;
set(h,'color','none') ;
title('error') ;
drawnow ;


text:
function [phi, Q, sigma_sq, mu0, sigma0] = mstep_lds_update(y, M, xt_n, Pt_n, Ptt1_n)

% Adapted from code written by Hrishi Deshpande

% M-step as described in Shumway, R.H. & Stoffer, D.S. (1982), "An approach to time series smoothing and forecasting using the EM algorithm",
% Journal of Time Series Analysis, 3, 253-264

q = size(y, 1);
n = size(y, 2)-1;
p = size(M{2}, 2);

%- - - - - - - - - -
% Eqns (9)-(11)
A = zeros(p); B = zeros(p); C = zeros(p);
for t=1:n
   tIdx = t+1;
   A = A + Pt_n{tIdx-1} + xt_n(:,tIdx-1) * xt_n(:,tIdx-1)';
   B = B + Ptt1_n{tIdx} + xt_n(:,tIdx) * xt_n(:,tIdx-1)';
   C = C + Pt_n{tIdx} + xt_n(:,tIdx) * xt_n(:,tIdx)';
end

%- - - - - - - - - -
% Eqns (12)-(14)
invA = inv(A);
phi = B * invA;
Q = (C - B*invA*B') / n;
R = zeros(q);
for t=1:n
   tIdx = t+1;
   R = R + (y(:,tIdx)-M{tIdx}*xt_n(:,tIdx)) * (y(:,tIdx)-M{tIdx}*xt_n(:,tIdx))' + ...
      M{tIdx}*Pt_n{tIdx}*M{tIdx}';
end
R = R / n;
R = eye(q)*mean(diag(R));
sigma_sq = R(1,1);
mu0 = xt_n(:, 1);            % (Eq 22, Ghahramani 1996)
sigma0 = Pt_n{1};            % (Eq 24, Ghahramani 1996)

author: TadasBaltrusaitis

repo: OpenFace


text:
function [p,pss] = propagate_bounds_from_qualities(p)
pss = [];
if isempty(p.nonlinear) && p.feasible
    if p.K.f>0
        Aeq = -p.F_struc(1:p.K.f,2:end);
        beq = p.F_struc(1:p.K.f,1);
        A = [Aeq;-Aeq];
        b = [beq;-beq];
        [p.lb,p.ub,redundant,pss] = tightenbounds(A,b,p.lb,p.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1));
    end
    pss=[];
    if p.K.l>0
        A = -p.F_struc(1+p.K.f:p.K.f+p.K.l,2:end);
        b = p.F_struc(1+p.K.f:p.K.f+p.K.l,1);
        [p.lb,p.ub,redundant,pss] = tightenbounds(A,b,p.lb,p.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1));
        if length(redundant)>0
            pss.AL0A(redundant,:)=[];
            pss.AG0A(redundant,:)=[];
            p.F_struc(p.K.f+redundant,:)=[];
            p.K.l = p.K.l - length(redundant);
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [seq, ground_truth] = load_video_info(video_path, video_name)

ground_truth = dlmread([video_path '/groundtruth_rect.txt']);
% ground_truth = dlmread([video_path, '/', video_name, '_gt', '.txt']);

seq.format = 'otb';
seq.len = size(ground_truth, 1);
seq.init_rect = ground_truth(1,:);

img_path = [video_path '/img/'];

if exist([img_path num2str(1, '%04i.png')], 'file'),
    img_files = num2str((1:seq.len)', [img_path '%04i.png']);
elseif exist([img_path num2str(1, '%04i.jpg')], 'file'),
    img_files = num2str((1:seq.len)', [img_path '%04i.jpg']);
elseif exist([img_path num2str(1, '%04i.bmp')], 'file'),
    img_files = num2str((1:seq.len)', [img_path '%04i.bmp']);
else
    error('No image files to load.')
end

seq.s_frames = cellstr(img_files);

end

repo: UDT


text:
%% ============================Notation============================ %%
% X_sub_super
% q_ToFrom
% p_ofWhat_expressedInWhatFrame

%% =============================Setup============================== %%
clear;
close all;
clc;
addpath('utils');

tic
dataDir = '../datasets';

% fileName = 'dataset3';
% fileName = 'dataset3_fresh_10noisy';
% fileName = 'dataset3_fresh_10lessnoisy';
% fileName = 'dataset3_fresh_20lessnoisy';
% fileName = 'dataset3_fresh_40lessnoisy';
% fileName = 'dataset3_fresh_60lessnoisy';
% fileName = 'dataset3_fresh_80lessnoisy';
% fileName = 'dataset3_fresh_100lessnoisy';
% fileName = 'dataset3_fresh_500lessnoisy';
% fileName = '2011_09_26_drive_0035_sync_KLT';
% fileName = '2011_09_26_drive_0005_sync_KLT';
% fileName = '2011_09_30_drive_0020_sync_KLT';
% fileName = '2011_09_26_drive_0027_sync_KLT';
% fileName = '2011_09_30_drive_0020_sync_KLT'; kStart = 2; kEnd = 900;

% Good KITTI runs
% fileName = '2011_09_26_drive_0001_sync_KLT'; kStart = 2; kEnd = 98;
fileName = '2011_09_26_drive_0036_sync_KLT'; kStart = 2; kEnd = 239;
% fileName = '2011_09_26_drive_0051_sync_KLT'; kStart = 2; kEnd = 114;
% fileName = '2011_09_26_drive_0095_sync_KLT'; kStart = 2; kEnd = 139;

load(sprintf('%s/%s.mat',dataDir,fileName));

% r_i_vk_i = p_vi_i;

%Dataset window bounds
% kStart = 2; kEnd = 177;
% kStart = 1215; kEnd = 1715;

%Set constant
numLandmarks = size(y_k_j,3);

%Set up the camera parameters
camera.c_u      = cu;                   % Principal point [u pixels]
camera.c_v      = cv;                   % Principal point [v pixels]
camera.f_u      = fu;                   % Focal length [u pixels]
camera.f_v      = fv;                   % Focal length [v pixels]
camera.b        = b;                    % Stereo baseline [m]
camera.q_CI     = rotMatToQuat(C_c_v);  % 4x1 IMU-to-Camera rotation quaternion
camera.p_C_I    = rho_v_c_v;            % 3x1 Camera position in IMU frame

%Set up the noise parameters
y_var = 11^2 * ones(1,4);               % pixel coord var
noiseParams.u_var_prime = y_var(1)/camera.f_u^2;
noiseParams.v_var_prime = y_var(2)/camera.f_v^2;

w_var = 4e-2 * ones(1,3);               % rot vel var
v_var = 4e-2 * ones(1,3);               % lin vel var
dbg_var = 1e-6 * ones(1,3);            % gyro bias change var
dbv_var = 1e-6 * ones(1,3);            % vel bias change var
noiseParams.Q_imu = diag([w_var, dbg_var, v_var, dbv_var]);

q_var_init = 1e-6 * ones(1,3);         % init rot var
p_var_init = 1e-6 * ones(1,3);         % init pos var
bg_var_init = 1e-6 * ones(1,3);        % init gyro bias var
bv_var_init = 1e-6 * ones(1,3);        % init vel bias var
noiseParams.initialIMUCovar = diag([q_var_init, bg_var_init, bv_var_init, p_var_init]);

% MSCKF parameters
msckfParams.minTrackLength = 10;        % Set to inf to dead-reckon only
msckfParams.maxTrackLength = Inf;      % Set to inf to wait for features to go out of view
msckfParams.maxGNCostNorm  = 1e-2;     % Set to inf to allow any triangulation, no matter how bad
msckfParams.minRCOND       = 1e-12;
msckfParams.doNullSpaceTrick = true;
msckfParams.doQRdecomp = true;

% IMU state for plotting etc. Structures indexed in a cell array
imuStates = cell(1,numel(t));
prunedStates = {};

% imuStates{k}.q_IG         4x1 Global to IMU rotation quaternion
% imuStates{k}.p_I_G        3x1 IMU Position in the Global frame
% imuStates{k}.b_g          3x1 Gyro bias
% imuStates{k}.b_v          3x1 Velocity bias
% imuStates{k}.covar        12x12 IMU state covariance

% We don't really need these outside of msckfState, do we?
% camState = cell(1,numel(t));
% camStates{k}.q_CG        4x1 Global to camera rotation quaternion
% camStates{k}.p_C_G       3x1 Camera Position in the Global frame
% camStates{k}.trackedFeatureIds  1xM List of feature ids that are currently being tracked from that camera state

%msckfState.imuState
%msckfState.imuCovar
%msckfState.camCovar
%msckfState.imuCamCovar
%msckfState.camStates

% Measurements as structures all indexed in a cell array
dT = [0, diff(t)];
measurements = cell(1,numel(t));
% groundTruthStates = cell(1,numel(t));
% groundTruthMap = rho_i_pj_i;

% Important: Because we're idealizing our pixel measurements and the
% idealized measurements could legitimately be -1, replace our invalid
% measurement flag with NaN
y_k_j(y_k_j == -1) = NaN;

for state_k = kStart:kEnd
    measurements{state_k}.dT    = dT(state_k);                      % sampling times
    measurements{state_k}.y     = squeeze(y_k_j(1:2,state_k,:));    % left camera only
    measurements{state_k}.omega = w_vk_vk_i(:,state_k);             % ang vel
    measurements{state_k}.v     = v_vk_vk_i(:,state_k);             % lin vel

    %Idealize measurements
    validMeas = ~isnan(measurements{state_k}.y(1,:));
    measurements{state_k}.y(1,validMeas) = (measurements{state_k}.y(1,validMeas) - camera.c_u)/camera.f_u;
    measurements{state_k}.y(2,validMeas) = (measurements{state_k}.y(2,validMeas) - camera.c_v)/camera.f_v;

    %Ground Truth
    q_IG = rotMatToQuat(axisAngleToRotMat(theta_vk_i(:,state_k)));
    p_I_G = r_i_vk_i(:,state_k);

    groundTruthStates{state_k}.imuState.q_IG = q_IG;
    groundTruthStates{state_k}.imuState.p_I_G = p_I_G;

    % Compute camera pose from current IMU pose
    C_IG = quatToRotMat(q_IG);
    q_CG = quatLeftComp(camera.q_CI) * q_IG;
    p_C_G = p_I_G + C_IG' * camera.p_C_I;

    groundTruthStates{state_k}.camState.q_CG = q_CG;
    groundTruthStates{state_k}.camState.p_C_G = p_C_G;

end

%Struct used to keep track of features
featureTracks = {};
trackedFeatureIds = [];

% featureTracks = {track1, track2, ...}
% track.featureId
% track.observations

%% ==========================Initial State======================== %%
%Use ground truth for first state and initialize feature tracks with
%feature observations
%Use ground truth for the first state

firstImuState.q_IG = rotMatToQuat(axisAngleToRotMat(theta_vk_i(:,kStart)));
firstImuState.p_I_G = r_i_vk_i(:,kStart);
% firstImuState.q_IG = [0;0;0;1];
% firstImuState.q_IG = rotMatToQuat(rotx(90));
% firstImuState.p_I_G = [0;0;0];

[msckfState, featureTracks, trackedFeatureIds] = initializeMSCKF(firstImuState, measurements{kStart}, camera, kStart, noiseParams);
imuStates = updateStateHistory(imuStates, msckfState, camera, kStart);
msckfState_imuOnly{kStart} = msckfState;

%% ============================MAIN LOOP========================== %%

numFeatureTracksResidualized = 0;
map = [];

for state_k = kStart:(kEnd-1)
    fprintf('state_k = %4d\n', state_k);

    %% ==========================STATE PROPAGATION======================== %%

    %Propagate state and covariance
    msckfState = propagateMsckfStateAndCovar(msckfState, measurements{state_k}, noiseParams);
    msckfState_imuOnly{state_k+1} = propagateMsckfStateAndCovar(msckfState_imuOnly{state_k}, measurements{state_k}, noiseParams);

    %Add camera pose to msckfState
    msckfState = augmentState(msckfState, camera, state_k+1);

    %% ==========================FEATURE TRACKING======================== %%
    % Add observations to the feature tracks, or initialize a new one
    % If an observation is -1, add the track to featureTracksToResidualize
    featureTracksToResidualize = {};

    for featureId = 1:numLandmarks
        %IMPORTANT: state_k + 1 not state_k
        meas_k = measurements{state_k+1}.y(:, featureId);

        outOfView = isnan(meas_k(1,1));

        if ismember(featureId, trackedFeatureIds)

            if ~outOfView
                %Append observation and append id to cam states
                featureTracks{trackedFeatureIds == featureId}.observations(:, end+1) = meas_k;

                %Add observation to current camera
                msckfState.camStates{end}.trackedFeatureIds(end+1) = featureId;
            end

            track = featureTracks{trackedFeatureIds == featureId};

            if outOfView ...
                    || size(track.observations, 2) >= msckfParams.maxTrackLength ...
                    || state_k+1 == kEnd

                %Feature is not in view, remove from the tracked features
                [msckfState, camStates, camStateIndices] = removeTrackedFeature(msckfState, featureId);

                %Add the track, with all of its camStates, to the
                %residualized list
                if length(camStates) >= msckfParams.minTrackLength
                    track.camStates = camStates;
                    track.camStateIndices = camStateIndices;
                    featureTracksToResidualize{end+1} = track;
                end

                %Remove the track
                featureTracks = featureTracks(trackedFeatureIds ~= featureId);
                trackedFeatureIds(trackedFeatureIds == featureId) = [];
            end

        elseif ~outOfView && state_k+1 < kEnd % && ~ismember(featureId, trackedFeatureIds)
            %Track new feature
            track.featureId = featureId;
            track.observations = meas_k;
            featureTracks{end+1} = track;
            trackedFeatureIds(end+1) = featureId;

            %Add observation to current camera
            msckfState.camStates{end}.trackedFeatureIds(end+1) = featureId;
        end
    end

    %% ==========================FEATURE RESIDUAL CORRECTIONS======================== %%
    if ~isempty(featureTracksToResidualize)
        H_o = [];
        r_o = [];
        R_o = [];

        for f_i = 1:length(featureTracksToResidualize)

            track = featureTracksToResidualize{f_i};

            %Estimate feature 3D location through Gauss Newton inverse depth
            %optimization
            [p_f_G, Jcost, RCOND] = calcGNPosEst(track.camStates, track.observations, noiseParams);
            % Uncomment to use ground truth map instead
%              p_f_G = groundTruthMap(:, track.featureId); Jcost = 0; RCOND = 1;
%              p_f_C = triangulate(squeeze(y_k_j(:, track.camStates{1}.state_k, track.featureId)), camera); Jcost = 0; RCOND = 1;

            nObs = size(track.observations,2);
            JcostNorm = Jcost / nObs^2;
            fprintf('Jcost = %f | JcostNorm = %f | RCOND = %f\n',...
                Jcost, JcostNorm,RCOND);

            if JcostNorm > msckfParams.maxGNCostNorm ...
                    || RCOND < msckfParams.minRCOND
%                     || norm(p_f_G) > 50

                break;
            else
                map(:,end+1) = p_f_G;
                numFeatureTracksResidualized = numFeatureTracksResidualized + 1;
                fprintf('Using new feature track with %d observations. Total track count = %d.\n',...
                    nObs, numFeatureTracksResidualized);
            end

            %Calculate residual and Hoj
            [r_j] = calcResidual(p_f_G, track.camStates, track.observations);
            R_j = diag(repmat([noiseParams.u_var_prime, noiseParams.v_var_prime], [1, numel(r_j)/2]));
            [H_o_j, A_j, H_x_j] = calcHoj(p_f_G, msckfState, track.camStateIndices);

            % Stacked residuals and friends
            if msckfParams.doNullSpaceTrick
                H_o = [H_o; H_o_j];

                if ~isempty(A_j)
                    r_o_j = A_j' * r_j;
                    r_o = [r_o ; r_o_j];

                    R_o_j = A_j' * R_j * A_j;
                    R_o(end+1 : end+size(R_o_j,1), end+1 : end+size(R_o_j,2)) = R_o_j;
                end

            else
                H_o = [H_o; H_x_j];
                r_o = [r_o; r_j];
                R_o(end+1 : end+size(R_j,1), end+1 : end+size(R_j,2)) = R_j;
            end
        end

        if ~isempty(r_o)
            % Put residuals into their final update-worthy form
            if msckfParams.doQRdecomp
                [T_H, Q_1] = calcTH(H_o);
                r_n = Q_1' * r_o;
                R_n = Q_1' * R_o * Q_1;
            else
                T_H = H_o;
                r_n = r_o;
                R_n = R_o;
            end

            % Build MSCKF covariance matrix
            P = [msckfState.imuCovar, msckfState.imuCamCovar;
                   msckfState.imuCamCovar', msckfState.camCovar];

            % Calculate Kalman gain
            K = (P*T_H') / ( T_H*P*T_H' + R_n ); % == (P*T_H') * inv( T_H*P*T_H' + R_n )

            % State correction
            deltaX = K * r_n;
            msckfState = updateState(msckfState, deltaX);

            % Covariance correction
            tempMat = (eye(12 + 6*size(msckfState.camStates,2)) - K*T_H);
%             tempMat = (eye(12 + 6*size(msckfState.camStates,2)) - K*H_o);

            P_corrected = tempMat * P * tempMat' + K * R_n * K';

            msckfState.imuCovar = P_corrected(1:12,1:12);
            msckfState.camCovar = P_corrected(13:end,13:end);
            msckfState.imuCamCovar = P_corrected(1:12, 13:end);

%             figure(1); clf; imagesc(deltaX); axis equal; axis ij; colorbar;
%             drawnow;

        end

    end

        %% ==========================STATE HISTORY======================== %%
        imuStates = updateStateHistory(imuStates, msckfState, camera, state_k+1);

        %% ==========================STATE PRUNING======================== %%
        %Remove any camera states with no tracked features
        [msckfState, deletedCamStates] = pruneStates(msckfState);

        if ~isempty(deletedCamStates)
            prunedStates(end+1:end+length(deletedCamStates)) = deletedCamStates;
        end

%         if max(max(msckfState.imuCovar(1:12,1:12))) > 1
%             disp('omgbroken');
%         end

        plot_traj;
%     figure(1); imagesc(msckfState.imuCovar(1:12,1:12)); axis equal; axis ij; colorbar;
%     drawnow;
end %for state_K = ...

toc

%% ==========================PLOT ERRORS======================== %%
kNum = length(prunedStates);
p_C_G_est = NaN(3, kNum);
p_I_G_imu = NaN(3, kNum);
p_C_G_imu = NaN(3, kNum);
p_C_G_GT = NaN(3, kNum);
theta_CG_err = NaN(3,kNum);
theta_CG_err_imu = NaN(3,kNum);
err_sigma = NaN(6,kNum); % cam state is ordered as [rot, trans]
err_sigma_imu = NaN(6,kNum);
%
tPlot = NaN(1, kNum);
%
for k = 1:kNum
    state_k = prunedStates{k}.state_k;

    p_C_G_GT(:,k) = groundTruthStates{state_k}.camState.p_C_G;
    p_C_G_est(:,k) = prunedStates{k}.p_C_G;
    q_CG_est  = prunedStates{k}.q_CG;

    theta_CG_err(:,k) = crossMatToVec( eye(3) ...
                    - quatToRotMat(q_CG_est) ...
                        * ( C_c_v * axisAngleToRotMat(theta_vk_i(:,kStart+k-1)) )' );

    err_sigma(:,k) = prunedStates{k}.sigma;
    imusig = sqrt(diag(msckfState_imuOnly{state_k}.imuCovar));
    err_sigma_imu(:,k) = imusig([1:3,10:12]);

    p_I_G_imu(:,k) = msckfState_imuOnly{state_k}.imuState.p_I_G;
    C_CG_est_imu = C_CI * quatToRotMat(msckfState_imuOnly{state_k}.imuState.q_IG);
    theta_CG_err_imu(:,k) = crossMatToVec( eye(3) ...
                    - C_CG_est_imu ...
                        * ( C_CI * axisAngleToRotMat(theta_vk_i(:,kStart+k-1)) )' );

    tPlot(k) = t(state_k);
end

% p_I_G_GT = p_vi_i(:,kStart:kEnd);
p_I_G_GT = r_i_vk_i(:,kStart:kEnd);
p_C_G_GT = p_I_G_GT + repmat(rho_v_c_v,[1,size(p_I_G_GT,2)]);
p_C_G_imu = p_I_G_imu + repmat(rho_v_c_v,[1,size(p_I_G_imu,2)]);

rotLim = [-0.5 0.5];
transLim = [-0.5 0.5];

% Save estimates
msckf_trans_err = p_C_G_est - p_C_G_GT;
msckf_rot_err = theta_CG_err;
imu_trans_err = p_C_G_imu - p_C_G_GT;
imu_rot_err = theta_CG_err_imu;
save(sprintf('../KITTI Trials/msckf_%s', fileName));

armse_trans_msckf = mean(sqrt(sum(msckf_trans_err.^2, 1)/3));
armse_rot_msckf = mean(sqrt(sum(msckf_rot_err.^2, 1)/3));
final_trans_err_msckf = norm(msckf_trans_err(:,end));

armse_trans_imu = mean(sqrt(sum(imu_trans_err.^2, 1)/3));
armse_rot_imu = mean(sqrt(sum(imu_rot_err.^2, 1)/3));
final_trans_err_imu = norm(imu_trans_err(:,end));

fprintf('Trans ARMSE: IMU %f, MSCKF %f\n',armse_trans_imu, armse_trans_msckf);
fprintf('Rot ARMSE: IMU %f, MSCKF %f\n',armse_rot_imu, armse_rot_msckf);
fprintf('Final Trans Err: IMU %f, MSCKF %f\n',final_trans_err_imu, final_trans_err_msckf);

% Translation Errors
figure
subplot(3,1,1)
plot(tPlot, p_C_G_est(1,:) - p_C_G_GT(1,:), 'LineWidth', 2)
hold on
plot(tPlot, 3*err_sigma(4,:), '--r')
plot(tPlot, -3*err_sigma(4,:), '--r')
% ylim(transLim)
xlim([tPlot(1) tPlot(end)])
title('Translational Error')
ylabel('\delta r_x')

subplot(3,1,2)
plot(tPlot, p_C_G_est(2,:) - p_C_G_GT(2,:), 'LineWidth', 2)
hold on
plot(tPlot, 3*err_sigma(5,:), '--r')
plot(tPlot, -3*err_sigma(5,:), '--r')
% ylim(transLim)
xlim([tPlot(1) tPlot(end)])
ylabel('\delta r_y')

subplot(3,1,3)
plot(tPlot, p_C_G_est(3,:) - p_C_G_GT(3,:), 'LineWidth', 2)
hold on
plot(tPlot, 3*err_sigma(6,:), '--r')
plot(tPlot, -3*err_sigma(6,:), '--r')
% ylim(transLim)
xlim([tPlot(1) tPlot(end)])
ylabel('\delta r_z')
xlabel('t_k')

% Rotation Errors
figure
subplot(3,1,1)
plot(tPlot, theta_CG_err(1,:), 'LineWidth', 2)
hold on
plot(tPlot, 3*err_sigma(1,:), '--r')
plot(tPlot, -3*err_sigma(1,:), '--r')
ylim(rotLim)
xlim([tPlot(1) tPlot(end)])
title('Rotational Error')
ylabel('\delta \theta_x')

subplot(3,1,2)
plot(tPlot, theta_CG_err(2,:), 'LineWidth', 2)
hold on
plot(tPlot, 3*err_sigma(2,:), '--r')
plot(tPlot, -3*err_sigma(2,:), '--r')
ylim(rotLim)
xlim([tPlot(1) tPlot(end)])
ylabel('\delta \theta_y')

subplot(3,1,3)
plot(tPlot, theta_CG_err(3,:), 'LineWidth', 2)
hold on
plot(tPlot, 3*err_sigma(3,:), '--r')
plot(tPlot, -3*err_sigma(3,:), '--r')
ylim(rotLim)
xlim([tPlot(1) tPlot(end)])
ylabel('\delta \theta_z')
xlabel('t_k')

author: utiasSTARS

repo: msckf-swf-comparison


text:
function [] = show_plot(gen, pop, varargin)
%   This function plots the front, stolen from KeLi and Zhanghu's code.
%   gen = generation count
%   pop = population at generation gen
%   optional parameters:
%       is_text = will show values for each point on the plot, boolean
%       xlist = list of variables to be plotted for the pareto set
%       flist = list of objectives to plot, when nobj > 3

    global nreal ;
    global nobj ;

    p = inputParser;
    addRequired(p, 'gen', @isnumeric);
    addRequired(p, 'pop', @ismatrix);
    addOptional(p, 'is_text', false, @islogical);
    addOptional(p, 'xlist', [], @(x) isvector(x) || isempty(x));
    addOptional(p, 'flist', [], @(x) isvector(x) || isempty(x));
    addOptional(p, 'xlim_', [], @(x) isvector(x) || isempty(x));
    addOptional(p, 'ylim_', [], @(x) isvector(x) || isempty(x));
    addOptional(p, 'zlim_', [], @(x) isvector(x) || isempty(x));
    parse(p, gen, pop, varargin{:});

    if(nobj == 2)
        f = p.Results.pop(:, nreal + [1 2]);
    elseif(isempty(p.Results.flist))
        f = p.Results.pop(:, nreal + [1 2 3]);
        flist = [1 2 3] ;
    else
        f = p.Results.pop(:, nreal + p.Results.flist);
        flist = p.Results.flist ;
    end

    if(~isempty(p.Results.xlist))
        xlist = p.Results.xlist ;
        x = p.Results.pop(:, nreal + xlist);
    end

    fstr = sprintf('objective space, gen = %d', p.Results.gen);
    xstr = sprintf('variable space, gen = %d', p.Results.gen);

    if(p.Results.gen == 1)
        if(~isempty(p.Results.xlist))
            cpos = get(gcf, 'Position');
            set(gcf,'Position',[cpos(1) cpos(2) 1000 400]);
        end
    end

    hold off;
    if(~isempty(p.Results.xlist))
        subplot(1,2,1);
    end

    % plot objectives
    if nobj == 2
        plot(f(:,1), f(:,2), 'ro', 'MarkerSize', 4);
        xlabel('f1', 'FontSize', 6);
        ylabel('f2', 'FontSize', 6);
        if(~isempty(p.Results.xlim_)); xlim(p.Results.xlim_); end;
        if(~isempty(p.Results.ylim_)); ylim(p.Results.ylim_); end;
        if(p.Results.is_text)
            strValues = strtrim(cellstr(num2str([f(:,1) f(:,2)],'(%.4f,%.4f)')));
            text(f(:,1), f(:,2), strValues, 'VerticalAlignment', 'bottom');
        end
    else
        plot3(f(:,1), f(:,2), f(:,3), ...
                            'ro', 'MarkerSize', 4);
        xlabel(sprintf('f%d', flist(1)), 'FontSize', 6);
        ylabel(sprintf('f%d', flist(2)), 'FontSize', 6);
        zlabel(sprintf('f%d', flist(3)), 'FontSize', 6);
        if(~isempty(p.Results.xlim_)); xlim(p.Results.xlim_); end;
        if(~isempty(p.Results.ylim_)); ylim(p.Results.ylim_); end;
        if(~isempty(p.Results.zlim_)); zlim(p.Results.zlim_); end;
        if(p.Results.is_text)
            strValues = strtrim(cellstr(num2str([f(:,1) f(:,2) f(:,3)],...
                                            '(%.4f,%.4f,%.4f)')));
            text(f(:,1), f(:,2), f(:,3), ...
                            strValues, 'VerticalAlignment', 'bottom');
        end
    end
    title(fstr, 'FontSize', 8);
    box on;
    drawnow;

    % plot variables
    if(~isempty(p.Results.xlist))
        subplot(1,2,2);
        if(length(xlist) == 2)
            plot(x(:,1), x(:,2), 'ro', 'MarkerSize', 4);
            xlabel('x1', 'FontSize', 6);
            ylabel('x2', 'FontSize', 6);
            if(p.Results.is_text)
                strValues = strtrim(cellstr(num2str([x(:,1) x(:,2)],...
                                            '(%.4f,%.4f)')));
                text(x(:,1), x(:,2), strValues, ...
                                'VerticalAlignment', 'bottom');
            end
        else
            plot3(x(:,1), x(:,2), x(:,3), ...
                            'ro', 'MarkerSize', 4);
            xlabel(sprintf('x%d', xlist(1)), 'FontSize', 6);
            ylabel(sprintf('x%d', xlist(2)), 'FontSize', 6);
            zlabel(sprintf('x%d', xlist(3)), 'FontSize', 6);
            if(p.Results.is_text)
                strValues = strtrim(cellstr(...
                        num2str([x(:,1) x(:,2) x(:,3)],...
                                '(%.4f,%.4f,%.4f)')));
                text(x(:,1), x(:,2), x(:,3), ...
                                strValues, 'VerticalAlignment', 'bottom');
            end
        end
        title(xstr, 'FontSize', 8);
        box on;
        drawnow;
    end

author: chudur-budur


text:
function kern = invcmpndKernExpandParam(kern, params)

% INVCMPNDKERNEXPANDPARAM Create kernel structure from INVCMPND kernel's parameters.
% FORMAT
% DESC returns a inv-compound kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : invcmpndKernParamInit, invcmpndKernExtractParam, cmpndkernExpandParam
%
% COPYRIGHT : Andreas C. Damianou, 2012

% KERN

kern = cmpndKernExpandParam(kern, params);


text:
function [ S, f, Serr ]= mtspectrumc_unequal_length_trials( data, movingwin, params, sMarkers )

% This routine computes the multi-taper spectrum for a given set of unequal length segments. It is
% based on modifications to the Chronux routines. The segments are continuously structured in the
% data matrix, with the segment boundaries given by markers. Below,
% movingwin is used in a non-overlaping way to partition each segment into
% various windows. Th spectrum is evaluated for each window, and then the
% window spectrum estimates averaged. Further averaging is conducted by
% repeating the process for each segment.
%
% Inputs:
%
%   data = data( samples, channels )- here segments must be stacked
%   as explained in the email
%   movingwin = [window winstep] i.e length of moving
%              window and step size. Note that units here have
%              to be consistent with units of Fs. If Fs=1 (ie normalized)
%              then [window winstep]should be in samples, or else if Fs is
%              unnormalized then they should be in time (secs).
%   sMarkers = N x 2 array of segment start & stop marks. sMarkers(n, 1) = start
%           sample index; sMarkers(n,2) = stop sample index for the nth segment
%   params = see Chronux help on mtspecgramc
%
% Output:
%
%       S       frequency x channels
%       f       frequencies x 1
%       Serr    (error bars) only for err(1)>=1
%
%

iwAvg = 1; % 0=no weighted average, 1=weighted average
debug = 1; % will display intermediate calcs.

if nargin < 2; error('avgSpectrum:: Need data and window parameters'); end;
if nargin < 3; params=[]; end;
if isempty( sMarkers ), error( 'avgSpectrum:: Need Markers...' ); end
[ tapers, pad, Fs, fpass, err, trialave, params ] = getparams( params );
if nargout > 2 && err(1)==0;
%   Cannot compute error bars with err(1)=0. change params and run again.
    error('avgSpectrum:: When Serr is desired, err(1) has to be non-zero.');
end;

% Set moving window parameters to no-overlapping
if abs(movingwin(2) - movingwin(1)) >= 1e-6, disp( 'avgSpectrum:: Warming: Window parameters for averaging should be non-overlapping. Set movingwin(2) = movingwin(1).' ); end

wLength = round( Fs * movingwin(1) ); % number of samples in window
wStep = round( movingwin(2) * Fs ); % number of samples to step through

% Check whether window lengths satify segment length > NW/2
if ( wLength < 2*tapers(1) ), error( 'avgSpectrum:: movingwin(1) > 2*tapers(1)' ); end

% Left align segment markers for easier coding
sM = ones( size( sMarkers, 1 ), 2 );
sM( :, 2 ) = sMarkers( :, 2 ) - sMarkers( :, 1 ) + 1;

% min-max segments
Nmax = max( sM(:,2) ); Nmin = min( sM(:,2) );
if ( Nmin < 2*tapers(1) ), error( 'avgSpectrum:: Smallest segment length > 2*tapers(1). Change taper settings' ); end

% max time-sample length will be the window length.
nfft = 2^( nextpow2( wLength ) + pad );
[ f, findx ] = getfgrid( Fs, nfft, fpass);

% Precompute all the tapers
sTapers = tapers;
sTapers = dpsschk( sTapers, wLength, Fs ); % compute tapers for window length

nChannels = size( data, 2 );
nSegments = size( sMarkers, 1 );

if debug
    disp( ['Window Length = ' num2str(wLength)] );
    disp( ['Window Step = ' num2str(wStep)] );
    disp( ' ' );
end

s = zeros( length(f), nChannels );
serr = zeros( 2, length(f), nChannels );
S = zeros( length(f), nChannels );
Serr = zeros( 2, length(f), nChannels );
nWins = 0;
for sg = 1 : nSegments
    % Window lengths & steps fixed above
    % For the given segment, compute the positions & number of windows
    N = sM(sg,2);
    wStartPos = 1 : wStep : ( N - wLength + 1 );
    nWindows = length( wStartPos );
    if nWindows
        nWins = nWins + nWindows; % for averaging purposes

        w=zeros(nWindows,2);
        for n = 1 : nWindows
            w(n,:) = [ wStartPos(n), (wStartPos(n) + wLength - 1) ]; % nWindows x 2. just like segment end points
        end

        % Shift window limits back to original sample-stamps
        w(:, 1) = w(:,1) + (sMarkers( sg, 1 ) - 1);
        w(:, 2) = w(:,2) + (sMarkers( sg, 1 ) - 1);

        if debug
            disp( ['Segment Start/Stop = ' num2str( w(1,1) ) ' ' num2str( w(end,2) ) ] );
            disp( ['Min / Max Window Positions = ' num2str( min(w(:,1)) ) ' ' num2str( max(w(:,1)) ) ] );
            disp( ['Total Number of Windows = ' num2str(nWindows) ]);
            disp( ' ' );
        end

        % Pile up window segments similar to segment pileup
        wData = zeros( wLength, nChannels, nWindows ); %initialize to avoid fragmentation
        for n = 1:nWindows
            %wData( :, :, n ) = detrend( data( w(n,1):w(n,2), : ), 'constant' );
            wData( :, :, n ) = detrend( data( w(n,1):w(n,2), : ) );
        end

        % J1 = frequency x taper x nWindows
        % J2 = frequency x taper x nWindows x nChannels
        J2 = zeros( length(f), tapers(2), nWindows, nChannels ); J2 = complex( J2, J2 );
        for c = 1 : nChannels
            J1 = mtfftc( squeeze(wData( :, c, : )), sTapers, nfft, Fs ); % FFT for the tapered data
            J2( :, :, :, c ) = J1(findx,:,:);
        end
        % J2 = frequency x taper x nWindows x nChannels
        % Inner mean = Average over tapers => frequency x nWindows x nChannels
        % Outer mean = Average over windows => frequency x nChannels
        dim1 = [length(f), nWindows, nChannels];
        dim2 = [length(f), nChannels];
        % s = frequency x nChannels
        s = reshape( squeeze( mean( reshape( squeeze( mean( conj(J2).*J2, 2 ) ), dim1), 2 ) ), dim2 );

        % Now treat the various "windowed data" as "trials"
        % serr = 2 x frequency x channels. Output from specerr = 2 x frequency x 1
        for c = 1 : nChannels
            serr( :, :, c ) = specerr( squeeze( s(:, c ) ), squeeze( J2(:,:,:, c ) ), err, 1 );
        end

        if iwAvg
            % Segment Weighted error estimates.
            S = S + nWindows*s;
            Serr = Serr + nWindows*serr;
        else
            S = S + s;
            Serr = Serr + serr;
        end

    else
        if debug, disp(['avgSpectrum:: Zero windows for segment: ' num2str(sg) ]); end
    end
end

% Segment Weighted error estimates.
% Only over those that had non-zero windows
if nWins && iwAvg
    S=S/nWins; Serr=Serr/nWins;
end
if ~nWins
    if debug, disp(['avgCoherence:: No segment long enough with movingwin parameters found. Reduce movingwin.' ]); end
end

author: VisLab

repo: EEG-Clean-Tools


text:
classdef (Abstract) extension < handle

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

    methods
        function task_class = task_class(obj, task_class, mpopt)
        end

        function dmc_class = dm_converter_class(obj, dmc_class, fmt, mpopt)
        end

        function dm_class = data_model_class(obj, dm_class, task_tag, mpopt)
        end

        function nm_class = network_model_class(obj, nm_class, task_tag, mpopt)
        end

        function mm_class = math_model_class(obj, mm_class, task_tag, mpopt)
        end

        function dmc_elements = dmc_element_classes(obj, dmc_class, fmt, mpopt)
            dmc_elements = {};      %% no modifications
        end

        function dm_elements = dm_element_classes(obj, dm_class, task_tag, mpopt)
            dm_elements = {};       %% no modifications
        end

        function nm_elements = nm_element_classes(obj, nm_class, task_tag, mpopt)
            nm_elements = {};       %% no modifications
        end

        function mm_elements = mm_element_classes(obj, mm_class, task_tag, mpopt)
            mm_elements = {};       %% no modifications
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function cm = getColorMap(num_colors)
%GETCOLORMAP    Return default k-Wave color map.
%
% DESCRIPTION:
%       getColorMap returns the default color map used for display and
%       visualisation across the k-Wave Toolbox. Zero values are displayed
%       as white, positive values are displayed as yellow through red to
%       black, and negative values are displayed as light to dark
%       blue-greys. If no value for num_colors is provided, cm will have
%       256 colors.
%
% USAGE:
%       cm = getColorMap()
%       cm = getColorMap(num_colors)
%
% OPTIONAL INPUTS:
%       num_colors  - number of colors in the color map (default = 256)
%
% OUTPUTS:
%       cm          - three column color map matrix which can be applied
%                     using colormap
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 3rd July 2009
%       last update - 17th July 2009
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also colormap

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% set literals
if nargin == 0
    neg_pad = 48;
    num_colors = 256;
else
    neg_pad = round(48*num_colors/256);
end

% define colour spectrums
neg = bone(num_colors/2 + neg_pad);
neg = neg(1 + neg_pad:end, :);
pos = flipud(hot(num_colors/2));

% create custom colour map
cm = [neg; pos];

repo: Aero-matlab


text:
function rowMatrix = computeRowsOfToeplitzHankelMatrix(rowNumber,...
        nColumns, crossCorrelationVectors, hankelMatrixIsAdded, dcIsIncluded)
%     if rowNumber == 1
%         toeplitzRows = crossCorrelationVectors(1:nColumns,:);
%     else
%         toeplitzRows = ...
%             [flip(crossCorrelationVectors(2:rowNumber,:),1);...
%             crossCorrelationVectors(1:nColumns-rowNumber+1,:)];
%     end
% modified by Liming Shi for fast computation in MATLAB
    toeplitzRows=crossCorrelationVectors(rowNumber:-1:1,:);

    if dcIsIncluded && hankelMatrixIsAdded
        hankelOffset = 1;
    else
        hankelOffset = 3;
    end

    hankelRows = crossCorrelationVectors((0:nColumns-1)+...
                                         hankelOffset+rowNumber-1,:);

    if hankelMatrixIsAdded
        rowMatrix = toeplitzRows + hankelRows;
    else
        rowMatrix = toeplitzRows - hankelRows;
    end

end

author: LimingShi


text:
% Set the path to your trackers
path_to_your_trackers = 'E:/Database/vot-toolkit-master_HCF/tracker';

if isempty(path_to_your_trackers)
    error('Set the path to your trackers!');
end

tracker_label = 'DeepSTRCF';
tracker_command = generate_matlab_command('vot_wrapper(''DeepSTRCF'', ''DeepSTRCF_VOT_setting'')', {[path_to_your_trackers '/STRCF/']});
tracker_interpreter = 'matlab';
tracker_trax = false;

repo: STRCF


text:
% Author: Jai Juneja
% Date: 12/02/2013
%
% Function transforms a point p_r from robot's local frame to the global
% frame.
%
% Inputs:
%   r = [x y alpha]'	:   Robot frame
%   p_r                 :	Point in robot's frame (2*n matrix for n points)
%
% Outputs:
%   p_g                 :   Point in global frame (2*n matrix)
%   Optional:
%   PG_r                :   Jacobian of p_g wrt. robot frame
%   PG_pr               :   Jacobian of p_g wrt. p_r
function [p_g, PG_r, PG_pr] = transToGlobal(r, p_r)

    r_pos = r(1:2);     % Robot's position
    r_ang = r(3);       % Robot's orientation

    % Transformation matrix:
    R = [cos(r_ang), -sin(r_ang) ; ...
        sin(r_ang), cos(r_ang)];

    % Hence point p_r tranformed to global co-ordinate frame is:
    p_g = R * p_r;
    p_g = [p_g(1, :) + r_pos(1); p_g(2, :) + r_pos(2)];

    if nargout > 1      % Compute Jacobians
        pr_x = p_r(1); pr_y = p_r(2);

        PG_r = [...
            1,  0,  - pr_y*cos(r_ang) - pr_x*sin(r_ang); ...
            0,  1,  pr_x*cos(r_ang) - pr_y*sin(r_ang)];

        PG_pr = R;
    end
end

author: jaijuneja

repo: ekf-slam-matlab