text:
classdef ParadigmSpecCSP < ParadigmDataflowSimplified
    % Advanced paradigm for oscillatory processes via the Spectrally weighted CSP algorithm.
    %
    % The Spec-CSP paradigm [1] is a more advanced variant of CSP, developed for the Berlin
    % Brain-Computer Interface (BBCI); the primary focus was motor imagery BCI, but the algorithm was
    % designed from the outset to be applicable for a wider range of applications. The implementation
    % closely follows the TR [2].
    %
    % The paradigm is applicable to the majority of oscillatory processes, and is the most advanced
    % spatio-spectrally adaptive method that is currently provided in the toolbox. Whenever the exact
    % frequency and location of some (conjectured) oscillatory process is not known exactly, Spec-CSP
    % can be used, and typically gives better results than CSP with an appropriately unrestricted (e.g.,
    % broad-band) spectral filter. Several other methods exist to adapt the spectrum to a process of
    % interest, among others Common Spatio-Spectral Patterns [3], Common Sparse Spectral Spatial Pattern
    % [4], r^2-based heuristics [5], automated parameter search, and manual selection based on visual
    % inspection. Several of these methods have been shown to give approx. comparable results [2]. An
    % alternative and competitive method, especially when there complex interactions between frequency
    % bands and time periods are to be modeled is the Dual-Augmented Lagrange paradigm
    % (para_dal/para_dal_hf).
    %
    % The method iteratively optimizes spatial and spectral filters in alternation and extracts
    % log-variance features from the resulting (filtered) signal. These features are subsequently
    % processed by a (typically simple) machine learning component, by default LDA. Learning is
    % therefore significantly slower than CSP. An option is to impose custom prior knowledge on the
    % relevant data spectrum, for example by placing a focus on the alpha rhythm, without ruling out
    % other frequencies. Note that there are parameters which constrain the spectrum: one is the
    % frequency prior and the other is the spectral filter that is applied before running the alorithm;
    % both need to be adapted when the considered spectrum shall be extended (e.g. to high-gamma
    % oscillations). Other parameters which are frequently adapted are the time window of interest and
    % the learner component (e.g., logistic regression is a good alternative choice).
    %
    % Some application areas include detection of major brain rhythm modulations (e.g. theta, alpha,
    % beta, gamma), for example related to relaxation/stress, aspects of workload, emotion,
    % sensori-motor imagery, and in general cortical idle oscillations in various modalities.
    %
    % Example: Consider the goal of predicting the emotion felt by a person at a given time. A possible
    % calibration data set for this task would contain a sequence of blocks in each of which the subject
    % is one out of several possible emotions, indicated by events 'e1','e2','e3','e4' covering these
    % blocks at regular rate. The data might for example be induced via guided imagery [6].
    %
    %   calib = io_loadset('data sets/bruce/emotions.eeg')
    %   myapproach = {'SpecCSP' 'SignalProcessing',{'EpochExtraction',[-2.5 2.5]}};
    %   [loss,model,stats] = bci_train('Data',calib,'Approach',myapproach, 'TargetMarkers',{'e1','e2','e3','e4'});
    %
    %
    % References:
    %  [1] Tomioka, R., Dornhege, G., Aihara, K., and Mueller, K.-R. "An iterative algorithm for spatio-temporal filter optimization."
    %      In Proceedings of the 3rd International Brain-Computer Interface Workshop and Training Course 2006.
    %  [2] Ryota Tomioka, Guido Dornhege, Guido Nolte, Benjamin Blankertz, Kazuyuki Aihara, and Klaus-Robert Mueller
    %      "Spectrally Weighted Common Spatial Pattern Algorithm for Single Trial EEG Classification",
    %      Mathematical Engineering Technical Reports (METR-2006-40), July 2006.
    %  [3] Steven Lemm, Benjamin Blankertz, Gabriel Curio, and Klaus-Robert M�ller.
    %      "Spatio-spectral filters for improving classification of single trial EEG."
    %      IEEE Trans Biomed Eng, 52(9):1541-1548, 2005.
    %  [4] G. Dornhege, B. Blankertz, M. Krauledat, F. Losch, G. Curio, and K.-R. M�ller,
    %      "Combined optimization of spatial and temporal filters for improving brain-computer interfacing,"
    %      IEEE Transactions on Biomedical Engineering, vol. 53, no. 11, pp. 2274?2281, 2006.
    %  [5] Benjamin Blankertz, Ryota Tomioka, Steven Lemm, Motoaki Kawanabe, and Klaus-Robert Mueller.
    %      "Optimizing spatial filters for robust EEG single-trial analysis."
    %      IEEE Signal Process Mag, 25(1):41-56, January 2008
    %  [6] Onton J & Makeig S. "Broadband high-frequency EEG dynamics during emotion imagination."
    %      Frontiers in Human Neuroscience, 2009.
    %
    % Name:
    %   Spectrally Weighted CSP
    %
    %                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                           2010-04-29

    methods

        function defaults = preprocessing_defaults(self)
            defaults = {'FIRFilter',{'Frequencies',[6 7 33 34],'Type','minimum-phase'}, 'EpochExtraction',[0.5 3.5], 'Resampling',100};
        end

        function model = feature_adapt(self,varargin)
            args = arg_define(varargin, ...
                arg_norep('signal'), ...
                arg({'patterns','PatternPairs'},3,uint32([1 1 64 1000]),'Number of CSP patterns (times two).','cat','Feature Extraction'), ...
                arg({'pp','ParameterP'},0,[-1 1],'Regularization parameter p''. Can be searched over -1:0.5:1.','cat','Feature Extraction','guru',true), ...
                arg({'qp','ParameterQ'},1,[0 4],'Regularization parameter q''. Can be searched over 0:0.5:4.','cat','Feature Extraction','guru',true), ...
                arg({'prior','SpectralPrior'},'@(f) f>=7 & f<=30',[],'Prior frequency weighting function.','cat','Feature Extraction', 'type','expression'), ...
                arg({'steps','MaxIterations'},3,uint32([1 3 10 50]),'Number of iterations. A step is spatial optimization, followed by spectral optimization.','cat','Feature Extraction'));

            [signal,n_of,pp,qp,prior,steps] = deal(args.signal,args.patterns,args.pp,args.qp,args.prior,args.steps);
            if signal.nbchan == 1
                error('Spec-CSP does intrinsically not support single-channel data (it is a spatial filter).'); end
            if signal.nbchan < args.patterns
                error('Spec-CSP prefers to work on at least as many channels as you request output patterns. Please reduce the number of pattern pairs.'); end

            % read a few parameters from the options (and re-parameterize the hyper-parameters p' and q' into p and q)
            p = pp+qp;
            q = qp;
            if isnumeric(prior) && length(prior) == 2
                prior = @(f) f >= prior(1) & f <= prior(2); end
            % number of C=Channels, S=Samples and T=Trials #ok<NASGU>
            [C,S,dum] = size(signal.data); %#ok<NASGU>
            % build a frequency table (one per DFT bin)
            freqs = (0:S-1)*signal.srate/S;
            % evaluate the prior I
            I = prior(freqs);
            % and find table indices that are supported by the prior
            bands = find(I);

            % preprocessing
            for c=1:2
                % compute the per-class epoched data X and its Fourier transform (along time), Xfft
                X{c} = exp_eval_optimized(set_picktrials(signal,'rank',c));
                [C,S,T] = size(X{c}.data);
                Xfft{c} = fft(X{c}.data,[],2);
                % the full spectrum F of covariance matrices per every DFT bin and trial of the data
                F{c} = single(zeros(C,C,max(bands),T));
                for k=bands
                    for t=1:T
                        F{c}(:,:,k,t) = 2*real(Xfft{c}(:,k,t)*Xfft{c}(:,k,t)'); end
                end
                % compute the cross-spectrum V as an average over trials
                V{c} = mean(F{c},4);
            end

            % 1. initialize the filter set alpha and the number of filters J
            J = 1; alpha{J}(bands) = 1;
            % 2. for each step
            for step=1:steps
                % 3. for each set of spectral coefficients alpha{j} (j=1,...,J)
                for j=1:J
                    % 4. calculate sensor covariance matrices for each class from alpha{j}
                    for c = 1:2
                        Sigma{c} = zeros(C);
                        for b=bands
                            Sigma{c} = Sigma{c} + alpha{j}(b)*V{c}(:,:,b); end
                    end
                    % 5. solve the generalized eigenvalue problem Eq. (2)
                    [VV,DD] = eig(Sigma{1},Sigma{1}+Sigma{2});
                    % and retain n_of top eigenvectors at both ends of the eigenvalue spectrum...
                    W{j} = {VV(:,1:n_of), VV(:,end-n_of+1:end)};
                    iVV = inv(VV)'; P{j} = {iVV(:,1:n_of), iVV(:,end-n_of+1:end)};
                    % as well as the top eigenvalue for each class
                    lambda(j,:) = [DD(1), DD(end)];
                end
                % 7. set W{c} from all W{j}{c} such that lambda(j,c) is minimal/maximal over j
                W = {W{argmin(lambda(:,1))}{1}, W{argmax(lambda(:,2))}{2}};
                P = {P{argmin(lambda(:,1))}{1}, P{argmax(lambda(:,2))}{2}};
                % 8. for each projection w in the concatenated [W{1},W{2}]...
                Wcat = [W{1} W{2}]; J = 2*n_of;
                Pcat = [P{1} P{2}];
                for j=1:J
                    w = Wcat(:,j);
                    % 9. calcualate (across trials within each class) mean and variance of the w-projected cross-spectrum components
                    for c=1:2
                        % part of Eq. (3)
                        s{c} = zeros(size(F{c},4),max(bands));
                        for k=bands
                            for t = 1:size(s{c},1)
                                s{c}(t,k) = w'*F{c}(:,:,k,t)*w; end
                        end
                        mu_s{c} = mean(s{c},1);
                        var_s{c} = var(s{c},0,1);
                    end
                    % 10. update alpha{j} according to Eqs. (4) and (5)
                    for c=1:2
                        for k=bands
                            % Eq. (4)
                            alpha_opt{c}(k) = max(0, (mu_s{c}(k)-mu_s{3-c}(k)) / (var_s{1}(k) + var_s{2}(k)) );
                            % Eq. (5), with prior from Eq. (6)
                            alpha_tmp{c}(k) = alpha_opt{c}(k).^q * (I(k) * (mu_s{1}(k) + mu_s{2}(k))/2).^p;
                        end
                    end
                    % ... as the maximum for both classes
                    alpha{j} = max(alpha_tmp{1},alpha_tmp{2});
                    % and normalize alpha{j} so that it sums to unity
                    alpha{j} = alpha{j} / sum(alpha{j});
                end
            end
            alpha = [vertcat(alpha{:})'; zeros(S-length(alpha{1}),length(alpha))];
            model = struct('W',{Wcat},'P',{Pcat},'alpha',{alpha},'freqs',{freqs},'bands',{bands},'chanlocs',{signal.chanlocs});
        end

        function features = feature_extract(self,signal,featuremodel)
            features = zeros(size(signal.data,3),size(featuremodel.W,2));
            for t=1:size(signal.data,3)
                features(t,:) = log(var(2*real(ifft(featuremodel.alpha.*fft(signal.data(:,:,t)'*featuremodel.W))))); end
        end

        function visualize_model(self,varargin) %#ok<*INUSD>
            args = arg_define([0 3],varargin, ...
                arg_norep({'myparent','Parent'},[],[],'Parent figure.'), ...
                arg_norep({'featuremodel','FeatureModel'},[],[],'Feature model. This is the part of the model that describes the feature extraction.'), ...
                arg_norep({'predictivemodel','PredictiveModel'},[],[],'Predictive model. This is the part of the model that describes the predictive mapping.'), ...
                arg({'patterns','PlotPatterns'},true,[],'Plot patterns instead of filters. Whether to plot spatial patterns (forward projections) rather than spatial filters.'), ...
                arg({'paper','PaperFigure'},false,[],'Use paper-style font sizes. Whether to generate a plot with font sizes etc. adjusted for paper.'), ...
                arg_nogui({'nosedir_override','NoseDirectionOverride'},'',{'','+X','+Y','-X','-Y'},'Override nose direction.'));
            arg_toworkspace(args);

            % no parent: create new figure
            if isempty(myparent)
                myparent = figure('Name','Common Spatial Patterns'); end
            % determine nose direction for EEGLAB graphics
            try
                nosedir = args.fmodel.signal.info.chaninfo.nosedir;
            catch
                disp_once('Nose direction for plotting not store in model; assuming +X');
                nosedir = '+X';
            end
            if ~isempty(nosedir_override)
                nosedir = nosedir_override; end
            % number of pairs, and index of pattern per subplot
            np = size(featuremodel.W,2)/2; idxp = [1:np np+(2*np:-1:np+1)]; idxf = [np+(1:np) 2*np+(2*np:-1:np+1)];
            % for each CSP pattern...
            for p=1:np*2
                subplot(4,np,idxp(p),'Parent',myparent);
                if args.patterns
                    plotdata = featuremodel.P(:,p);
                else
                    plotdata = featuremodel.W(:,p);
                end
                topoplot(plotdata,featuremodel.chanlocs,'nosedir',nosedir);
                subplot(4,np,idxf(p),'Parent',myparent);
                alpha = featuremodel.alpha(:,p);
                range = 1:max(find(alpha)); %#ok<MXFND>
                pl=plot(featuremodel.freqs(range),featuremodel.alpha(range,p));
                xlim([min(featuremodel.freqs(range)) max(featuremodel.freqs(range))]);
                l1 = xlabel('Frequency in Hz');
                l2 = ylabel('Weight');
                t=title(['Spec-CSP Pattern ' num2str(p)]);
                if args.paper
                    set([gca,t,l1,l2],'FontUnits','normalized');
                    set([gca,t,l1,l2],'FontSize',0.2);
                    set(pl,'LineWidth',2);
                end
            end
            try set(gcf,'Color',[1 1 1]); end
        end

        function layout = dialog_layout_defaults(self)
            layout = {'SignalProcessing.Resampling.SamplingRate', 'SignalProcessing.FIRFilter.Frequencies', ...
                'SignalProcessing.FIRFilter.Type', 'SignalProcessing.EpochExtraction', '', ...
                'Prediction.FeatureExtraction', '', ...
                'Prediction.MachineLearning.Learner'};
        end

        function tf = needs_voting(self)
            tf = true;
        end

    end
end

repo: REST


text:
function ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid, visBool)
% matRad function to import dicom ct data
%
% call
%   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool)
%   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid)
%   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, visBool)
%   ct = matRad_importDicomCt(ctList, resolution, dicomMetaBool, grid, visBool)
%
% input
%   ctList:         list of dicom ct files
%   resolution:   	resolution of the imported ct cube, i.e. this function
%                   will interpolate to a different resolution if desired
%   dicomMetaBool:  store complete dicom information if true
%   grid:           optional: a priori grid specified for interpolation
%   visBool:        optional: turn on/off visualization
%
% output
%   ct:             matRad ct struct. Note that this 3D matlab array
%                   contains water euqivalent electron denisities.
%                   Hounsfield units are converted using a standard lookup
%                   table in matRad_calcWaterEqD
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matRad_cfg = MatRad_Config.instance();

matRad_cfg.dispInfo('\nimporting ct-cube...');

%% processing input variables
if ~exist('visBool','var')
  visBool = 0;
end

% creation of ctInfo list
numOfSlices = size(ctList,1);
matRad_cfg.dispInfo('\ncreating info...')

sliceThicknessStandard = true;
for i = 1:numOfSlices

    if verLessThan('matlab','9')
        tmpDicomInfo = dicominfo(ctList{i,1});
    else
        tmpDicomInfo = dicominfo(ctList{i,1},'UseDictionaryVR',true);
    end

    % remember relevant dicom info - do not record everything as some tags
    % might not been defined for individual files
    ctInfo(i).PixelSpacing            = tmpDicomInfo.PixelSpacing;
    ctInfo(i).ImagePositionPatient    = tmpDicomInfo.ImagePositionPatient;
    ctInfo(i).SliceThickness          = tmpDicomInfo.SliceThickness;
    ctInfo(i).ImageOrientationPatient = tmpDicomInfo.ImageOrientationPatient;
    ctInfo(i).PatientPosition         = tmpDicomInfo.PatientPosition;
    ctInfo(i).Rows                    = tmpDicomInfo.Rows;
    ctInfo(i).Columns                 = tmpDicomInfo.Columns;
    ctInfo(i).Width                   = tmpDicomInfo.Width;
    ctInfo(i).Height                  = tmpDicomInfo.Height;
    ctInfo(i).RescaleSlope            = tmpDicomInfo.RescaleSlope;
    ctInfo(i).RescaleIntercept        = tmpDicomInfo.RescaleIntercept;

    %Problem due to some CT files using non-standard SpacingBetweenSlices

    if isempty(ctInfo(i).SliceThickness)
        %Print warning ocne
        if sliceThicknessStandard
            matRad_cfg.dispWarning('Non-standard use of SliceThickness Attribute (empty), trying to overwrite with SpacingBetweenSlices');
            sliceThicknessStandard = false;
        end
        ctInfo(i).SliceThickness = tmpDicomInfo.SpacingBetweenSlices;
    end

    if i == 1
        completeDicom = tmpDicomInfo;
    end

    matRad_progress(i,numOfSlices);
end

% adjusting sequence of slices (filenames may not be ordered propperly....
% e.g. CT1.dcm, CT10.dcm, CT100zCoordList = [ctInfo.ImagePositionPatient(1,3)]';.dcm, CT101.dcm,...
CoordList = [ctInfo.ImagePositionPatient]';
[~, indexing] = sort(CoordList(:,3)); % get sortation from z-coordinates

ctList = ctList(indexing);
ctInfo = ctInfo(indexing);

%% check data set for consistency
if size(unique([ctInfo.PixelSpacing]','rows'),1) > 1
    matRad_cfg.dispError('Different pixel size in different CT slices');
end

coordsOfFirstPixel = [ctInfo.ImagePositionPatient];
if numel(unique(coordsOfFirstPixel(1,:))) > 1 || numel(unique(coordsOfFirstPixel(2,:))) > 1
    matRad_cfg.dispError('Ct slices are not aligned');
end
if sum(diff(coordsOfFirstPixel(3,:))<=0) > 0
    matRad_cfg.dispError('Ct slices not monotonically increasing');
end
if numel(unique([ctInfo.Rows])) > 1 || numel(unique([ctInfo.Columns])) > 1
    matRad_cfg.dispError('Ct slice sizes inconsistent');
end

%% checking the patient position
% As of now, the matRad treatment planning system is only valid for
% patients in a supine position. Other orientations (e.g. prone, decubitus
% left/right) are not supported.
% Defined Terms:
% HFP     Head First-Prone                  (not supported)
% HFS     Head First-Supine                 (supported)
% HFDR    Head First-Decubitus Right        (not supported)
% HFDL    Head First-Decubitus Left         (not supported)
% FFDR    Feet First-Decubitus Right        (not supported)
% FFDL    Feet First-Decubitus Left         (not supported)
% FFP     Feet First-Prone                  (not supported)
% FFS     Feet First-Supine                 (supported)

if isempty(regexp(ctInfo(1).PatientPosition,{'S','P'}, 'once'))
    matRad_cfg.dispError(['This Patient Position is not supported by matRad.'...
        ' As of now only ''HFS'' (Head First-Supine), ''FFS'''...
        ' (Feet First-Supine), '...
        '''HFP'' (Head First-Prone), and ''FFP'''...
        ' (Feet First-Prone) can be processed.'])
end

%% creation of ct-cube
matRad_cfg.dispInfo('reading slices...')
origCt = zeros(ctInfo(1).Height, ctInfo(1).Width, numOfSlices);
for i = 1:numOfSlices
    currentFilename = ctList{i};
    [currentImage, map] = dicomread(currentFilename);
    origCt(:,:,i) = currentImage(:,:); % creation of the ct cube

    % draw current ct-slice
    if visBool
        if ~isempty(map)
            image(ind2rgb(uint8(63*currentImage/max(currentImage(:))),map));
            xlabel('x [voxelnumber]')
            ylabel('y [voxelnumber]')
            title(['Slice # ' int2str(i) ' of ' int2str(numOfSlices)])
        else
            image(ind2rgb(uint8(63*currentImage/max(currentImage(:))),bone));
            xlabel('x [voxelnumber]')
            ylabel('y [voxelnumber]')
            title(['Slice # ' int2str(i) ' of ' int2str(numOfSlices)])
        end
        axis equal tight;
        pause(0.1);
    end
    matRad_progress(i,numOfSlices);
end

%% correction if not lps-coordinate-system
% when using the physical coordinates (ctInfo.ImagePositionPatient) to
% arrange the  slices in z-direction, there is no more need for mirroring
% in the z-direction
matRad_cfg.dispInfo('\nz-coordinates taken from ImagePositionPatient\n')

% The x- & y-direction in lps-coordinates are specified in:
% ImageOrientationPatient
xDir = ctInfo(1).ImageOrientationPatient(1:3); % lps: [1;0;0]
yDir = ctInfo(1).ImageOrientationPatient(4:6); % lps: [0;1;0]
nonStandardDirection = false;

% correct x- & y-direction
%
% if xDir(1) == 1 && xDir(2) == 0 && xDir(3) == 0
%     matRad_cfg.dispInfo('x-direction OK\n')
% elseif xDir(1) == -1 && xDir(2) == 0 && xDir(3) == 0
%     matRad_cfg.dispInfo('\nMirroring x-direction...')
%     origCt = flip(origCt,1);
%     matRad_cfg.dispInfo('finished!\n')
% else
%     nonStandardDirection = true;
% end
%
% if yDir(1) == 0 && yDir(2) == 1 && yDir(3) == 0
%     matRad_cfg.dispInfo('y-direction OK\n')
% elseif yDir(1) == 0 && yDir(2) == -1 && yDir(3) == 0
%     matRad_cfg.dispInfo('\nMirroring y-direction...')
%     origCt = flip(origCt,2);
%     matRad_cfg.dispInfo('finished!\n')
% else
%     nonStandardDirection = true;
% end

if nonStandardDirection
    matRad_cfg.dispInfo(['Non-standard patient orientation.\n'...
        'CT might not fit to contoured structures\n'])
end

%% interpolate cube
matRad_cfg.dispInfo('\nInterpolating CT cube...');
if exist('grid','var')
    ct = matRad_interpDicomCtCube(origCt, ctInfo, resolution, grid);
else
    ct = matRad_interpDicomCtCube(origCt, ctInfo, resolution);
end
matRad_cfg.dispInfo('finished!\n');

%% remember some parameters of original dicom
ct.dicomInfo.PixelSpacing            = ctInfo(1).PixelSpacing;
                                       tmp = [ctInfo.ImagePositionPatient];
ct.dicomInfo.SlicePositions          = tmp(3,:);
ct.dicomInfo.SliceThickness          = [ctInfo.SliceThickness];
ct.dicomInfo.ImagePositionPatient    = ctInfo(1).ImagePositionPatient;
ct.dicomInfo.ImageOrientationPatient = ctInfo(1).ImageOrientationPatient;
ct.dicomInfo.PatientPosition         = ctInfo(1).PatientPosition;
ct.dicomInfo.Width                   = ctInfo(1).Width;
ct.dicomInfo.Height                  = ctInfo(1).Height;
ct.dicomInfo.RescaleSlope            = ctInfo(1).RescaleSlope;
ct.dicomInfo.RescaleIntercept        = ctInfo(1).RescaleIntercept;
if isfield(completeDicom, 'Manufacturer')
ct.dicomInfo.Manufacturer            = completeDicom.Manufacturer;
end
if isfield(completeDicom, 'ManufacturerModelName')
ct.dicomInfo.ManufacturerModelName   = completeDicom.ManufacturerModelName;
end
if isfield(completeDicom, 'ConvolutionKernel')
ct.dicomInfo.ConvolutionKernel       = completeDicom.ConvolutionKernel;
end

% store patientName only if user wants to
if isfield(completeDicom,'PatientName') && dicomMetaBool == true
    ct.dicomInfo.PatientName         = completeDicom.PatientName;
end
if dicomMetaBool == true
    ct.dicomMeta                     = completeDicom;
end

ct.timeStamp = datestr(clock);

% convert to Hounsfield units
matRad_cfg.dispInfo('\nconversion of ct-Cube to Hounsfield units...');
ct = matRad_calcHU(ct);
matRad_cfg.dispInfo('finished!\n');

end

repo: matRad


text:
function M = shapefitfactory(VJt)
% Linear manifold structure for optimization over the ShapeFit search space
%
% function M = shapefitfactory(VJt)
%
% Input: VJt is a matrix of size dxn, such that VJt * ones(n, 1) = 0.
%
% Returns M, a structure describing the Euclidean space of d-by-n matrices
% equipped with the standard Frobenius distance and associated trace inner
% product, as a manifold for Manopt. Matrices on M, denoted by T, have size
% dxn and obey T*ones(n, 1) = 0 (centered columns) and <VJt, T> = 1, where
% <A, B> = Trace(A' * B).
%
% See this paper: http://arxiv.org/abs/1506.01437
% ShapeFit: Exact location recovery from corrupted pairwise directions, 2015
% Paul Hand, Choongbum Lee, Vladislav Voroninski
%
% See also: shapefit_smoothed

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 18, 2015.
% Contributors:
% Change log:
%
%   Jan. 25, 2017 (NB):
%       M.tangent = M.proj now, instead of being identity. This is notably
%       necessary so that checkgradient will pick up on gradients that do
%       not lie in the appropriate tangent space.
%
%   Jan. 4, 2021 (NB):
%       Changes for compatibility with Octave 6.1.0.

    [d, n] = size(VJt);

    M.name = @() sprintf('ShapeFit space of size %d x %d', d, n);

    M.dim = @() d*n - d - 1;

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d, 'fro');

    M.dist = @(x, y) norm(x-y, 'fro');

    M.typicaldist = @() sqrt(d*n);

    VJt_normed = VJt / norm(VJt, 'fro');
    M.proj = @(T, U) projection(U, VJt_normed);
    function PU = projection(U, VJt_normed)
        % Center the columns
        PU = bsxfun(@minus, U, mean(U, 2));
        % Remove component along VJt
        % Note: these two actions can be executed separately, without
        % interference, owing to VJt having centered columns itself.
        PU = PU - (VJt_normed(:)'*U(:))*VJt_normed;
    end

    M.egrad2rgrad = M.proj;

    M.ehess2rhess = @(x, eg, eh, d) projection(eh, VJt_normed);

    M.tangent = M.proj;

    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end

    M.retr = M.exp;

    M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.randvec = @(x) randvec(VJt, VJt_normed);
    function u = randvec(VJt, VJt_normed)
        u = projection(randn(size(VJt)), VJt_normed);
        u = u / norm(u, 'fro');
    end

    % We exploit the fact that VJt_normed belongs to the manifold
    M.rand = @() VJt_normed + randn(1) * randvec(VJt, VJt_normed);

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(d, n);

    M.transp = @(x1, x2, d) d;

    M.pairmean = @(x1, x2) .5*(x1+x2);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [d, n]);
    M.vecmatareisometries = @() true;

end

author: NicolasBoumal

repo: manopt


text:
function draw_dot(adj,varargin);
%DRAW_DOT   Draw a graph.
% DRAW_DOT(ADJ) plots the graph ADJ in the current figure window, using
% 'neato' to optimize the layout.
%
% Optional arguments can be passed as name/value pairs: [default]
%
% 'isbox'     - a vector specifying which nodes should be boxed [0]
% 'rotate'    - rotate the graph so that nodes are vertically aligned [1]
% 'tolerance' - alignment tolerance for 'rotate' [0.001]
% 'start'     - a random seed (to select different solutions)
% 'options'   - a string of command-line options for 'neato' ['']
% All of the optional arguments to graph_to_dot are also supported, such as
% 'node_label'.
%
% See also GRAPH_TO_DOT.
%
% Example:
% size=15; Adj = rand(size) > .8;
% Adj2 = triu(Adj,1)+ triu(Adj,1)' + diag(zeros(size,1));
% draw_dot(Adj2)

% Original: Leon Peshkin
% Modified by Tom Minka

% minka
N = size(adj,1);
unique_labels = cellstr(num2str((1:N)','%-1d'));
labels = unique_labels;
isbox = zeros(N,1);
rotate_flag = 1;
tolerance = 0.001;
options = '';
for i = 1:2:length(varargin)
  switch varargin{i}
    case 'node_label', labels = varargin{i+1};
      % replace with unique labels
      varargin{i+1} = unique_labels;
    case 'isbox', isbox = varargin{i+1};
    case 'rotate', rotate_flag = varargin{i+1};
    case 'tolerance', tolerance = varargin{i+1};
    case 'start', start = varargin{i+1};
      options = [options ' -Gstart=' num2str(start)];
    case 'options', options = [options ' ' varargin{i+1}];
  end
end

if ispc, shell = 'dos'; else, shell = 'unix'; end  %  Which OS ?

cmdline = strcat(shell,'(''neato -V'')');
status = eval(cmdline);
%[status, result] = dos('neato -V');  % request version to check NEATO
if status == 1,  fprintf('Complaining \n'); exit, end

tmpDOTfile = '_GtDout.dot';            % to be platform independant no use of directories
tmpLAYOUT  = '_LAYout.dot';
graph_to_dot(adj > 0, 'filename', tmpDOTfile, 'node_label', unique_labels, varargin{:});  % save in file

cmdline = strcat([shell '(''neato -Tdot ' tmpDOTfile options ' -o ' tmpLAYOUT ''')']); % preserve trailing spaces
status = eval(cmdline);         %  get NEATO todo layout

[adj, permuted_labels, x, y] = dot_to_graph(tmpLAYOUT);  %  load layout
delete(tmpLAYOUT); delete(tmpDOTfile);     % clean up temporary files

% permute the original arguments to match permuted_labels.
order = [];
for i = 1:length(permuted_labels)
  j = strmatch(permuted_labels{i},unique_labels,'exact');
  order(i) = j(1);
end
labels = labels(order);
isbox = isbox(order);
if rotate_flag
  [x,y] = best_rotation(x,y,tolerance);
end

figure(1); clf; axis square      %  now plot
[x, y, h] = draw_graph(adj>0, labels, isbox, x, y, varargin{:});

function [x,y] = best_rotation(x,y,h)
% Rotate the points to maximize the horizontal and vertical alignment.
% Written by Tom Minka.

xm = mean(x);
ym = mean(y);
xr = max(x)-min(x);
yr = max(y)-min(y);
x = (x-xm)/xr;
y = (y-ym)/yr;

xy = [x(:) y(:)];
if 1
  angle = fminbnd(@rotation_cost,-pi/4,pi/4,[],xy,h);
else
  angles = linspace(-pi/4,pi/4,40);
  e = [];
  for i = 1:length(angles)
    e(i) = rotation_cost(angles(i),xy,h);
  end
  %figure(2)
  %plot(angles*180/pi,e)
  angle = angles(argmin(e));
end
%angle*180/pi
c = cos(angle); s = sin(angle);
xy = xy*[c s; -s c];

x = xy(:,1)*xr+xm;
y = xy(:,2)*yr+ym;

function e = rotation_cost(angle,xy,h)
% xy is 2-column matrix.
% e is small if many x's and y's are aligned.

c = cos(angle); s = sin(angle);
xy = xy*[c s; -s c];
dx = sqdist(xy(:,1)',xy(:,1)');
dy = sqdist(xy(:,2)',xy(:,2)');
dx = setdiag(dx,Inf);
dy = setdiag(dy,Inf);
e = sum(exp(-dx(:)/h))+sum(exp(-dy(:)/h));
e = -e;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function L = computeL(xK,K)
% COMPUTEL computes L-function from ripleys K function (L=sqrt(K/pi)-xK)
% L = computeL(xK,K)
% xK - position where K function was estimated
% K - K-function
L=bsxfun(@minus,sqrt(K/pi),xK);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef DCTBasis < spx.dict.Operator

properties(SetAccess=private)
    % Dimensions
    N
end

methods

function self = DCTBasis(N)
    if nargin < 1
        error('Basis dimensions must be specified');
    end
    self.N = N;
end

function [mm, nn]  = get_size(self)
    mm = self.N;
    nn = self.N;
end

function result = apply(self, vectors)
    if size(vectors, 1) ~= self.N
        error('Dimensions mismatch');
    end
    result = idct(vectors);
end

function result = apply_ctranspose(self, vectors)
    if size(vectors, 1) ~= self.N
        error('Dimensions mismatch');
    end
    result = dct(vectors);
end

function result = double(self)
    % Converts the operator into a MATRIX
    n = self.N;
    result = zeros(n);
    for i=1:n
        x = zeros(n, 1);
        x(i) = 1;
        result(:, i) = self.apply(x);
    end
end

function result = norm(self)
    result = 1;
end

end

end

author: indigits

repo: sparse-plex


text:
%run_computeScores_localModel applies the local model to the bounding-box proposals to compute their scores

% SETUP THESE PATHS TO RUN THE CODE
dataPath = 'data/HollywoodHeads';
resultPath = 'results/HollywoodHeads';

% network to evaluate
netFile = fullfile( 'models', 'local.mat' );

% file to store the scores
resultFile = fullfile( resultPath, 'local', 'localModel-scores-test.mat' );

% Choose subset of data to compute the scores.
%   To run the evaluation of the local model you need the test test (3).
%   To run the training of the pairwise model you need all the subsets ([1,2,3]). It might be a good idea to do this in batches on multiple GPUs.
scoreSubset = 3; % 1 - train subset, 2 - validation, 3 - test; can do [1,2,3] to compute scores on all the subsets

%% setup paths
filePath = fileparts( mfilename('fullpath') );
run( fullfile( fileparts( filePath ), 'setup.m' ) );

%% parameters
opts_cnn = struct;
opts_cnn.dataPath = dataPath;
opts_cnn.dataset.trainingSetFile = fullfile('Splits', 'train.txt');
opts_cnn.dataset.validationSetFile = fullfile('Splits', 'val.txt');
opts_cnn.dataset.testSetFile = fullfile('Splits', 'test.txt');
opts_cnn.dataset.groundTruthLocalPrefix = 'Annotations';
opts_cnn.dataset.imageLocalPrefix = 'JPEGImages';
opts_cnn.dataset.candidateLocalPrefix = 'Candidates';

opts_cnn.expDir =  resultPath;
opts_cnn.maxGpuImagesEvaluation = 3000;
opts_cnn.gpuBatchSize = 128;
opts_cnn.scoreMode = 'beforeSoftMax'; % 'beforeSoftMax' or 'afterSoftMax' or 'scoreDifference';

opts_cnn.evaluation.iouThreshold = 0.5;
opts_cnn.evaluation.nmsMaxCandidateNumber = inf;
opts_cnn.evaluation.useDifficultImages = false;

%% load dataset
opts_cnn.imdbPath = fullfile(opts_cnn.expDir, 'imdb.mat');
if exist(opts_cnn.imdbPath, 'file')
    fprintf('Reading imdb file %s\n', opts_cnn.imdbPath);
    imdb = load(opts_cnn.imdbPath);
    if isfield(imdb, 'opts') && isequal( imdb.opts, opts_cnn.dataset )
        fprintf('imdb.opts is compatible with the opt_cnn.dataset\n')
    else
        warning('opts_cnn.dataset parameters are not compatible with the provided imdb file. Be careful!');
    end
else
    warning('imdb file is not found. Making it will take some time.');
    fprintf('Generating imdb file %s\n', opts_cnn.imdbPath);
    imdb = cnn_prepareData_localModel( opts_cnn.dataset, 'dataPath', opts_cnn.dataPath );
    imdb.opts = opts_cnn.dataset;
    mkdir(opts_cnn.expDir);
    save(opts_cnn.imdbPath, '-struct', 'imdb', '-v7.3') ;
end

%% load network
net = load( netFile, '-mat');
if ~isfield(net, 'layers') && isfield(net, 'net')
    net = net.net;
end
if ~isfield(net, 'layers')
    error('Could not load the network!');
end
net = vl_simplenn_move(net, 'gpu');

%% select the set of images to run evaluation
imageSetToDoPr = false( numel( imdb.images.set ), 1 );
for iSubset = 1 : numel( scoreSubset )
    imageSetToDoPr = imageSetToDoPr | ( imdb.images.set(:) == scoreSubset( iSubset ) );
end
imageSetToDoPr = find( imageSetToDoPr );

%% start the evaluation
opts_cnn.getBatchEvaluation = struct;
opts_cnn.getBatchEvaluation.cropMode = 'warp';
opts_cnn.getBatchEvaluation.cropPad = [18, 18, 18, 18];
opts_cnn.getBatchEvaluation.meanImage = net.normalization.averageImage;
opts_cnn.getBatchEvaluation.numPatchesPerImage = inf;
opts_cnn.getBatchEvaluation.maxPositives = inf;
opts_cnn.getBatchEvaluation.maxGpuImages = opts_cnn.maxGpuImagesEvaluation;
opts_cnn.getBatchEvaluation.iouPositiveThreshold = opts_cnn.evaluation.iouThreshold;
opts_cnn.getBatchEvaluation.iouNegativeThreshold = opts_cnn.evaluation.iouThreshold;
opts_cnn.getBatchEvaluation.dataPath = opts_cnn.dataPath;
opts_cnn.getBatchEvaluation.jitterStd = 0;

batchWrapperEvaluation = @(imdb, batch) cnn_getBatch_localModel(imdb, batch, ...
    opts_cnn.getBatchEvaluation) ;

[scores, candidateIds] = cnn_computeScores_localModel( net, imdb, batchWrapperEvaluation, ...
    'imageSet', imageSetToDoPr, ...
    'gpuBatchSize', opts_cnn.gpuBatchSize, ...
    'scoreMode', opts_cnn.scoreMode );

if exist(resultFile, 'file')
    warning('The results file already exists. Overwriting!');
end
if ~exist(fileparts(resultFile), 'dir')
    mkdir(fileparts(resultFile));
end
save( resultFile, 'scores', 'candidateIds', '-v7.3' );

%% save detections to files
detSavePath = fullfile( resultPath, 'local', 'dets');
if ~exist(detSavePath, 'dir')
    mkdir(detSavePath)
end
detSaveFormat = fullfile(detSavePath, '%s.mat');
disp('Saving detection files');
for i=1:length(imageSetToDoPr)
    imgIdx = imageSetToDoPr(i);
    cand_path = fullfile(dataPath, imdb.candidateFiles{imgIdx});
    load(cand_path);
    %load candidate
    BB = [convertBb_X1Y1X2Y2_to_X1Y1WH(boxes(candidateIds{imgIdx},[2 1 4 3])) scores{imgIdx}];
    [~, im_name, ~] = fileparts(imdb.imageFiles{imgIdx});
    savePath = sprintf(detSaveFormat, im_name);
    save(savePath, 'BB');
end

author: aosokin


text:
function [W, S, latent] = PCA(X, pcaDims)
%% [W, S, latent] = PCA(X, pcaDims)
%
% Function for learning a Principle Component Analysis (PCA) subspace.
%
% Inputs:
%   X: training data, with rows corresponding to samples and columns features.
%   pcaDims: optional parameter specifying the output PCA subspace
%       dimensions.
%
% Outputs:
%   W: the learned PCA subspace projection matrix.
%   S: projected coefficients of X by W.
%   latent: latent values of the covariant matrix
%
% Version: 1.0
% Date: 2014-07-22
%
% Author: Shengcai Liao
% Institute: National Laboratory of Pattern Recognition,
%   Institute of Automation, Chinese Academy of Sciences
% Email: scliao@nlpr.ia.ac.cn

[n, d] = size(X);

if nargin == 1
    pcaDims = min(d,n) - 1;
end

if pcaDims <= 0 || pcaDims >= min(d,n)
    error('Incorrect pcaDims.');
end

fprintf('Start to train PCA...');
t0 = tic;

mu = mean(X);
X = bsxfun(@minus, X, mu);

if n >= d
    C = X' * X; %[d,d]
else
    C = X * X'; %[n,n]
end

[W, D] = eig(C);
latent = diag(D);
[latent, index] = sort(latent, 'descend');
W = W(:, index);

W = W(:, 1:pcaDims);
latent = latent(1:pcaDims);

if n < d
    W = X' * W * diag(1 ./ sqrt(latent));
end

if nargout >= 2 && ~isempty(S)
    S = X * W;
end

fprintf('%.3g seconds.\n', toc(t0));

author: happynear

repo: FaceVerification


text:
% ir_mex_build
%
% Some of the functions in this toolbox use calls to mex files.
% The source code for *some*, but not all, of those mex files
% is included with the toolbox, and this script generates the corresponding mex files.

% The full MIRT toolbox with complete mex files can be found at:
% http://web.eecs.umich.edu/~fessler/code/index.html

% This script and the scripts it calls have been tested
% for compiling on Mac, PC, and Linux.

ir_mex_build_mri
ir_mex_build_table

% todo: ir_mex_build_penal

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
#!/usr/bin/octave -q
% Copyright 2018 Marc René Schädler
%
% This file is part of the mobile hearing aid prototype project
% The the mobile hearing aid prototype project is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
%
% The mobile hearing aid prototype project is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with the mobile hearing aid prototype project. If not, see http://www.gnu.org/licenses/.

close all
clear
clc

[signal, fs] = audioread('/tmp/inout.wav');

analysis_length = 100; % ms

analysis_samples = floor(fs.*analysis_length./1000);
analysis_frames = size(signal,1)./analysis_samples;

level_in = -inf*ones(analysis_frames,2);
level_out = -inf*ones(analysis_frames,2);

freq_in = zeros(analysis_frames,2);
freq_out = zeros(analysis_frames,2);

analysis_window = hanning(analysis_samples);
analysis_window = analysis_window ./ sqrt(mean(analysis_window.^2));

for i=1:analysis_frames
  in1 = signal(1+(i-1).*analysis_samples:i.*analysis_samples,1) .* analysis_window;
  in2 = signal(1+(i-1).*analysis_samples:i.*analysis_samples,2) .* analysis_window;
  out1 = signal(1+(i-1).*analysis_samples:i.*analysis_samples,3) .* analysis_window;
  out2 = signal(1+(i-1).*analysis_samples:i.*analysis_samples,4) .* analysis_window;

  level_in(i,1) = 10*log10(mean(in1.^2));
  level_in(i,2) = 10*log10(mean(in2.^2));
  level_out(i,1) = 10*log10(mean(out1.^2));
  level_out(i,2) = 10*log10(mean(out2.^2));

  [~, maxidx_in1] = max(abs(fft(in1)));
  [~, maxidx_in2] = max(abs(fft(in2)));
  [~, maxidx_out1] = max(abs(fft(out1)));
  [~, maxidx_out2] = max(abs(fft(out2)));

  freq_in(i,1) = (maxidx_in1-1)./analysis_samples.*fs;
  freq_in(i,2) = (maxidx_in2-1)./analysis_samples.*fs;
  freq_out(i,1) = (maxidx_out1-1)./analysis_samples.*fs;
  freq_out(i,2) = (maxidx_out2-1)./analysis_samples.*fs;
end

freq_bins = logspace(log10(88.4),log10(16000),16);
freq_color = hsv(length(freq_bins)).*0.8;

figure('Position',[0 0 800 800],'Visible','on');
sample_levels = [-80 -60 -40 -20 0];

for i=1:2
  subplot(2,2,i);
  plot([0 0],[-100 10],'k');
  hold on;
  plot([-100  10],[  0   0],'k');
  plot([-100  10],[-100 10],'k');
  plot([-100 -10],[ -80 10],'k--');
  plot([-100 -30],[ -60 10],'k--');
  plot([-100 -50],[ -40 10],'k--');
  plot([-100 -70],[ -20 10],'k--');
  plot([-100 -90],[   0 10],'k--');

  out_mean = nan(length(freq_bins),length(sample_levels));

  h = [];
  for j=1:length(freq_bins)
    if j == 1
      freq_mask = freq_in(:,i) >= freq_bins(j)  & freq_in(:,i) < freq_bins(j+1);
    elseif j==length(freq_bins)
      freq_mask = freq_in(:,i) >= freq_bins(j-1) & freq_in(:,i) < freq_bins(j);
    else
      freq_mask = freq_in(:,i) >= freq_bins(j-1) & freq_in(:,i) < freq_bins(j+1);
    end
    x = level_in(freq_mask,i);
    y = level_out(freq_mask,i);
    validmask = x > -100 & y > -100 & ~isnan(x) & ~isnan(y) & ~isinf(x) & ~isinf(y);
    x = x(validmask);
    y = y(validmask);
    c = freq_color(j,:);
    if length(x) > 6
      p = polyfit(x,y,4);
      out_mean(j,:) = polyval(p,sample_levels);
    end
    scatter(x,y,20,c);
    if mod(j-1,2) == 0
      h(j) = plot(-90:10,polyval(p,-90:10),'Color',c);
    end
  end
  legend(h,num2str(round(freq_bins(2:2:end).')),'location','southeast');

  xlim([-100 10]);
  ylim([-100 10]);
  xlabel('Input / dB FS');
  ylabel('Output / dB FS');

  subplot(2,2,i+2);

  plot(log([50 16000]),[0 0],'k');
  hold on;

  level_color = hsv(length(sample_levels)).*0.8;

  h = [];
  for j=1:length(sample_levels)
    c = level_color(j,:);
    plot(log(freq_bins(1:end-1)),sample_levels(j).*ones(length(freq_bins)-1,1),'--','Color',c);
    h(j) = plot(log(freq_bins),out_mean(:,j),'-','Color',c);
  end

  xlim(log([50 16000]));
  ylim([-100 10]);
  xlabel('Frequency / Hz');
  ylabel('Levels / dB');
  set(gca,'xtick',log(freq_bins(2:2:end-1)));
  set(gca,'xticklabel',round(freq_bins(2:2:end-1)));
  legend(h,num2str(sample_levels.'),'location','southeast');
end

set(gcf,'PaperUnits','inches','PaperPosition',[0 0 8 8].*1.4);
print('-depsc2','-r300','inout.eps');

author: m-r-s

repo: hearingaid-prototype


text:
function result = graph_adj_is_node_connected ( adj, nnode )

%*****************************************************************************80
%
%% GRAPH_ADJ_IS_NODE_CONNECTED determines if a graph is nodewise connected.
%
%  Definition:
%
%    A graph is nodewise connected if, from every node, there is a path
%    to any other node.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ADJ(NNODE,NNODE), the adjacency matrix for the
%    graph.  ADJ(I,J) is nonzero if there is an edge from node I to node J.
%
%    Input, integer NNODE, the number of nodes.
%
%    Output, integer RESULT.
%    0, the graph is not nodewise connected.
%    1, the graph is nodewise connected.
%

%
%  FOUND(I) is 1 if node I has been reached.
%  LIST(I) contains a list of the nodes as they are reached.
%
  list(1:nnode) = 0;
  found(1:nnode) = 0;
%
%  Start at node 1.
%
  found(1) = 1;
  list(1) = 1;
  ilo = 1;
  ihi = 1;
%
%  From the batch of nodes found last time, LIST(ILO:IHI),
%  look for unfound neighbors, and store their indices in LIST(JLO:JHI).
%
  while ( 1 )

    jlo = ihi + 1;
    jhi = ihi;

    for ii = ilo : ihi

      i = list(ii);

      for j = 1 : nnode

        if ( adj(i,j) ~= 0 || adj(j,i) ~= 0 )

          if ( found(j) == 0 )
            jhi = jhi + 1;
            list(jhi) = j;
            found(j) = 1;
          end

        end

      end

    end
%
%  If no neighbors were found, exit.
%
    if ( jhi < jlo )
      break
    end
%
%  If neighbors were found, then go back and find THEIR neighbors.
%
    ilo = jlo;
    ihi = jhi;

  end
%
%  No more neighbors were found.  Have we reached all nodes?
%
  if ( ihi == nnode )
    result = 1;
  else
    result = 0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = process_evt_detect_eog( varargin )
% PROCESS_EVT_DETECT_EOG: Detect eye blinks in a continuous file, and create set of events called "blink"
%
% USAGE:  OutputFiles = process_evt_detect_eog('Run', sProcess, sInputs)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, Elizabeth Bock, 2011-2016

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Detect eye blinks';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 44;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/ArtifactsDetect#Detection:_Blinks';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data'};
    sProcess.OutputTypes = {'raw', 'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Channel name
    sProcess.options.channelname.Comment = 'Channel name: ';
    sProcess.options.channelname.Type    = 'channelname';
    sProcess.options.channelname.Value   = '';
    % Channel name comment
    sProcess.options.channelhelp.Comment = 'You can use the montage syntax here: "ch1, -ch2"';
    sProcess.options.channelhelp.Type    = 'label';
    % Time window
    sProcess.options.timewindow.Comment = 'Time window:';
    sProcess.options.timewindow.Type    = 'timewindow';
    sProcess.options.timewindow.Value   = [];
    % Event name
    sProcess.options.eventname.Comment = 'Event name: ';
    sProcess.options.eventname.Type    = 'text';
    sProcess.options.eventname.Value   = 'blink';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Pre-defined set of properties
    %sProcess.options.eventname.Value = 'blink';
    sProcess.options.threshold.Value = {2, 'x std', 2};
    sProcess.options.bandpass.Value  = {[1.5 15], 'Hz', 2};
    sProcess.options.blanking.Value  = {0.8, 'ms', []};
    sProcess.options.isnoisecheck.Value = 1;
    sProcess.options.isclassify.Value   = 1;
    % Call the generic version of the event detection
    OutputFiles = process_evt_detect('Run', sProcess, sInputs);
end

author: brainstorm-tools


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function qu = unit(q)
%Quaternion.unit Unitize a quaternion
%
% QU = Q.unit() is a quaternion which is a unitized version of  Q

  qu = q / norm(q);
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function p13_demo ( iteration_max, h )

%*****************************************************************************80
%
%% P13_DEMO runs the 2D demo problem #13, with mesh size H.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2006
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer ITERATION_MAX, the maximum number of iterations that DISTMESH
%    should take.  (The program might take fewer iterations if it detects convergence.)
%
%    Input, real H, the mesh spacing parameter.
%
  if ( nargin < 1 )
    iteration_max = 200;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P13_DEMO - Note:\n' );
    fprintf ( 1, '  No value of ITERATION_MAX was supplied.\n' );
    fprintf ( 1, '  The default value ITERATION_MAX = %d will be used.\n', ...
      iteration_max );
  end

  if ( nargin < 2 )
    h = 0.025;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P13_DEMO - Note:\n' );
    fprintf ( 1, '  No value of H was supplied.\n' );
    fprintf ( 1, '  The default value H = %f will be used.\n', h );
  end
%
%  Put the random number generator into a fixed initial state.
%
  rand ( 'state', 111 );
%
%  Set the rendering method for the current figure to Z-buffering.
%
  set ( gcf, 'rend', 'z' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 13:\n' );
  fprintf ( 1, '  The Sandia Fork, h = %f\n', h )

  fd = @p13_fd;
  fh = @p13_fh;
  box = [0.0, 0.0; 1.0, 1.0];

  fixed = [ ...
   0.10000, 0.00000; ...
   0.20000, 0.00000; ...
   0.80000, 0.00000; ...
   0.90000, 0.00000; ...
   0.55000, 0.39686; ...
   0.55000, 0.90000; ...
   0.45000, 0.90000; ...
   0.45000, 0.39686 ];

  [ p, t ] = distmesh_2d ( fd, fh, h, box, iteration_max, fixed );

  post_2d ( p, t, fh )
%
%  Write a PostScript image of the triangulation.
%
  [ node_num, junk ] = size ( p );
  [ tri_num, junk ] = size ( t );
  p = p';
  t = t';
  node_show = 0;
  triangle_show = 1;

  triangulation_order3_plot ( 'p13_mesh.eps', node_num, p, tri_num, ...
    t, node_show, triangle_show );
%
%  Write a text file containing the nodes.
%
  r8mat_write ( 'p13_nodes.txt', 2, node_num, p );
%
%  Write a text file containing the triangles.
%
  i4mat_write ( 'p13_elements.txt', 3, tri_num, t );

  return
end


text:
function positions = GetPositions(varargin)

%GetPositions - Get position samples.
%
%  USAGE
%
%    positions = GetPositions(<options>)
%
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'mode'         'all' gets all position samples, 'clean' discards bad
%                    position samples (lights undetected, out of boundaries,
%                    too close, etc. - see SETTINGS) (default 'clean')
%     'coordinates'  'video' gets position samples in video frame coordinates
%                    (e.g., [0..368]x[0..284]), 'normalized' uses normalized
%                    values (i.e. [0..1]x[0..1]), 'real' converts pixels to
%                    centimeters (default 'normalized')
%     'pixel'        size of the video pixel in cm (no default value)
%     'discard'      discard position samples when one or more LED is missing
%                    ('partial', default) or only when no LEDs are detected
%                    ('none', coordinates for missing LEDs are set to NaN)
%     'distances'    min and max distances [m M] allowed between LEDs
%                    (in pixels, default [0 Inf])
%    =========================================================================
%
%  EXAMPLES
%
%    p = GetPositions;
%    p = GetPositions('mode','all');
%
%  CUSTOM DEFAULTS
%
%    Properties 'mode', 'coordinates', 'pixel' and 'distances' can have custom
%    default values (type 'help CustomDefaults' for details).

% Copyright (C) 2004-2014 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

global DATA SETTINGS;
if isempty(DATA),
	error('No session defined (did you forget to call SetCurrentSession? Type ''help Data'' for details).');
end

% Default values (customizable defaults must be empty at this point)
mode = '';
coordinates = '';
pixel = [];
distances = [];
discard = 'partial';

if mod(length(varargin),2) ~= 0,
	error('Incorrect number of parameters (type ''help GetPositions'' for details).');
end

% Parse options
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+1) ' is not a property (type ''help GetPositions'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'mode',
			mode = lower(varargin{i+1});
			if ~isstring_FMAT(mode,'clean','all'),
				error('Incorrect value for property ''mode'' (type ''help GetPositions'' for details).');
			end
		case 'coordinates',
			coordinates = lower(varargin{i+1});
			if ~isstring_FMAT(coordinates,'video','normalized','real'),
				error('Incorrect value for property ''coordinates'' (type ''help GetPositions'' for details).');
			end
		case 'pixel',
			pixel = varargin{i+1};
			if ~isdscalar(pixel,'>0'),
				error('Incorrect value for property ''pixel'' (type ''help GetPositions'' for details).');
			end
		case 'discard',
			discard = lower(varargin{i+1});
			if ~isstring_FMAT(discard,'partial','none'),
				error('Incorrect value for property ''discard'' (type ''help GetPositions'' for details).');
			end
		case 'distances',
			distances = lower(varargin{i+1});
			if ~isdvector(distances,'#2','>=0'),
				error('Incorrect value for property ''distances'' (type ''help GetPositions'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help GetPositions'' for details).']);
	end
end

% Customizable defaults (if not already set)
mode = GetCustomDefaults(mode,'mode','clean');
coordinates = GetCustomDefaults(coordinates,'coordinates','normalized');
pixel = GetCustomDefaults(pixel,'pixel',[]);
distances = GetCustomDefaults(distances,'distances',[0 Inf]);

minDistance = distances(1);
maxDistance = distances(2);

if strcmp(coordinates,'real') & isempty(pixel),
	error(['Missing ''pixel'' property-value pair (type ''help GetPositions'' for details).']);
end

positions = DATA.positions;
if isempty(positions), return; end

% Discard positions with incorrect distance between lights
if strcmp(mode,'clean'),
	if size(positions,2) == 5,
		% Two head lamps
		distance = sqrt((positions(:,4)-positions(:,2)).^2+(positions(:,5)-positions(:,3)).^2);
		selected = ( ...
				distance >= minDistance ...
				& distance <= maxDistance ...
			);
		positions = positions(selected,:);
	end
end

% Mark undetects (by convention, values of -1)
if strcmp(discard,'none'),
	undetected = all(positions(:,2:end) == -1,2);
	partiallyDetected = positions == -1;
	positions(partiallyDetected) = NaN;
else
    undetected = any(positions(:,2:end) == -1,2);

end
%  positions(undetected,2:end) = -1; % Make sure all values are -1
good = ~undetected;

% Clip if necessary (but issue warning)
m = min(positions(good,2:end));
M = max(positions(good,2:end));
%  if any(m<0)||any(M(1:2:end)>DATA.maxX|M(2:2:end)>DATA.maxY),
%  	warning(['Position data should stay within [0,' num2str(DATA.maxX) ']x[0,' num2str(DATA.maxY) ']. The data will now be clipped accordingly.']);
%  	positions(good,2) = Clip(positions(good,2),0,DATA.maxX);
%  	positions(good,3) = Clip(positions(good,3),0,DATA.maxY);
%  end

if strcmp(coordinates,'normalized'),
	% Normalize coordinates
	maxima = M(1:2);%[DATA.maxX DATA.maxY];
	positions(good,2:end) = positions(good,2:end) ./ repmat(maxima,sum(good),(size(positions,2)-1)/2);
elseif strcmp(coordinates,'real'),
	% Convert to cm
	positions(good,2:end) = positions(good,2:end) * pixel;
end

% Discard undetects
if strcmp(mode,'clean'),
	positions(undetected,:) = [];
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function r = subsasgn(r, s, v)

	if s(1).type  ~= '.'
		error('only .field supported')
	end
	switch s(1).subs,
	case 'links'
		r.links = v;
	case 'lineopt',
		r.lineopt = v;
	case 'shadowopt',
		r.shadowopt = v;
	case 'offset',
		L = r.links;
		for i=1:r.n,
			L(i).offset = v(i);
		end
	case 'qlim',
		if numrows(v) ~= r.n,
			error('insufficient rows in joint limit matrix');
		end
		L = r.links;
		for i=1:r.n,
			L(i).qlim = v(i,:);
		end
	case 'q',
		r.q = v;
	case 'handle',
		r.handle = v;
	case 'plotopt',
		r.plotopt = v;
	case 'gravity',
		r.gravity = v;
	case 'tool',
		if ~ishomog(v)
			error('base must be a homogeneous transform');
		end
		r.tool = v;
	case 'base',
		if ~ishomog(v)
			error('base must be a homogeneous transform');
		end
		r.base = v;
	case 'name',
		r.name = v;
	case 'manufacturer',
		r.manufacturer = v;
	case 'comment',
		r.comment = v;
	otherwise, error('Unknown method in subsasgn')
	end


text:
function xf_detect = spike_detection_filter(x, par)
%this function filter the signal, using the detection filter. Is used in the
%readInData class.

sr = par.sr;
fmin_detect = par.detect_fmin;
fmax_detect = par.detect_fmax;

% HIGH-PASS FILTER OF THE DATA
if par.detect_order>0
    if exist('ellip','file')                         %Checks for the signal processing toolbox
        [b,a] = ellip(par.detect_order,0.1,40,[fmin_detect fmax_detect]*2/sr);
        if exist('FiltFiltM','file')
            xf_detect = FiltFiltM(b, a, x);
        else
            xf_detect = filtfilt(b, a, x);
        end
    else
        xf_detect = fix_filter(x);                   %Does a bandpass filtering between [300 3000] without the toolbox.
    end
else
    xf_detect = x;
end

author: csn-le


text:
function [strctident,signreslabels]=loadcorrelres(strctident,endo,names,startdate,enddate,lags,n,IRFt,favar,pref)

signreslabels=strctident.signreslabels;
signreslabels_shocksindex=strctident.signreslabels_shocksindex;

if strctident.CorrelInstrument=="" && strctident.CorrelShock=="" %if  CorrelInstrument CorrelShock are empty then skip this part completly
    strctident.checkCorrelInstrumentShock=0; %there is nothing to check
    strctident.hbartext_CorrelInstrumentShock='';
else
    % however if we have no sign, relmagn FEVD restrictions at all, but we have correl restrictions then generate them
    if strctident.signres==0 && strctident.zerores==0 && strctident.magnres==0 && strctident.favar_signres==0 && strctident.favar_zerores==0 && strctident.favar_magnres==0 && strctident.relmagnres==0 && strctident.favar_relmagnres==0 && strctident.FEVDres==0 && strctident.favar_FEVDres==0
        % loop over endogenous (columns)
        for ii=1:n
            % give the generic name 'shock ii'
            if IRFt==6 && ii==1 % special case for the IV shock in IRFt==6
                signreslabels{ii,1}=strcat('IV Shock (',strctident.Instrument,')');
                signreslabels_shocksindex=[signreslabels_shocksindex; ii];
            else
                signreslabels{ii,1}=['shock ' num2str(ii)];
                if strctident.signrestableempty(ii)==0 % restrictions, however no label is found, count this column
                    signreslabels_shocksindex=[signreslabels_shocksindex; ii];
                end
            end
        end
        % save the shock index to later determine the number of identified shocks
        strctident.signreslabels_shocksindex=unique(signreslabels_shocksindex);
        strctident.signreslabels_shocks=signreslabels(strctident.signreslabels_shocksindex);
        strctident.signreslabels=signreslabels;
        strctident.periods=[0;0]; % also provide some periods here to generate irfs
    end

    % check if we have IV correlation restrictions
    [IVcorrel,txtcorrel]=xlsread(pref.excelFile,'IV');
    checkCorrelInstrument_index=ismember(txtcorrel(1,2:end),strctident.CorrelInstrument);
    checkCorrelShock_index=ismember(strctident.signreslabels,strctident.CorrelShock); %%%%% evt. change to strcmp loop

    %% we might have a correl shock only (no other restrictions), identify
    if sum(checkCorrelShock_index)==0
        % load the data from Excel
        % sign restrictions values
        if strctident.signres==1
            [~,~,strngs1]=xlsread(pref.excelFile,'sign res values');
        elseif strctident.relmagnres==1
            [~,~,strngs1]=xlsread(pref.excelFile,'relmagn res values');
        elseif strctident.FEVDres==1
            [~,~,strngs1]=xlsread(pref.excelFile,'FEVD res values');
        else
            [~,~,strngs1]=xlsread(pref.excelFile,'sign res values');
        end
        % replace NaN entries by blanks
        strngs1(cellfun(@(x) any(isnan(x)),strngs1))={[]};
        % convert all numeric entries into strings
        strngs1(cellfun(@isnumeric,strngs1))=cellfun(@num2str,strngs1(cellfun(@isnumeric,strngs1)),'UniformOutput',0);
        % identify the non-empty entries (pairs of rows and columns), changed the
        % routine here
        % empty strngs1 columns and rows, to make sure empty rows (with empty strings) are dismissed
        for ii=1:size(strngs1,1) %rows
            strngs1emptyrows(ii,1)=isempty(cat(2,strngs1{ii,:}));
        end
        strngs1emptyrows_index=find(strngs1emptyrows==0);

        for ii=1:size(strngs1,2) %columns
            strngs1emptycolumns(1,ii)=isempty(cat(2,strngs1{:,ii}));
        end
        strngs1emptycolumns_index=find(strngs1emptycolumns==0);
        %
        strngs1=strngs1(strngs1emptyrows_index,strngs1emptycolumns_index);

        [nerows1,neclmns1]=find(~cellfun('isempty',strngs1));
        % count the number of such entries
        neentries1=size(nerows1,1);
        % all these entries contrain strings: fix them to correct potential user formatting errors
        % loop over entries (value table)
        for ii=1:neentries1
            strngs1{nerows1(ii,1),neclmns1(ii,1)}=bear.utils.fixstring(strngs1{nerows1(ii,1),neclmns1(ii,1)});
        end

        % loop over endogenous variables
        for ii=1:n
            % for each variable, there should be two entries in the table corresponding to its name
            % one is the column lable, the other is the row label
            [r,c]=find(strcmp(strngs1,endo{ii,1}));
            % if it is not possible to find two entries, return an error
            if size(r,1)<2
                message=['Endogenous variable ' endo{ii,1} ' cannot be found in both rows and columns of the table. Please verify that the ''sign res values'' sheet of the Excel data file is properly filled.'];
                msgbox(message,'Sign restriction error');
                error('programme termination: sign restriction error');
            end
            % otherwise, the greatest number in r corresponds to the row of the column labels: record it
            rows(ii,1)=max(r);
            % the greatest number in c corresponds to the column of the row labels: record it
            clmns(ii,1)=max(c);
        end

        % find the shock, %%here a user specified label is provided in  the res
        % sheets in the excel file
        [r2,c2]=find(strcmp(strngs1,strctident.CorrelShock));

        if ~isempty(r2)==1 && ~isempty(c2)==1
            c2index=find(c2==clmns);
            signreslabels{c2index,1}=strctident.CorrelShock;
            signreslabels_shocksindex=[signreslabels_shocksindex;c2index];

            %save in structure
            strctident.signreslabels=signreslabels;
            strctident.signreslabels_shocksindex=sort(signreslabels_shocksindex);

            %again, check if we have IV correlation restrictions
            [IVcorrel,txtcorrel]=xlsread(pref.excelFile,'IV');
            checkCorrelInstrument_index=ismember(txtcorrel(1,2:end),strctident.CorrelInstrument);
            checkCorrelShock_index=ismember(strctident.signreslabels,strctident.CorrelShock); %%%%% evt. change to strcmp loop
        elseif isempty(r2)==1 && isempty(c2)==1
            % do nothing in this case
        end

    end
    % check for correlation restrictions
    % when there is no shock with an extra instrument
    if sum(checkCorrelShock_index)==0 && sum(checkCorrelInstrument_index)==0
        strctident.checkCorrelInstrumentShock=0; %there is nothing to check
        strctident.hbartext_CorrelInstrumentShock='';
        % warning if we found a CorrelShock only
    elseif sum(checkCorrelShock_index)~=0 && sum(checkCorrelInstrument_index)==0
        strctident.checkCorrelInstrumentShock=0;
        strctident.hbartext_CorrelInstrumentShock='';
        message=['Found CorrelShock ' strctident.CorrelShock ', but CorrelInstrument ' strctident.CorrelInstrument ' cannot be found in the "IV" excel sheet. Correlation restrictions are ignored.'];
        msgbox(message,'Correlation restriction warning');
        % warning if we found a CorrelShock only
    elseif sum(checkCorrelShock_index)==0 && sum(checkCorrelInstrument_index)~=0 && strcmp(strctident.CorrelShock,'CorrelShock')==0
        strctident.checkCorrelInstrumentShock=0;
        strctident.hbartext_CorrelInstrumentShock='';
        message=['Found CorrelInstrument ' strctident.CorrelInstrument ', but CorrelShock ' strctident.CorrelShock ' cannot be found. Correlation restrictions are ignored.'];
        msgbox(message,'Correlation restriction warning');
        % we found a CorrelInstrument and a CorrelShock
    elseif sum(checkCorrelInstrument_index)~=0 && sum(checkCorrelShock_index)~=0 | strcmp(strctident.CorrelShock,'CorrelShock')==1
        % find index for the Shock with correlation restriction, %%here no
        % label is provided in the res sheets in the excel file
        strctident.CorrelShock_index=find(checkCorrelShock_index); %save index
        strctident.checkCorrelInstrumentShock=1;
        Index = strcmp(txtcorrel(1,:),strctident.CorrelInstrument);           %find the instrument in the IV sheet
        IVnum = find(Index==1, 1, 'first')-1;
        IVcorrel = IVcorrel(:, IVnum);
        IVcorrel = IVcorrel(~isnan(IVcorrel));
        txtcorrel = txtcorrel(2:length(IVcorrel)+1,1);              % drop IV names from txt
        date = names(2:end,1);                                   %get the datevector of the VAR
        startlocationY_in_Y=find(strcmp(date,startdate));        %location of sample startdate in Y datevector
        endlocationY_in_Y=find(strcmp(date,enddate));            %location of sample enddate in Y datevector
        date = date(startlocationY_in_Y+lags:endlocationY_in_Y,:);  %cut datevector of Y such that it corresponds to the time dates used in the VAR
        strctident.OverlapIVcorrelinY = ismember(date,txtcorrel);           %Use this to cut EPS
        strctident.OverlapYinIVcorrel = ismember(txtcorrel,date);           %Use this to cut IV
        strctident.IVcorrel = IVcorrel(strctident.OverlapYinIVcorrel,:);               %cut all the entries from IV that are not in the sample
        strctident.hbartext_CorrelInstrumentShock='correlation, ';

        % check for empty other res shocks and non-empty correl res shocks
        correlshocktest=0;
        for ii=1:size(strctident.signreslabels,1)
            if correlshocktest==1
                break
            end
            if strctident.signrestableempty(ii,1)==1 && strctident.favar_signrestableempty(ii,1)==1 && strctident.relmagnrestableempty(ii,1)==1 && strctident.favar_relmagnrestableempty(ii,1)==1 && strctident.FEVDrestableempty(ii,1)==1 && strctident.favar_FEVDrestableempty(ii,1)==1 && checkCorrelShock_index(ii,1)==0 && strcmp(strctident.CorrelShock,'CorrelShock')==1
                % add them to the list of identified shocks in this case
                strctident.signreslabels{ii,1}=['CorrelShock (' strctident.CorrelInstrument ')'];    % or generate a label
                strctident.signreslabels_shocksindex=[strctident.signreslabels_shocksindex; ii]; % add it in the shocksindex
                checkCorrelShock_index(ii,1)=1;
                correlshocktest=1;
            end
        end

        % update the output in strctident
        strctident.signreslabels_shocks=strctident.signreslabels(strctident.signreslabels_shocksindex);
        signreslabels=strctident.signreslabels;
        strctident.CorrelShock_index=find(checkCorrelShock_index);

        %finally check if there are sign restrictions on the shock of interest. If
        %not we can use the flipped entry of the rotation matrix aswell
        if ~isempty(strctident.CorrelShock)
            if isempty(strctident.Scell{1,strctident.CorrelShock_index})
                strctident.FlipCorrel=1;
            else
                strctident.FlipCorrel=0;
            end
        end
    end
end

if favar.FAVAR==1
    % check here if the shocks of interest actually exist
    if favar.IRF.npltXshck > size(strctident.signreslabels_shocks,1)
        % error if no shock to plot is found, otherwise code crashes at a later stage
        message=['Error: Shock(s) (' favar.IRFplotXshock ') cannot be found.'];
        error('bear:loadcorrelres:ShockNotFound',message)
    end
end

author: european-central-bank

repo: BEAR-toolbox