text:
function tests = tapas_physio_examples_test()
% Tests all PhysIO examples both as matlab scripts and via batch editor and
% compares relevant parts of physio-output structure and
% multiple_regressors file to reference results, saved by
% physio_update_examples when deploying
%
%    tests = tapas_physio_examples_test()
%
% NOTE: In order to run these tests, the corresponding example data
%       for PhysIO has be downloaded via tapas_download_example_data
%       (in misc/ subfolder of tapas release)
%
% IN
%
% OUT
%
% EXAMPLE
%   tapas_physio_examples_test
%
%   See also tapas_download_example_data physio_update_examples

% Author:   Lars Kasper
% Created:  2019-08-12
% Copyright (C) 2019 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under
% the terms of the GNU General Public License (GPL), version 3. You can
% redistribute it and/or modify it under the terms of the GPL (either
% version 3 or, at your option, any later version). For further details,
% see the file COPYING or <http://www.gnu.org/licenses/>.

tests = functiontests(localfunctions);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Auxiliary functions for automation of code folder structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% path to examples, needed for all test cases
function setupOnce(testCase)
% run GE example and extract physio
testCase.TestData.pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..', '..');
% TODO: Make generic!
testCase.TestData.pathExamples =  fullfile(testCase.TestData.pathPhysioPublic, '..', 'examples');
end

% close all created figures from examples after each test
function teardown(testCase)
close(testCase.TestData.createdFigHandles);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Matlab-only tests (w/o SPM) start here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function test_bids_cpulse3t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of BIDS/CPULSE3T example using matlab only
dirExample = 'BIDS/CPULSE3T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_bids_ppu3t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of BIDS/PPU3T example using matlab only
dirExample = 'BIDS/PPU3T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_bids_ppu3t_separate_files_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of BIDS/PPU3T_Separate_Files example using matlab only
dirExample = 'BIDS/PPU3T_Separate_Files';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_biopac_txt_ppu3t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of BioPac_txt/PPU3T example using matlab only
dirExample = 'BioPac_txt/PPU3T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_ge_ppu3t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of GE/PPU3T example using matlab only
dirExample = 'GE/PPU3T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_philips_ecg3t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Philips/ECG3T example using matlab only
dirExample = 'Philips/ECG3T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_philips_ecg3t_v2_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Philips/ECG3T_V2 example using matlab only
dirExample = 'Philips/ECG3T_V2';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_philips_ecg7t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Philips/ECG7T example using matlab only
dirExample = 'Philips/ECG7T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_philips_ppu3t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Philips/PPU3T example using matlab only
dirExample = 'Philips/PPU3T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_hcp_ppu3t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_HCP/PPU3T example using matlab only
dirExample = 'Siemens_HCP/PPU3T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vb_ecg3t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VB/ECG3T example using matlab only
dirExample = 'Siemens_VB/ECG3T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vb_ppu3t_sync_first_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VB/PPU3T example using matlab only
% synced to first DICOM volume of run
dirExample = 'Siemens_VB/PPU3T_Sync_First';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vb_ppu3t_sync_last_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VB/PPU3T example using matlab only
% synced to last DICOM volume of run
dirExample = 'Siemens_VB/PPU3T_Sync_Last';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vd_ppu3t_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VD/PPU3T example using matlab only
dirExample = 'Siemens_VD/PPU3T';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vd_ppu3t_for_bids_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VD/PPU3T_For_BIDS example using matlab only
dirExample = 'Siemens_VD/PPU3T_For_BIDS';
doUseSpm = false;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vd_ppu3t_for_bids_vs_bids_converted_matlab_only(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% from same Siemens VD data externally converted to BIDS
% to current output of re-run of Siemens_VD/PPU3T_For_BIDS example using matlab only
% !!!TODO: FIX Results, not equivalent to external BIDS converter results
dirExample = 'Siemens_VD/PPU3T_For_BIDS';
dirRefResults = 'BIDS/PPU3T_Separate_Files';
doUseSpm = false;
idxTests = []; % empty for now, should be 1:5, or at least 4:5
run_example_and_compare_reference(testCase, dirExample, doUseSpm, ...
    dirRefResults, idxTests)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SPM-requiring tests start here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function test_bids_cpulse3t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of BIDS/CPULSE3T example using SPM Batch Editor
dirExample = 'BIDS/CPULSE3T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_bids_ppu3t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of BIDS/PPU3T example using SPM Batch Editor
dirExample = 'BIDS/PPU3T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_bids_ppu3t_separate_files_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of BIDS/PPU3T_Separate_Files example using SPM Batch Editor
dirExample = 'BIDS/PPU3T_Separate_Files';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_biopac_txt_ppu3t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of BioPac_txt/PPU3T example using SPM Batch Editor
dirExample = 'BioPac_txt/PPU3T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_ge_ppu3t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of GE/PPU3T example using SPM Batch Editor
dirExample = 'GE/PPU3T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_philips_ecg3t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Philips/ECG3T example using SPM Batch Editor
dirExample = 'Philips/ECG3T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_philips_ecg3t_v2_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Philips/ECG3T_V2 example using SPM Batch Editor
dirExample = 'Philips/ECG3T_V2';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_philips_ecg7t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Philips/ECG7T example using SPM Batch Editor
dirExample = 'Philips/ECG7T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_philips_ppu3t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Philips/PPU3T example using SPM Batch Editor
dirExample = 'Philips/PPU3T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_hcp_ppu3t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_HCP/PPU3T example using SPM Batch Editor
dirExample = 'Siemens_HCP/PPU3T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vb_ecg3t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VB/ECG3T example using SPM Batch Editor
dirExample = 'Siemens_VB/ECG3T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vb_ppu3t_sync_first_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VB/PPU3T example
% (sync to first DICOM volume time stamp) using SPM Batch Editor
dirExample = 'Siemens_VB/PPU3T_Sync_First';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vb_ppu3t_sync_last_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VB/PPU3T example
% (sync to last DICOM volume time stamp) using SPM Batch Editor
dirExample = 'Siemens_VB/PPU3T_Sync_Last';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vd_ppu3t_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VD/PPU3T example using SPM Batch Editor
dirExample = 'Siemens_VD/PPU3T';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vd_ppu3t_for_bids_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of Siemens_VD/PPU3T_For_BIDS example using SPM Batch Editor
dirExample = 'Siemens_VD/PPU3T_For_BIDS';
doUseSpm = true;
run_example_and_compare_reference(testCase, dirExample, doUseSpm)
end

function test_siemens_vd_ppu3t_for_bids_vs_bids_converted_with_spm(testCase)
%% Compares previously saved physio-structure and multiple regressors file
% from same Siemens VD data externally converted to BIDS
% to current output of re-run of Siemens_VD/PPU3T_For_BIDS example using SPM Batch Editor
% !!!TODO: FIX Results, not equivalent to external BIDS converter results
dirExample = 'Siemens_VD/PPU3T_For_BIDS';
dirRefResults = 'BIDS/PPU3T_Separate_Files';
doUseSpm = true;
idxTests = []; % empty for now, should be 1:5, or at least 4:5
run_example_and_compare_reference(testCase, dirExample, doUseSpm, ...
    dirRefResults, idxTests)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Auxiliary Functions for automation and code re-use
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function run_example_and_compare_reference(testCase, dirExample, doUseSpm, ...
    dirRefResults, idxTests)
%% Compares previously saved physio-structure and multiple regressors file
% to current output of re-run of example in specified example sub-folder
% Note: both SPM or matlab-script based execution is possible
% (check parameter doUseSpm below!)
%
% IN
%   testCase    TestCase object, required by Matlab in this test function
%   dirExample  subfolder of examples where example data is stored
%               e.g., BIDS/PPU3T
%   doUseSpm    if true, _spm_job.mat versions of example scripts are used, using batch editor
%               if false, _matlab_script.m versions of example scripts are
%               used
%   dirRefResults
%               sub-folder of examples/TestReferenceResults that is used to
%               load expected reference solutions for all test cases
%               default: dirExample is used
%               if you want to cross-check two different integrations based
%               on the same data, this might be useful to cross-reference
%               expected test results (e.g., converted-to-BIDS data vs native
%               vendor read-in)
%   idxTests    Tests to be included in verification
%               sometimes it is necessary to skip some tests, if no
%               equivalency is expected (e.g., using the same data window,
%               but from a shorter logfile, ons_secs.raw will differ)
%               default: [1 2 3 4 5] (all)
%
% OUT
%

if nargin < 5
    idxTests = 1:5;
end

if nargin < 4
    dirRefResults = dirExample;
end

% hard-coded relative tolerance
relTol = 0.01; % 0.01 means 1 percent deviation from expected value allowed

%% Generic settings
% methods for recursively comparing structures, see
% https://ch.mathworks.com/help/matlab/ref/matlab.unittest.constraints.structcomparator-class.html
import matlab.unittest.constraints.IsEqualTo
import matlab.unittest.constraints.AbsoluteTolerance
import matlab.unittest.constraints.RelativeTolerance
import matlab.unittest.constraints.StructComparator
import matlab.unittest.constraints.NumericComparator
import matlab.unittest.constraints.StringComparator

pathExamples = testCase.TestData.pathExamples;
pathCurrentExample = fullfile(pathExamples, dirExample);

%% Actual run of example, via batch editor or as matlab script
if doUseSpm
    pathNow = pwd;
    cd(pathCurrentExample); % for prepending absolute paths correctly

    fileJobMat = [regexprep(lower(dirExample), '/', '_') '_spm_job.mat'];
    fileExample = fullfile(pathCurrentExample, fileJobMat);
    load(fileExample, 'matlabbatch');

    % remove unnecessary verbosity for test
    matlabbatch{1}.spm.tools.physio.verbose.level = 0;

    spm_jobman('run', matlabbatch);
    cd(pathNow);

    dirExampleOutput =  matlabbatch{1}.spm.tools.physio.save_dir{1};
    fileExampleOutputPhysio = matlabbatch{1}.spm.tools.physio.model.output_physio;
    fileExampleOutputTxt = matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors;

else % has verbosity...cannot switch it off

    fileJobMScript = [regexprep(lower(dirExample), '/', '_') '_matlab_script.m'];
    fileExample = fullfile(pathExamples, dirExample, fileJobMScript);

    % runs example Matlab script
    % will output a PhysIO-struct, from which we can parse output files
    run(fileExample);

    % retrieve output files, remove preprending path
    [~, dirExampleOutput] = fileparts(physio.save_dir);
    [~, fn,ext] = fileparts(physio.model.output_physio);
    fileExampleOutputPhysio = [fn ext];
    [~, fn, ext] = fileparts(physio.model.output_multiple_regressors);
    fileExampleOutputTxt = [fn ext];
end

%% Retrieve current results from file
pathExampleOutput = fullfile(pathCurrentExample, dirExampleOutput);

load(fullfile(pathExampleOutput, fileExampleOutputPhysio), 'physio');
R = load(fullfile(pathExampleOutput,fileExampleOutputTxt));
actPhysio = physio;
actRegressorsFromTxt = R;

% for later closing
testCase.TestData.createdFigHandles = physio.verbose.fig_handles;

%% Load reference data and compare to actual run for certain subfields

% Load physio structure from reference data
fileReferenceData = fullfile(pathExamples, 'TestReferenceResults', 'examples', ...
    dirRefResults, fileExampleOutputPhysio);
load(fileReferenceData, 'physio');
expPhysio = physio;

% Compare all numeric sub-fields of physio with some tolerance
% ons_secs has all the computed preprocessed physiological and scan timing
% sync data, from which .model derives the physiological regressors later
% on
doTestOnsSecsRaw            = ismember(1, idxTests);
doTestSpulse                = ismember(2, idxTests);
doTestOnsSecs               = ismember(3, idxTests);
doTestMultipleRegressorsMat = ismember(4, idxTests);
doTestMultipleRegressorsTxt = ismember(5, idxTests);

% 1. Test ons_secs.raw only to check whether data read-in and basic
% filtering before cropping at least worked!
if doTestOnsSecsRaw
    testCase.verifyThat(actPhysio.ons_secs.raw, ...
        IsEqualTo(expPhysio.ons_secs.raw,  ...
        'Using', StructComparator(NumericComparator, 'Recursively', true), ...
        'Within', RelativeTolerance(relTol), ...
        'IgnoringFields',  {'spulse_per_vol'}...
        ), 'Comparing all numeric subfields of ons_secs.raw to check read-in and basic filtering of phys recordings');
end

% 2. Check some crucial timing parameters more vigorously
if doTestSpulse
    verifyEqual(testCase, actPhysio.ons_secs.raw.spulse, expPhysio.ons_secs.raw.spulse, ...
        'RelTol', relTol/10, ...
        'Comparing spulse (onset time of slice pulse (scan) events)');
end

% 3. Test other fields of ons_secs populated during preprocessing and some modeling steps
% Note: spulse_per_vol cannot be compared, because cell!
if doTestOnsSecs
    testCase.verifyThat(actPhysio.ons_secs, ...
        IsEqualTo(expPhysio.ons_secs,  ...
        'Using', StructComparator(NumericComparator, 'Recursively', true), ...
        'Within', RelativeTolerance(relTol), ...
        'IgnoringFields',  {'spulse_per_vol', 'raw'}...
        ), 'Comparing all numeric subfields of ons_secs to check full preprocessing of phys recordings');
end

% recursive with string
% testCase.verifyThat(actPhysio, ...
%     IsEqualTo(expPhysio, 'Using', ...
%     StructComparator(StringComparator, 'Recursively', true), ...
%     'IgnoringCase', true, ...
%     'IgnoringWhitespace', true ...
%     ));

% 4. Compare final multiple regressor matrix in physio.mat structure
% Check Multiple_Regressors output
if doTestMultipleRegressorsMat
    actSolution = actPhysio.model.R;
    expSolution = expPhysio.model.R;

    verifyEqual(testCase, actSolution, expSolution, ...
        'RelTol', relTol, ...
        'Comparing multiple regressors in physio.model.R');
end

%5. Test: Load reference data from multiple regressors txt file and test as
%well

if doTestMultipleRegressorsTxt
    fileReferenceData = fullfile(pathExamples, 'TestReferenceResults', 'examples', ...
        dirRefResults, fileExampleOutputTxt);
    R = load(fileReferenceData);
    expRegressorsFromTxt = R;

    verifyEqual(testCase, actRegressorsFromTxt, expRegressorsFromTxt, ...
        'RelTol', relTol, ...
        'Comparing multiple regressors in txt-files');
end

end

author: translationalneuromodeling

repo: tapas


text:
function G = setupgraph(G,b,s)

if s==1
    for i=1 : size(G,1)
        for j=1 :size(G,1)
            if G(i,j)==0
                G(i,j)=b;
            end
        end
    end
end
if s==2
    for i=1 : size(G,1)
        for j=1 : size(G,1)
            if G(i,j)==b
                G(i,j)=0;
            end
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function vObjRef=relVecAdd(vObsFrame,vObjInFrame)
%%RELVECADD  Special relativistic addition of velocity vectors. In the
%            inertial reference coordinate system, an observer moves with
%            constant velocity vObsFrame. In the coordinate system of the
%            observer, an object moves with constant velocity vObjInFrame.
%            This computes the velocity of the observed object in the
%            inertial reference frame. Under special relativity, it is not
%            just vObsFrame+vObjInFrame as it is in Newtonian mechanics.
%
%INPUTS:  vObsFrame The 3XN set of N velocity vectors in meters per second
%                   of the observer with respect to the inertial reference
%                   coordinate system. The magnitude of the velocity must
%                   be less than the speed of light.
%       vObjInFrame The 3XN set of N velocity vectors in meters per second
%                   of the object with respect to the observer's coordinate
%                   system. The magnitude of the velocity must be less than
%                   or equal to the speed of light.
%
%OUTPUTS: vObjRef  The 3XN set of velocity vectors in vObjInFrame
%                  transformed into the inertial reference coordinate
%                  system.
%
%The formulae for special relativistic velocity addition is derived in
%Chapter 1.4 of [1]. The magnitudes of vObsFrame and vObjInFrame must both
%be less than the speed of light (Constants.speedOfLight).
%
%REFERENCES:
%[1] G. Ludyk, Einstein in Matrix Form: Exact Derivation of the Theory of
%    Special and General Relativity without Tensors. Heidelberg: Springer,
%    2013.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

c=Constants.speedOfLight;

v=vObsFrame;
u=vObjInFrame;

%The magnitudes of all of the v vectors.
vMag=sqrt(sum(v.*v,1));

gamma=1./sqrt(1-vMag.^2/c^2);

%v^T*u for each of the velocity vectors.
uv=sum(v.*u,1);

Num=v+u+bsxfun(@times,(1./gamma-1),(u-bsxfun(@times,(uv./vMag.^2),v)));
Denom=1+uv/c^2;

vObjRef=bsxfun(@rdivide,Num,Denom);

%If vMag=0 for any of the vectors, then NaNs will appear. In such an
%instance, the correct solution is vObjInFrame, because the frame is not
%moving.
colSel=any(~isfinite(vObjRef),1);
vObjRef(:,colSel)=u(:,colSel);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function bernoulli_number3_test ( )

%*****************************************************************************80
%
%% BERNOULLI_NUMBER3_TEST tests BERNOULLI_NUMBER3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BERNOULLI_NUMBER3_TEST\n' );
  fprintf ( 1, '  BERNOULLI_NUMBER3 computes Bernoulli numbers.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I     Exact           BERNOULLI_NUMBER3\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, c0 ] = bernoulli_number_values ( n_data );

    if ( n_data == 0 )
      break
    end

    c1 = bernoulli_number3 ( n );

    fprintf ( 1, '  %2d  %14e  %14e\n', n, c0, c1 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function res = ifft2shift(x)

warning('Deprecated! Consider using ifftnshift(x) instead');
res = ifftnshift(x,1:2);

end

author: thomaskuestner


text:
function Y = VBA_quantile(X,p)
% Y = VBA_QUANTILE(X,P) returns quantiles of the values in the vector X.
% VBA alternative to the quantile function from the statistics toolbox.
% P is a scalar or a vector of cumulative probability values.
% Y is the same size as P, and Y(i) contains the P(i)-th quantile.

% check inputs
if isempty(X)
    Y = nan(numel(p));
    return
end

X(isnan(X)) = [];
X = VBA_vec(X);

% compute quantiels
n = numel(X);

if n==1
    Y = X*ones(1,numel(p));
    return
end
q = [0 (0.5:1:(n-0.5))/n 1];
SX = sort(X);
Y = interp1(q,SX([1 1:n n]),p,'linear');

author: MBB-team

repo: VBA-toolbox


text:
% VL_TEST_IMWBACKWARDMX  Test: imwbackwardmx

I = vl_test_pattern(102) ;

figure(1) ; clf ;
imagesc(I) ;
colormap(gray(256)) ;
axis equal ; axis off ;

[M,N] = size(I) ;
ur    = linspace(-N/2, N/2, N) ;
vr    = linspace(-M/2, M/2, M) ;
[u,v] = meshgrid(ur,vr) ;

for s=.75*(1+cos(linspace(0,2*pi,100)))/2+.25
  up = (cos(s-1)*u - sin(s-1)*v) / s ;
  vp = (sin(s-1)*u + cos(s-1)*v) / s ;
  J = vl_imwbackward(ur,vr,vl_imsmooth(I,1/s*.5),up,vp) ;
  imagesc(J) ; drawnow ;
end

author: jianxiongxiao

repo: ProfXkit


text:
function [maxBSL,BKG,Thresh,equiVol,mPT,bsl] = BSLestimate(PT,maskRTS,voxVol)
%"BSLestimate"
%   Wrapper for Background Subtracted Lesion Estimation -- Returns BSL
%   estimate, Background estimate, Effective Thershold esimate,
%   Equivalent Volume estmate, and WS mask used for estimate.
%
% CRS 05/20/13
%
%Usage:
%   [maxBSL,BKG,Thresh,equiVol,mPT] = BSLestimate(PT,maskRTS,voxVol)
%       PT      = PET images in SUV
%       maskRTS = VOI mask for restricting the BSL estimate
%       voxVol  = volume of a voxel in ml
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.%
%
%% preprocess slices for WS
regPad = [2 2];
s0 = regionprops(maskRTS, {'Centroid','BoundingBox'});
x0(1) = s0.BoundingBox(2); x0(2) = s0.BoundingBox(2+3);
y0(1) = s0.BoundingBox(1); y0(2) = s0.BoundingBox(1+3);
z0(1) = s0.BoundingBox(3); z0(2) = s0.BoundingBox(3+3);
X = floor(x0(1) + 1 - regPad(1):x0(1) + x0(2) + regPad(1));
Y = floor(y0(1) + 1 - regPad(2):y0(1) + y0(2) + regPad(2));
Z = floor(z0(1) + 1:z0(1) + z0(2));
regionPT = PT(X,Y,Z);
maskPT = maskRTS(X,Y,Z);

%% Compute Water Sheds
[Sheds,nShedsT] = BSLwatershed(regionPT,maskPT);
Sheds = maskPT.*Sheds;

Wsheds = struct( ...
    'Shed',Sheds, ...
    'PET',maskPT.*regionPT, ...
    'voxVol',voxVol ...
    );

%% Calculate BSL
[bsl,bkg] = BSLcalcShedsVis(Wsheds,nShedsT,0);

%% Find Max BSL & Background
sBSL = smooth(bsl, max(ceil(numel(bsl)/15),8));
% figure,plot(bsl,'-k'), hold on
% plot(sBSL,'-r'), hold off
[maxBSL, imaxBSL] = max(sBSL);
BKG = bkg(imaxBSL);

% Get WS Region Mask
mPT = zeros(size(PT));
mPTregion = BSLshedMaskVis(Wsheds,imaxBSL);
mPT(X,Y,Z) = mPTregion;

% Get equivalent volume and threshold
[Thresh, equiVol] = BSLthresh(Wsheds,maxBSL);

return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function exactness_test09 ( )

%*****************************************************************************80
%
%% EXACTNESS_TEST09 tests Gauss-Chebyshev2 rules for the Chebyshev2 integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 May 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXACTNESS_TEST09\n' );
  fprintf ( 1, '  Test Gauss-Chebyshev2 rules for the Chebyshev2 integral.\n' );
  fprintf ( 1, '  Density function rho(x) = sqrt(1-x^2).\n' );
  fprintf ( 1, '  Region: -1 <= x <= +1.\n' );
  fprintf ( 1, '  Exactness: 2*N-1.\n' );

  for n = 1 : 5

    [ x, w ] = chebyshev2_set ( n );
    p_max = 2 * n;
    chebyshev2_exactness ( n, x, w, p_max );

  end

  return
end


text:
function success =  runSegForDicomInCondaEnv(inputDicomPath,outputDicomPath,...
    sessionPath,algorithm,condaEnvList,wrapperFunctionList,skipMaskExport)
% function success =  runSegForDicomInCondaEnv(inputDicomPath,outputDicomPath,...
%     sessionPath,algorithm,condaEnvList,wrapperFunctionList,skipMaskExport)
%
% This function serves as a wrapper for different types of segmentations.
%---------------------------------------------------------------------------------------
% INPUT:
% inputDicomPath - path to input DICOM directory which needs to be segmented.
% outputDicomPath - path to write DICOM RTSTRUCT for resulting segmentation.
% sessionPath - path to write temporary segmentation metadata.
% algorithm - string which specifies segmentation algorithm
% --- Optional---
% skipMaskExport - flag to skip export of prior mask
%%---------------------------------------------------------------------------------------
% Following directories are created within the session directory:
% --- ctCERR: contains CERR file/s of input DICOM.
% --- segmentedOrigCERR: CERR file with resulting segmentation fused with
% original CERR file.
% --- segResultCERR: CERR file with segmentation. Note that CERR file can
% be cropped based on initial segmentation.
%
% EXAMPLE:
% inputDicomPath = /path/to/input/dicom/folder;
% outputDicomPath = /path/to/output/dicom/folder;
% sessionPath = '/path/to/directory/for/temp/session';
% algorithm = 'CT_name_of_model';
% condaEnvList = '/location/of/conda/env';
% wrapperFunctionList = '/location/of/inference/wrapper/function'
% success =  runSegForDicomInCondaEnv(inputDicomPath,outputDicomPath,...
%     sessionPath,algorithm,condaEnvList,wrapperFunctionList)
%
% ------------------------------------------------------------------------------------
% APA, AI 04/28/2022

if ~exist('wrapperFunctionList','var') || isempty(skipMaskExport)
    wrapperFunctionList = '';
end

if ~exist('skipMaskExport','var') || isempty(skipMaskExport)
    skipMaskExport = true;
end

%% Create session directory to write segmentation metadata
if inputDicomPath(end) == filesep
    [~,folderNam] = fileparts(inputDicomPath(1:end-1));
else
    [~,folderNam] = fileparts(inputDicomPath);
end
dateTimeV = clock;
randStr = sprintf('%6.3f',rand*1000);
sessionDir = ['session',folderNam,num2str(dateTimeV(4)), num2str(dateTimeV(5)),...
    num2str(dateTimeV(6)), randStr];
fullSessionPath = fullfile(sessionPath,sessionDir);
while exist(fullSessionPath,'dir')
    randStr = sprintf('%6.3f',rand*1000);
    sessionDir = [sessionDir, randStr];
    fullSessionPath = fullfile(sessionPath,sessionDir);
end

%% Create sub-directories
%-For  CERR files
mkdir(fullSessionPath)
cerrPath = fullfile(fullSessionPath,'dataCERR');
mkdir(cerrPath)
outputCERRPath = fullfile(fullSessionPath,'segmentedOrigCERR');
mkdir(outputCERRPath)
segResultCERRPath = fullfile(fullSessionPath,'segResultCERR');
mkdir(segResultCERRPath)
%-For structname-to-label map
labelPath = fullfile(fullSessionPath,'outputLabelMap');
mkdir(labelPath);

%% Import DICOM to CERR
tic
recursiveFlag = true;
importDICOM(inputDicomPath,cerrPath,recursiveFlag);
toc

%% Build DICOM export options

% Parse algorithm and convert to cell arrray
if iscell(algorithm)
    algorithmC = algorithm;
else
    algorithmC = strsplit(algorithm,'^');
end

%% Read user configurations
userOptS = struct();
userOptS(1) = [];
allLabelNamesC = cell(1,length(algorithmC));
dcmExportOptS = struct();
dcmExportOptS(1) = [];
for nAlg = 1:length(algorithmC)

    configFilePath = fullfile(getCERRPath,'ModelImplementationLibrary',...
        'SegmentationModels','ModelConfigurations',...
        [algorithmC{nAlg}, '_config.json']);
    if isempty(userOptS)
        userOptS = readDLConfigFile(configFilePath);
    else
        userOptS = dissimilarInsert(userOptS,readDLConfigFile(configFilePath));
    end

    %Get list of label names
    if ischar(userOptS(nAlg).strNameToLabelMap)
        labelDatS = readDLConfigFile(fullfile(labelPath,...
            userOptS(nAlg).strNameToLabelMap));
        labelMapS = labelDatS.strNameToLabelMap;
    else
        labelMapS = userOptS(nAlg).strNameToLabelMap;
    end
    allLabelNamesC{nAlg} = {labelMapS.structureName};

    % Get DICOM export settings
    if isfield(userOptS(nAlg), 'dicomExportOptS')
        if isempty(dcmExportOptS)
            dcmExportOptS = userOptS(nAlg).dicomExportOptS;
        else
            dcmExportOptS = dissimilarInsert(dcmExportOptS,...
                userOptS(nAlg).dicomExportOptS);
        end
    end
end

%% Load planC
planCfiles = dir(fullfile(cerrPath,'*.mat'));
for p=1:length(planCfiles)

    % Load planC
    planCfiles(p).name
    fileNam = fullfile(planCfiles(p).folder,planCfiles(p).name);
    planC = loadPlanC(fileNam, tempdir);
    planC = quality_assure_planC(fileNam,planC);

    % Run Segmentation
    batchSize = [];
    planC = runSegForPlanCInCondaEnv(planC,sessionPath,algorithm,...
    condaEnvList,wrapperFunctionList,batchSize,skipMaskExport);

    % Save planC to cerrPath
    planC = save_planC(planC,[],'passed',fileNam);

    % Write RTSTRUCT to DICOM
    [~,dcmFileName,~] = fileparts(fileNam);
    for nAlg = 1:length(algorithmC)
        % Get scan index to associate final segmentation
        identifierS = userOptS(nAlg).structAssocScan.identifier;
        if ~isempty(fieldnames(userOptS(nAlg).structAssocScan.identifier))
            origScanNum = getScanNumFromIdentifiers(identifierS,planC);
        else
            origScanNum = 1; %Assoc with first scan by default
        end
        if isempty(dcmExportOptS)
            algExportS = struct();
        else
            algExportS = dcmExportOptS(nAlg);
        end
        exportAISegToDICOM(planC,origScanNum,outputDicomPath,...
            dcmFileName,algExportS(nAlg),allLabelNamesC{nAlg})
    end

end

success = 0;


text:
%VL_IMREADJPEG Load and transform images asynchronously
%   IMAGES = VL_IMREADJPEG(FILES) reads the specified cell array
%   FILES of JPEG files and returns at cell array of images IMAGES.
%
%   IMAGES = VL_IMREADJPEG(FILES, 'NumThreads', T) uses T parallel
%   threads to accelerate the operation. Note that this is independent
%   of the number of computational threads used by MATLAB.
%
%   VL_IMREADJPEG(FILES, 'Prefetch') starts reading the specified
%   images but returns immediately to MATLAB. Reading happens
%   concurrently with MATLAB in one or more separated threads.  A
%   subsequent call IMAGES=VL_IMREADJPEG(FILES) *specifying exactly
%   the same files in the same order* will then return the loaded
%   images. This can be sued to quickly load a batch of JPEG images as
%   MATLAB is busy doing something else.
%
%   The function can transforms the images on the fly in various
%   ways. Transformatiosn are applied as follows:
%
%   1) An (H,W) image is loaded from disk.
%
%   2) A rectangular subset of the image is cropped and resized. The
%      geometry of the crop is determined as follows:
%
%      1) First, the shape (Ho,Wo) of the output image (i.e. the
%         resized crop) is determined. This, as determined by the
%         `Resize` option, can be either the same as the input image
%         (H,W) or one or both of height and witdth can be set to an
%         arbitrary value.
%

%      2) Given the output shape (Ho,Wo) of the crop, the shape
%         (Hi,Wi) and location of the crop in the input image is
%         determined.  First, an anisotropy ratio (change in aspect
%         ratio) is selected according to `CropAnisotropy`. Given
%         that, the input crop rectangle is scaled to fill a certain
%         percentage of the input image according to `CropSize`.
%         Finally, the crop is extracted either from the middle of the
%         input image or at a random location according to
%         `CropLocation`.
%
%   3) The cropedd and resized image undergoes color post
%      post-processing, including mean subtraction
%      (`SubtractAverage`), random color shift (`Brightness`), and
%      random changes in saturation (`Saturation`), and contrats
%      (`Contrast`).
%
%   The function takes the following options:
%
%   `Prefetch`:: not specified
%     If specified, run without blocking (see above).
%
%   `Verbose`:: not specified
%     If specified, increase the verbosity level.
%
%   `Pack`:: not specified
%     If specified, pack all K images in a single H x W x 3 x K
%     array. This requires using the option `'Resize',[H W]` in order
%     to set the output size to be the same for all images. Note
%     furthermore that gray-scale images are automatically extended to
%     three color channels.
%
%   `GPU`:: not specified
%     If specified, return GPU arrays instead of standard arrays. Note
%     that, using the `Prefetch` option, the data is copied to the GPU
%     in parallel with computations.
%
%   `NumThreads`:: `1`
%     Specify the number of threads used to read images. This number
%     must be at least 1. Note that it does not make sense to specify
%     a number larger than the number of available CPU cores, and
%     often fewer threads are sufficient as reading images is memory
%     access bound rather than CPU bound.
%
%   `Resize`:: not specified
%     If specified, turn on image resizing. The argument can either
%     specify the desired height and width `[H, W]` or be simply a
%     scalar `S`. In the latter case, the image is resized
%     isotropically so that the shorter side is equal to `S`.
%
%     Resizing uses bilinear interpolation. When shrinking, a simple
%     form of antialiasing is used by stretching the bilinear filter
%     over several input pixels to average them. The method is the
%     same as MATLAB `imresize` function (the two functions are
%     numerically equivalent).
%
%   `CropAnisotropy`:: `[1 1]`
%     Specify the minimum and maximum value of the anisotropy ratio
%     for the selected crops. The default value [1 1] means that the
%     crops are extracted isotropically, i.e. with the same amount of
%     stretch along horizontal and vertical directions. A value of `[0
%     0]` stretches the crop to fit the input image. A value `[A1 A2]`
%     means that an anisotropy ratio A = Wcrop/Hcrop is sampled
%     uniformly at random for each image from the specified interval.
%
%   `CropSize`:: `[1 1]`
%     Specifies the maximum and minimum size of a crop, as a
%     percentage of the maximum size it can take to fill the input
%     image (the latter is determined *after* fixing the
%     anisotropy). Using `[1 1]` means that the crop is always sized
%     to fill up as much area as possible.
%
%   `CropLocation`:: `'center'`
%     If set to `'random'`, the crop location in the input image is
%     sampled uniformly at random among all possible shifts such that
%     the crop is fully contained in the image.
%
%   `Flip`:: not specified
%     If specified, randomly flips a crop horizontally with 50%
%     probability.
%
%   `SubtractAverage`:: `[0 0 0]`
%     Subtract the specified RGB triplet from each pixel.
%
%   `Brightness`:: `zeros(3)`
%     Standard deviation `B` of the brightness shift. Each image is
%     modified by adding the color shift `B*w`, where `w` is a 3D
%     vector sampled from the standard Normal distribution.
%
%   `Contrast`:: `0`
%     Deviation `C` of the contrast shift. It must be between 0 and
%     1. A value of 0 amounts to no contrast shift. Otherwise, a
%     contrast shift is sampled at random, using `C` as an upper
%     bound.
%
%   `Saturation`:: `0`
%     The same as `Contrast`, but for saturation.
%
%   `Interpolation`:: `'bilinear'`
%     The interpolation method; one of `box`, `bilinear`, `bicubic`,
%     `lanczos2`, and `lanczos3`. The function uses interpolators
%     equivalent to MATLAB's `imresize`.
%
%   Further details on the processing performed by the function can be
%   found in the PDF manual.

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: layumi

repo: Image-Text-Embedding


text:
function menucros() % autogenerated function wrapper
    %  menucros allows the user to establish parameters for zmaps
    %
    % turned into function by Celso G Reyes 2017
    error('replaced (enveloped) by magrcro');
    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun(mfilename('fullpath'));

    %issue a warning if no grid exist
    if ~exist('cumuall', 'var')
        ZmapMessageCenter.set_message('Warning', ' You have to create a grid first!')
        return
    end

    % This is the info window text
    %
    ttlStr='The various z-value maps             ';
    hlpStr1= ...
        ['                                                          '
        ' This window allows you to calculate various z-value maps:'
        ' Once you have created a grid you can apply               '
        ' several different statistics to it. See the              '
        ' User Guide for additional information                    '
        '                                                          '
        ' Maxz-LTA: Uses the LTA statistic and calculates a LTA    '
        ' z-value map for each time-step. Only the maximum value   '
        ' for all time at each grid point is finally plotted in the'
        ' map. Thsi allows to find the most outstanding decreases  '
        ' in the entire catalog.                                   '
        ' ANIMATION: Will produca a movie of several z-value maps  '
        ' at differnt times. Either the LTA or AS statistic will be'
        ' used. ! Uses large amounts of memory!                    '
        ' TIMECUT: Will display one z-value map for a specified    '
        ' time, using LTA, AS, Rubberband or Percent.              '
        '                                                          '
        ' ALARM: Will calculate an Alarm cube using LTA or         '
        ' Rubberband. If any z-value (all time and space) exceeds a'
        ' threshold an alarm is declared and plotted as a "o" in   '
        ' cube representing time and space. Will first calculate   '
        ' Histogram of all z-value                                 '
        '                                                          '];

    hlpStr2= ...
        ['                                                '
        '                                                '];

    % make the interface
    %

    figure(mess);
    report_this_filefun(mfilename('fullpath'));
    clf
    in2 ='input';
    set(gcf,'Name','ZMAP Parameters');
    set(gca,'visible','off');
    set(gcf,'Units','points','pos',[200 200 220 350 ])

    te = text(0.1,0.90,'You can choose from the \newlinefollowing options:');
    set(te,'FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold')

    te3 = text(0.15,0.38,'Timecuts') ;
    set(te3,'FontSize',ZmapGlobal.Data.fontsz.m,'Color','r','FontWeight','bold')

    LTA2_button=uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
        'Position',[.10 .35 .45 .04 ],...
        'Units','normalized',...
        'callback',@(~,~)show_cro('lta',in2),...
        'String','LTA ');

    AS2_button=uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
        'Position',[.10 .28 .45 .04 ],...
        'Units','normalized',...
        'callback',@(~,~)show_cro('ast',in2),...
        'String','AS ');

    per_button=uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
        'Position',[.10 .20 .45 .04 ],...
        'Units','normalized',...
        'callback',@(~,~)show_cro('per',in2),...
        'String','Percent ');

    rub_button=uicontrol('BackGroundColor',[0.8 0.8 0.8]','Style','Pushbutton',...
        'Position',[.10 .12 .45 .04 ],...
        'Units','normalized',...
        'callback',@(~,~)show_cro('rub',in2),...
        'String','Rubberband');

    uicontrol('Style','Pushbutton',...
        'Position',[.70 .01 .25 .08 ],...
        'Units','normalized','Callback',@(~,~)ZmapMessageCenter(),'String','Cancel');

    uicontrol('Units','normal',...
        'Position',[.1 .01 .25 .08],'String','Info ',...
        'callback',@callbackfun_009)

    watchoff

    function callbackfun_009(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        zmaphelp(ttlStr,hlpStr1,hlpStr2);
    end

    function callbackfun_010(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ZmapMessageCenter();
        drawnow;
        incube;
    end

end

author: CelsoReyes


text:
clear; close all; clc;

% add required libraries to the path
addpath(genpath('baseline'));
addpath(genpath('LFCC'));
addpath(genpath('CQCC_v1.0'));
addpath(genpath('GMM'));
addpath(genpath('bosaris_toolkit'));
addpath(genpath('tDCF_v1'));

% set here the experiment to run (access and feature type)
access_type = 'LA'; % LA for logical or PA for physical
feature_type = 'CQCC'; % LFCC or CQCC

% set paths to the wave files and protocols

% TODO: in this code we assume that the data follows the directory structure:
%
% ASVspoof_root/
%   |- data_physical
%      |- ASVspoof2019_LA_dev_asv_scores_v1.txt
% 	   |- ASVspoof2019_LA_dev_v1/
% 	   |- ASVspoof2019_LA_protocols_v1/
% 	   |- ASVspoof2019_LA_train_v1/
%   |- data_physical
%      |- ASVspoof2019_PA_dev_asv_scores_v1.txt
%      |- ASVspoof2019_PA_dev_v1/
%      |- ASVspoof2019_PA_protocols_v1/
%      |- ASVspoof2019_PA_train_v1/

% pathToASVspoof2019Data = '/path/to/ASVspoof_root/';
pathToASVspoof2019Data = './';
if strcmp(access_type, 'LA')
    pathToDatabase = fullfile(pathToASVspoof2019Data, 'data_logical');
else
    pathToDatabase = fullfile(pathToASVspoof2019Data, 'data_physical');
end

trainProtocolFile = fullfile(pathToDatabase, horzcat('ASVspoof2019_', access_type, '_protocols'), horzcat('ASVspoof2019.', access_type, '.cm.train.trn.txt'));
devProtocolFile = fullfile(pathToDatabase, horzcat('ASVspoof2019_', access_type, '_protocols'), horzcat('ASVspoof2019.', access_type, '.cm.dev.trl.txt'));

%% Feature extraction for training data

% read train protocol
fileID = fopen(trainProtocolFile);
protocol = textscan(fileID, '%s%s%s%s%s');
fclose(fileID);

% get file and label lists
filename = protocol{2};
sys_id = protocol{4};
key = protocol{5};

% get indices of genuine and spoof files
% bonafideIdx = find(strcmp(key,'bonafide'));
% spoofIdx = find(strcmp(key,'spoof'));

disp('Extracting features for training data...');
dset_name = 'train';
cache_file_name = fullfile(['cache_' dset_name '_' access_type '_cqcc.mat']);

parfor i=1:length(filename)
    % extract CQCC feature
    filePath = fullfile(pathToDatabase,['ASVspoof2019_' access_type '_train/flac'],[filename{i} '.flac']);
    [x,fs] = audioread(filePath);

    % padding to make data-length = 64000(4s)
    x_len = size(x,1);
    max_len = 64000;
    if x_len >= max_len
        x = x(1:max_len);
    else % need to pad
        num_repeats = floor(max_len / x_len)+1;
        x_repeat = repmat(x,num_repeats,1);
        x= x_repeat(1:max_len);
    end

    if strcmp(feature_type,'LFCC')
        [stat,delta,double_delta] = extract_lfcc(x,fs,20,512,20);
        trainFeatureCell{i} = [stat delta double_delta]';
    elseif strcmp(feature_type,'CQCC')
        trainFeatureCell{i} = cqcc(x, fs, 96, fs/2, fs/2^10, 16, 29, 'ZsdD');
    end

    % convert sys_id
    switch sys_id{i}
        case '-'
            sys_id{i} = 0;
        case 'SS_1'
            sys_id{i} = 1;
        case 'SS_2'
            sys_id{i} = 2;
        case 'SS_4'
            sys_id{i} = 3;
        case 'US_1'
            sys_id{i} = 4;
        case 'VC_1'
            sys_id{i} = 5;
        case 'VC_4'
            sys_id{i} = 6;
        %For PA:
        case 'AA'
            sys_id{i} = 7;
        case 'AB'
            sys_id{i} = 8;
        case 'AC'
            sys_id{i} = 9;
        case 'BA'
            sys_id{i} = 10;
        case 'BB'
            sys_id{i} = 11;
        case 'BC'
            sys_id{i} = 12;
        case 'CA'
            sys_id{i} = 13;
        case 'CB'
            sys_id{i} = 14;
        case 'CC'
            sys_id{i} = 15;
        otherwise
            disp('error converting system id!')
    end
end

% convert key
data_y = strcmp(key,'bonafide');
data_x = transpose(trainFeatureCell);
save(cache_file_name, 'filename', 'data_x', 'data_y', 'sys_id', '-v7.3')

%% Feature extraction and scoring of development data

% read dev protocol
fileID = fopen(devProtocolFile);
protocol = textscan(fileID, '%s%s%s%s%s');
fclose(fileID);

% get file and label lists
filename = protocol{2};
sys_id = protocol{4};
key = protocol{5};

% get indices of genuine and spoof files
% bonafideIdx = find(strcmp(key,'bonafide'));
% spoofIdx = find(strcmp(key,'spoof'));

disp('Extracting features for development data...');
dset_name = 'dev';
cache_file_name = fullfile(['cache_' dset_name '_' access_type '_cqcc.mat']);

parfor i=1:length(filename)
    % extract CQCC feature
    filePath = fullfile(pathToDatabase,['ASVspoof2019_' access_type '_dev/flac'],[filename{i} '.flac']);
    [x,fs] = audioread(filePath);

    % padding to make data-length = 64000(4s)
    x_len = size(x,1);
    max_len = 64000;
    if x_len >= max_len
        x = x(1:max_len);
    else % need to pad
        num_repeats = floor(max_len / x_len) + 1;
        x_repeat = repmat(x,num_repeats,1);
        x= x_repeat(1:max_len);
    end

    if strcmp(feature_type,'LFCC')
        [stat,delta,double_delta] = extract_lfcc(x,fs,20,512,20);
        devFeatureCell{i} = [stat delta double_delta]';
    elseif strcmp(feature_type,'CQCC')
        devFeatureCell{i} = cqcc(x, fs, 96, fs/2, fs/2^10, 16, 29, 'ZsdD');
    end

    % convert sys_id
    switch sys_id{i}
        case '-'
            sys_id{i} = 0;
        case 'SS_1'
            sys_id{i} = 1;
        case 'SS_2'
            sys_id{i} = 2;
        case 'SS_4'
            sys_id{i} = 3;
        case 'US_1'
            sys_id{i} = 4;
        case 'VC_1'
            sys_id{i} = 5;
        case 'VC_4'
            sys_id{i} = 6;
         %For PA:
        case 'AA'
            sys_id{i} = 7;
        case 'AB'
            sys_id{i} = 8;
        case 'AC'
            sys_id{i} = 9;
        case 'BA'
            sys_id{i} = 10;
        case 'BB'
            sys_id{i} = 11;
        case 'BC'
            sys_id{i} = 12;
        case 'CA'
            sys_id{i} = 13;
        case 'CB'
            sys_id{i} = 14;
        case 'CC'
            sys_id{i} = 15;
        otherwise
            disp('error converting system id!')
    end
end
% convert key
data_y = strcmp(key,'bonafide');
data_x = transpose(devFeatureCell);
save(cache_file_name, 'filename', 'data_x', 'data_y', 'sys_id', '-v7.3')

author: nesl


text:
function [M,Minit,output] = cp_apr(X, R, varargin)
%CP_APR Compute nonnegative CP with alternating Poisson regression.
%
%   M = CP_APR(X, R) computes an estimate of the best rank-R
%   CP model of a tensor X using an alternating Poisson regression.
%   The input X can be a tensor, sptensor, ktensor, or ttensor. The
%   result P is a ktensor.
%
%   M = CP_APR(X, R, 'param', value, ...) specifies optional parameters and
%   values. Valid parameters and their default values are:
%      'tol' - Tolerance on the inner KKT violation {1.0e-4}
%      'maxiters' - Maximum number of iterations {1000}
%      'maxinneriters' = Maximum number of inner iterations {10}
%      'init' - Initial guess [{'random'}|ktensor]
%      'epsilon' - parameter to avoid divide by zero {100*eps}
%      'kappatol' - tolerance on complementary slackness {100*eps}
%      'kappa' - offset to fix complementary slackness {10*eps}
%      'printitn' - Print every n outer iterations; 0 for no printing {1}
%      'printinneritn' - Print every n inner iterations {0}
%
%   [M,M0] = CP_APR(...) also returns the initial guess.
%
%   [M,M0,out] = CP_APR(...) also returns additional output.
%      out.kktViolations - maximum kkt violation per iteration
%      out.nInnerIters   - number of inner iterations per iteration
%      out.nViolations   - number of factor matrices needing complementary
%                          slackness adjustment per iteration
%      out.nTotalIters   - total number of inner iterations
%
%   REFERENCE: E. C. Chi and T. G. Kolda. On Tensors, Sparsity, and
%   Nonnegative Factorizations, arXiv:1112.2414 [math.NA], December 2011,
%   URL: http://arxiv.org/abs/1112.2414. Submitted for publication.
%
%   See also CP_ALS, KTENSOR, TENSOR, SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

%% Extract dimensions of X and number of dimensions of X.
N = ndims(X);

%% Set algorithm parameters from input or by using defaults.
params = inputParser;
params.addParamValue('epsilon',1e-10,@isscalar);
params.addParamValue('tol',1e-4,@isscalar);
params.addParamValue('maxiters',1000,@(x) isscalar(x) & x > 0);
params.addParamValue('init','random',@(x) (isa(x,'ktensor') || ismember(x,{'random'})));
params.addParamValue('printitn',1,@isscalar);
params.addParamValue('kappa',1e-2,@isscalar);
params.addParamValue('kappatol',1e-10,@isscalar);
params.addParamValue('maxinneriters',10,@isscalar);
params.addParamValue('printinneritn',0,@isscalar);
params.parse(varargin{:});

%% Copy from params object.
epsilon = params.Results.epsilon;
tol = params.Results.tol;
maxOuterIters = params.Results.maxiters;
Minit = params.Results.init;
kappa = params.Results.kappa;
kappaTol = params.Results.kappatol;
maxInnerIters = params.Results.maxinneriters;
printOuterItn = params.Results.printitn;
printInnerItn = params.Results.printinneritn;
kktViolations = -ones(maxOuterIters,1);
nInnerIters = zeros(maxOuterIters,1);

%% Set up and error checking on initial guess for U.
if isa(Minit,'ktensor')
    if ndims(Minit) ~= N
        error('Initial guess does not have the right number of dimensions');
    end

    if ncomponents(Minit) ~= R
        error('Initial guess does not have the right number of components');
    end

    for n = 1:N
        if size(Minit,n) ~= size(X,n)
            error('Dimension %d of the initial guess is the wrong size',n);
        end
    end
elseif strcmp(Minit,'random')
    F = cell(N,1);
    for n = 1:N
        F{n} = rand(size(X,n),R);
    end
    Minit = ktensor(F);
else
    error('The selected initialization method is not supported');
end

%% Set up for iterations - initializing M and Phi.
M = normalize(Minit,[],1);
Phi = cell(N,1);
kktModeViolations = zeros(N,1);

if printOuterItn > 0
  fprintf('\nCP_APR:\n');
end

nViolations = zeros(maxOuterIters,1);

%% Main Loop: Iterate until convergence.
for iter = 1:maxOuterIters

    isConverged = true;
    for n = 1:N

        % Make adjustments to entries of M{n} that are violating
        % complementary slackness conditions.
        if (iter > 1)
            V = (Phi{n} > 1) & (M{n} < kappaTol);
            if any(V(:))
                nViolations(iter) = nViolations(iter) + 1;
                M{n}(V>0) = M{n}(V>0) + kappa;
            end
        end

        % Shift the weight from lambda to mode n
        M = redistribute(M,n);

        % Calculate product of all matrices but the n-th
        % (In sparse case, only calcuates entries corresponding to nonzeros in X.)
        Pi = calculatePi(X, M, R, n, N);

        % Do the multiplicative updates
        for i = 1:maxInnerIters

            % Count the inner iterations
            nInnerIters(iter) = nInnerIters(iter) + 1;

            % Calculate matrix for multiplicative update
            Phi{n} = calculatePhi(X, M, R, n, Pi, epsilon);

            % Check for convergence
            kktModeViolations(n) = max(abs(vec(min(M.U{n},1-Phi{n}))));
            if (kktModeViolations(n) < tol)
                break;
            else
                isConverged = false;
            end

            % Do the multiplicative update
            M{n} = M{n} .* Phi{n};

            % Print status
             if mod(i, printInnerItn)==0
                 fprintf('    Mode = %1d, Inner Iter = %2d, KKT violation = %.6e\n', n, i, kktModeViolations(n));
             end
        end

        % Shift weight from mode n back to lambda
        M = normalize(M,[],1,n);

    end

    kktViolations(iter) = max(kktModeViolations);

    if (mod(iter,printOuterItn)==0)
        fprintf(' Iter %4d: Inner Its = %2d KKT violation = %.6e, nViolations = %2d\n', ...
        iter, nInnerIters(iter), kktViolations(iter), nViolations(iter));
    end

    % Check for convergence
    if (isConverged)
        break;
    end
end

%% Clean up final result
M = normalize(M,'sort',1);

if printOuterItn>0
    normX = norm(X);
    normresidual = sqrt( normX^2 + norm(M)^2 - 2 * innerprod(X,M) );
    fit = 1 - (normresidual / normX); %fraction explained by model
    fprintf('===========================================\n');
    fprintf(' Final log-likelihood = %e \n', tt_loglikelihood(X,M));
    fprintf(' Final least squares fit = %e \n', fit);
    fprintf(' Final KKT violation = %7.7e\n', kktViolations(iter));
    fprintf(' Total inner iterations = %d\n', sum(nInnerIters));
end

output = struct;
output.params = params.Results;
output.kktViolations = kktViolations(1:iter);
output.nInnerIters = nInnerIters(1:iter);
output.nViolations = nViolations(1:iter);
output.nTotalIters = sum(nInnerIters);

end

function Pi = calculatePi(X, M, R, n, N)

if (isa(X,'sptensor'))
    Pi = ones(nnz(X), R);
    for nn = [1:n-1,n+1:N]
        Pi = M{nn}(X.subs(:,nn),:).*Pi;
    end
else
    U = M.U;
    Pi = khatrirao(U{[1:n-1,n+1:N]},'r');
end

end

function Phi = calculatePhi(X, M, R, n, Pi, epsilon)

if (isa(X,'sptensor'))
    Phi = -ones(size(X,n),R);
    xsubs = X.subs(:,n);
    v = sum(M.U{n}(xsubs,:).*Pi,2);
    wvals = X.vals ./ max(v, epsilon);
    for r = 1:R
        Yr = accumarray(xsubs, wvals .* Pi(:,r), [size(X,n) 1]);
        Phi(:,r) = Yr;
    end
else
    Xn = double(tenmat(X,n));
    V = M.U{n}*Pi';
    W = Xn ./ max(V, epsilon);
    Y = W * Pi;
    Phi = Y;
end

end

function y = vec(x)
y = x(:);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [CTUniform3D, CTUniformInfoS] = getUniformizedCTScan(native, scanNum, planC)
%function CTUniform3D = getUniformizedCTScan
%Return the uniformized (constant slice spacing) CT scan.
%CTUniformA is a 3D double matrix.
%JOD, 17 Oct 03.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC','var')
    global planC
end

indexS = planC{end};

CTUniformInfoS = planC{indexS.scan}(scanNum).uniformScanInfo;

scanData = planC{indexS.scan}(scanNum);
scanInfo = scanData.scanInfo(1);
clear scanData;

scanArraySup = planC{indexS.scan}(scanNum).scanArraySuperior;
scanArrayInf = planC{indexS.scan}(scanNum).scanArrayInferior;
uniformScanInfo = planC{indexS.scan}(scanNum).uniformScanInfo;
sliceNumSup  = uniformScanInfo.sliceNumSup;
sliceNumInf  = uniformScanInfo.sliceNumInf;

%creation options:
origOpts = planC{indexS.CERROptions};
% this code uses the superior and inferior matrices as well as the original CT scan to graph the CT image.
if isfield(origOpts,'uniformizedDataType')
    if strcmpi(planC{indexS.CERROptions}.uniformizedDataType,'uint8')
        CTMin = planC{indexS.scan}(scanNum).uniformScanInfo.minCTValue;
        CTMax = planC{indexS.scan}(scanNum).uniformScanInfo.maxCTValue;
        CTScale = 255 / (CTMax - CTMin);
    elseif strcmpi(planC{indexS.CERROptions}.uniformizedDataType,'uint16')
        CTMin = planC{indexS.scan}(scanNum).uniformScanInfo.minCTValue;
        CTMax = planC{indexS.scan}(scanNum).uniformScanInfo.maxCTValue;
        CTScale = 65535 / (CTMax - CTMin);
    end
else %scaling wasn't used to create the arrays.  This is here to retain backwards compatibility.
    CTScale = 1;
    CTMin = 0;
end

emptySup = isempty(scanArraySup);
emptyInf = isempty(scanArrayInf);

if exist('native') & native == 1
    if emptySup & emptyInf
        CTUniform3D = single(planC{indexS.scan}(scanNum).scanArray(:, :, sliceNumSup:sliceNumInf));
    elseif emptySup & ~emptyInf
        infM = single((double(scanArrayInf) / CTScale) + CTMin);
        CTUniform3D = cat(3, single(planC{indexS.scan}(scanNum).scanArray(:, :, sliceNumSup:sliceNumInf)),infM);
    elseif ~emptySup & emptyInf
        supM = single((double(scanArraySup) / CTScale) + CTMin);
        CTUniform3D = cat(3, supM, single(planC{indexS.scan}(scanNum).scanArray(:, :, sliceNumSup:sliceNumInf)));
    else % ~emptySup & ~emptyInf
        infM = single((double(scanArrayInf) / CTScale) + CTMin);
        supM = single((double(scanArraySup) / CTScale) + CTMin);
        CTUniform3D = cat(3, supM, single(planC{indexS.scan}(scanNum).scanArray(:, :, sliceNumSup:sliceNumInf)), ...
            infM);
    end
else
    if emptySup & emptyInf
        CTUniform3D = planC{indexS.scan}(scanNum).scanArray(:, :, sliceNumSup:sliceNumInf);
    elseif emptySup & ~emptyInf
        infM = ((double(scanArrayInf) / CTScale) + CTMin);
        CTUniform3D = cat(3, planC{indexS.scan}(scanNum).scanArray(:, :, sliceNumSup:sliceNumInf), uint16(infM));
    elseif ~emptySup & emptyInf
        supM = ((double(scanArraySup) / CTScale) + CTMin);
        CTUniform3D = cat(3, uint16(supM), planC{indexS.scan}(scanNum).scanArray(:, :, sliceNumSup:sliceNumInf));
    else % ~emptySup & ~emptyInf
        infM = (double(scanArrayInf) / CTScale) + CTMin;
        supM = (double(scanArraySup) / CTScale) + CTMin;
        clear scanArrayInf scanArraySup;
        CTUniform3D = cat(3, uint16(supM), planC{indexS.scan}(scanNum).scanArray(:, :, sliceNumSup:sliceNumInf), ...
            uint16(infM));
    end
end

CTUniformInfoS.size = size(CTUniform3D);


text:
function value = p23_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P23_F evaluates the integrand for problem 23.
%
%  Dimension:
%
%    DIM_NUM is arbitrary.
%
%  Region:
%
%    The interior of the unit simplex, for which all X's are nonnegative,
%    and sum ( X(1:N) ) <= 1.
%
%  Integral Parameters:
%
%    C defaults to 1.0.
%    Call P23_R8 to get or set this value.
%
%    E(1:DIM_NUM) defaults to (/ 2, 2, ..., 2 /).
%    Call P23_I4VEC to get or set this value.
%
%  Integrand:
%
%    F(X) = C * X1^E1 * X2^E2 * ... * Xn^En
%
%    C is real, all exponents E are nonnegative integers.
%
%  Exact Integral:
%
%    C * Gamma(E1+1) * Gamma(E2+1) * ... * Gamma(En+1) / Gamma(E1+E2+...+En+1)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  c = 0.0;
  c = p23_r8 ( 'G', 'C', c );

  e = [];
  e = p23_i4vec ( 'G', 'E', dim_num, e );

  value(1:point_num) = c;

  for point = 1 : point_num

    for dim = 1 : dim_num
      value(point) = value * x(dim,point)^e(dim);
    end

  end

  p23_i4 ( 'I', '#', point_num );

  return
end


text:
function this = rmfield(this, vals);

    this.EEG = rmfield(this.EEG, vals);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function a = transpose(a)
%TRANSPOSE    Implements  a.'  for gradients
%
%   c = a.'
%

% written  10/16/98     S.M. Rump
% modified 11/03/03     S.M. Rump  improved performance
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  [m n] = size(a.x);
  a.x = a.x.';
  if m*n~=1
    index = reshape( 1:(m*n) , m , n )';
    a.dx = a.dx( index , : );
  end

repo: OpenMAS


text:
function str = printor(strs, quote)

% PRINTOR prints a single or mutiple strings as "x1, x2, x3 or x4". If there is
% only one string, that string is returned without additional formatting.
%
% See also PRINTAND

if nargin<2
  quote = false;
end

if quote
  if numel(strs)>2
    str = sprintf('''%s'', ', strs{1:(end-2)});
    str = sprintf('%s''%s'' or ''%s''', str, strs{end-1}, strs{end});
  elseif numel(strs)==2
    str = sprintf('''%s'' or ''%s''', strs{end-1}, strs{end});
  else
    str = strs{1};
  end
else
  if numel(strs)>1
    str = sprintf('%s, ', strs{1:(end-2)});
    str = sprintf('%s%s or %s', str, strs{end-1}, strs{end});
  else
    str = strs{1};
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function txt = event_datacursor_txt(~,event_obj)
    % show events off a map, filtering for x,y,z on main catalog
    % capable of showing multiple events if they
    % event_obj   Object containing event data structure
    % output_txt  Data cursor text
    %
    % currently selected event(s) can be found on desktop as datacursor_catalog

    %event_obj.Target : handle of object data cursor is referencing
    %event_obj.Position : array specifying x,y,z

    % redefine data cursor Updatefcn
    % set(dcm_obj, 'UpdateFcn',@myupdatefcn)
    %ZG=ZmapGlobal.Data;
    event_obj.Target
    event_obj.Position

    switch lower(event_obj.Target.DisplayName)
        case {'stations'}
            txt=parse_stations(event_obj);
        case {'volcanoes'}
            txt=parse_volcanoes(event_obj);
        case {'main faultline'}
            txt=parse_faults(event_obj);
        case {'plate boundaries'}
            txt=parse_plates(event_obj);

        otherwise

            txt=parse_quakes(event_obj);
    end

end

%% These parsing routines should probably be assigned to the features themselves.

function txt = parse_quakes(event_obj)
    disp(event_obj.Target)
    latTol = 0.001;
    lonTol = 0.001;
    depthTol = 0.1;
    OUTPUT_NAME='datacursor_catalog';

    assignin('base','clickedevent',event_obj)
    evpos=event_obj.Position;
    ZG=ZmapGlobal.Data;
    globalcatalog = ZG.primeCatalog;
      evNum=find(...
         abs(evpos(2) - globalcatalog.Y) < latTol &...
         abs(evpos(1) - globalcatalog.X) < lonTol &...
         abs(evpos(3) - globalcatalog.Z < depthTol));

     minicat=globalcatalog.subset(evNum);
     minicat.Name=sprintf('datacursor:%s', char(datetime(),'uuuuMMdd''T''hhmmss'));
     assignin('base',OUTPUT_NAME,minicat);

     txt={};
     if minicat.Count > 1
         txt={sprintf('%d events within position tolerance',minicat.Count)};
     end

     if minicat.Count > 3
         try
             for i =1:minicat.Count
                 if i > 5 && i<minicat.Count
                     continue
                 end
                 if i== 5 && i < minicat.Count
                     txt = [txt, {sprintf('  ... skipping %d events ...',minicat.Count-5)}];
                     continue
                 end

                 lattxt=lat_text(event_obj,minicat,i);
                 lontxt=lon_text(event_obj,minicat,i);

                 mt = minicat.MagnitudeType{i};
                 if isempty(mt)
                     mt='magnitude';
                 end
                 txt=[txt,{[sprintf('[# %5d] ',evNum(i)),...
                     sprintf(' %s',char(minicat.Date(i),'uuuu-MM-dd HH:mm:ss ; ')),...
                     sprintf('( %8s , %8s ) at %5.2f km ; ',lattxt,lontxt,evpos(3)),...
                     sprintf('%s %.1f',mt, minicat.Magnitude(i))]};
                     ];
             end

         catch ME
             warning(ME.message)
         end

     else
         try
             for i =1:minicat.Count
                 lattxt=lat_text(event_obj,minicat,i);
                 lontxt=lon_text(event_obj,minicat,i);
                 mt = minicat.MagnitudeType{i};
                 if isempty(mt)
                     mt='magnitude';
                 end
                 txt=[txt,{sprintf('- Event # %d -',evNum(i)),...
                     sprintf('%s',char(minicat.Date(i),'uuuu-MM-dd HH:mm:ss')),...
                     sprintf('( %s , %s )',lattxt,lontxt),...
                     sprintf(' %.2f km depth',minicat.Depth(i)),...
                     sprintf('%s : %.1f',mt, minicat.Magnitude(i))};
                     ];
                 if i < minicat.Count
                     txt=[txt,{''}];
                 end
             end

         catch ME
             warning(ME.message)
         end
     end
end

function txt = parse_stations(event_obj)
    % get table information

    latTol = 0.01;
    lonTol = 0.01;

    ZG=ZmapGlobal.Data;
    f=ZG.features('stations');
    evpos=event_obj.Position;
    idx=find(...
        abs(evpos(2) - f.Latitude) < latTol &...
        abs(evpos(1) - f.Longitude) < lonTol);
    txt={sprintf('%d station(s):',numel(idx))};
    for i=1:numel(idx)
        lattxt=lat_text(event_obj,f,i);
        lontxt=lon_text(event_obj,f,i);

        txt=[txt,f.Names(idx(i)),...
            sprintf('( %s , %s )',lattxt,lontxt),...
            sprintf(' %.2f km elev',-f.Depth(i))];
    end
end

function txt = parse_faults(event_obj)
    txt='fault'
    %f=ZG.features('faults');
end
function txt = parse_plates(event_obj)
    txt='plate boundary'
    %f=ZG.features('plates');
end

function txt = parse_volcanoes(event_obj)    % get table information
    latTol = 0.001;
    lonTol = 0.001;

    ZG=ZmapGlobal.Data;
    f=ZG.features('volcanoes');
    evpos=event_obj.Position;
    idx=find(...
        abs(evpos(2) - f.Latitude) < latTol &...
        abs(evpos(1) - f.Longitude) < lonTol)
    txt={sprintf('%d volcanic feature(s):',numel(idx))};
    for i=1:numel(idx)
        lattxt=lat_text(event_obj,f,i);
        lontxt=lon_text(event_obj,f,i);

        txt=[txt,f.Names(idx(i)),...
            sprintf('( %s , %s )',lattxt,lontxt),...
            sprintf(' %.2f km elev',-f.Depth(i))];
    end

end

function lattxt=lat_text(event_obj,obj,idx)
    if ~exist('idx','var')
        idx=1;
    end
    evpos=event_obj.Position;
    if evpos(1)>0
        lattxt=[num2str(obj.Latitude(idx)) ' N'];
    else
        lattxt=[num2str(obj.Latitude(idx)) ' S'];
    end
end

function lontxt=lon_text(event_obj,obj,idx)
    if ~exist('idx','var')
        idx=1;
    end
    evpos=event_obj.Position;
    if evpos(1)>0
        lontxt=[num2str(obj.Longitude(idx)) ' E'];
    else
        lontxt=[num2str(obj.Longitude(idx)) ' W'];
    end
end


text:
classdef BatchNorm < dagnn.ElementWise
  properties
    numChannels
    epsilon = 1e-5
    opts = {'NoCuDNN'} % ours seems slightly faster
    noise_param_idx
    noise_cache_size
    bnorm_moment_type_trn
    bnorm_moment_type_tst
  end

  properties (Transient)
    moments
  end

  methods
    function outputs = forward(obj, inputs, params)

	  if strcmp(obj.net.mode, 'test')
        outputs{1} = vl_nnbnorm(inputs{1}, params{1}, params{2}, ...
                                'moments', params{3}, ...
                                'epsilon', obj.epsilon, ...
                                obj.opts{:}) ;
      else

      	% Guosheng: changed for RefineNet
      	bnorm_moment_type=obj.bnorm_moment_type_trn;

        if strcmp(bnorm_moment_type, 'batch')

	        [outputs{1},obj.moments] = ...
	            vl_nnbnorm(inputs{1}, params{1}, params{2}, ...
	                       'epsilon', obj.epsilon, ...
	                       obj.opts{:}) ;

		elseif strcmp(bnorm_moment_type, 'global')
			outputs{1} = vl_nnbnorm(inputs{1}, params{1}, params{2}, ...
                                'moments', params{3}, ...
                                'epsilon', obj.epsilon, ...
                                obj.opts{:}) ;
		else
            error('not support!');
        end

      end

    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)

		% Guosheng: changed for RefineNet

    	bnorm_moment_type=obj.bnorm_moment_type_trn;

    	if strcmp(bnorm_moment_type, 'batch')

	      [derInputs{1}, derParams{1}, derParams{2}, derParams{3}] = ...
	        vl_nnbnorm(inputs{1}, params{1}, params{2}, derOutputs{1}, ...
	                   'epsilon', obj.epsilon, ...
	                   'moments', obj.moments, ...
	                   obj.opts{:}) ;

  		elseif strcmp(bnorm_moment_type, 'global')

  			[derInputs{1}, derParams{1}, derParams{2}, derParams{3}] = ...
	        	vl_nnbnorm(inputs{1}, params{1}, params{2}, derOutputs{1}, ...
	        		'moments', params{3}, ...
	                   'epsilon', obj.epsilon, ...
	                   obj.opts{:}) ;

	    else
            error('not support!');
  		end

  		obj.moments = [] ;

      % multiply the moments update by the number of images in the batch
      % this is required to make the update additive for subbatches
      % and will eventually be normalized away
      derParams{3} = derParams{3} * size(inputs{1},4) ;
    end

    % ---------------------------------------------------------------------
    function obj = BatchNorm(varargin)
      obj.load(varargin{:}) ;
    end

    function params = initParams(obj)
      params{1} = ones(obj.numChannels,1,'single') ;
      params{2} = zeros(obj.numChannels,1,'single') ;
      params{3} = zeros(obj.numChannels,2,'single') ;
    end

    function attach(obj, net, index)
      attach@dagnn.ElementWise(obj, net, index) ;
      p = net.getParamIndex(net.layers(index).params{3}) ;
      net.params(p).trainMethod = 'average' ;
      net.params(p).learningRate = 0.1 ;
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function value = i4_btest ( i4, pos )

%*****************************************************************************80
%
%% I4_BTEST returns TRUE if the POS-th bit of an I4 is 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Military Standard 1753,
%    FORTRAN, DoD Supplement To American National Standard X3.9-1978,
%    9 November 1978.
%
%  Parameters:
%
%    Input, integer I4, the integer to be tested.
%
%    Input, integer POS, the bit position, between 0 and 31.
%
%    Output, logical VALUE, is TRUE if the POS-th bit of I4 is 1.
%
  i4_huge = 2147483647;

  if ( pos < 0 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_BTEST - Fatal error!\n' );
    fprintf ( 1, '  POS < 0.\n' );
    error ( 'I4_BTEST - Fatal error!' );

  elseif ( pos < 31 )

    if ( 0 <= i4 )
      j = floor ( i4 );
    else
      j = floor ( ( i4_huge + i4 ) + 1 );
    end

    for k = 1 : pos
      j = floor ( j / 2 );
    end

    if ( mod ( j, 2 ) == 0 )
      value = 0;
    else
      value = 1;
    end

  elseif ( pos == 31 )

    if ( i4 < 0 )
      value = 1;
    else
      value = 0;
    end

  elseif ( 31 < pos )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_BTEST - Fatal error!\n' );
    fprintf ( 1, '  31 < POS.\n' );
    error ( 'I4_BTEST - Fatal error!' );

  end

  return
end


text:
function bz_meta_extracted = bz_database_extract_meta(bz_datasets,j)
% Submit dataset to the Buzsakilab metadata-database
% v1.1
%
% INPUT
% bz_datasets
% structure with subfields named after the mySQL columns. The subfields has
% to have the right format for submitting to the database.
%
%
% By Peter Petersen
% petersen.peter@gmail.com

bz_meta_extracted = bz_datasets;
if ~exist('j')
    j = size(bz_meta_extracted,2);
end
Session = bz_datasets(j).SessionPath;
fprintf(['\n ', num2str(j), '. Animal ' bz_datasets(j).Animal, ', session ', bz_datasets(j).Session, ': '])

% LFP
if exist([Session,'.eeg'])
    bz_meta_extracted(j).LFP = 'Yes';
elseif exist([Session,'.lfp'])
    bz_meta_extracted(j).LFP = 'Yes';
else
    bz_meta_extracted(j).LFP = 'No';
end

% Raw data
if exist(([Session,'.dat']))
    bz_meta_extracted(j).RawData = 'Raw data';
elseif exist(([Session,'.fil']))
    bz_meta_extracted(j).RawData = 'High-pass filtered';
else
    bz_meta_extracted(j).RawData = 'No';
end
if ~strcmp(bz_meta_extracted(j).RawData,'No')
    fprintf(['Data format: ' char(bz_meta_extracted(j).RawData) '. '])
end

% Sampling rate and channel count
if exist(([Session,'.xml']))
    xml = LoadXml(([Session,'.xml']));
    bz_meta_extracted(j).SampleRate = xml.SampleRate;
    bz_meta_extracted(j).lfpSampleRate = xml.lfpSampleRate;
    bz_meta_extracted(j).nChannels = xml.nChannels;
    bz_meta_extracted(j).nElecGps = xml.nElecGps;
    bz_meta_extracted(j).Unitsformat = 'NA';
end

% Number of Cells
bz_meta_extracted(j).Cells = 0;
for i = 1:bz_meta_extracted(j).nElecGps
    if exist(([Session,'.clu.' num2str(i)]))
        cluster_index = load(([Session,'.clu.' num2str(i)]));
        bz_meta_extracted(j).Cells = bz_meta_extracted(j).Cells + sum(unique(cluster_index(2:end))>1);
        bz_meta_extracted(j).Unitsformat = 'Neurosuite';
    end
end
if bz_meta_extracted(j).Cells>0
    fprintf([num2str(bz_meta_extracted(j).Cells) ' units detected. '])
else
    fprintf('No units detected. ')
end

% Duration
if exist(([Session,'.eeg']))
    temp_ = dir(([Session,'.eeg']));
    bz_meta_extracted(j).Duration = temp_.bytes/bz_meta_extracted(j).nChannels/bz_meta_extracted(j).lfpSampleRate/2;
elseif exist(([Session,'.lfp']))
    temp_ = dir(([Session,'.lfp']));
    bz_meta_extracted(j).Duration = temp_.bytes/bz_meta_extracted(j).nChannels/bz_meta_extracted(j).lfpSampleRate/2;
else
    bz_meta_extracted(j).Duration = 0;
end

% Sleep scored
if exist(([Session,'.states.REM'])) | exist(([Session,'.states.SWS'])) | exist(([Session,'.states.Wake']))
    bz_meta_extracted(j).SleepScored = 'Yes';
elseif exist(([Session,'-states.mat']))
    bz_meta_extracted(j).SleepScored = 'Yes';
elseif exist(([Session,'.sts.REM'])) | exist(([Session,'.sts.SWS'])) | exist(([Session,'.sts.Wake']))
    bz_meta_extracted(j).SleepScored = 'Yes';
elseif exist(([Session,'_SleepScore.mat']))
    bz_meta_extracted(j).SleepScored = 'Yes';
elseif exist(([Session,'.SleepState.states.mat']))
    bz_meta_extracted(j).SleepScored = 'Yes';
else
    bz_meta_extracted(j).SleepScored = 'No';
end
if strcmp(bz_meta_extracted(j).SleepScored, 'Yes')
    fprintf('SleepScoring detected. ')
else
    fprintf('No sleepScoring. ')
end

% Ripples
if exist(([Session,'.ripplesALL.event.mat']))
    fprintf('Ripple analysis detected. ')
    bz_meta_extracted(j).RipplesDetected = 'Yes';
else
    fprintf('No ripple analysis detected. ')
    bz_meta_extracted(j).RipplesDetected = 'No';
end

% Tracking
if exist(([Session,'.whl']))
    fprintf('whl tracking files detected. ')
    bz_meta_extracted(j).Tracking = 'Camera';
    bz_meta_extracted(j).TrackingFormat = 'whl';
    bz_meta_extracted(j).TrackingFrameRate = 39.06;
else
    bz_meta_extracted(j).Tracking = 'No';
    bz_meta_extracted(j).TrackingFormat = nan;
end

% Buzcode
if exist(([Session,'.sessionInfo.mat']))
    fprintf('Buzcode sessionInfo file detected. ')
    sessionInfo = load([Session,'.sessionInfo.mat']);
    if isfield(sessionInfo.sessionInfo,'region') && ~isempty(sessionInfo.sessionInfo.region)
        bz_meta_extracted(j).BrainRegion = strjoin(unique(sessionInfo.sessionInfo.region),',');
    end
    if isfield(sessionInfo.sessionInfo,'depth') && ~isempty(sessionInfo.sessionInfo.region)
        bz_meta_extracted(j).Depth = char(sessionInfo.sessionInfo.depth);
    end
    if isfield(sessionInfo.sessionInfo,'Date')
        bz_meta_extracted(j).Date = sessionInfo.sessionInfo.Date;
    end
    bz_meta_extracted(j).FileFormat = 'Buzcode';
    clear sessionInfo
end

if exist(([Session,'.behavior.mat']))
    fprintf('Buzcode behavior file detected. ')
    behavior = load([Session,'.behavior.mat']);
    bz_meta_extracted(j).TrackingFormat = 'buzcode';
    if isfield(behavior.behavior,'acquisitionsystem') && ~isempty(behavior.behavior.acquisitionsystem)
        bz_meta_extracted(j).Tracking = behavior.behavior.acquisitionsystem;
    end
    if isfield(behavior.behavior,'samplingRate') && ~isempty(behavior.behavior.samplingRate)
        bz_meta_extracted(j).TrackingFrameRate = behavior.behavior.samplingRate;
    end
    if isfield(behavior.behavior,'events') && isfield(behavior.behavior.events,'trials') && ~isempty(behavior.behavior.events.trials)
    bz_meta_extracted(j).BehavioralTrials = size(behavior.behavior.events.trials,2);
    end
    clear behavior
end

% Brendon's standard
if exist(([Session,'_BasicMetaData.mat']))
    fprintf('_BasicMetaData file detected (Brendons format). ')
    BasicMetaData = load([Session,'_BasicMetaData.mat']);
    if isfield(BasicMetaData,'CortexRegion')
        bz_meta_extracted(j).BrainRegion = strjoin(unique(BasicMetaData.CortexRegion),',');
    end
    if isfield(BasicMetaData,'bmd') && isfield(BasicMetaData.bmd,'Par') && isfield(BasicMetaData.bmd.Par,'Date')
        bz_meta_extracted(j).Date = BasicMetaData.bmd.Par.Date;
    end
    bz_meta_extracted(j).FileFormat = 'Brendon_metadata';
    clear bmd CortexRegion
end

% Neurodata Without Borders

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [RMS_position, RMS_orientation, NEES_pose, NEES_orientation] = RocEKF_plot_rms_nees( estimation_results, data, do_vis )
% plot rms and nees for R-EKF
N  = size(estimation_results, 2);

T = 2:N;
RMS_position=[];
RMS_orientation=[];
for i = T
    position = estimation_results{i}.position;
    ap = data.poses.position(:,i);

    RMS_position = [RMS_position norm(position-ap)];
    RMS_orientation=  [RMS_orientation norm(so3_log(estimation_results{i}.orientation*(data.poses.orientation(3*i-2:3*i,1:3))'))];

end

fprintf( 'mean(RMS:Position) = %f\n', mean(RMS_position) );
fprintf( 'mean(RMS:Orientation) = %f\n', mean(RMS_orientation) );

if do_vis == 1
    figure;

    subplot(2,2,1);plot(T, RMS_position);
    title('RMS:position(meter)');xlim([0,N]);
    subplot(2,2,2); plot(T,RMS_orientation);
    title('RMS:orientation(radius)');xlim([0,N]);
end

NEES_pose=[];
NEES_orientation=[];
for i = T
    position = estimation_results{i}.position;
    ap = data.poses.position(:,i);

    dw=so3_log(estimation_results{i}.orientation*(data.poses.orientation(3*i-2:3*i,1:3))');
    dv=(jaco_r(-dw))\( (estimation_results{i}.position-estimation_results{i}.orientation*data.poses.orientation(3*i-2:3*i,1:3)'*data.poses.position(:,i)));

     cov_o=estimation_results{i}.cov(1:3,1:3);
     cov_pose=estimation_results{i}.cov(1:6,1:6);

     invcov_o=eye(3)/cov_o;
     invcov_pose=eye(6)/cov_pose;
    dP=[dw;dv];

    NEES_orientation = [NEES_orientation dw'*invcov_o*dw/3];
    NEES_pose=  [NEES_pose dP'*invcov_pose*dP/6];

end

fprintf( 'mean(NEES:Position) = %f\n', mean(NEES_pose(2:end)) );
fprintf( 'mean(NEES:Orientation) = %f\n', mean(NEES_orientation(2:end)) );

if do_vis == 1
    subplot(2,2,3);
    plot(T,NEES_orientation);
    title('NEES:orientation');xlim([0,N]);
    subplot(2,2,4);
    plot(T,NEES_pose);
    title('NEES:pose');xlim([0,N]);
end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
%FINDESSENTIALMAT  Calculates an essential matrix from the corresponding points in two images
%
%     E = cv.findEssentialMat(points1, points2)
%     [E, mask] = cv.findEssentialMat(...)
%     [...] = cv.findEssentialMat(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __points1__ Cell array of N (N>=5) 2D points from the first image, or numeric array
%   Nx2/Nx1x2/1xNx2. The point coordinates should be floating-point (single or
%   double precision).
% * __points2__ Cell array or numeric array of the second image points of the
%   same size and format as `points1`.
%
% ## Output
% * __E__ Essential matrix, 3x3.
% * __mask__ Output vector of N elements, every element of which is set to 0
%   for outliers and to 1 for the other points (inliers). The array is
%   computed only in the RANSAC and LMedS robust methods.
%
% ## Options
% * __CameraMatrix__ Camera matrix `K = [fx 0 cx; 0 fy cy; 0 0 1]`. Note that
%   this function assumes that `points1` and `points2` are feature points from
%   cameras with the same camera matrix. default `eye(3)`.
% * __Method__ Method for computing an essential matrix. One of:
%   * __Ransac__ for the RANSAC algorithm. (default)
%   * __LMedS__ for the LMedS algorithm.
% * __Confidence__ Parameter used for the RANSAC or LMedS methods only. It
%   specifies a desirable level of confidence (probability) that the estimated
%   matrix is correct. In the range 0..1 exclusive. default 0.999
% * __Threshold__ Parameter used for RANSAC. It is the maximum distance from a
%   point to an epipolar line in pixels, beyond which the point is considered
%   an outlier and is not used for computing the final essential matrix. It
%   can be set to something like 1-3, depending on the accuracy of the point
%   localization, image resolution, and the image noise. default 1.0
%
% This function estimates essential matrix based on the five-point algorithm
% solver in [Nister03]. [SteweniusCFS] is also a related. The epipolar
% geometry is described by the following equation:
%
%     [p2;1]' * inv(K)' * E * inv(K) * [p1;1] = 0
%
% where `E` is an essential matrix, `p1` and `p2` are corresponding points in
% the first and the second images, respectively. The result of this function
% may be passed further to cv.decomposeEssentialMat or cv.recoverPose to
% recover the relative pose between cameras.
%
% `K` is the camera matrix with focal length `fx` and `fy` and principal point
% `[cx,cy]`:
%
%     K = [fx  0 cx;
%          0  fy cy;
%          0   0  1]
%
% ## Example
% Estimation of essential matrix using the RANSAC algorithm:
%
%     % initialize the points here
%     points1 = {[1,1],[3,1],[5,1],...}
%     points2 = {[2,3],[4,3],[6,3],...}
%     % estimate essential matrix
%     [E, mask] = cv.findEssentialMat(points1, points2, 'Method','Ransac');
%
% ## References
% [Nister03]:
% > David Nister. "An efficient solution to the five-point relative pose
% > problem". Pattern Analysis and Machine Intelligence, IEEE Transactions on,
% > 26(6):756-770, 2004.
%
% [SteweniusCFS]:
% > Henrik Stewenius. "Calibrated fivepoint solver".
%
% See also: cv.findFundamentalMat, cv.recoverPose, estimateEssentialMatrix
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% -------------------------------------------------------------------------
function net = cnn_init_cafferef(net,opts)
% -------------------------------------------------------------------------

drop6p = find(cellfun(@(a) strcmp(a.name, 'dropout6'), net.layers)==1);
drop7p = find(cellfun(@(a) strcmp(a.name, 'dropout7'), net.layers)==1);

if ~isempty(drop6p)
  assert(~isempty(drop7p));
  net.layers{drop6p}.rate = opts.DropOutRate;
  net.layers{drop7p}.rate = opts.DropOutRate;
else
  relu6p = find(cellfun(@(a) strcmp(a.name, 'relu6'), net.layers)==1);
  relu7p = find(cellfun(@(a) strcmp(a.name, 'relu7'), net.layers)==1);

  drop6 = struct('type','dropout','rate', opts.DropOutRate,'name','dropout6') ;
  drop7 = struct('type','dropout','rate', opts.DropOutRate,'name','dropout7') ;
  net.layers = [net.layers(1:relu6p) drop6 net.layers(relu6p+1:relu7p) drop7 net.layers(relu7p+1:end)];
end

% replace fc8
fc8l = cellfun(@(a) strcmp(a.name, 'fc8'), net.layers)==1;

nCls = opts.nCls ;
% nCls = 101;
sizeW = size(net.layers{fc8l}.weights{1});

if sizeW(4)~=nCls
  net.layers{fc8l}.weights = {zeros(sizeW(1),sizeW(2),sizeW(3),nCls,'single'), ...
    zeros(1, nCls, 'single')};
end

% change loss
% net.layers(end) = [];
net.layers{end} = struct('name','loss', 'type','softmaxloss') ;

% convert to dagnn
net = dagnn.DagNN.fromSimpleNN(net, 'canonicalNames', true) ;

poolLyr1 = find(arrayfun(@(a) strcmp(a.name, opts.pool1Layer), net.layers)==1);
assert(~isempty(poolLyr1));
% configure appr-rank-pool
switch opts.pool1Type
  case 'arpool'
    if strcmp(opts.pool1Layer,'conv1')
      net.addLayer('arpool',AppRankPooling('scale',1),{net.layers(poolLyr1).inputs{1},'VideoId1'},'DynImgN');
      net.addLayer('l2normalize',L2Normalize('scale',6000,'clip',[-128 128]),...
        'DynImgN','DynImg');
    else
      net.addLayer('arpool',AppRankPooling('scale',0.1),{net.layers(poolLyr1).inputs{1},'VideoId1'},'DynImgN');
      net.addLayer('reluP',dagnn.ReLU(),...
      {'DynImgN'},'DynImg');
    end
    net.setLayerInputs(opts.pool1Layer,{'DynImg'}) ;
  case 'ppool1'
    if strcmp(opts.pool1Layer,'conv1')
      net.addLayer('parampool',LinComb('pad',[1 1 10 1]),...
      {net.layers(poolLyr1).inputs{1},'VideoId1'},'DynImg',{'conv0f','conv0b'});
    else
      net.addLayer('parampool',LinComb('pad',[1 1 10 1]),...
      {net.layers(poolLyr1).inputs{1},'VideoId1'},'DynImgN',{'conv0f','conv0b'});
    net.addLayer('reluP',dagnn.ReLU(),...
      {'DynImgN'},'DynImg');
    end

    net.layers(poolLyr1).inputs{1} = 'DynImg' ;
%     net.params(end-1).value = 0.01 * randn(1,1,10,1,'single');
    net.params(end-1).value = 0.1 * ones(1,1,10,1,'single');
    net.params(end).value = zeros(1,1,'single');

    net.params(end-1).learningRate = 0.1 ;
    net.params(end).learningRate = 0.2 ;
  case 'ppool2'
    if strcmp(opts.pool1Layer,'conv1')
      net.addLayer('parampool',LinComb('pad',[1 1 10 1]),...
      {net.layers(poolLyr1).inputs{1},'VideoId1'},'DynImg',{'conv0f','conv0b'});
    else
      net.addLayer('parampool',LinComb('pad',[1 1 10 1]),...
      {net.layers(poolLyr1).inputs{1},'VideoId1'},'DynImgN',{'conv0f','conv0b'});
    net.addLayer('reluP',dagnn.ReLU(),...
      {'DynImgN'},'DynImg');
    end

    net.layers(poolLyr1).inputs{1} = 'DynImg' ;
%     net.params(end-1).value = 0.01 * randn(1,1,10,1,'single');
    net.params(end-1).value = 0.1 * ones(1,1,10,1,'single');
    net.params(end).value = zeros(1,1,'single');

    net.params(end-1).learningRate = 0.1 ;
    net.params(end).learningRate = 0.2 ;
  case 'none'

  otherwise
    error('Unknown pool type %s', opts.pool1Type) ;
end

% second pool layer (max pooling)
poolLyr2 = find(arrayfun(@(a) strcmp(a.name, opts.pool2Layer), net.layers)==1);
net.addLayer('tempPoolMax',TemporalPooling('method','max'),...
  {net.layers(poolLyr2(1)).inputs{1},'VideoId2'},'tempPoolMax');

net.layers(poolLyr2).inputs{1} = 'tempPoolMax';

% add multi-class error
net.addLayer('errMC',ErrorMultiClass(),{'prediction','label'},'mcerr');

net_ = net.saveobj ;
net = dagnn.DagNN.loadobj(net_) ;

net.removeLayer('loss') ;
net.addLayer('loss', ...
             LossNormalized('loss', 'softmaxlog') ,...
             {'prediction', 'label'}, ...
             'objective') ;

% replace standard matconvnet bnorm with my version
bns = find(arrayfun(@(a) strcmp(class(a.block), 'dagnn.BatchNorm'), net.layers)==1);
for i=1:numel(bns)
  bb = net.layers(bns(i)).block ;
  net.layers(bns(i)).block = BatchNormN('numChannels',bb.numChannels,...
  'epsilon',bb.epsilon,...
  'opts',bb.opts) ;
end

author: hbilen

repo: dynamic-image-nets


text:
function runSimSCFDE()

SP.FFTsize = 512;
SP.CPsize = 20;

SP.SNR = [0:2:30];
SP.numRun = 10^4;

% TS 25.104
pedAchannel = [1 10^(-9.7/20) 10^(-22.8/20)];
pedAchannel = pedAchannel/sqrt(sum(pedAchannel.^2));
vehAchannel = [1 0 10^(-1/20) 0 10^(-9/20) 10^(-10/20) 0 0 0 10^(-15/20) 0 0 0 10^(-20/20)];
vehAchannel = vehAchannel/sqrt(sum(vehAchannel.^2));
idenChannel = 1;

SP.channel = idenChannel;
%SP.channel = pedAchannel;
%SP.channel = vehAchannel;

SP.equalizerType ='ZERO';
%SP.equalizerType ='MMSE';

SER_scfde = scfde(SP);
SER_ofdm = ofdm(SP);
save scfde_awgn

SP.channel = pedAchannel;
SP.equalizerType ='ZERO';
SER_scfde = scfde(SP);
SER_ofdm = ofdm(SP);
save scfde_pedA_zero

SP.channel = pedAchannel;
SP.equalizerType ='MMSE';
SER_scfde = scfde(SP);
SER_ofdm = ofdm(SP);
save scfde_pedA_mmse

SP.channel = vehAchannel;
SP.equalizerType ='ZERO';
SER_scfde = scfde(SP);
SER_ofdm = ofdm(SP);
save scfde_vehA_zero

SP.channel = vehAchannel;
SP.equalizerType ='MMSE';
SER_scfde = scfde(SP);
SER_ofdm = ofdm(SP);
save scfde_vehA_mmse

SP.CPsize = 0;

SP.channel = pedAchannel;
SP.equalizerType ='MMSE';
SER_scfde = scfde(SP);
SER_ofdm = ofdm(SP);
save scfde_pedA_noCP_mmse

SP.channel = vehAchannel;
SP.equalizerType ='MMSE';
SER_scfde = scfde(SP);
SER_ofdm = ofdm(SP);
save scfde_vehA_noCP_mmse

author: Sable