text:
function varargout = cvx_version( varargin )

% CVX_VERSION   Returns version and environment information for CVX.
%
%    When called with no arguments, CVX_VERSION prints out version and
%    platform information that is needed when submitting CVX bug reports.
%
%    This function is also used internally to return useful variables that
%    allows CVX to adjust its settings to the current environment.

global cvx___

args = varargin;
compile = false;
quick = nargout > 0;
if nargin
    if ~ischar( args{1} ),
        quick = true;
    else
        tt = strcmp( args, '-quick' );
        quick = any( tt );
        if quick, args(tt) = []; end
        tt = strcmp( args, '-compile' );
        compile = any( tt );
        if compile, quick = false; args(tt) = []; end
    end
end

if isfield( cvx___, 'loaded' ),

    if quick, return; end
    fs = cvx___.fs;
    mpath = cvx___.where;
    isoctave = cvx___.isoctave;
    msub = cvx___.msub;

else

    % Matlab / Octave flag
    isoctave = exist( 'OCTAVE_VERSION', 'builtin' );

    % File and path separators, MEX extension
    if isoctave,
        comp = octave_config_info('canonical_host_type');
        mext = 'mex';
        ispc = false;
        ismac = false;
        if octave_config_info('mac'),
            msub = 'mac';
            ismac = true;
        elseif octave_config_info('windows'),
            msub = 'win';
            ispc = true;
        elseif octave_config_info('unix') && any(strfind(comp,'linux')),
            msub = 'lin';
        else
            msub = 'unknown';
        end
        if ~isempty( msub ),
            msub = [ 'o_', msub ];
            if strncmp( comp, 'x86_64', 6 ),
                msub = [ msub, '64' ];
            else
                msub = [ msub, '32' ];
            end
        end
    else
        comp = computer;
        ispc = strncmp( comp, 'PC', 2  );
        ismac = strncmp( comp, 'MAC', 3 );
        mext = mexext;
        msub = '';
    end
    if ispc,
        fs = '\';
        fsre = '\\';
        ps = ';';
        cs = false;
    else
        fs = '/';
        fsre = '/';
        ps = ':';
        cs = ~ismac;
    end

    % Install location
    mpath = mfilename('fullpath');
    temp = strfind( mpath, fs );
    mpath = mpath( 1 : temp(end) - 1 );

    % Numeric version
    nver = version;
    nver(nver=='.') = ' ';
    nver = sscanf(nver,'%d');
    nver = nver(1) + 0.01 * ( nver(2) + 0.01 * nver(3) );

    if isoctave || ~usejava('jvm'),
        jver = 0;
    else
        jver = char(java.lang.System.getProperty('java.version'));
        try
            ndxs = strfind( jver, '.' );
            jver = str2double( jver(1:ndxs(2)-1) );
        catch
            jver = 0;
        end
    end

    cvx___.where = mpath;
    cvx___.isoctave = isoctave;
    cvx___.nver = nver;
    cvx___.jver = jver;
    cvx___.comp = comp;
    cvx___.mext = mext;
    cvx___.msub = msub;
    cvx___.fs = fs;
    cvx___.fsre = fsre;
    cvx___.ps = ps;
    cvx___.cs = cs;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quick exit for non-verbose output %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if quick,
    if nargout,
        varargout = { fs, cvx___.ps, mpath, cvx___.mext };
    end
    cvx_load_prefs( false );
    cvx___.loaded = true;
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Verbose output (cvx_setup, cvx_version plain) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cvx_ver = '3.0beta';
cvx_bld = '1119';
cvx_bdate = 'Tue Jun 17 17:13:54 2014';
cvx_bcomm = 'c3c52c8';
line = '---------------------------------------------------------------------------';
fprintf( '\n%s\n', line );
fprintf( 'CVX: Software for Disciplined Convex Programming       (c)2014 CVX Research\n' );
fprintf( 'Version %7s, Build %4s (%7s)%38s\n', cvx_ver, cvx_bld, cvx_bcomm, cvx_bdate );
fprintf( '%s\n', line );
fprintf( 'Installation info:\n    Path: %s\n', cvx___.where );
if isoctave,
    fprintf( '    GNU Octave %s on %s\n', version, cvx___.comp );
else
    verd = ver('MATLAB');
    fprintf( '    MATLAB version: %s %s\n', verd.Version, verd.Release );
    if usejava( 'jvm' ),
        os_name = char(java.lang.System.getProperty('os.name'));
        os_arch = char(java.lang.System.getProperty('os.arch'));
        os_version = char(java.lang.System.getProperty('os.version'));
        java_str = char(java.lang.System.getProperty('java.version'));
        fprintf('    OS: %s %s version %s\n', os_name, os_arch, os_version );
        fprintf('    Java version: %s\n', java_str );
    else
        fprintf( '    Architecture: %s\n', cvx___.comp );
        fprintf( '    Java version: disabled\n' );
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check for valid version %
%%%%%%%%%%%%%%%%%%%%%%%%%%%

issue = false;
isoctave = cvx___.isoctave;
nver = cvx___.nver;
if isoctave,
    if nver <= 3.08,
        fprintf( '%s\nCVX is not yet supported on Octave.\n(Please do not waste further time trying: changes to Octave are required.\nBut they are coming! Stay tuned.)\n%s\n', line, line );
        issue = true;
    end
elseif nver < 7.08 && strcmp( cvx___.comp(end-1:end), '64' ),
    fprintf( '%s\nCVX requires MATLAB 7.8 or later (7.5 or later on 32-bit platforms).\n' , line, line );
    issue = true;
elseif nver < 7.05,
    fprintf( '%s\nCVX requires MATLAB 7.5 or later (7.8 or later on 64-bit platforms).\n' , line, line );
    issue = true;
end

%%%%%%%%%%%%%%%%%%%%%%%%
% Verify file contents %
%%%%%%%%%%%%%%%%%%%%%%%%

fid = fopen( [ mpath, fs, 'MANIFEST' ], 'r' );
if fid > 0,
    fprintf( 'Verfying CVX directory contents:' );
    manifest = textscan( fid, '%s' );
    manifest = manifest{1};
    fclose( fid );
    newman = get_manifest( mpath, fs );
    if ~isequal( manifest, newman ),
        missing = setdiff( manifest, newman );
        additional = setdiff( newman, manifest );
        if ~isempty( missing ) || ~isempty( additional ),
            if ispc,
                missing = strrep( missing, '/', fs );
                additional = strrep( additional, '/', fs );
            end
            if ~isempty( missing ),
                fprintf( '\n    WARNING: The following files/directories are missing:\n' );
                isdir = cellfun(@(x)x(end)==fs,missing);
                missing_d = missing(isdir);
                missing_f = missing(~isdir);
                while ~isempty( missing_d ),
                    mdir = missing_d{1};
                    ss = strncmp( missing_d, mdir, length(mdir) );
                    tt = strncmp( missing_f, mdir, length(mdir) );
                    fprintf( '        %s%s%s + %d files, %d subdirectories\n', mpath, fs, mdir, nnz(tt), nnz(ss) - 1 );
                    missing_d(ss) = [];
                    missing_f(tt) = [];
                end
                for k = 1 : min(length(missing_f),10),
                    fprintf( '        %s%s%s\n', mpath, fs, missing_f{k} );
                end
                if length(missing_f) > 10,
                    fprintf( '        (and %d more files)\n', length(missing_f) - 10 );
                end
                fprintf( '    These omissions may prevent CVX from operating properly.\n'  );
            end
            if ~isempty( additional ),
                if isempty( missing ), fprintf( '\n' ); end
                fprintf( '    WARNING: The following extra files/directories were found:\n' );
                isdir = cellfun(@(x)x(end)==fs,additional);
                issedumi = cellfun(@any,regexp( additional, [ '^sedumi.*[.]', mexext, '$' ] ));
                additional_d = additional(isdir&~issedumi);
                additional_f = additional(~isdir&~issedumi);
                additional_s = additional(issedumi);
                while ~isempty( additional_d ),
                    mdir = additional_d{1};
                    ss = strncmp( additional_d, mdir, length(mdir) );
                    tt = strncmp( additional_f, mdir, length(mdir) );
                    fprintf( '        %s%s%s + %d files, %d subdirectories\n', mpath, fs, mdir, nnz(tt), nnz(ss) - 1 );
                    additional_d(ss) = [];
                    additional_f(tt) = [];
                end
                for k = 1 : min(length(additional_f),10),
                    fprintf( '        %s%s%s\n', mpath, fs, additional_f{k} );
                end
                if length(additional_f) > 10,
                    fprintf( '        (and %d more files)\n', length(additional_f) - 10 );
                end
                fprintf( '    These files may alter the behavior of CVX in unsupported ways.\n' );
                if ~isempty( additional_s ),
                    fprintf( '    ERROR: obsolete versions of SeDuMi MEX files were found:\n' );
                    for k = 1 : length(additional_s),
                        fprintf( '        %s%s%s\n', mpath, fs, additional_f{k} );
                    end
                    fprintf( '    These files are now obsolete, and must be removed to ensure\n' );
                    fprintf( '    that SeDuMi operates properly and produces sound results.\n' );
                    if ~issue,
                        fprintf( '    Please remove these files and re-run CVX_SETUP.\n' );
                        issue = true;
                    end
                end
            end
        else
            fprintf( '\n    No missing files.\n' );
        end
    else
        fprintf( '\n    No missing files.\n' );
    end
else
    fprintf( 'Manifest missing; cannot verify file structure.\n' ) ;
end
if ~compile,
    mexpath = [ mpath, fs, 'lib', fs ];
    mext = cvx___.mext;
    if ( ~exist( [ mexpath, 'cvx_eliminate_mex.', mext ], 'file' ) || ...
         ~exist( [ mexpath, 'cvx_bcompress_mex.', mext ], 'file' ) ) && ~issue,
        issue = true;
        if ~isempty( msub ),
          mexpath = [ mexpath, msub, fs ];
          issue = ~exist( [ mexpath, 'cvx_eliminate_mex.mex' ], 'file' ) || ...
                         ~exist( [ mexpath, 'cvx_bcompress_mex.mex' ], 'file' );
        end
        if issue,
          fprintf( '    ERROR: one or more MEX files for this platform are missing.\n' );
          fprintf( '    These files end in the suffix ".%s". CVX will not operate\n', mext );
          fprintf( '    without these files. Please visit\n' );
          fprintf( '        http://cvxr.com/cvx/download\n' );
          fprintf( '    And download a distribution targeted for your platform.\n' );
        end
    end
end

%%%%%%%%%%%%%%%
% Preferences %
%%%%%%%%%%%%%%%

cvx_load_prefs( true );

%%%%%%%%%%%%%%%%
% License file %
%%%%%%%%%%%%%%%%

if isoctave,
    if ~isempty( cvx___.license ),
        fprintf( 'CVX Professional is not supported with Octave.\n' );
    end
elseif cvx___.jver < 1.6,
    fprintf('       WARNING: full support for CVX Professional licenses\n' );
    fprintf('       requres Java version 1.6.0 or later. Please upgrade.\n' );
elseif exist( 'cvx_license', 'file' ),
    cvx_license( args{:} );
end

%%%%%%%%%%%%%%%
% Wrapping up %
%%%%%%%%%%%%%%%

if ~issue,
    cvx___.loaded = true;
end
clear fs;
fprintf( '%s\n', line );
if length(dbstack) <= 1,
    fprintf( '\n' );
end

%%%%%%%%%%%%%%%%%%%%%%
% Preference loading %
%%%%%%%%%%%%%%%%%%%%%%

function cvx_load_prefs( verbose )

global cvx___
fs = cvx___.fs;
isoctave = cvx___.isoctave;
errmsg = '';
if isoctave,
    pfile = [ prefdir, fs, '.cvx_prefs.mat' ];
else
    pfile = [ regexprep( prefdir(1), [ cvx___.fsre, 'R\d\d\d\d\w$' ], '' ), fs, 'cvx_prefs.mat' ];
end
if verbose,
	fprintf( 'Preferences:\n    Path: %s\n', pfile );
end
outp = [];
if exist( pfile, 'file' ),
	try
		outp = load( pfile );
    catch errmsg
        errmsg = strcat( { '    ' }, cvx_error( errmsg, 67, true ) );
        header = 'An unexpected error occurred loading the preferences:';
    	errmsg = [ { header }, errmsg ];
	end
else
	errmsg = 'No preference file was found.';
end
if ~isempty( outp ),
    try
        cvx___.path       = outp.path;
        cvx___.license    = outp.license;
        cvx___.solvers    = outp.solvers;
        cvx___.broadcast  = outp.broadcast;
        cvx___.expert     = outp.expert;
        cvx___.precision  = outp.precision;
        cvx___.precflag   = outp.precflag;
        cvx___.quiet      = outp.quiet;
        cvx___.profile    = outp.profile;
    catch
        outp = [];
        errmsg = 'Preferences are outdated.';
    end
end
if isempty( outp )
    cvx___.expert     = false;
    cvx___.precision  = [eps^0.5,eps^0.5,eps^0.25];
    cvx___.precflag   = 'default';
    cvx___.broadcast  = cvx___.isoctave;
    cvx___.quiet      = false;
    cvx___.profile    = false;
    cvx___.path       = [];
    cvx___.solvers    = [];
    cvx___.license    = [];
end
cvx___.pfile = pfile;
if ~isempty( errmsg ),
    if ischar( errmsg ),
        errmsg = { errmsg };
    end
	if verbose,
    	errmsg{end+1} = 'Default preferences will be used.';
		fprintf( '    %s\n', errmsg{:} );
    else
        errmsg{end+1} = [ 'Path: ', pfile ];
    	errmsg{end+1} = 'Please re-run CVX_SETUP.';
        temp = sprintf( '%s\n', errmsg{:} );
	    cvx_throw( temp(1:end-1) );
	end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Recursive manifest building function %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function newman = get_manifest( mpath, fs )
dirs   = {};
files  = {};
nfiles = dir( mpath );
ndir   = '';
dndx   = 0;
pat2   = '^\.|~$|';
pat    = '^\.|~$|^cvx_license.[md]at$|^doc$|^examples$';
while true,
    isdir  = [ nfiles.isdir ];
    nfiles = { nfiles.name };
    tt     = cellfun( @isempty, regexp( nfiles, pat ) ); pat = pat2;
    isdir  = isdir(tt);
    nfiles = nfiles(tt);
    ndirs  = nfiles(isdir);
    if ~isempty(ndirs),
        dirs = horzcat( dirs, strcat(strcat(ndir,ndirs), fs ) ); %#ok
    end
    nfiles = nfiles(~isdir);
    if ~isempty(nfiles),
        files = horzcat( files, strcat(ndir,nfiles) ); %#ok
    end
    if length( dirs ) == dndx, break; end
    dndx = dndx + 1;
    ndir = dirs{dndx};
    nfiles = dir( [ mpath, fs, ndir ] );
end
[tmp,ndxs1] = sort(upper(dirs)); %#ok
[tmp,ndxs2] = sort(upper(files)); %#ok
newman = horzcat( dirs(ndxs1), files(ndxs2) );
if ispc,
    newman = strrep( newman, fs, '/' );
end
newman = newman(:);

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
classdef nme_bus3p_acp < mp.nme_bus3p & mp.form_acp

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

%     properties
%         name = 'bus3p';
%     end

    methods
        function obj = add_vvars(obj, nm, dm, idx)
            dme = obj.data_model_element(dm);
            nb = obj.nk;
            p = idx{1};

            %% prepare angle bounds for ref buses
            ref = dme.type == mp.NODE_TYPE.REF;
            va_lb = -Inf(nb, 1);
            va_ub =  Inf(nb, 1);
            vm_start = dme.(sprintf('vm%d_start', p));
            va_start = dme.(sprintf('va%d_start', p));
            va1_lb(ref) = va_start(ref);
            va1_ub(ref) = va_start(ref);

            if p == 1
                nm.init_indexed_name('va', 'Va3', {obj.nn});
                nm.init_indexed_name('vm', 'Vm3', {obj.nn});
            end
            nm.add_var('va', 'Va3', {p}, nb, va_start, va_lb, va_ub);
            nm.add_var('vm', 'Vm3', {p}, nb, vm_start, 0, 2);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% The COBRAToolbox: testSubSystemModification.m
%
% Purpose:
%     - test Manipulation of subSystems and retrieval by subSystems in
%     addSubSystemsToReactions, setRxnSubSystem, removeSubSystemFromReaction,
%     findRxnsFromSubSystem, getModelSubSystems, isReactionInSubSystem
%
% Authors:
%     Thomas Pfau - Nov 2017

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testSubSystemModification.m'));
cd(fileDir);

model = getDistributedModel('ecoli_core_model.mat');
%This is a correct model. Now lets do some subSystem Manipulation
originalSubSystems = {'Exchange','Anaplerotic reactions','Citric Acid Cycle','Glutamate Metabolism','Glycolysis/Gluconeogenesis','Inorganic Ion Transport and Metabolism','Oxidative Phosphorylation','Pentose Phosphate Pathway','Pyruvate Metabolism','Transport, Extracellular'};
subSysToManipulate = 'Citric Acid Cycle';
originalPositions = [4;5;8;15;46;59;64;90];

%First get all reactions with 'Citric Acid Cycle' subSystem annotation
[reactionNames,reactionPos] = findRxnsFromSubSystem(model,subSysToManipulate);
assert(isempty(setxor(reactionPos,[4;5;8;15;46;59;64;90]))); %These and only these
assert(isequal(reactionNames,model.rxns(reactionPos))); %Same order

%Now, get all subSystems in the model
subsInModel = getModelSubSystems(model);
assert(isempty(setxor(subsInModel,originalSubSystems)));

%Now, add Citric Acid Cycle to the first three reactions
originalSubsystems = model.subSystems;
modelNew = addSubSystemsToReactions(model,model.rxns(1:3),subSysToManipulate);
[reactionNames,reactionPos] = findRxnsFromSubSystem(modelNew,subSysToManipulate);
assert(isempty(setxor(union(originalPositions,[1:3]),reactionPos))); %Added only to those and not changed elsewhere.
assert(all(~cellfun(@(x) isequal(x,{''}),originalSubsystems(1:3)))); %The first three subSystems are non Empty, i.e. this empty is not removed
assert(all(cellfun(@(x,y) all(ismember(y,x)),modelNew.subSystems,originalSubsystems))); %Old subSystems are retained.

%Add it to an empty reactions
modelWithEmptyReplaced = addSubSystemsToReactions(model,13,subSysToManipulate); %Biomass gets assigned to Citric Acid Cycle
assert(isequal(modelWithEmptyReplaced.subSystems{13},{subSysToManipulate})); %Assert that the '' was removed and only the Citric Acid cycle is there.

%Remove the Citric Acid Cycle from all reactions
rxnToDo = true(numel(model.rxns),1);
modelWOTCA = removeSubSystemsFromRxns(modelNew,rxnToDo,subSysToManipulate); %Remove it from model with the added SubSystems.
assert(all(cellfun(@(x) isequal(x,{''}),modelWOTCA.subSystems(originalPositions)))); %The former TCA reactions should have no subSystem left.
rxnToDo(originalPositions) = false; %Ignore them
assert(all(cellfun(@(x,y) isempty(setxor(x,y)),modelWOTCA.subSystems(rxnToDo),model.subSystems(rxnToDo)))); %All others should have their original

%Now, set the TCA as subSystem for a batch of reactions:
reactionsToSet = [3 4 7 8 13];
modelWithSet = setRxnSubSystems(model,reactionsToSet,subSysToManipulate);
assert(all(cellfun(@(x) isequal(x,{subSysToManipulate}),modelWithSet.subSystems(reactionsToSet))));
notChanged = setdiff(1:numel(model.rxns),reactionsToSet);
assert(all(cellfun(@(x,y) isequal(x,y), modelWithSet.subSystems(notChanged),model.subSystems(notChanged))));

%Test for multiple values
SomeSubSys = 'testtSubSys';
modelTest = addSubSystemsToReactions(model,1:4,SomeSubSys);
newSubs = getModelSubSystems(modelTest);
assert(isempty(setxor(newSubs,union(originalSubSystems,SomeSubSys))));

%finally test isReactionInSubSystem
assert(~isReactionInSubSystem(model,model.rxns{9},'Citric Acid Cycle'))
assert(all(isReactionInSubSystem(model,model.rxns(4:5),'Citric Acid Cycle')));
assert(isequal(isReactionInSubSystem(model,21:-1:18,'Exchange'),[true; true; false; false]));

%Return to old path
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [ll,ei] = minlen2(pp,tt)
%MINLEN2 return the minimum length edge for each triangle in
%a two-dimensional triangulation.
%   [ELEN,IMIN] = MINLEN2(VERT,TRIA) returns the minimum le-
%   ngth ELEN and local edge index IMIN for all triangles in
%   the triangulation {VERT,TRIA}.

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 16/01/2017

%---------------------------------------------- basic checks
    if (~isnumeric(pp) || ~isnumeric(tt) )
        error('minlen2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(pp) ~= +2 || ndims(tt) ~= +2 )
        error('minlen2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(pp,2)~= +2 || size(tt,2) < +3 )
        error('minlen2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    nnod = size(pp,1) ;

%---------------------------------------------- basic checks
    if (min(min(tt(:,1:3))) < +1 || ...
            max(max(tt(:,1:3))) > nnod )
        error('minlen2:invalidInputs', ...
            'Invalid TRIA input array.') ;
    end

%------------------------------------------ compute edge-len
    l1 = sum((pp(tt(:,2),:) ...
             -pp(tt(:,1),:)).^2,2);
    l2 = sum((pp(tt(:,3),:) ...
             -pp(tt(:,2),:)).^2,2);
    l3 = sum((pp(tt(:,1),:) ...
             -pp(tt(:,3),:)).^2,2);

%------------------------------------------ compute min.-len
   [ll,ei] = min([l1,l2,l3],[],+2);

end

author: CHLNDDEV


text:
function feats =concatenate_feat(frame_feat, nprev, nnext)
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
% concatenate nprev + self + nnext frames together

    [nframes, ndim] = size(frame_feat);
    nDim= (nprev+nnext+1)*ndim;

    % duplicate the beginning and the end of frame
    begin_frame = repmat(frame_feat (1, :), nprev, 1);
    end_frame = repmat(frame_feat ( end, : ), nnext, 1) ;

    frame_feat2 = [begin_frame; frame_feat; end_frame];

    % original frame list
    ndx= nprev+1: nprev+nframes;

    % frame list + nprev + nnext frames
    ndx_merge= repmat( ndx, nprev+nnext+1, 1) + ...
               repmat( (0:nprev+nnext)'-nprev, 1, nframes) ;

    ndx_merge=ndx_merge(:);

    % merge them together
    feats=reshape( frame_feat2(ndx_merge,:)', nDim, nframes)';
end

author: posenhuang

repo: deeplearningsourceseparation


text:
function x = i4vec_decrement ( n, x )

%*****************************************************************************80
%
%% I4VEC_DECREMENT decrements an I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input/output, integer X(N), the array to be decremented.
%
  x(1:n) = x(1:n) - 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Some example for running filterbanks and computers such as stft computer
%
% This filterbank library is the implementation of paper:
% "Learning Filter Banks from Raw Speech for Phone Recognition":
% https://arxiv.org/abs/1711.01161
%
% Please find the path management in the startup.m script in the root directory
% of this repository. Note that by starting matlab in the root directory, this
% script should automatically run. If it is not the case, you can also move to the
% root directory and run this script manually.
%
%
% Copyright (c) 2018 Department of Computer Science,
%                    University of Toronto, Canada,
%                    Vector Institute, Canada
%
% License
% This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Yingxue Wang <yingxue@cs.toronto.edu>
%  Sean Robertson <sdrobert@cs.toronto.edu>
%

clear all;

%% Using Scaling Functions
% Mel Scaling Function
% Please use the Voicebox
% (http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html)
% to convert between Hz, Mel, Erb and MIDI frequency scales
disp('Example - Using Mel Scaling Function')
scaling_function = MelScaling();
hertz = 200;
scale = scaling_function.hertz_to_scale(hertz);
disp('done')

%% Using Window Functions
% Hann Window, Hamming Window, Gamma Window
disp('Example - Using Window Function: Impulse Response of a Gamma Window')
window_size = 1000;
peak_ratio = 0.75;
order = 2;
window = GammaWindow(order, peak_ratio).get_impulse_response(window_size);
disp('done')

%% Using Compute Method to Compute FilterBanks
% Example: Compute a Standard GaborFilterBank
disp('Example - Compute a Standard GaborFilterBank')

bank = GaborFilterBank(MelScaling());
frame_length_ms = 25;
frame_shift_ms = 10;
frame_style = 'causal'; %{'causal', 'centered'}
include_energy = true;  %{true, false}
pad_to_nearest_power_of_two = true; % {true, false}
use_log = true; %{true, false}
use_power = true; %{true, false}

% Create a random signal, for example
signal_len = 2 ^ 10;
signal = rand(signal_len,1);    % Example Lengths: 'empty buffer',
                                % 'length 1 buffer',
                                % 'medium buffer'(2^8),
                                % 'large buffer'(2^10)
% disp(signal)

% Create Computer
computer = ShortTimeFourierTransformFrameComputer(...
                            bank, ...
                            frame_length_ms, ...
                            frame_shift_ms, ...
                            frame_style, ...
                            include_energy, ...
                            pad_to_nearest_power_of_two, ...
                            use_log, ...
                            use_power);

% Below are different methods to compute the filterbank
    % Use computer.compute_full to obtain the features
    feats_full = computer.compute_full(signal);
    % disp(feats_full)

    % Use frame_by_frame_calculation to obtain the features
    feats_framewise = Util.frame_by_frame_calculation(computer, signal);
    % disp(feats_framewise)

    % Compute a piece of signal
    feats_chunk = computer.compute_chunk(signal(1:ceil(length(signal)/4)));
    feats_chunk = computer.finalize();
    % disp(feats_chunk)

disp('done')

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function varargout = midiOut(varargin)

% midiOut -- connect to a MIDI output device and send notes
%
% >> S = midiOut('L')
% Returns a list of MIDI devices as a structure array. Fields are 'index',
% 'name', 'input' and 'output'. The last two are flags (1 or 0) that determine
% whether the device can be used as an input, or output, or both.
%
% >> midiOut('O', index)
% Opens a MIDI device by its index (as returned by aforementioned 'L' call).
% If the device is already opened, nothing will happen. If another device is
% opened, that one will be closed first.
%
% >> midiOut('C')
% Closes the current connection.
%
% >> midiOut('+', channel, notes, velocities)
% Send one or more 'note on' messages at the given channel (1..16).
% 'notes' and 'velocities' must be given as double-precision arrays with
% equal number of elements. Messages will be transmitted in order.
%
% >> midiOut('-', channel, notes, velocities)
% Send one or more 'note off' messages, same parameters as for 'note on'.
%
% >> midiOut('.', channel)
% Send an 'all notes off' message at the given channel.
%
% >> midiOut('P', channel, program)
% Send a 'program change' to the given channel.
%
% >> midiOut(msg)
% Send raw data. 'msg' must be of type 'uint8' and have a multiple of 3 elements.
%
% For example, play the note C4 (261.63 Hz) on the piano
% >> midiOut('O', 2)
% >> midiOut('+', 1, 60, 127)
%
% See also MIDIIN

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

error('the %s mex file is not available for your platform (%s)', mfilename, mexext);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function logProbs = makeLogProbs(counts,varargin)
%function logProbs = makeLogProbs(counts,varargin)
%each row is a set of counts. convert to normalized log probabilities.
[min_log_prob] = process_options(varargin,'min-log-prob',-200);
logProbs = log(full(normalize_rows(counts)));
logProbs(logProbs < min_log_prob) = min_log_prob;
logProbs = log(normalize_rows(exp(logProbs)));
end

author: jacobeisenstein

repo: SAGE


text:
function  [im_out, par] = SC_Sigma_1AG(par)
im_out    =   par.nim;
% parameters for noisy image
[h,  w, ch]      =  size(im_out);
par.h = h;
par.w = w;
par.ch = ch;
par = SearchNeighborIndex( par );
for ite  =  1 : par.outerIter
    %     % iterative regularization
    im_out = im_out+par.delta*(par.nim - im_out);
    % image to patches and estimate local noise variance
    Y = Image2PatchNew( im_out, par );
    % estimation of noise variance
    dif = mean( mean( mean( (par.nim-im_out).^2 ) ) );
    par.sigma = sqrt( abs( par.nSig^2 - dif ) );
    % estimation of noise variance
    if mod(ite-1, par.innerIter)==0
        par.nlsp = par.nlsp - par.nlspgap;
        % searching  non-local patches
        blk_arr = Block_Matching( Y, par );
    end
    % Weighted Sparse Coding
    Y_hat = zeros(par.ps2ch, par.maxrc, 'single');
    W_hat = zeros(par.ps2ch, par.maxrc, 'single');
    for i = 1:par.lenrc
        index = blk_arr(:, i);
        nlY = Y( : , index );
        DC = mean(nlY, 2);
        nDCnlY = bsxfun(@minus, nlY, DC);
        % update D and S
        [D, ~, ~] = svd( full(nDCnlY), 'econ' );
        % update C by soft thresholding
        B = D' * nDCnlY;
        C = sign(B) .* max( abs(B) -  par.lambda*par.sigma^2, 0 );
        % update Y
        nDCnlYhat = D * C;
        % add back DC components
        nlYhat = bsxfun(@plus, nDCnlYhat, DC);
        % aggregation
        Y_hat(:, index) = Y_hat(:, index) + nlYhat;
        W_hat(:, index) = W_hat(:, index) + ones(par.ps2ch, par.nlsp);
    end
    % Reconstruction
    im_out = PGs2Image(Y_hat, W_hat, par);
    % calculate the PSNR and SSIM
    PSNR =   csnr( im_out*255, par.I*255, 0, 0 );
    SSIM      =  cal_ssim( im_out*255, par.I*255, 0, 0 );
    fprintf('Iter %d : PSNR = %2.4f, SSIM = %2.4f\n', ite, PSNR, SSIM);
    par.PSNR(ite, par.image) = PSNR;
    par.SSIM(ite, par.image) = SSIM;
end
return;

author: csjunxu


text:
% Script demonstrating usage of the cbpdnjnt function.
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2016-07-01
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'Copyright' and 'License' files
% distributed with the library.

% Load dictionary
load([sporco_path '/Data/ConvDict.mat']);
dmap = containers.Map(ConvDict.Label, ConvDict.Dict);
D = dmap('12x12x36');

% Load test image
s = single(stdimage('lena'))/255;
if isempty(s),
  error('Data required for demo scripts has not been installed.');
end

% Highpass filter test image
npd = 16;
fltlmbd = 5;
[sl, sh] = lowpass(s, fltlmbd, npd);

% Compute representation
lambda = 1e-2;
mu = 1e-3;
opt = [];
opt.Verbose = 1;
opt.MaxMainIter = 500;
opt.rho = 100*lambda + 1;
opt.RelStopTol = 1e-3;
opt.AuxVarObj = 0;
opt.HighMemSolve = 1;
[X, optinf] = cbpdnjnt(D, sh, lambda, mu, opt);

% Compute reconstruction
DX = squeeze(ifft2(sum(bsxfun(@times, fft2(D, size(X,1), ...
             size(X, 2)), fft2(X)),3), 'symmetric'));

figure;
subplot(1,3,1);
plot(optinf.itstat(:,2));
xlabel('Iterations');
ylabel('Functional value');
subplot(1,3,2);
semilogy(optinf.itstat(:,6));
xlabel('Iterations');
ylabel('Primal residual');
subplot(1,3,3);
semilogy(optinf.itstat(:,7));
xlabel('Iterations');
ylabel('Dual residual');

figure;
imagesc(squeeze(sum(abs(X),3)));
title('Sum of absolute value of coefficient maps');

figure;
subplot(1,3,1);
imdisp(s);
title('Original image');
subplot(1,3,2);
imdisp(DX + sl);
title(sprintf('Reconstructed image (SNR: %.2fdB)', snr(s, DX + sl)));
subplot(1,3,3);
sdf = DX + sl - s;
sdf = sdf - min(sdf(:));
sdf = sdf / max(sdf(:));
imagesc(sdf);
axis image; axis off;
title('Difference');

author: thfylsty


text:
function j_polynomial_plot ( n_vec, alpha_vec, beta_vec, filename )

%*****************************************************************************80
%
%% J_POLYNOMIAL_PLOT plots Jacobi polynomials J(n,a,b,x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N_VEC(*), the orders of 1 or more polynomials
%    to be plotted together.
%
%    Input, integer ALPHA_VEC(*), BETA_VEC(*), the alpha and beta values
%    for each polynomial.
%
%    Input, string FILENAME, the name into which the graphics information is
%    to be stored.  Note that the PNG format will be used.
%
  a = -1.0;
  b = +1.0;
  m = 501;
  x = linspace ( a, b, m );
  x = x';
  vec_num = length ( n_vec );

  clf
  hold on
  for i = 1 : vec_num
    n = n_vec(i);
    alpha = alpha_vec(i);
    beta = beta_vec(i);
    y = j_polynomial ( m, n, alpha, beta, x );
    plot ( x, y(:,n+1), 'LineWidth', 2 );
  end
  grid on
  xlabel ( '<--- X --->' )
  ylabel ( '<--- J(n,a,b,x) --->' )
  title ( 'Jacobi polynomials J(n,a,b,x)' )
  hold off
  print ( '-dpng', filename )

  return
end


text:
function mfilelist = mfilegrep(searchstring,startpath,recurse,wholewords,casesens,showline)

% MFILELIST = MFILEGREP(SEARCHSTRING, STARTPATH, RECURSE, WHOLEWORDS, CASESENS, SHOWLINE)
% returns in a cell array the names of all m-files within the startpath
% (searching recursively or non-recursively), that contain the text searchstring.
%
% The function calls the embedded subfunction 'filesearch' to locate the files.
%
% Takes from 1 to 6 input arguments; at a minimum, the user must specify the desired searchstring.
% STARTPATH and RECURSE (= 0 or 1) may also be inputs; otherwise, a gui interface will prompt the
%       user for these values.
% WHOLEWORDS and CASESENS (= 0 or 1) may be indicated as well; if not, the default values of 0 and 0,
%       respectively, will be used. Wholewords == 1 specifies that only whole-word matches will be returned;
%       wholewords == 0 will also return partial matches. Casesens == 1 returns case-sensitive matches only;
%       casesens == 0 ignores case.
% The final (optional) input is SHOWLINE (= 0 or 1; default = 0). If showline == 1, the output also
%       shows the text of all lines found in the file(s) containing the search string.
%
% EXAMPLES:
%
% FILELIST = MFILEGREP('domct') prompts the user via uicontrols for starting path and inclusion or
%             exclusion of subdirectories, then searches all mfiles on the indicated path for 'domct'.
% FILELIST = MFILEGREP('domct', 'C:\Mfiles') searches specified directory for m-files; prompts the user for:
%             inclusion or exclusion of subdirectories.
% FILELIST = MFILEGREP('domct','C:\Mfiles',1) searches specified directory recursively for m-files.
% FILELIST = MFILEGREP('domct','C:\Mfiles',1, 0, 1) searches specified directory recursively for m-files;
%             in this case, the search returns partial word matches (e.g., 'examp' matches 'example'),
%             and the search is case sensitive.
% FILELIST = MFILEGREP('domct','C:\Mfiles',0, 0, 1, 1) searches specified directory non-recursively for m-files;
%             returns partial and case-sensitive word matches, and includes all lines containing the search string
%             in the output cell array.
%
% Copyright Brett Shoelson, Ph.D.
% Shoelson Consulting
% brett.shoelson@joslin.harvard.edu
% V1: 6/22/01. Based extensively on the new version of collectcode.m.
% V2. 9/30/01. Modified to include "showline" option.
% V3. 10/25/01. Modified to match default behavior to program description.

if ~nargin
	error('Requires at least one argument indicating desired searchstring.');
elseif nargin == 1
	% Get starting path
	[filename,startpath]=uigetfile('*.m','Select any m-file in the desired starting pathname.');
	%Remove trailing '\' character
	if strcmp(startpath(end),'\'),startpath=startpath(1:end-1);end
	% Query for inclusion of subdirectories
	recurse=questdlg('Include subdirectories?','','YES','No','YES');
	if strcmp(recurse,'YES'),recurse=1;else;recurse=0;end
	wholewords = 0; casesens = 0; showline = 0;
elseif nargin == 2
	recurse=questdlg('Include subdirectories?','','YES','No','YES');
	if strcmp(recurse,'YES'),recurse=1;else;recurse=0;end
	wholewords = 0; casesens = 0; showline = 0;
elseif nargin == 3
	wholewords = 0; casesens = 0; showline = 0;
elseif nargin == 4
	casesens = 0; showline = 0;
elseif nargin == 5
    showline = 0;
elseif nargin > 6
	error('Too many input arguments.');
end

% Check for appropriateness of inputs:
if ~ischar(searchstring)
	error('searchstring must be a character string.');
end
if recurse ~= 1 & recurse ~= 0
	error('recurse must be 1 or 0.');
end
if wholewords ~= 1 & wholewords ~= 0
	error('wholewords must be 1 or 0.');
end
if casesens ~= 1 & casesens ~= 0
	error('casesens must be 1 or 0.');
end
if showline ~= 1 & showline ~= 0
	error('showline must be 1 or 0.');
end

currpath = cd;

% Compile filenames (calls subfunction filesearch)
files = filesearch('m', startpath, recurse);

cd(startpath);

line = 0;
mfilelist = {};

if ~ casesens
	searchstring = lower(searchstring);
end

if wholewords
	alphachars = ['A' ; 'B' ; 'C' ; 'D' ; 'E' ; 'F' ; 'G' ; 'H' ; 'I' ; 'J' ; 'K' ; 'L' ; ...
			'M' ; 'N' ; 'O' ; 'P' ; 'Q' ; 'R' ; 'S' ; 'T' ; 'U' ; 'V' ; 'W' ; 'X' ; ...
			'Y' ; 'Z'];
	alphachars = [alphachars; lower(alphachars); '1' ; '2' ; '3' ; '4' ; '5' ; '6' ; '7' ; '8' ; '9' ; '0'];
end

% Loop through all mfiles on search path

fprintf('Searching directory(ies)....\n');

for i = 1 : size(files, 1)
	filename = files{i};
	%fprintf('Searching %s....\n', filename);
	% Open mfile(i) and read contents line by line; write lines to codetext.txt
	fid3 = fopen(filename, 'rt');
	while line ~=-1
		line = fgets(fid3);
		if line ~= -1
			if ~ casesens
				line = lower(line);
			end
			if ~isempty(findstr(line, searchstring))
				% findstr finds the instance of the shorter string in the longer; the following if/then ensures that
				% the function won't return instances of subelements (eg., findstr('y','classify') is non-empty).
				if length(line) >= length(searchstring)
					if wholewords
						wholefound = iswhole(searchstring, line, alphachars);
						if wholefound
							if isempty(intersect(mfilelist,filename))
								mfilelist{size(mfilelist, 1) + 1, 1} = filename;
							end
							if showline
								mfilelist{size(mfilelist, 1) + 1, 1} = line;
							else
								break
							end
						end
					else
						if isempty(intersect(mfilelist,filename))
							mfilelist{size(mfilelist, 1) + 1, 1} = filename;
						end
						if showline
							mfilelist{size(mfilelist, 1) + 1, 1} = line;
						else
							break
						end
					end % if wholewords
				end % if length(line) >= length(searchstring)
			end % if ~isempty(findstr(line, searchstring))
		end % if line ~=- 1
	end % while line ~=- 1
	line = 0;
	% Close mfile(i) and continue
	fclose(fid3);
end % for i = 1 : size(files, 1)

cd(currpath);

function files = filesearch(extension,startpath,recurse)
%
%  Search a specified path (recursively or non-recursively) for all instances of files with a specified extension.
%  Results are returned in a cell array ('files').
%
%  FILES = FILESEARCH(EXTENSION,STARTPATH,1) searches for all instances of *.extension in the starting directory OR
%          in any subdirectory beneath the starting directory.
%  FILES = FILESEARCH(EXTENSION,STARTPATH,0) searches for all instances of *.extension in the starting directory ONLY,
%          and excludes files in subdirectories.
%  Examples:
%
%  files = filesearch('txt','C:\WinNT',1)
%  mfiles = filesearch('m','D:\Mfiles',0)

% Copyright Brett Shoelson, Ph.D., 5/29/2001
% brett.shoelson@joslin.harvard.edu

if recurse
	files = {};
	%Create string of recursive directories/subdirectories
	paths = genpath(startpath);
	%Find instances of the path separator
	seplocs = findstr(paths,pathsep);
	%Parse paths into individual (vertically catenated) directories
	if ~isempty(seplocs)
		directories = paths(1:seplocs(1)-1);
		for i = 1:length(seplocs)-1
			directories = strvcat(directories,paths(seplocs(i)+1:seplocs(i+1)-1));
		end
		%Search each directory for instances of *.extension, appending to current list
		for i = 1:size(directories,1)
			%Compile located files as vertically catenated strings
			tmp = dir([deblank(directories(i,:)) '\*.' extension]);
			if ~isempty(tmp),tmp = char(tmp.name);end
			%Update files to reflect newly detected files. (Omit trailing blanks.)
			for j = 1:size(tmp,1)
				files{size(files,1)+1,1} = [deblank(directories(i,:)) '\' deblank(tmp(j,:))];
			end
		end
	end
else %Search non-recursively
	files = {};
	tmp = dir([startpath '\*.' extension]);
	if ~isempty(tmp),tmp = char(tmp.name);end
	for j = 1:size(tmp,1)
		files{size(files,1)+1,1} = [startpath '\' deblank(tmp(j,:))];
	end
end

function wholefound = iswhole(searchstring, teststring, alphachars)
charlocs = findstr(searchstring, teststring);
prevlocs = charlocs - 1;
postlocs = charlocs + length(searchstring);
delprevposns = [find(prevlocs == 0)];
delpostposns = [find(postlocs >= length(teststring))];
prevlocs(delprevposns) = '';
postlocs(delpostposns) = '';
prevmbrs = ismember(teststring(prevlocs)', alphachars);
postmbrs = ismember(teststring(postlocs)', alphachars);
if ~ isempty(delprevposns), prevmbrs = [0; prevmbrs]; end
if ~ isempty(delpostposns), postmbrs = [postmbrs; 0]; end
wholefound = any(~ or (prevmbrs, postmbrs) );

author: Sable

repo: mcbench-benchmarks


text:
%% Example
% real time control of the KUKA iiwa 7 R 800
% Moving first joint of the robot, using a sinisoidal function

% An example script, it is used to show how to use the different
% functions of the KUKA Sunrise matlab toolbox

% First start the server on the KUKA iiwa controller
% Then run this script using Matlab

% Copyright: Mohammad SAFEEA, 15th of June 2017

% Important: Be careful when runnning the script, be sure that no human, nor obstacles
% are around the robot

close all,clear all;clc;

warning('off')

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  disp('Connection could not be establised, script aborted');
  return;
end

%% Go to initial position

jPos={0,0,0,0,0,0,0};

setBlueOn(t_Kuka); % turn on blue light

relVel=0.5;
movePTPJointSpace( t_Kuka , jPos, relVel); % move to initial configuration

%% Start direct servo in joint space
realTime_startDirectServoJoints(t_Kuka);
%% Some variables
w=2.0; % motion constants, frequency rad/sec
A=pi/6; % motion constants, amplitude of motion
dt=0;
counter=0;
%% Initiate timing
tic;
%% Control loop
try
    daCount=0;
    while(dt<(80*pi/w))
        if daCount==0
            t0=toc; % calculate initial time
            t_0=toc;
            daCount=1;
        end
     %% perform trajectory calculation here
      time=toc;
      dt=time-t0;
      jPos{1}=A*(1-cos(w*dt));
      %% Send joint positions to robot
        if(toc-t_0>0.003)
            counter=counter+1;
            sendJointsPositionsf( t_Kuka ,jPos);
            t_0=toc;
        end
    end
    tstart=t0;
    tend=time;
    rate=counter/(tend-tstart);
    %% Stop the direct servo motion
    realTime_stopDirectServoJoints( t_Kuka );
    fprintf('\nTotal execution time is %f: \n',tend-t0 );
    fprintf('\nThe rate of joint nagles update per second is: \n');
    disp(rate);
    fprintf('\n')
    pause(2);
    %% turn off light
    setBlueOff(t_Kuka);
    net_turnOffServer( t_Kuka )
    disp('Direct servo motion completed successfully')
    warning('on')
    return;
catch
    % in case of error turn off the server
    net_turnOffServer( t_Kuka );
    disp('Error during execution the direct servo motion')
    fclose(t_Kuka);
    warning('on')
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function model = ResNet50_for_DetectTrack_ILSVRCvid(model)
% ResNet 50 layers D&T with OHEM training (finetuned from res3a)

model.solver_def_file        = fullfile(pwd, 'models', 'rfcn_prototxts', 'ResNet-50L_ILSVRCvid_corr', 'solver_160k240k_lr1_4.prototxt');
model.test_net_def_file      = fullfile(pwd, 'models', 'rfcn_prototxts', 'ResNet-50L_ILSVRCvid_corr', 'test_track.prototxt');

model.net_file               = fullfile(pwd, 'models', 'pre_trained_models', 'ResNet-50L', 'ResNet-50-D-ilsvrc-vid.caffemodel');
model.mean_image             = fullfile(pwd, 'models', 'pre_trained_models', 'ResNet-50L', 'mean_image');

end

author: feichtenhofer

repo: Detect-Track


text:
function Iout=scorfilt(I)
% SPATIAL CORRELATION FILTER
%
% IN
% I:        Image to filter, (RGB color or grayscale)
% %
% OUT
% Iout:     Filtered image (double)
%
% EXAMPLE:
%
% %read image
% imfile = fullfile(matlabroot,...
%          'toolbox','images','imdemos','greens.jpg');
% Io=im2double(imread(imfile));
% %add noise
% sigma_n=20; %noise std. dev.
% randn('seed', 0);
% I=Io+sigma_n/255*randn(size(Io));
% figure('Name', 'input'), imshow(I);
% %filter image
% Iout = scorfilt(I);
% figure('Name', 'output'), imshow(Iout);
%
% Author:   Carlos Estrada
% Date:     May 25, 2012
% Version   1.0
%
%% Configuration
r_max=5;    %max. half-size window
c=1/5;      %c 'neighborhood' constant default value

%% labels
I=im2double(I);
[m n ch]=size(I);

if ch==1           %gray scale
   maxI=max(I(:));
   minI=min(I(:));
   bins=255;   %quantization bins
   Ieq=uint8((bins-1)*histeq((I-minI)./(maxI-minI), bins))+1;

   labels=zeros(m,n,'uint8');
   N=0;
   bin_mark=zeros(bins,1,'uint8');

   for i=1:m
       for j=1:n
          id=Ieq(i,j);
         if bin_mark(id)==0
            N=N+1;
            labels(i,j)=N;
            bin_mark(id)=N;
         else
            labels(i,j)=bin_mark(id);
         end
       end
   end

elseif ch==3       %color
%    maxI=max(I(:));
%    minI=min(I(:));
%    Ifit=(I-minI)/(maxI-minI);
   bins=255;   %quantization bins
   [labels, centers] = rgb2ind(I,bins,'nodither');
   N=size(centers,1);
   labels=labels+1;
end

%% adjacency matrix
Acum=zeros(N,'uint32');
for i=1:m
    for j=1:n
       a=labels(i,j);
       if i>1
         e=labels(i-1,j);
         Acum(a,e)=Acum(a,e)+1;
       end
       if j>1
         e=labels(i,j-1);
         Acum(a,e)=Acum(a,e)+1;
       end
    end
end
A=double(Acum+Acum');

h4=sum(A,2);
P=A./repmat(h4,1,N); %transition matrix
Q=((P+P')./2)^2;
%% neighborhood size
teta=diag(P)'*h4./sum(h4);
w=sqrt((1-teta)/teta);
sigma_d=c*w;
r=max(1,min(r_max,round(sigma_d*2.5)));
%fprintf('\tN: %d\tw: %.4f\tsigma_d: %.4f', N,w,sigma_d);
%% compute Iout
Iout=zeros(m,n,ch);
[x,y] = meshgrid(-r:r,-r:r);
F=exp(-(x.^2 + y.^2)./(2*sigma_d^2));

for i=1:m
    iMin = max(i-r,1);
    iMax = min(i+r,m);
    di=iMin:iMax;
    for j=1:n
       jMin = max(j-r,1);
       jMax = min(j+r,n);
       dj=jMin:jMax;

       a=labels(i,j);
       l = labels(di,dj);
       f=F(di-i+r+1,dj-j+r+1);

       w=Q(l(:),a).*f(:);
       w=w./sum(w);

       y=I(di,dj,:);

       for cc=1:ch
         yc=y(:,:,cc);
         Iout(i,j,cc)=w'*yc(:);
       end
    end
end


text:
function opt=varargin2struct(varargin)
%
% opt=varargin2struct('param1',value1,'param2',value2,...)
%   or
% opt=varargin2struct(...,optstruct,...)
%
% convert a series of input parameters into a structure
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2012/12/22
%
% input:
%      'param', value: the input parameters should be pairs of a string and a value
%       optstruct: if a parameter is a struct, the fields will be merged to the output struct
%
% output:
%      opt: a struct where opt.param1=value1, opt.param2=value2 ...
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

len=length(varargin);
opt=struct;
if(len==0) return; end
i=1;
while(i<=len)
    if(isstruct(varargin{i}))
        opt=mergestruct(opt,varargin{i});
    elseif(ischar(varargin{i}) && i<len)
        opt=setfield(opt,lower(varargin{i}),varargin{i+1});
        i=i+1;
    else
        error('input must be in the form of ...,''name'',value,... pairs or structs');
    end
    i=i+1;
end


text:
function stats = signtest_matrix(dat)
% :Usage:
% ::
%
%     stats = signtest_matrix(dat)
%
% This is a matric-ized version of Matlab 2010's signtest.m
% it returns identical P-values to Matlab's function
%
% ..
%    copyright 2011, tor wager
% ..

n = size(dat, 1); % n: how many obs in sign test
n_omit = sum(dat == 0 | isnan(dat));
n = (n - n_omit)';

isvalid = n > 0;

dat(:, ~isvalid) = [];
n(~isvalid) = [];

npos = sum(dat > 0);
nneg = sum(dat < 0);
sgn = min(npos, nneg)';
mydirection = double(npos > nneg)';
mydirection(mydirection == 0) = -1;

if max(n) < 100
    method = 'exact';
else
    method = 'approximate';
end

if isequal(method,'exact')
    p = min(1, 2 .* binocdf(sgn, n, 0.5));  % p > 1 means center value double-counted
    zval = NaN .* ones(size(p));
else
    % Do a continuity correction, keeping in mind the right direction
    z = (npos-nneg - sign(npos-nneg)) ./ (n .^ .5);
    p = 2 .* normcdf(-abs(z), 0, 1);
    stats.zval = z;
    stats.p = p;
end

% FDR correction
% p = p(stats.isvalid); % already done
pthr = FDR(p, .05);
if isempty(pthr), pthr = Inf; end

stats.pthr_fdr05 = pthr;
sig = p < .05;
fdr_sig = p < pthr;

% Re-insert removed CASES (rows) and fill with zeros
mydirection = zeroinsert(~isvalid, mydirection);
n = zeroinsert(~isvalid, n);
p = zeroinsert(~isvalid, p);
zval = zeroinsert(~isvalid, zval);
sig = zeroinsert(~isvalid, sig);
fdr_sig = zeroinsert(~isvalid, fdr_sig);

stats = struct('method', method, 'isvalid', isvalid, 'n', n, ...
'direction', mydirection, 'p', p, 'zval', zval, 'pthr_fdr05', pthr, ...
'sig', sig, 'fdr_sig', fdr_sig);

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function blockplay(f)
%BLOCKPLAY Schedules block to be played
%   Usage: blockplay(L)
%
%   Input parameters:
%      f    : Samples.
%
%   Function schedules samples in *f* to be played. Since playrec handles
%   playing and recording in a single command, the actual relay of samples
%   to playrec is done in the next call of |blockread|.
%   In case no audio output is expected (in the rec only mode),
%   the function does nothing.

complainif_notenoughargs(nargin,1,'BLOCKPLAY');

source = block_interface('getSource');

if ( iscell(source) && strcmp(source{1},'rec')) || ...
   strcmp(source,'rec')
   % Do nothing in rec only mode.
   return;
   % error('%s: Blocks cannot be played in the rec only mode.',upper(mfilename));
end

% Reformat f if necessary
f = comp_sigreshape_pre(f,'BLOCKPLAY',0);

block_interface('setToPlay',f);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [] = aps_gacos_files(orderflag_GACOS_website)
% script that runs and checks which GACOS data files are needed based on
% the satellite pass time.
% INPUTS:
%
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% By David Bekaert - November 2017
% Modifications:

stdargin = nargin ;
if stdargin<1
    orderflag_GACOS_website = 0;
end

% getting the variables from the parms_aps file
workdir = pwd;
stamps_processed = getparm_aps('stamps_processed');
UTC_sat =  getparm_aps('UTC_sat');
gacos_datapath = getparm_aps('gacos_datapath');
datestructure = 'yyyymmdd';                               % assumed date structure for era
if isempty(gacos_datapath)
    error('please specify gacos_datapath')
end

% loading the data
if strcmp(stamps_processed,'y')
    ll_matfile = getparm_aps('ll_matfile');
    ps = load(ll_matfile);
    dates = ps.day;
else
    ifgday_matfile = getparm_aps('ifgday_matfile');
    ifgs_dates = load(ifgday_matfile);
    ifgs_dates = ifgs_dates.ifgday;
    dates = reshape(ifgs_dates,[],1);
    dates = unique(dates);
    dates = datenum(num2str(dates),'yyyymmdd');
end

% getting the dates
n_dates = length(dates);

% find two closest times with respect the 1hr spaced ERA5 data
time = str2num(UTC_sat(1:2)) + str2num(UTC_sat(4:5))/60;
fprintf(['Satellite pass is ' num2str(time) ' UTC \n'])
datelist =  datestr(dates,datestructure);
filelist = [datelist repmat('.ztd',size(datelist,1),1) ];

if orderflag_GACOS_website==0
    % outputing this information to a file
    fid = fopen('GACOS_files.txt','w');
    fprintf(['Required GACOS files for all interferograms \n'])
    for k=1:size(filelist,1)
        fprintf([filelist(k,:) '\n']);
        fprintf(fid,[filelist(k,:) '\n']);
    end
    fclose(fid);

    % getting the BBOX needed for GACO, in case user wants to download manually
    crop_range_in = 0.5; % increasing extent of weather data region by this value (degree)
    region_lat_range = getparm_aps('region_lat_range');
    region_lon_range = getparm_aps('region_lon_range');
    if isempty(region_lat_range) == 1
        error('Specify the region for the weather model data')
    end
    fprintf('increasing crop area by %s deg in each direction \n',num2str(crop_range_in))
    S = num2str(min(round(region_lat_range)) - crop_range_in);
    N = num2str(max(round(region_lat_range)) + crop_range_in);
    W = num2str(min(round(region_lon_range)) - crop_range_in);
    E = num2str(max(round(region_lon_range)) + crop_range_in);       % DB fixed typo min to max
    weatherregstr = [N,'/',W,'/',S,'/',E];   % N/W/S/E
    fprintf('weather model region N/W/S/E %s \n',weatherregstr);

    fid2 = fopen('GACOS_download_info.txt','w');
    % give the metadata needed
    fprintf(fid2,['UTC: ' UTC_sat '\n']);
    fprintf(fid2,['South: ' S '\n']);
    fprintf(fid2,['North: ' N '\n']);
    fprintf(fid2,['West: ' W '\n']);
    fprintf(fid2,['East: ' E '\n\n']);
    fprintf(fid2,['dates:\n']);
    for k=1:size(datelist,1)
        fprintf(fid2,[datelist(k,1:end) '\n']);
    end
    fclose(fid2);
end

%% Below is specific for the ECMWF data website
crop_range_in = 2; % increasing extent of weather data region by this value (degree)
overwrite_flag=-1;
if orderflag_GACOS_website==1

    error('GACOS website does not support api call')

    % weather model region
    fprintf('getting the data from the GACOS website ...')
    region_lat_range = getparm_aps('region_lat_range');
    region_lon_range = getparm_aps('region_lon_range');
    if isempty(region_lat_range) == 1
        error('Specify the region for the weather model data')
    end
    fprintf('increasing crop area by %s deg in each direction \n',num2str(crop_range_in))
    S = num2str(min(round(region_lat_range)) - crop_range_in);
    N = num2str(max(round(region_lat_range)) + crop_range_in);
    W = num2str(min(round(region_lon_range)) - crop_range_in);
    E = num2str(max(round(region_lon_range)) + crop_range_in);       % DB fixed typo min to max
    weatherregstr = [N,'/',W,'/',S,'/',E];   % N/W/S/E
    fprintf('weather model region N/W/S/E %s \n',weatherregstr);

    %%%%% CODE BELOW IS LEGACY CODE FROM ECMWF DOWNLOAD PROGRAM
    %%%%% TB CLEANED

    for l = 1:size(filelist,1)
        subdirpath = [gacos_datapath,'/',filelist(l,5:12),'/'];
        if exist(subdirpath,'dir') == 0
            fprintf('creating directory %s \n',subdirpath);
            mkdir(subdirpath)
        end
        if exist([subdirpath,'ggap', filelist(l,5:16) '.nc'],'file') == 0
            cd(subdirpath)
            fprintf('Order and downloading %s \n',filelist(l,5:16))
            aps_era5_ECMWF_Python(filelist(l,5:16),weatherregstr) %write python donwload file
            python_str = ['python ',filelist(l,5:16),'.py > ',filelist(l,5:16),'down.log &'];
            [a,b] = system(python_str); % start python script
            clear a b
            cd ..
        else
            if overwrite_flag==-1
                str='';
                while ~strcmpi(str,'y') && ~strcmpi(str,'n')
                    fprintf(['Do you want to overwrite existing files? \n'])
                    str = input('[y: for yes, n: no] \n','s');
                end
                if strcmpi(str,'n')
                    overwrite_flag=0;
                else
                    overwrite_flag=1;
                end
            end
            % check if the files need to be overwritten
            if overwrite_flag==1
                cd(subdirpath)
                delete([subdirpath,'ggap', filelist(l,5:16) '.nc'])
                fprintf('Order and downloading %s \n',filelist(l,5:16))
                aps_era5_ECMWF_Python(filelist(l,5:16),weatherregstr) %write python donwload file
                python_str = ['python ',filelist(l,5:16),'.py > ',filelist(l,5:16),'down.log &'];
                [a,b] = system(python_str); % start python script
                clear a b
                cd ..
            elseif overwrite_flag==0
                fprintf('File %s has already been downloaded \n',filelist(l,:))
            end
        end
    end
    cd(workdir)
end

% check if the files needs to be organized into a YYYYMMDD folder structure
try
    cd(gacos_datapath)
    for k=1:size(filelist,1)
       file = filelist(k,:);
        date = file(1:8);
        if exist(date,'dir')~=7
            mkdir(date)
        end
        if exist(file,'file')==2
            movefile([file '*'],[date filesep '.'])
        end

    end
    cd(workdir)

catch
    cd(workdir)
end

cd(workdir)     % back to work directory

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function varargout = fevalm(f, r, lambda, theta)
%FEVAL   Evaluate a BALLFUNV
%   FEVAL(F, R, L, T) evaluates a BALLFUNV F at the points (R,L,T) in spherical coordinates
%   at a tensor-product grid R x L x T.
%
% See also SUBSREF.

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

if isempty( f )
    varargout = {};
    return
end

F = f.comp;

% Evaluate at the tensor grid
valsX = fevalm(F{1}, r, lambda, theta);
valsY = fevalm(F{2}, r, lambda, theta);
valsZ = fevalm(F{3}, r, lambda, theta);

varargout = {valsX, valsY, valsZ};
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% function XF = fftBandpass(X, FS, FS1, FP1, FP2, FS2)
%
% Bandpass filter for the signal X (time x trials). An acuasal fft
% algorithm is applied (i.e. no phase shift). The filter functions is
% constructed from a Hamming window.
%
% Fs : sampling frequency
%
% The passbands (Fp1 Fp2) and stop bands (Fs1 Fs2) are defined as
%                 -----------
%                /           \
%               /             \
%              /               \
%             /                 \
%   ----------                   -----------------
%           Fs1  Fp1       Fp2  Fs2
%
% If no output arguments are assigned the filter function H(f) and
% impulse response are plotted.
%
% NOTE: for long data traces the filter is very slow.
%
%------------------------------------------------------------------------
% Ole Jensen, Brain Resarch Unit, Low Temperature Laboratory,
% Helsinki University of Technology, 02015 HUT, Finland,
% Report bugs to ojensen@neuro.hut.fi
%------------------------------------------------------------------------

%    Copyright (C) 2000 by Ole Jensen
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You can find a copy of the GNU General Public License
%    along with this package (4DToolbox); if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
function xf = fftBandpass(x,Fs,Fs1,Fp1,Fp2,Fs2)
    if size(x,1) == 1
        x = x';
    end
    % Make x even
    Norig = size(x,1);
    if rem(Norig,2)
        x = [x' zeros(size(x,2),1)]';
    end

    % Normalize frequencies
    Ns1 = Fs1/(Fs/2);
    Ns2 = Fs2/(Fs/2);
    Np1 = Fp1/(Fs/2);
    Np2 = Fp2/(Fs/2);

    % Construct the filter function H(f)
    N = size(x,1);
    Nh = N/2;

    B = fir2(N-1,[0 Ns1 Np1 Np2 Ns2 1],[0 0 1 1 0 0]);
    H = abs(fft(B));  % Make zero-phase filter function
    IPR = real(ifft(H));
    if nargout == 0
        subplot(2,1,1)
        f = Fs*(0:Nh-1)/(N);
        plot(f,H(1:Nh));
        xlim([0 2*Fs2])
        ylim([0 1]);
        title('Filter function H(f)')
        xlabel('Frequency (Hz)')
        subplot(2,1,2)
        plot((1:Nh)/Fs,IPR(1:Nh))
        xlim([0 2/Fp1])
        xlabel('Time (sec)')
        ylim([min(IPR) max(IPR)])
        title('Impulse response')
    end

    if size(x,2) > 1
        for k=1:size(x,2)
            xf(:,k) = real(ifft(fft(x(:,k)) .* H'));
        end
        xf = xf(1:Norig,:);
    else
        xf = real(ifft(fft(x') .* H));
        xf = xf(1:Norig);
    end
end

author: kavli-ntnu


text:
function spy( prob, reduce )
if nargin < 2 || ~reduce,
    A = extract( prob );
else
    A = eliminate( prob );
end
spy( A' );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function R = fit_rotation(S)
  % FIT_ROTATION find rotation for a given covariance matrix
  %
  % R = fit_rotation(S)
  %
  % Inputs:
  %   S  n by n covariance matrix
  % Outputs:
  %   R  n by n rotation matrix closest to S
  %

  % svd
  [su,ss,sv]=svd(S);
  R = sv*su';
  % if reflection then flip last column
  if( det(R) < 0 )
    su(:,end) = -su(:,end);
    R = sv*su';
  end
  % should definitely be rotation now
  %assert( det(R) >= 0 );
end

author: alecjacobson

repo: gptoolbox


text:
function test_suite=test_moxunit_util_escape_xml
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_moxunit_util_escape_xml_basics
    aeq=@(a,b)assertEqual(moxunit_util_escape_xml(a),b);

    aeq('asdk','asdk');
    aeq('&','&');
    aeq('"','"');
    aeq('''',''');
    aeq('<','<');
    aeq('>','>');

    aeq('apos&&&apos','apos&&&apos');

    aeq('&amp&','&amp&');
    aeq('&quot','&quot');
    aeq('&apos','&apos');
    aeq('&lt','&lt');
    aeq('&gt','&gt');

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function mnt = mnt_getNIRSMontage(varargin)
%MNT_GETNIRSMONTAGE - Get montage with NIRS channel positions
%
%Description:
% get montage with NIRS channel positions for a given set of sources and
% detectors. NIRS channels are placed half-way between a source and a
% detector on a spherical head model.
%
%Synopsis:
% MNT = nirs_getMontage(CLABSOURCE,CLABDETECTOR,<OPT>)
%       If you give channels labels of the source and detectors (should
%       correspond to EEG channel labels), montages for source and
%       detectors are produced, as well as montages for the
%       source-detector combinations.
%
% MNT = nirs_getMontage(MNT,<OPT>)
%       You can provide a montage struct wherein source and detector
%       montages are already given in .source and .detector fields.
%       In this case, only the source-detector combinations are determined.
%
%Input
% OPT: struct or property/value list of optional properties:
%  .File       - montage file used for channel positions
%  .ClabPolicy - specifies how the source-detector combinations are
%                labeled. 'Label' (default) concatenates the labels of
%                source and detector (eg. FzCz), 'Number' concatenates
%                the channel numbers (eg. 112)
%  .Connector  - specify how source and detectors labels are connected
%                (default ''). To have no connector, set to ''.
%  .Projection - determine projection of 3D coordinates to 2D
%                'orthogonal' or 'euclidean' (default)
%
%Output:
% MNT: NIRS montage
%  .x           - x coordinate of channel positions (2d projection)
%  .y           - y coordinate of channel positions (2d projection)
%  .pos_3d      - 3d positions on spherical head model
%  .clab        - channel labels
%  .angulardist - angular distances between sources and detectors on the
%                 spherical model head (in radians)
%  .sd          - source and detector index to which eg the channel belongs [1 12]
%  .source      - struct with montage for the sources
%  .detector    - struct with montage for the detectors
%
%See also: setElectrodeMontage, nirs_reduceMontage

% matthias.treder@tu-berlin.de 2011
% Markus Wenzel 2013 (adapted it to the new toolbox)
% Jan Mehnert February 2014 (ready for public BBCI toolbox) (jan@mehnert.org)

props= {'ClabPolicy'	'label'             'CHAR'
        'Projection'	'euclidean'         'CHAR'
        'PositionFcn'   @mntutil_posExt55   'FUNC'
        'Connector'     ''                  'CHAR'
       };

if nargin==0,
    mnt= props; return
end

if isstruct(varargin{1})
    mnt = varargin{1};
    opt= opt_proplistToStruct(varargin{2:end});
    misc_checkType(mnt, 'STRUCT'); %(redundant)
else
    clabSource = varargin{1};
    clabDetector = varargin{2};
    opt= opt_proplistToStruct(varargin{3:end});
    misc_checkType(clabSource, 'CHAR|CELL{CHAR}');
    misc_checkType(clabDetector, 'CHAR|CELL{CHAR}');
end

[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

if ~exist('mnt','var')
  % Determine source and detector montages
  mnt = struct();
  mnt.source = mnt_setElectrodePositions(clabSource, 'PositionFcn', opt.PositionFcn);
  mnt.detector = mnt_setElectrodePositions(clabDetector, 'PositionFcn', opt.PositionFcn);
end

if any(isnan(mnt.source.pos_3d(:))) || any(isnan(mnt.detector.pos_3d(:)))
  warning(['Some source/detector positions are NaN. Note that the NIRS channels '...
    'involving these sources/detectors will also have NaN positions.\n '])
end

% Determine positions of source-detector combinations
mnt.x = [];
mnt.y = [];
mnt.pos_3d = [];
mnt.clab = {};
mnt.angulardist = [];

chanIdx = 1;
for ss=1:numel(mnt.source.clab)
  for dd=1:numel(mnt.detector.clab)
    ps = mnt.source.pos_3d(:,ss);
    pd = mnt.detector.pos_3d(:,dd);
    % 3d position of new channel (half-way between source and detector)
    % Set to norm=1 so that it's located on the head surface
    newpos = (ps + pd) / norm(ps+pd);
    % Angular distance in radians = angle between the two vectors1
    costheta = dot(ps,pd);
    theta = acos(costheta);
    % label
    if strcmp(opt.ClabPolicy,'label')
      label = [mnt.source.clab{ss} opt.Connector mnt.detector.clab{dd}];
    elseif strcmp(opt.ClabPolicy,'Number')
      label = [num2str(ss) opt.Connector num2str(dd)];
    end
    % Add to struct
    mnt.pos_3d(:,chanIdx) = newpos;
    mnt.angulardist = [mnt.angulardist theta];
    mnt.clab = {mnt.clab{:} label};
    mnt.sd(chanIdx,:) = [ss,dd];
    chanIdx = chanIdx + 1;
  end
end
mnt = mntutil_projectNIRSChannels3dTo2d(mnt,opt.Projection);

author: bbci


text:
function [coeff,variance]=OLSPriorTheta(y,timetrend,f)
% OLSPriorTheta(y,timetrend,f) computes the prior for coefficients of the
% deterministic trends in the TVE-VAR approach. The prior is normally
% distributed with mean "coeff" and variance "variance". The prior
% parameters are computed regressing the data y on timetrend. Then, coeff is
% the OLS estimate, and variance is the asymptotic variance of the OLS,
% scaled up of a factor f.

%% Prepare data: make sure the dimensions of the inputs are correct
if size(y,1)>size(y,2)
    y=y';
end
if size(timetrend,1)>size(timetrend,2)
    timetrend=timetrend';
end

%% Run OLS
T=size(y,2); % number of data available
k=size(timetrend,1); % number of regressors

P=y*timetrend'*inv(timetrend*timetrend'); % the OLS coefficients
resid=y-P*timetrend; % the residuals from the regression
O=1/(T-k)*resid*resid';
Q=timetrend*timetrend';
varpar=kron(O,Q^-1); % the asymptotic variance

%% The output
variance=varpar*f;
coeff=P';

end

author: european-central-bank

repo: BEAR-toolbox


text:
%########################################################################
%
%	- PPGI Toolbox -
%   A MATLAB toolbox for Photoplethysmography Imaging (PPGI)
%
% Author   : Christian S. Pilz
% Company  : The Nature of Space of Time
% Date     : 07.05.2019
%
% Contact  : cpi@partofthestars.com
% Web Page : www.partofthestars.com
%
% Version  : beta0.1
%
%########################################################################
%
%	test_spherical_mean.m:
%
% Description:
%
%   test of the spherical mean feature extraction on given sample rgb data
%
%
clear all;
close all;

load('./../media/data/example_data.mat');

if ~exist('skin_pixels')
    disp('error: no skin pixels available. execute test_skin.m first!');
    return;
end

spm=spherical_mean();

for f=1:size(skin_pixels,2)
    f
    [signal(f,:) ssr]=spm.get(skin_pixels{f});
end

fs=25;
low_frequency=0.5;
high_frequency=2.5;
bpf=bandpass_filter(fs,low_frequency,high_frequency);
signal_filtered=bpf.get(signal);

[pearson, rmse, snr, snr_var, bpm] = ground_truth_stats.get(ppg,signal_filtered(:,1),fs);

author: partofthestars

repo: PPGI-Toolbox


text:
function [candidates, candScores, alignInfo, otherInfo] = lstmDecoder(models, data, params)
% Decode from an LSTM model.
%   stackSize: the maximum number of translations we want to get.
% Output:
%   - candidates: list of candidates
%   - candScores: score of the corresponding candidates (stackSize * batchSize)
%
% Thang Luong @ 2015, <lmthang@stanford.edu>
%   With help from Hieu Pham.

  % backward compatibility: not a cell, single model, put into a cell format
  if ~iscell(models)
    tmpModels = cell(1, 1);
    tmpModels{1} = models;
    tmpModels{1}.params = params;
    models = tmpModels;
  end

  beamSize = params.beamSize;
  stackSize = params.stackSize;
  batchSize = size(data.srcInput, 1);

  srcMaxLen = data.srcMaxLen;
  params.srcMaxLen = srcMaxLen;

  %% init
  minLen = floor(data.srcMinLen*params.minLenRatio);
  if params.forceDecoder
    maxLen = data.tgtMaxLen;
  else
    maxLen = floor(srcMaxLen*params.maxLenRatio);
  end
  fprintf(2, '# Decoding batch of %d sents, minLen=%d, maxLen=%d, tgtEos=%d, %s\n', batchSize, minLen, maxLen, params.tgtEos, datestr(now));
  fprintf(params.logId, '# Decoding batch of %d sents, minLen=%d, maxLen=%d, tgtEos=%d, %s\n', batchSize, minLen, maxLen, params.tgtEos, datestr(now));

  startTime = clock;

  %%%%%%%%%%%%
  %% encode %%
  %%%%%%%%%%%%
  %% multiple models
  numModels = length(models);
  modelData = cell(numModels, 1);
  zeroStates = cell(numModels, 1);
  for mm=1:numModels
    models{mm}.params.curBatchSize = batchSize;
    models{mm}.params.srcMaxLen = srcMaxLen;
    [models{mm}.params] = setAttnParams(models{mm}.params);

    [zeroStates{mm}] = createZeroState(models{mm}.params);
    modelData{mm} = data;

    % compatible
    [models{mm}.params] = backwardCompatible(models{mm}.params, {'unkDiscount'}, 0);
  end

  % encoder
  prevStates = cell(numModels, 1);
  for mm=1:numModels
    [~, prevStates{mm}, ~, modelData{mm}, ~] = encoderLayerForward(models{mm}, zeroStates{mm}, modelData{mm}, models{mm}.params, 1);
  end

  %%%%%%%%%%%%
  %% decode %%
  %%%%%%%%%%%%
  [candidates, candScores, alignInfo, otherInfo] = rnnDecoder(models, params, prevStates, minLen, maxLen, beamSize, stackSize, batchSize, ...
  modelData, data, params.tgtEos, 0);

  % char: now generate words for <unk>
  if params.charTgtGen % assume all models are char-level models
    assert(stackSize == 1); % TODO: remove this, then char_initEmb needs to be initialized differently.
    ss = 1;

    % prepare starting states
    char_initEmb = cell(numModels, 1);
    for mm=1:numModels
      char_initEmb{mm} = zeroMatrix([params.lstmSize, batchSize * maxLen], params.isGPU, params.dataType);
    end

    numRareWords = 0;
    otherInfo.rarePositions = cell(batchSize, 1);
    for sentId=1:batchSize
      flags = candidates{sentId}{ss} == params.tgtUnk;
      positions = find(flags);

      for mm=1:numModels
        char_initEmb{mm}(:, numRareWords+1:numRareWords+length(positions)) = otherInfo.transStates{ss, mm}(:, positions, sentId);
        % assert(isempty(find(otherInfo.transStates{ss, mm}(:, length(flags), sentId) ~= 0, 1)));
      end

      numRareWords = numRareWords + length(positions);
      otherInfo.rarePositions{sentId} = positions;
    end

    for mm=1:numModels
      char_initEmb{mm}(:, numRareWords+1:end) = [];
    end

    otherInfo.numRareWords = numRareWords;
    otherInfo.rareWords = cell(batchSize, 1);
    if numRareWords
      % char params
      charParams = params;
      charParams.curBatchSize = numRareWords;

      % model
      charModels = cell(numModels, 1);
      for mm=1:numModels
        charModels{mm}.W_tgt = models{mm}.W_tgt_char;
        charModels{mm}.W_emb_tgt = models{mm}.W_emb_tgt_char;
        charModels{mm}.W_soft = models{mm}.W_soft_char;

        charParams.numLayers = models{mm}.params.charNumLayers;
        charParams.tgtVocab = models{mm}.params.tgtCharVocab;

        charModels{mm}.params = charParams;
      end

      % prev states
      char_minLen = 2;
      char_maxLen = 20;
      fprintf(2, '  # Decoding char , minLen=%d, maxLen=%d, charEos=%d, %s\n', char_minLen, char_maxLen, params.tgtCharEos, datestr(now));
      fprintf(params.logId, '  # Decoding char , minLen=%d, maxLen=%d, charEos=%d, %s\n', char_minLen, char_maxLen, params.tgtCharEos, datestr(now));

      char_modelData = cell(numModels, 1);
      char_prevStates = cell(numModels, 1);
      zeroBatch = zeroMatrix([params.lstmSize, charParams.curBatchSize], params.isGPU, params.dataType);

      for mm=1:numModels
        char_prevStates{mm} = cell(charModels{mm}.params.numLayers, 1);
        for ll=1:charModels{mm}.params.numLayers % layer
          if ll == 1
            char_prevStates{mm}{ll}.h_t = char_initEmb{mm}(:, :);
          else
            char_prevStates{mm}{ll}.h_t = zeroBatch;
          end
          char_prevStates{mm}{ll}.c_t = zeroBatch;
        end
      end

      [char_candidates, ~, ~, ~] = rnnDecoder(charModels, charParams, char_prevStates, char_minLen, ...
        char_maxLen, beamSize, stackSize, numRareWords, char_modelData, [], params.tgtCharEos, 1);

      % get words
      count = 0;
      for sentId=1:batchSize
        numCharWords = length(otherInfo.rarePositions{sentId});
        otherInfo.rareWords{sentId} = cell(numCharWords, 1);
        for jj=1:numCharWords
          otherInfo.rareWords{sentId}{jj} = [params.tgtCharVocab{char_candidates{count+jj}{1}(1:end-1)}];
        end
        count = count + numCharWords;
      end
    end
  end

  endTime = clock;
  timeElapsed = etime(endTime, startTime);
  fprintf(2, '  Done, minLen=%d, maxLen=%d, speed %f sents/s, time %.0fs, %s\n', minLen, maxLen, batchSize/timeElapsed, timeElapsed, datestr(now));
  fprintf(params.logId, '  Done, minLen=%d, maxLen=%d, speed %f sents/s, time %.0fs, %s\n', minLen, maxLen, batchSize/timeElapsed, timeElapsed, datestr(now));
end

%%
% Take the initial states, starting symbol, and then decode!
% This method should be reusable.
%%
function [candidates, candScores, alignInfo, otherInfo] = rnnDecoder(models, params, prevStates, minLen, maxLen, beamSize, stackSize, ...
  batchSize, modelData, data, tgtEos, isChar)
  numModels = length(models);

  % first decoder timestep
  attnInfos = cell(numModels, 1);
  firstAlignIdx = [];
  for mm=1:numModels
    % char
    if isChar
      isFeedInput = 0;
      attnFunc = 0;
    else
      isFeedInput = models{mm}.params.feedInput;
      attnFunc = models{mm}.params.attnFunc;
    end
    initEmb = models{mm}.W_emb_tgt(:, repmat(models{mm}.params.tgtSos, batchSize, 1));
    decRnnFlags = struct('decode', 1, 'test', 1, 'attn', attnFunc, 'feedInput', isFeedInput);
    [prevStates{mm}, attnInfos{mm}] = rnnStepLayerForward(models{mm}.W_tgt, initEmb, ...
      prevStates{mm}, ones(batchSize, 1), models{mm}.params, decRnnFlags, modelData{mm}, models{mm});
  end

  % output alignment
  if params.align && isChar == 0
    [~, firstAlignIdx] = getAlignWeights(attnInfos, data.srcLens, models, params);
  end

  if isChar
    sentIndices = 1:batchSize;
  else
    sentIndices = data.startId:(data.startId+batchSize-1);
  end
  [candidates, candScores, alignInfo, otherInfo] = decodeBatch(models, params, prevStates, attnInfos, minLen, maxLen, beamSize, stackSize, batchSize, ...
    sentIndices, modelData, firstAlignIdx, data, tgtEos, isChar);
end

%%
% Beam decoder from an LSTM model, works for multiple sentences
% Input:
%   - encoded vector of the source sentences
%   - maximum length willing to go
%   - beamSize
%   - stackSize: maximum number of translations collected for one example
%%
function [candidates, candScores, alignInfo, otherInfo] = decodeBatch(models, params, prevStates, attnInfos, minLen, maxLen, beamSize, stackSize, batchSize, ...
originalSentIndices, modelData, firstAlignIdx, data, tgtEos, isChar)
  numElements = batchSize*beamSize;

  candidates = cell(batchSize, 1);
  candScores = -1e10*oneMatrix([stackSize, batchSize], params.isGPU, params.dataType); % set to a very small value
  numDecoded = zeros(batchSize, 1);
  for ii=1:batchSize
    candidates{ii} = cell(stackSize, 1);
  end

  numModels = length(models);
  otherInfo = [];

  % align
  if params.align && isChar == 0
    alignInfo = cell(batchSize, 1);
    for ii=1:batchSize
      alignInfo{ii} = cell(stackSize, 1);
    end
  else
    alignInfo = [];
  end

  %% first prediction
  % scores, words: beamSize * batchSize
  if params.forceDecoder
    assert(isChar == 0);
    [scores, words, otherData] = nextBeamStep(models, prevStates, beamSize, 0, isChar, params.tgtUnk, params.unkDiscount, data.tgtOutput(:, 1));
    otherInfo.forceDecodeOutputs = zeroMatrix([maxLen, numElements], params.isGPU, params.dataType); % maxLen * (numElements)
    otherInfo.forceDecodeOutputs(1, :) = otherData.maxWords;
  else
    [scores, words] = nextBeamStep(models, prevStates, beamSize, tgtEos, isChar, params.tgtUnk, params.unkDiscount);
  end

  % TODO: by right, we should filter out words == tgtEos, but I
  % think for good models, we don't have to worry :)

  %% matrix dimension note
  % note that we order matrices in the following dimension: n * numElements
  % columns correspond to the 1st sent go first, then the 2nd one, until the batchSize-th sent.
  sentIndices = repmat(1:batchSize, [beamSize, 1]);
  sentIndices = sentIndices(:)'; % 1 ... 1, 2 ... 2, ...., batchSize ... batchSize . 1 * (beamSize*batchSize)

  %% init beam
  beamScores = scores(:)'; % 1 * numElements
  beamHistory = zeroMatrix([maxLen, numElements], params.isGPU, params.dataType); % maxLen * (numElements)
  beamHistory(1, :) = words(:); % words for sent 1 go together, then sent 2, ...
  if params.align && isChar == 0
    alignHistory = zeroMatrix([maxLen, numElements], params.isGPU, params.dataType); % maxLen * (numElements)
    for ii=1:batchSize
      alignHistory(1,(ii-1)*beamSize+1:ii*beamSize) = firstAlignIdx(ii);
    end
  end

  %% replicate
  beamStates = cell(numModels, 1);
  for mm=1:numModels % model
    beamStates{mm} = cell(models{mm}.params.numLayers, 1);
    for ll=1:models{mm}.params.numLayers % lstmSize * numElements
      beamStates{mm}{ll}.c_t = reshape(repmat(prevStates{mm}{ll}.c_t, beamSize, 1),  models{mm}.params.lstmSize, numElements);
      beamStates{mm}{ll}.h_t = reshape(repmat(prevStates{mm}{ll}.h_t, beamSize, 1),  models{mm}.params.lstmSize, numElements);
    end
    beamStates{mm}{end}.softmax_h = reshape(repmat(prevStates{mm}{end}.softmax_h, beamSize, 1),  models{mm}.params.lstmSize, numElements);
  end

  %% char
  if params.charTgtGen && isChar == 0
    % assert(stackSize == 1);
    beamHistTopStates = cell(numModels, 1);
    otherInfo.transStates = cell(stackSize, numModels);
    for ss=1:stackSize
      for mm=1:numModels
        % sent 1, ..., sent2, ..., sent_batchSize
        otherInfo.transStates{ss, mm} = zeroMatrix([models{mm}.params.lstmSize, maxLen-1, batchSize], params.isGPU, params.dataType);
      end
    end
    for mm=1:numModels
      assert(models{mm}.params.charTgtGen == 1);
      beamHistTopStates{mm} = zeroMatrix([models{mm}.params.lstmSize, numElements, maxLen], params.isGPU, params.dataType);
      if models{mm}.params.charFeedOpt
        tgtWordStates = hiddenLayerForward(models{mm}.W_h_char, attnInfos{mm}.input, params.nonlinear_f);
        beamHistTopStates{mm}(:, :, 1) = reshape(repmat(tgtWordStates, beamSize, 1),  models{mm}.params.lstmSize, numElements);
      else
        beamHistTopStates{mm}(:, :, 1) = beamStates{mm}{end}.softmax_h;
      end
    end
  end

  oneMask = ones(1, numElements);

  if isChar == 0
    data.srcLens = reshape(repmat(data.srcLens, beamSize, 1), 1, []);

    % attentional / positional models
    for mm=1:numModels % model
      if models{mm}.params.attnFunc
        models{mm}.params.curBatchSize = numElements;

        modelData{mm}.curBatchSize = numElements;
        modelData{mm}.srcLens = reshape(repmat(modelData{mm}.srcLens, beamSize, 1), 1, []);

        % duplicate srcHidVecs
        modelData{mm}.srcHidVecsOrig = duplicateSrcHidVecs(modelData{mm}.srcHidVecsOrig, batchSize, beamSize);

        %% replicate srcMask
        % first transpose to srcMaxLen * batchSize
        % replicate (srcMaxLen*beamSize)*batchSize
        % reshape maxLen*(beamSize*batchSize)
        % then transpose back to (beamSize*batchSize)*maxLen
        srcMaxLen = size(modelData{mm}.srcMask, 2);
        modelData{mm}.srcMask = reshape(repmat(modelData{mm}.srcMask', beamSize, 1), srcMaxLen, beamSize*batchSize)';
      end
    end
  end

  % decodeCompleteCount = 0; % count how many sentences we have completed collecting the translations
  beamWords = zeroMatrix([1, numElements], params.isGPU, params.dataType);
  beamIndices = zeroMatrix([1, numElements], params.isGPU, params.dataType);
  % align
  if params.align && isChar == 0
    beamAlignIds = zeroMatrix([1, numElements], params.isGPU, params.dataType);
  end

  if beamSize == 1 % useful for force decoding
    doneFlags = zeros(1, batchSize); % mark if we have finished decoding a sentence
  end

  attnInfos = cell(numModels, 1);
  for sentPos=1:(maxLen-1)
    %% Description:
    % At this point, hypotheses of length sentPos are completed.
    % This loop will prepare hypotheses of length(sentPos+1) by:
    %   (a) first, finding out what the top beamSize^2 nextWords are.
    %   (b) among these nextWords, those which are equal to <eos>, will signal
    %      a complete translations and we will collect. These complete
    %      translations have length (sentPos+1), inclusive of <eos>.
    %   (c) we keep beamSize non-eos nextWords to build hypotheses of length (sentPos+1).
    tgtPos = sentPos+1;

    %% compute next lstm hidden states
    for mm=1:numModels
      % char
      if isChar
        isFeedInput = 0;
        attnFunc = 0;
      else
        isFeedInput = models{mm}.params.feedInput;
        attnFunc = models{mm}.params.attnFunc;
      end

      decRnnFlags = struct('decode', 1, 'test', 1, 'attn', attnFunc, 'feedInput', isFeedInput);
      [beamStates{mm}, attnInfos{mm}] = rnnStepLayerForward(models{mm}.W_tgt, models{mm}.W_emb_tgt(:, beamHistory(sentPos, :)), beamStates{mm}, ...
        oneMask, models{mm}.params, decRnnFlags, modelData{mm}, models{mm});
    end

    %% output alignment
    if params.align && isChar == 0
      [~, alignIdx] = getAlignWeights(attnInfos, data.srcLens, models, params);

      % we want to mimic the structure of allBestWords
      % size (beamSize * beamSize) x batchSize
      alignIdx = reshape(repmat(alignIdx, beamSize, 1), [], batchSize);
    end

    %% predict the next word
    % allBestScores, allBestWords should have size beamSize * (beamSize*batchSize)
    if params.forceDecoder
      assert(isChar == 0);
      [allBestScores, allBestWords, otherData] = nextBeamStep(models, beamStates, beamSize, 0, isChar, params.tgtUnk, params.unkDiscount, data.tgtOutput(:, tgtPos));
      otherInfo.forceDecodeOutputs(tgtPos, :) = otherData.maxWords;
    else
      [allBestScores, allBestWords] = nextBeamStep(models, beamStates, beamSize, 0, isChar, params.tgtUnk, params.unkDiscount);
    end

    % allBestWords, allBestScores should have size: (beamSize*beamSize) * batchSize
    [allBestScores, allBestWords, indices] = addSortScores(allBestScores, allBestWords, beamScores, beamSize, batchSize);

    %% build new beam
    for sentId=1:batchSize
      rowIndices = indices(:, sentId)';
      bestWords = allBestWords(rowIndices, sentId);

      % get candidates
      selectedIndices = find(bestWords~=tgtEos, beamSize);
      if ~isempty(selectedIndices)
        % update words
        startId = (sentId-1)*beamSize+1;
        endId = sentId*beamSize;
        sentWords = bestWords(selectedIndices);
        beamWords(startId:endId) = sentWords;

        % align
        if params.align && isChar == 0
          bestAlignIds = alignIdx(rowIndices, sentId);
          beamAlignIds(startId:endId) = bestAlignIds(selectedIndices);
        end

        % update beam
        sentBeamIndices = floor((rowIndices(selectedIndices)-1)/beamSize) + 1;
        beamIndices(startId:endId) = sentBeamIndices;

        % update scores
        beamScores(startId:endId) = allBestScores(selectedIndices, sentId);

        % get words that are eos and ranked before the last hypothesis in the next beam
        endIndices = find(bestWords(1:selectedIndices(end))==tgtEos);
      else % special case, happen for beamSize = 1, useful for force decoding
        assert(beamSize == 1);
        assert(bestWords == tgtEos);
        if doneFlags(sentId) == 1
          continue;
        else
          doneFlags(sentId) = 1;
          endIndices = 1;
        end
      end

      %% store translations
      if ~isempty(endIndices) && (sentPos+1)>=minLen % we don't want to start recording very short translations
        numTranslations = length(endIndices);
        eosBeamIndices = floor((rowIndices(endIndices)-1)/beamSize) + 1;
        histIndices = (sentId-1)*beamSize + eosBeamIndices;

        transScores = allBestScores(endIndices, sentId);
        for ii=1:numTranslations
          if numDecoded(sentId)<stackSize || transScores(ii) > candScores(stackSize, sentId) % haven't collected enough translations or better than the last
            if numDecoded(sentId)<stackSize
              numDecoded(sentId) = numDecoded(sentId) + 1;
            end

            candScores(numDecoded(sentId), sentId) = transScores(ii);
            candidates{sentId}{numDecoded(sentId)} = [beamHistory(1:sentPos, histIndices(ii)); tgtEos];

            % align
            if params.align && isChar == 0
              alignInfo{sentId}{numDecoded(sentId)} = [alignHistory(1:sentPos, histIndices(ii)); beamAlignIds(histIndices(ii))];
            end

            % char
            if params.charTgtGen && isChar == 0
              for mm=1:numModels % copy from existing states
                otherInfo.transStates{numDecoded(sentId), mm}(:, 1:sentPos, sentId) = beamHistTopStates{mm}(:, histIndices(ii), 1:sentPos);
              end
            end
          end
        end
      end
    end

    %% update history
    % overwrite previous history
    colIndices = (sentIndices-1)*beamSize + beamIndices;
    beamHistory(1:sentPos, :) = beamHistory(1:sentPos, colIndices);
    beamHistory(sentPos+1, :) = beamWords;
    % char
    if params.charTgtGen && isChar == 0
      for mm=1:numModels
        beamHistTopStates{mm}(:, :, 1:sentPos) = beamHistTopStates{mm}(:, colIndices, 1:sentPos);
        if models{mm}.params.charFeedOpt
          beamHistTopStates{mm}(:, :, sentPos+1) = hiddenLayerForward(models{mm}.W_h_char, attnInfos{mm}.input(:, colIndices), params.nonlinear_f);
        else
          beamHistTopStates{mm}(:, :, sentPos+1) = beamStates{mm}{end}.softmax_h(:, colIndices);
        end
      end
    end

    % align
    if params.align && isChar == 0
      alignHistory(1:sentPos, :) = alignHistory(1:sentPos, colIndices);
      alignHistory(sentPos+1, :) = beamAlignIds;
    end

    %% update lstm states
    for mm=1:numModels % model
      for ll=1:models{mm}.params.numLayers
        % lstmSize * (numElements): h_t and c_t vectors of each sent are arranged near each other
        beamStates{mm}{ll}.c_t = beamStates{mm}{ll}.c_t(:, colIndices);
        beamStates{mm}{ll}.h_t = beamStates{mm}{ll}.h_t(:, colIndices);
      end
      beamStates{mm}{end}.softmax_h = beamStates{mm}{end}.softmax_h(:, colIndices);
    end
  end % for sentPos

  for sentId=1:batchSize
    if numDecoded(sentId) == 0 % no translations found, use what we have in the beam
      fprintf(2, '  ! Sent %d: no translations end in %s, appending.\n', originalSentIndices(sentId), params.tgtVocab{tgtEos});
      for bb = 1:beamSize
        eosIndex = (sentId-1)*beamSize + bb;

        if numDecoded(sentId)<stackSize || beamScores(eosIndex) > candScores(stackSize, sentId) % haven't collected enough translations or better than the last
          if numDecoded(sentId)<stackSize
            numDecoded(sentId) = numDecoded(sentId) + 1;
          end

          candidates{sentId}{numDecoded(sentId)} = [beamHistory(1:maxLen, eosIndex); tgtEos]; % append eos at the end
          candScores(numDecoded(sentId), sentId) = beamScores(eosIndex);
          fprintf('  %.2f, num decoded %d:', candScores(numDecoded(sentId), sentId), numDecoded(sentId));
          fprintf(' %s', candidates{sentId}{numDecoded(sentId)});
          fprintf('\n');

          % align
          if params.align && isChar == 0
            if params.isReverse
              alignInfo{sentId}{numDecoded(sentId)} = [alignHistory(1:maxLen, eosIndex); 1];
            else
              alignInfo{sentId}{numDecoded(sentId)} = [alignHistory(1:maxLen, eosIndex); modelData{1}.srcLens(sentId)-1];
            end
          end

          % char
          if params.charTgtGen && isChar == 0
            for mm=1:numModels
              otherInfo.transStates{numDecoded(sentId), mm}(:, 1:maxLen, sentId) = beamHistTopStates{mm}(:, eosIndex, 1:maxLen);
            end
          end
        end
      end % for bb
    end
    candidates{sentId}(numDecoded(sentId)+1:end) = [];
  end
end

function [allBestScores, allBestWords, indices] = addSortScores(allBestScores, allBestWords, beamScores, beamSize, batchSize)
  % use previous beamScores, 1 * (beamSize*batchSize), update along the first dimentions
  allBestScores = bsxfun(@plus, allBestScores, beamScores);
  allBestScores = reshape(allBestScores, [beamSize*beamSize, batchSize]);
  allBestWords = reshape(allBestWords, [beamSize*beamSize, batchSize]);

  % for each sent, select the best beamSize candidates, out of beamSize*beamSize ones
  if beamSize > 1
    [allBestScores, indices] = sort(allBestScores, 'descend'); % (beamSize*beamSize) * batchSize
  else
    indices = ones(1, batchSize);
  end
end

%%
% return bestLogProbs, bestWords of sizes beamSize * curBatchSize
%%
function [bestLogProbs, bestWords, otherData] = nextBeamStep(models, lastDecStates, beamSize, ignoreSymbol, isChar, tgtUnk, unkDiscount, varargin)
  softmax_input = models{1}.W_soft*lastDecStates{1}{end}.softmax_h;
  if length(models)>1 % aggregate predictions from multiple models
    for ii=2:length(models)
      softmax_input = softmax_input + models{ii}.W_soft*lastDecStates{ii}{end}.softmax_h;
    end
    softmax_input = softmax_input./length(models);
  end
  [logProbs] = softmaxDecode(softmax_input);

  % sort
  if length(varargin)==1 % force decoder
    correctWords = reshape(varargin{1}, 1, []);
    bestWords = correctWords;
    bestLogProbs = logProbs(sub2ind(size(logProbs), correctWords, 1:size(logProbs, 2)));
    [~, otherData.maxWords] = max(logProbs, [], 1);
  else
    [sortedLogProbs, sortedWords] = sort(logProbs, 'descend');
    bestWords = sortedWords(1:beamSize, :);
    bestLogProbs = sortedLogProbs(1:beamSize, :);
    otherData = [];

    if ignoreSymbol>0
      flags = bestWords == ignoreSymbol;
      cols = find(any(flags)); % mark which column has the ignoreSymbol
      if ~isempty(cols)
        indices = find(flags);
        assert(length(indices) == length(cols)); % each col has a single ignoreSymbol
        bestWords(indices) = sortedWords(beamSize+1, cols);
        bestLogProbs(indices) = sortedLogProbs(beamSize+1, cols);
      end
    end

    if isChar==0 && unkDiscount>0
      flags = bestWords == tgtUnk;
      bestLogProbs(flags) = bestLogProbs(flags) - unkDiscount;
    end
  end
end

function [logProbs] = softmaxDecode(scores)
%% only compute logProbs
  mx = max(scores);
  scores = bsxfun(@minus, scores, mx); % subtract max elements
  logProbs = bsxfun(@minus, scores, log(sum(exp(scores))));
end

function [srcHidVecs] = duplicateSrcHidVecs(srcHidVecs, batchSize, beamSize)
  numElements = batchSize*beamSize;
  lstmSize = size(srcHidVecs, 1);
  numPositions = size(srcHidVecs, 3);
  srcHidVecs = permute(srcHidVecs, [1, 3, 2]); % lstmSize * numPositions * batchSize
  srcHidVecs = reshape(srcHidVecs, lstmSize*numPositions, batchSize);
  srcHidVecs = repmat(srcHidVecs, beamSize, 1);
  srcHidVecs = reshape(srcHidVecs, lstmSize, numPositions, numElements);
  srcHidVecs = permute(srcHidVecs, [1, 3, 2]); % lstmSize * numElements * numPositions
end

author: lmthang


text:
function [tt]=tt_qlaplace_dn(d)

% returns a rank-4,5...5 QTT decomposition of
% Delta_{1} \otimes \Id_{2} \ otimes \ldots \otimes \Id_{D} + \ldots +
%  + \Id_{1} \ otimes \ldots \otimes \Id_{D-1} \otimes Delta_{D},
% Delta_{k} being a discretization of Laplace operator on 2^d(k) points
% uniform grid,
% Dirichlet-Neumann boundary conditions being imposed
%
% D=size(d,2) must be >= 1
%
% September 3, 2010
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%
% Look for details in the Preprint No. 75, 2010 of
% Max-Planck Institute for Mathematics in the Sciences
% Vladimir A. Kazeev and Boris N. Khoromskij
% On explicit QTT representation of Laplace operator and its inverse
% http://www.mis.mpg.de/publications/preprints/2010/prepr2010-75.html

d=fliplr(d);
D=size(d,2);
tt=cell(sum(d),1);
I=eye(2);
J=zeros(2);
J(1,2)=1;
I2=zeros(2);
I2(2,2)=1;

if (D == 1)
	for key=1 : d
		if (key == 1)
			tt{key}=zeros(2,2,4);
			tt{key}(:,:,1)=2*I-J-J';
			tt{key}(:,:,2)=-J;
			tt{key}(:,:,3)=-J';
			tt{key}(:,:,4)=-I2;
		elseif (key == d)
			tt{key}=zeros(2,2,4);
			tt{key}(:,:,1)=I;
			tt{key}(:,:,2)=J';
			tt{key}(:,:,3)=J;
			tt{key}(:,:,4)=I2;
		else
			tt{key}=zeros(2,2,4,4);
			tt{key}(:,:,1,1)=I;
			tt{key}(:,:,2,2)=J;
			tt{key}(:,:,3,3)=J';
			tt{key}(:,:,4,4)=I2;
			tt{key}(:,:,2,1)=J';
			tt{key}(:,:,3,1)=J;
		end
	end
else

	key=0;
	for k=1 : D
		for kappa=1 : d(k)
			key=key+1;
			if (kappa == 1)
				if (k == 1)
					tt{key}=zeros(2,2,5);
					tt{key}(:,:,1)=2*I-J-J';
					tt{key}(:,:,2)=-J;
					tt{key}(:,:,3)=-J';
					tt{key}(:,:,4)=-I2;
					tt{key}(:,:,5)=I;
				elseif (k == D)
					tt{key}=zeros(2,2,2,4);
					tt{key}(:,:,1,1)=2*I-J-J';
					tt{key}(:,:,1,2)=-J;
					tt{key}(:,:,1,3)=-J';
					tt{key}(:,:,1,4)=-I2;
					tt{key}(:,:,2,1)=I;
				else
					tt{key}=zeros(2,2,2,5);
					tt{key}(:,:,1,1)=2*I-J-J';
					tt{key}(:,:,1,2)=-J;
					tt{key}(:,:,1,3)=-J';
					tt{key}(:,:,1,4)=-I2;
					tt{key}(:,:,1,5)=I;
					tt{key}(:,:,2,1)=I;
				end
			elseif (kappa == d(k))
				if (k == D)
					tt{key}=zeros(2,2,4);
					tt{key}(:,:,1)=I;
					tt{key}(:,:,2)=J';
					tt{key}(:,:,3)=J;
					tt{key}(:,:,4)=I2;
				else
					tt{key}=zeros(2,2,5,2);
					tt{key}(:,:,5,1)=I;
					tt{key}(:,:,1,2)=I;
					tt{key}(:,:,2,2)=J';
					tt{key}(:,:,3,2)=J;
					tt{key}(:,:,4,2)=I2;
				end
			else
				if (k == D)
					tt{key}=zeros(2,2,4,4);
					tt{key}(:,:,1,1)=I;
					tt{key}(:,:,2,2)=J;
					tt{key}(:,:,3,3)=J';
					tt{key}(:,:,4,4)=I2;
					tt{key}(:,:,2,1)=J';
					tt{key}(:,:,3,1)=J;
				else
					tt{key}=zeros(2,2,5,5);
					tt{key}(:,:,1,1)=I;
					tt{key}(:,:,2,2)=J;
					tt{key}(:,:,3,3)=J';
					tt{key}(:,:,4,4)=I2;
					tt{key}(:,:,2,1)=J';
					tt{key}(:,:,3,1)=J;
					tt{key}(:,:,5,5)=I;
				end
			end
		end
	end
end
tt=tt_matrix(tt); % @Bydlocode
return
end

author: oseledets

repo: TT-Toolbox


text:
function C = cell_add(A,B,a,b)

% cell_add - add two cell arrays
%
%   C = cell_add(A,B,a,b);
%
%   C{i} = a*A{i} + b*B{i};
%
%   Copyright (c) 2008 Gabriel Peyre

if nargin<3
    a = 1;
end
if nargin<4
    b = 1;
end

if iscell(A)
    if length(A)~=length(B)
        error('A and B must be of the same size');
    end
    for i=1:length(A)
        C{i} = cell_add(A{i},B{i},a,b);
    end
    return;
end

C = a*A + b*B;

author: gpeyre

repo: matlab-toolboxes