text:
%   CLASS File_Name_Processor
% =========================================================================
%
% DESCRIPTION
%   Class to process file names with keywords
%
% EXAMPLE
%   fnp = File_Name_Processor();
%   fnp.dateKeyRep('${WWWW}/igs${WWWWD}_${6H}.clk_30s', GPS_Time(datenum('2017/02/05 23:12:11')))
%
% FOR A LIST OF CONSTANTs and METHODS use doc FTP_Server
%
% REQUIRES:
%   goGPS settings;
%
% COMMENTS
%   Nothing to report

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti, Giulio Tagliaferro
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

classdef File_Name_Processor < handle

    properties (Constant)
        GPS_WEEK = '${WWWW}';
        GPS_WD = '${WWWWD}';
        GPS_DOW = '${D}';
        GPS_3H = '${3H}';
        GPS_6H = '${6H}';
        GPS_HH = '${HH}';
        GPS_YY = '${YY}';
        GPS_YYDOY = '${YYDOY}';
        GPS_YYYY = '${YYYY}';
        GPS_DOY = '${DOY}';
        GPS_FIRST_DOY_TRIMESTER = '${DOYT}';
        GPS_SESSION = '${S}';
        GPS_MM = '${MM}';
        GPS_DD = '${DD}';
        GPS_1D = '${1D}';
        GPS_QQ = '${QQ}';
        GPS_5M = '${5M}';
        VMF_RES = '${VMFR}';
        VMF_SOURCE = '${VMFS}';
    end

    properties (SetAccess = private, GetAccess = public)
        log = Logger.getInstance();
    end

    properties (SetAccess = protected, GetAccess = protected)
        addr;        % ip address (or alias)
        port = 21;   % connection port
        path = '/';  % path to files
    end

    methods
        function this = File_Name_Processor()
            % Constructor
        end

        function file_name_out = dateKeyRep(this, file_name, date, session,vmf_res,vmf_source)
            % substitute time placeholder with the proper format
            % SYNTAX: file_name = this.dateKeyRep(file_name, date, session)
            narginchk(3,6)
            if (nargin < 4)
                session = '0';
            end
            if (nargin < 5)
                vmf_res = '';
            end
            if (nargin < 6)
                vmf_source = '';
            end
            file_name_out = file_name;
            if any(file_name_out == '$') && ~isempty(regexp(file_name_out, '\$\{Y|(DOY)', 'once'))
                [year, doy] = date.getDOY();
                file_name_out = strrep(file_name_out, this.GPS_DOY, sprintf('%03d', doy));
                file_name_out = strrep(file_name_out, this.GPS_FIRST_DOY_TRIMESTER, sprintf('%03d', GPS_Time.doy2fdoytr(doy,year)));
                file_name_out = strrep(file_name_out, this.GPS_YY, sprintf('%02d', mod(year,100)));
                file_name_out = strrep(file_name_out, this.GPS_YYYY, sprintf('%04d', year));
                file_name_out = strrep(file_name_out, this.GPS_YYDOY, sprintf('%02d%03d', mod(year,100), doy));
            end
            if any(file_name_out == '$') && ~isempty(regexp(file_name_out, '\$\{(WWWW)|(WD)|(DOW)|(3H)|(6H)|(HH)|(QQ)|(5M)\}', 'once'))
                [gps_week, gps_sow, gps_dow] = date.getGpsWeek();
                file_name_out = strrep(file_name_out, this.GPS_WEEK, sprintf('%04d', gps_week(1)));
                if any(file_name_out == '$')
                    file_name_out = strrep(file_name_out, this.GPS_DOW, sprintf('%01d', gps_dow(1)));
                end
                if any(file_name_out == '$')
                    file_name_out = strrep(file_name_out, this.GPS_WD, sprintf('%04d%01d', gps_week(1), gps_dow(1)));
                    file_name_out = strrep(file_name_out, this.GPS_3H, sprintf('%02d', fix((gps_sow(1) - double(gps_dow(1)) * 86400)/(3*3600))*3));
                    file_name_out = strrep(file_name_out, this.GPS_6H, sprintf('%02d', fix((gps_sow(1) - double(gps_dow(1)) * 86400)/(6*3600))*6));
                    file_name_out = strrep(file_name_out, this.GPS_HH, sprintf('%02d', fix((gps_sow(1) - double(gps_dow(1)) * 86400)/(3600))));
                    file_name_out = strrep(file_name_out, this.GPS_QQ, sprintf('%02d', mod(15 * fix((gps_sow(1) - double(gps_dow(1)) * 86400)/(900)), 60)));
                    file_name_out = strrep(file_name_out, this.GPS_5M, sprintf('%02d', mod(5 * fix((gps_sow(1) - double(gps_dow(1)) * 86400)/(300)), 60)));
                end
            end
            if any(file_name_out == '$') && ~isempty(regexp(file_name_out, '\$\{(MM)|DD)|(1D)\}', 'once'))
                [date] = datevec(date.getMatlabTime());
                file_name_out = strrep(file_name_out, this.GPS_MM, sprintf('%02d', date(2)));
                file_name_out = strrep(file_name_out, this.GPS_DD, sprintf('%02d', date(3)));
                file_name_out = strrep(file_name_out, this.GPS_1D, sprintf('%d', date(3)));
            end
            if any(file_name_out == '$')
                file_name_out = strrep(file_name_out, this.GPS_SESSION, sprintf('%01d', session));
            end
            if any(file_name_out == '$')
                file_name_out = strrep(file_name_out, this.VMF_RES, char(vmf_res));
            end
            if any(file_name_out == '$')
                file_name_out = strrep(file_name_out, this.VMF_SOURCE, char(vmf_source));
            end
        end

        function step_sec = getStepSec(this, file_name)
            % get the shorter time keyword present in file_name [seconds]
            % SYNTEX: step_sec = this.getStepSec(file_name);
            % Check for GPS time placeholders
            step_sec = 0;
            if ~isempty(strfind(file_name, this.GPS_5M))
                step_sec = 300;
            elseif ~isempty(strfind(file_name, this.GPS_QQ))
                step_sec = 900;
            elseif ~isempty(strfind(file_name, this.GPS_HH))
                step_sec = 3600;
            elseif ~isempty(strfind(file_name, this.GPS_3H))
                step_sec = 3 * 3600;
            elseif ~isempty(strfind(file_name, this.GPS_6H))
                step_sec = 6 * 3600;
            elseif ~isempty((strfind(file_name, this.GPS_DOW))) || ~isempty(strfind(file_name, this.GPS_DD)) || ~isempty(strfind(file_name, this.GPS_WD)) || ~isempty(strfind(file_name, this.GPS_DOY)) || ~isempty(strfind(file_name, this.GPS_1D)) || ~isempty(strfind(file_name, this.GPS_YYDOY))
                step_sec = 24 * 3600;
            elseif ~isempty(strfind(file_name, this.GPS_WEEK))
                step_sec = 24 * 3600 * 7;
            elseif ~isempty(strfind(file_name, this.GPS_FIRST_DOY_TRIMESTER))
                step_sec = 24 * 3600 * 91;
            elseif ~isempty(strfind(file_name, this.GPS_YYYY)) || ~isempty(strfind(file_name, this.GPS_YY))
                step_sec = 24 * 3600 * 365;
            end
        end

        function [file_name_lst, date_list] = dateKeyRepBatch(this, file_name, date_start, date_stop, session_list, session_start, session_stop, vmf_res, vmf_source)
            % substitute time placeholder with the proper format
            % SYNTAX: file_name = this.dateKeyRepBatch(file_name, date_start, date_stop)
            % NOTE: I consider only two possible file formats:
            %         - dependend on UTC time (year, doy)
            %         - dependent on GPS time (week, day of week, h of the day)

            if nargin < 5
                session_list = '0';
                session_start = '0';
                session_stop = '0';
            end
            if nargin < 8
                vmf_res = '';
                vmf_source = '';
            end
            session = ~isempty(strfind(file_name, this.GPS_SESSION));
            sss_start = strfind(session_list, session_start);
            sss_stop = strfind(session_list, session_stop);

            % Check for GPS time placeholders
            step_sec = this.getStepSec(file_name);

            if (step_sec > 0)
                file_name_lst = {};
                date_list = date_start.getCopy();
                if not(abs(step_sec - (24 * 3600 * 91)) < 1e-3)
                    date0 = GPS_Time((floor(((date_start.getNominalTime(0.5).getMatlabTime() - GPS_Time.GPS_ZERO) * 86400) / step_sec) * step_sec) / 86400 + GPS_Time.GPS_ZERO);
                else
                    date0 = date_start.getCopy();
                end
                date0 = date0.getNominalTime(60);
                date1 = date_stop.getNominalTime(1);

                date_list.toUnixTime(); % keep an higher precision

                i = 1;
                % Find all the file in the interval of dates
                while (date0.getMatlabTime() <= date1.getMatlabTime())
                    if session
                        % run over session
                        for s = sss_start : length(session_list)
                            file_name_lst{i} = this.keyRep(file_name, this.GPS_SESSION, session_list(s)); %#ok<AGROW>
                            file_name_lst{i} = this.dateKeyRep(file_name_lst{i}, date0, '0',vmf_res,vmf_source); %#ok<AGROW>
                            i = i + 1;
                        end
                        sss_start = 1;
                    else
                        file_name_lst{i} = this.dateKeyRep(file_name, date0, '0', vmf_res,vmf_source); %#ok<AGROW>
                        i = i + 1;
                    end
                    date0.addIntSeconds(step_sec);
                    date0 = date0.getNominalTime(60);
                    if (date0.getMatlabTime() <= date1.getMatlabTime())
                        date_list.append(date0);
                    end
                end
            else
                date_list = GPS_Time();
                i = 1;
                % run over session
                if session
                    for s = sss_start : length(session_list)
                        file_name_lst{i} = this.keyRep(file_name, this.GPS_SESSION, session_list(s)); %#ok<AGROW>
                        i = i + 1;
                    end
                else
                    if iscell(file_name)
                        file_name_lst = file_name;
                    else
                        file_name_lst = {file_name};
                    end
                end
            end

            if session
                sss_clean = length(file_name_lst) - (length(session_list)-sss_stop) + 1 : length(file_name_lst);
                file_name_lst(sss_clean) = [];
            end
            if iscell(file_name_lst)
                file_name_lst = file_name_lst';
            end
        end
    end

    methods (Static)
        function dir_path_out = getInnerDir(dir_path_list)
            % Get recursive dir list
            %
            % SYNTAX
            %   dir_path_out = getInnerDir(dir_path_list)

            dir_path_out = cell(0);
            % make base_dir a cell
            if ~iscell(dir_path_list)
                dir_path_list = {dir_path_list};
            end

            % loop on base_dirs
            for d = 1 : numel(dir_path_list)
                dir_path_out = [dir_path_out dir_path_list(d)];
                dir_path = dir_path_list{d};
                fprintf(' Scanning %s for subdirs\n', dir_path);
                % Search fo subdirs
                tmp = dir(dir_path);
                tmp = tmp([tmp.isdir]);
                dir_list = {};
                for i = numel(tmp) : -1 : 1
                    % remove hidden or special dirs (log)
                    if not(tmp(i).name(1) == '.') && not(strcmp(tmp(i).name, 'log'))
                        dir_list = [dir_list {fullfile(dir_path, tmp(i).name)}];
                    end
                end

                if ~isempty(dir_list)
                    for c = 1 : numel(dir_list)
                        cur_dir = dir_list(c);
                        dir_path_out = [dir_path_out File_Name_Processor.getInnerDir(cur_dir)];
                    end
                end
            end
        end

        function file_list_out = findFiles(dir_path_list, prefix)
            % Get recursive file list
            %
            % SYNTAX
            %   file_list_out = findFiles(dir_path_list, prefix)
            file_list_out = cell(0);
            % make base_dir a cell
            if ~iscell(dir_path_list)
                dir_path_list = {dir_path_list};
            end

            % loop on base_dirs
            for d = 1 : numel(dir_path_list)
                dir_path = dir_path_list{d};
                fprintf(' Scanning %s for subdirs\n', dir_path);
                % Search fo subdirs
                tmp = dir(dir_path);
                tmp_dir_list = tmp([tmp.isdir]);
                tmp_file_list = tmp(~[tmp.isdir]);
                dir_list = {};
                for i = 1 : numel(tmp_dir_list)
                    % remove hidden or special dirs (log)
                    if not(tmp_dir_list(i).name(1) == '.') && not(strcmp(tmp_dir_list(i).name, 'log'))
                        dir_list = [dir_list {fullfile(dir_path, tmp_dir_list(i).name)}];
                    end
                end

                if ~isempty(dir_list)
                    for c = 1 : numel(dir_list)
                        cur_dir = dir_list(c);
                        file_list_out = [file_list_out File_Name_Processor.findFiles(cur_dir, prefix)];
                    end
                end

                for i = 1 : numel(tmp_file_list)
                    if length(tmp_file_list(i).name) > length(prefix) &&  strcmp(tmp_file_list(i).name(1:length(prefix)), prefix)
                        file_list_out = [file_list_out {fullfile(dir_path, tmp_file_list(i).name)}];
                    end
                end
            end
        end

        function dir_path = getFullDirPath(dir_path, dir_base, dir_fallback, empty_fallback)
            % Get the full path given the relative one and the relative dir_base
            % It changes the folder with dir_fallback only if the changed folder exist
            %
            % SYNTAX:
            %   dir_path = getFullDirPath(dir_path, <dir_base default = pwd>, fallback_dir_base);

            if nargin == 1
                dir_base = pwd;
            end
            if (nargin < 4)
                empty_fallback = dir_base;
            end
            if isempty(dir_base)
                dir_base = pwd;
            end

            dir_path_bk = dir_path;

            fnp = File_Name_Processor;
            dir_path = fnp.checkPath(dir_path);
            if isunix
                if ~isempty(dir_path)
                    if (dir_path(1) == '/')
                        dir_path = fnp.checkPath(dir_path);
                    else
                        if ~isempty(dir_base)
                            if (dir_base(1) ~= '/')
                                dir_base = fnp.getFullDirPath(dir_base, pwd);
                            end
                        end
                        dir_path = fnp.checkPath([dir_base filesep dir_path]);
                    end
                else
                    dir_path = fnp.checkPath(empty_fallback);
                end
            else
                if length(dir_path) > 1
                    if (dir_path(2) == ':')
                        dir_path = fnp.checkPath(dir_path);
                    else
                        if length(dir_base) > 1
                            if (dir_base(2) == ':')
                                dir_base = fnp.getFullDirPath(dir_base, pwd);
                            end
                        end
                        if ~(numel(dir_path) >= 2 && strcmp(dir_path(1:2), '\\'))
                            dir_path = fnp.checkPath([dir_base filesep dir_path]);
                        end
                    end
                else
                    dir_path = fnp.checkPath(empty_fallback);
                end
            end

            dir_path = File_Name_Processor.remDirModifiers(dir_path);

            % Fallback if not exist
            if (nargin >= 3) && ~isempty(dir_fallback)
                if ~exist(dir_path, 'file')
                    prefix = '';
                    if numel(dir_fallback) > 2
                        prefix = iif(strcmp(dir_fallback(1:2),'\\'), '\\', '');
                    end
                    tmp = fnp.getFullDirPath(dir_path_bk, dir_fallback);
                    % change dir_path only if the new path exist
                    if exist(tmp, 'file')
                        dir_path = tmp;
                    end
                    dir_path = [prefix dir_path];
                end
            end
        end

        function path = remDirModifiers(path)
            % Utility to remove from a path dir changes with "." and ".."
            %
            % SYNTAX
            %   path = File_Name_Processor.remDirModifiers(path)

            if ~isempty(path)
                % remove './'
                path = strrep(path, [filesep '.' filesep], filesep);

                prefix = '';
                if numel(path) > 2
                    prefix = iif(strcmp(path(1:2),'\\'), '\\', '');
                end

                % extract sub folder names
                list = regexp(path,['[^' iif(filesep == '\', '\\', filesep) ']*'],'match');

                % search for "../"
                dir_up = find(strcmp(list,'..'));
                offset = 0;
                for i = dir_up
                    list(max(1, (i - 1 : i) - offset)) = [];
                    offset = offset + 2;
                end

                % restore full path start
                if isunix
                    path = [prefix filesep strCell2Str(list, filesep)];
                else
                    path = [prefix strCell2Str(list, filesep)];
                end
            end
        end

        function dir_path = getRelDirPath(dir_path, dir_base)
            % Get the full path given the relative one and the relative dir_base
            % SYNTAX: dir_path = getRelativeDirPath(dir_path, <dir_base default = pwd>);

            if nargin == 1
                dir_base = pwd;
            end

            fnp = File_Name_Processor;
            dir_base = fnp.getFullDirPath(dir_base, pwd);
            list_base = regexp(dir_base, ['[^' iif(filesep == '\', '\\', filesep) ']*'], 'match');
            n_dir_base = numel(list_base);

            if iscell(dir_path)
                for j = 1 : numel(dir_path)
                    dir_path{j} = fnp.getFullDirPath(dir_path{j}, dir_base);
                    drive_path = regexp(dir_path{j},'.(?=\:)', 'match', 'once');
                    drive_base = regexp(dir_base,'.(?=\:)', 'match', 'once');
                    if (isempty(drive_path) || drive_path == drive_base)
                        list_path = regexp(dir_path{j}, ['[^' iif(filesep == '\', '\\', filesep) ']*'], 'match');
                        n_dir_path = numel(list_path);
                        i = 0;
                        while (i < n_dir_base) && (i < n_dir_path) && strcmp(list_base{i+1}, list_path{i+1})
                            i = i + 1;
                        end
                        if (n_dir_base -i) > 0
                            list_path = [{repmat(['..' filesep],1, n_dir_base - i)} list_path(i+1:end)];
                        else
                            list_path = list_path(i+1:end);
                        end

                        dir_path{j} = strrep(strCell2Str(list_path, filesep), [filesep filesep], filesep);
                    end
                end
            else
                % if dir_path is a network address && dir_base is not
                prefix = '';
                if (numel(dir_path) > 2 && strcmp(dir_path(1:2),'\\')) && ~((numel(dir_base) > 2 && strcmp(dir_base(1:2),'\\')))
                    n_dir_base = 0;
                    list_base = {};
                    prefix = '\\';
                end
                dir_path = fnp.getFullDirPath(dir_path, dir_base);
                drive_path = regexp(dir_path,'.(?=\:\\)', 'match', 'once');
                drive_base = regexp(dir_base,'.(?=\:\\)', 'match', 'once');
                if (isempty(drive_path) || drive_path == drive_base)
                    list_path = regexp(dir_path, ['[^' iif(filesep == '\', '\\', filesep) ']*'], 'match');
                    n_dir_path = numel(list_path);
                    i = 0;
                    while (i < n_dir_base) && (i < n_dir_path) && strcmp(list_base{i+1}, list_path{i+1})
                        i = i + 1;
                    end
                    if (n_dir_base -i) > 0
                        if i > 0
                            list_path = [{repmat(['..' filesep],1, n_dir_base - i)} list_path(i+1:end)];
                        elseif isempty(prefix)
                            prefix = filesep;
                        end

                    else
                        list_path = list_path(i+1:end);
                    end

                    dir_path = [prefix strrep(strCell2Str(list_path, filesep), [filesep filesep], filesep)];
                end
            end
        end

        function file_name = getFileName(file_name)
            % Get only the file name of a full path
            [~, name, extension] = fileparts(file_name);
            file_name = [name extension];
        end

        function path = getPath(file_name)
            % Get only the path up to the folder name from a full path
            [path, ~, ~] = fileparts(file_name);
        end

        function file_name = keyRep(file_name, key, substitution)
            % Substitute a key in the file_name with another value
            file_name = strrep(file_name, key, substitution);
        end

        function [str_cell] = toIniStringComment(str_cell)
            % Get the list of accepted substitution
            if (nargin == 0)
                str_cell = {};
            end
            str_cell = Ini_Manager.toIniStringComment('Special Keywords that can be used in file names:', str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s    4 char GPS week', File_Name_Processor.GPS_WEEK), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s 4+1 char GPS week + day of the week', File_Name_Processor.GPS_WD), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s       1 char day of the week', File_Name_Processor.GPS_DOW), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s      2 char GPS hour (00, 03, 06, 09, 12, 15, 18, 21)', File_Name_Processor.GPS_3H), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s      2 char GPS hour (00, 06, 12, 18)', File_Name_Processor.GPS_6H), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s      2 char GPS hour', File_Name_Processor.GPS_HH), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s      2 char GPS quarter of hour (00, 15, 30, 45)', File_Name_Processor.GPS_QQ), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s      2 char GPS five minutes (05, 10, ... , 55)', File_Name_Processor.GPS_5M), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s 2+3 char GPS year + day of year', File_Name_Processor.GPS_YYDOY), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s    4 char GPS year', File_Name_Processor.GPS_YYYY), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s      2 char GPS year', File_Name_Processor.GPS_YY), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s      2 char GPS month', File_Name_Processor.GPS_MM), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s      2 char GPS day', File_Name_Processor.GPS_DD), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s     3 char GPS day of the year', File_Name_Processor.GPS_DOY), str_cell);
            str_cell = Ini_Manager.toIniStringComment(sprintf(' - %s       1 char session', File_Name_Processor.GPS_SESSION), str_cell);
        end

        function full_path = getFullPath(dir_path, file_path)
            % Get the list of files combining dir_path + file_path
            % It works when file_path is a cell or a simple string
            % SYNTAX: full_path = IO_Settings.getFullPath(dir_path, file_path);

            fnp = File_Name_Processor();

            if iscell(file_path)
                full_path = cell(numel(file_path),1);
                for f = 1 : numel(file_path)
                    full_path{f} = fnp.checkPath(fullfile(dir_path, file_path{f}));
                end
            else
                full_path = fnp.checkPath(fullfile(dir_path, file_path));
            end
        end

        function [universal_path, is_valid] = checkPath(path, prj_home)
            % Conversion of path OS specific to a universal one: "/" or "\" are converted to filesep
            %
            % SYNTAX:
            %   universal_path = checkPath(path)
            %
            % INPUT:
            %   path
            %
            % OUTPUT:
            %   universal_path
            %   < is_valid >            optional, contains the status of existence
            %
            % DESCRIPTION:
            %   Conversion of path OS specific to a universal one: "/" or "\" are converted to filesep
            %   if the second parameter is present is_valid contains the status of existence
            %       2 => is a file
            %       7 => is a folder

            if instr(path, '${PRJ_HOME}')
                if nargin < 2
                %    state = Core.getState();
                %    prj_home = state.getHomeDir;
                    prj_home = '';
                end
                path = strrep(path, '${PRJ_HOME}', prj_home);
            end

            prefix = '';
            if numel(path) > 2
                prefix = iif(strcmp(path(1:2),'\\'), '\', '');
            end

            if not(isempty(path))
                if (iscell(path))
                    % for each line of the cell
                    universal_path = cell(size(path));
                    for c = 1 : length(path)
                        check_drive_win = strfind(path{c},':');
                        if (~isempty(check_drive_win))
                            path{c} = path{c}(max(1, check_drive_win - 1) : end);
                        end
                        universal_path{c} = regexprep(path{c}, '(\\(?![ ]))|(\/)', filesep);
                        universal_path{c} = regexprep(universal_path{c}, ['\' filesep '\' filesep], filesep);
                    end
                    if (nargout == 2)
                        is_valid = zeros(size(path));
                        for c = 1 : length(path)
                            is_valid(c) = exist(universal_path{c}, 'file'); % if it is a file is_valid contains 2, if it is a dir it contains 7
                        end
                    end
                else
                    if ischar(path)
                        check_drive_win = strfind(path,':');
                        if (~isempty(check_drive_win))
                            path = path(max(1, check_drive_win - 1) : end);
                        end
                        universal_path = regexprep(path, '(\\(?![ ]))|(\/)', filesep);
                        universal_path = regexprep(universal_path, ['\' filesep '\' filesep], filesep);
                        if (nargout == 2)
                            is_valid = exist(path, 'file'); % if it is a file is_valid contains 2, if it is a dir it contains 7
                        end
                    else
                        universal_path = '';
                        is_valid = 0;
                    end
                end
            else
                universal_path = '';
                is_valid = 0;
            end
            universal_path = [prefix universal_path];
        end

    end
end

author: goGPS-Project


text:
classdef snr
    %SNR helper functions for calculating SNR

    properties
    end

    methods(Static)
        function result = SNR(signal, noise)
            sigEnergy = sum(signal .* conj(signal));
            noiseEnergy = sum(noise .* conj(noise));
            result=10*log10(sigEnergy./noiseEnergy);
        end
        function result = recSNRdB(signal, reconstruction)
            result=20*log10(norm(signal)/norm(signal-reconstruction));
        end

        function result = energyDB(signal)
            result = 20*log10(norm(signal));
        end
        function result = recSNR(signal, reconstruction)
            noise = (signal - reconstruction);
            result = (signal' * signal) / (noise'*noise);
        end
        function result = normalizedErrorEnergy(signal, reconstruction)
            noise = (signal - reconstruction);
            result = (noise' * noise) / (signal'*signal);
        end

        function result = recSNRsdB(signals, reconstructions)
            ratios = spx.snr.recSNRs(signals, reconstructions);
            result= 10 * log10(ratios);
        end

        function result = recSNRs(signals, reconstructions)
            sigEngergies =  spx.snr.energies(signals);
            noiseEnergies =  spx.snr.energies(signals - reconstructions);
            result = sigEngergies ./ noiseEnergies;
        end

        function result = energies(signals)
            signals = abs(signals);
            signals = signals .^2;
            result = sum(signals, 1);
        end

        function result = energiesDB(signals)
            energies = spx.snr.energies(signals);
            result = 10* log10(energies);
        end
    end

end

author: indigits

repo: sparse-plex


text:
function addToPath(arg)

% Adds the MFE Toolbox to the current path and optionally saves the path
silent = false;
if nargin>0 && strcmpi(arg,'-silent')
    silent = true;
end

dirList = {'bootstrap',...
    'crosssection',...
    'distributions',...
    'GUI',...
    'multivariate',...
    'tests',...
    'timeseries',...
    'univariate',...
    'utility',...
    'realized',...
    };
optionalDirList = {'duplication'};
mexDir = 'dlls';

for i=1:length(dirList)
    addpath(fullfile(pwd,dirList{i}));
end

user_entry = '';
while ~ismember({user_entry},{'Y','N'}) && ~silent
    user_entry = input('Add MATLAB work-a-like functions? [Y/N(default)] ','s');
    user_entry = upper(user_entry);
    if isempty(user_entry)
        user_entry = 'N';
    end
end

if strcmpi(user_entry,'Y')
    for i=1:length(optionalDirList)
        addpath(fullfile(pwd,optionalDirList{i}));
    end
end

user_entry = '';
while ~ismember({user_entry},{'Y','N'}) && ~silent
    user_entry = input('Make path permanent? [Y/N (default)] ','s');
    user_entry = upper(user_entry);
    if isempty(user_entry)
        user_entry = 'N';
    end
end

if strcmpi(user_entry,'Y')
    try
        savepath
    catch ME
        warning('MFEToolbox:Path',['There was a problem saving the path.  The error was: \n' ME.message])
    end
end

if strcmpi(computer,'PCWIN64')
    user_entry = '';
    while ~ismember({user_entry},{'Y','N'}) && ~silent
        user_entry = input('Add MEX files to path? [Y(default)/N] ','s');
        user_entry = upper(user_entry);
        if isempty(user_entry)
            user_entry = 'Y';
        end
    end
    if strcmp(user_entry,'Y') || silent
        addpath(fullfile(pwd,mexDir));
    end
end

author: bashtage

repo: mfe-toolbox


text:
function [results, success, raw] = dcopf_solver(om, mpopt)
%DCOPF_SOLVER  Solves a DC optimal power flow.
%
%   [RESULTS, SUCCESS, RAW] = DCOPF_SOLVER(OM, MPOPT)
%
%   Inputs are an OPF model object and a MATPOWER options struct.
%
%   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.
%
%   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
%   branch, gen, gencost fields, along with the following additional
%   fields:
%       .order      see 'help ext2int' for details of this field
%       .x          final value of optimization variables (internal order)
%       .f          final objective function value
%       .mu         shadow prices on ...
%           .var
%               .l  lower bounds on variables
%               .u  upper bounds on variables
%           .lin
%               .l  lower bounds on linear constraints
%               .u  upper bounds on linear constraints
%
%   SUCCESS     1 if solver converged successfully, 0 otherwise
%
%   RAW         raw output in form returned by MINOS
%       .xr     final value of optimization variables
%       .pimul  constraint multipliers
%       .info   solver specific termination code
%       .output solver specific output information
%
%   See also OPF, OPT_MODEL/SOLVE.

%   MATPOWER
%   Copyright (c) 2000-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
%
%   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.

%%----- initialization -----
%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

%% unpack data
mpc = om.get_mpc();
[baseMVA, bus, gen, branch, gencost] = ...
    deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
cp = om.params_legacy_cost();
Bf = om.get_userdata('Bf');
Pfinj = om.get_userdata('Pfinj');
[vv, ll] = om.get_idx();

%% problem dimensions
nb = size(bus, 1);          %% number of buses
nl = size(branch, 1);       %% number of branches
ny = om.getN('var', 'y');   %% number of piece-wise linear costs

%% options
model = om.problem_type();
opt = mpopt2qpopt(mpopt, model);
if strcmp(opt.alg, 'OSQP')
    opt.x0 = [];    %% disable provided starting point for OSQP
end

%% try to select an interior initial point, unless requested not to
if mpopt.opf.start < 2 && ...
        (strcmp(opt.alg, 'MIPS') || strcmp(opt.alg, 'IPOPT'))
    [x0, xmin, xmax] = om.params_var();     %% init var & bounds
    s = 1;                      %% set init point inside bounds by s
    lb = xmin; ub = xmax;
    lb(xmin == -Inf) = -1e10;   %% replace Inf with numerical proxies
    ub(xmax ==  Inf) =  1e10;
    x0 = (lb + ub) / 2;         %% set x0 mid-way between bounds
    k = find(xmin == -Inf & xmax < Inf);    %% if only bounded above
    x0(k) = xmax(k) - s;                    %% set just below upper bound
    k = find(xmin > -Inf & xmax == Inf);    %% if only bounded below
    x0(k) = xmin(k) + s;                    %% set just above lower bound
    Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
    x0(vv.i1.Va:vv.iN.Va) = Varefs(1);  %% angles set to first reference angle
    if ny > 0
        ipwl = find(gencost(:, MODEL) == PW_LINEAR);
        c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));    %% largest y-value in CCV data
        x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
    end
    opt.x0 = x0;
end

%%-----  run opf  -----
[x, f, eflag, output, lambda] = om.solve(opt);
success = (eflag == 1);

%%-----  calculate return values  -----
if ~any(isnan(x))
    %% update solution data
    Va = x(vv.i1.Va:vv.iN.Va);
    Pg = x(vv.i1.Pg:vv.iN.Pg);

    %% update voltages & generator outputs
    bus(:, VM) = ones(nb, 1);
    bus(:, VA) = Va * 180/pi;
    gen(:, PG) = Pg * baseMVA;

    %% compute branch flows
    branch(:, [QF, QT]) = zeros(nl, 2);
    branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
    branch(:, PT) = -branch(:, PF);
end

%% package up results
mu_l = lambda.mu_l;
mu_u = lambda.mu_u;
muLB = lambda.lower;
muUB = lambda.upper;

%% update Lagrange multipliers
il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
bus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);
gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);
branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);
bus(:, LAM_P)       = (mu_u(ll.i1.Pmis:ll.iN.Pmis) - mu_l(ll.i1.Pmis:ll.iN.Pmis)) / baseMVA;
if ~isempty(il)
    branch(il, MU_SF)   = mu_u(ll.i1.Pf:ll.iN.Pf) / baseMVA;
    branch(il, MU_ST)   = mu_l(ll.i1.Pf:ll.iN.Pf) / baseMVA;
end
gen(:, MU_PMIN)     = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
gen(:, MU_PMAX)     = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
pimul = [
  mu_l - mu_u;
 -ones(ny>0, 1);    %% dummy entry corresponding to linear cost row in A (in MINOS)
  muLB - muUB
];

mu = struct( ...
  'var', struct('l', muLB, 'u', muUB), ...
  'lin', struct('l', mu_l, 'u', mu_u) );

results = mpc;
[results.bus, results.branch, results.gen, ...
    results.om, results.x, results.mu, results.f] = ...
        deal(bus, branch, gen, om, x, mu, f);

raw = struct('xr', x, 'pimul', pimul, 'info', eflag, 'output', output);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [varargout] = binotest_dependent(X, Po)
% This function runs several different types of tests on dependent binomial data.
%
% :Usage:
% ::
%
%     [varargout] = binotest_dependent(X, Po)
%
% Overall, it tests the number of "hits" for each subject (row in X) against a null-hypothesis
% proportion p, across all subjects using a Z-test (two-tailed).
% The second level null hypothesis should be approximated by a normal
% distribution with a mean of p.  This approach assumes that each subject
% has an equal number of independent Bernoulli trials (columns in X) and
% that the number of subjects exceeds n=20 the test will be more accurate as n -> infinity.
% Also, calculates tests for each separate trial (e.g., subject columns),
% and the difference between proportions (two proportion z-test).
%
% :Inputs:
%
%   **X:**
%        X is a matrix of "hits" and "misses", coded as 1s and 0s.
%        where rows = subjects and columns = observations within
%        subject
%
%   **Po:**
%        Po is the null hypothesis proportion of "hits", e.g., often p = 0.5
%
% :Outputs:
%
%   **RES [1:5]:**
%        a structure containing the output of the stats for the
%        z-test, includes the number of subject (N), number of overall
%        hits (hits), the overall proportion of hits (prop), the standard
%        deviation (SE), z-statistic (Z), and the two tailed p-value (pval)
%        trial across all subjects.  Assumes independence
%
%   **RES1:**
%        Independent Single Interval Test for Column 1 (Column 1 only against Po)
%
%   **RES2:**
%        Independent Single Interval Test for Column 2 (Column 2 only against Po)
%
%   **RES3:**
%        Two proportion dependent difference z-test (Column 1 minus Column 2 against 0)
%
%   **RES4:**
%        Dependent single-interval test (Mean of Column 1 and Column 2 against Po)
%
%   **RES5:**
%        Two proportion dependent addition z-test (Column 1 plus Column 2 against 2 * Po)
%              (Similar to mean, not sure what this will be used for)
%
% :Examples:
% ::
%
%    [RES1, RES2, RES3, RES4, RES5] = binotest_dependent([1,1,1,1,0; 1,0,1,0,1]',.5)
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2014  Luke Chang & Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..

% ..
%    Process Data
% ..

X = double(X); % just in case
[N, k] = size(X);  % number of subjects x dependent obs per person
A = X(:, 1);
B = X(:, 2);

% -------------------------------------------------------------------------
% Test if Data is in Correct Format
% -------------------------------------------------------------------------
if k > 2, error('NOT IMPLEMENTED FOR MORE THAN 2 DEP OBSERVATIONS PER UNIT'); end

u = unique(X(:));
if any(u ~= 1 & u ~= 0), error('X MUST HAVE VALUES OF 1 OR 0 ONLY.'); end

if N < 20
    warning(['Running dependent binomial test requires many subjects to approximate normal distribution.  You are only using ' num2str(length(unique(subject_id))) ' subjects.  Interpret results with caution.'])
end

% -------------------------------------------------------------------------
% Calculate Variance and Covariance
% -------------------------------------------------------------------------
% First, calculate variance of the quantity we want for a single trial,
% using Pa and Pb as estimates for the binomial parameters
% (This is not the SD of the SAMPLING distribution, but the SD of the
% distribution of interest (e.g., variance of sum, difference)
% Later, we will use this to get the SD for the SAMPLING distribution,
% which is used to construct confidence intervals and tests.

%Calculate Probabilities of each Subject Trial
Pa = sum(A) ./ N;  % P-hat for A
Pb = sum(B) ./ N;  % P-hat for B

% Calculate Variance for each Subject Trial
Va = Pa * (1 - Pa); % ./ N;
Vb = Pb * (1 - Pb); % ./ N;

% Calculate Covariance between Trials
% C = ( (sum(A .* B) ./ N) * (sum(~A .* ~B) ./ N) ) - ( (sum(A .* ~B) ./ N) * (sum(B .* ~A) ./ N) ); % % Agresti, p. 410 SAME as below
C = ( (sum(A .* B) ./ N) - (Pa * Pb) );   % Covariance - expected prop under independence

% Calculate Variance for Test Statistics adjusting for covariance
Vaplusb = Va + Vb - (2 * C);
Vameanb = ( Va + Vb - (2 * C) ) / 4;  % / 4 because V(aX) = a^2Var(X)
Vaminusb = Va + Vb + (2 * C);

% -------------------------------------------------------------------------
% Define Functions
% -------------------------------------------------------------------------

sdfunction = @(V, N) (V ./ N) .^ .5;  % Calculate Standard Deviation
zfunc = @(Px, Po, SD) (Px - Po) ./ SD; % Calculate Z-Value
pfunction_gauss = @(z) 2 .* min(normcdf(z, 0, 1), normcdf(z, 0, 1, 'upper')); %Calculate p-value for z-distribution

% pfunction_bino = @(N, Po, hits) 2 * min(binocdf(hits, N, Po), (1 - binocdf(hits - 1, N, Po))); %Calculate p-value for binomial distribution

% -------------------------------------------------------------------------
% Define Various Confidence Intervals - TO DO - See Agresti Book
% -------------------------------------------------------------------------

% switch  cimethod
%     case 'Wald'
%
%         cifunc = @(z)
%
%     case 'scoring'  % default
%
%     otherwise
%         error('Unknown CI method. Check inputs.');
% end

% Now, we will use this to get the SD for the SAMPLING distribution,
% which is used to construct confidence intervals and tests.
% e.g., SDa = (Va ./ N) .^ .5;

% -------------------------------------------------------------------------
% Calculate test statistics
% -------------------------------------------------------------------------
% Use z-test (large-sample approximation; ok for N = 20+)

names = {'P1' 'P2' 'P1minusP2' 'PavgP1P2' 'P1plusP2' };
est = [Pa Pb Pa-Pb mean([Pa Pb]) Pa+Pb];
Vs = [Va Vb Vaminusb Vameanb Vaplusb];
Pnull = [Po Po 0 Po 2*Po];

for i = 1:length(names)
    SDs(i) = sdfunction(Vs(i), N);
    Z(i) = zfunc(est(i), Pnull(i), SDs(i));
    pval(i) = pfunction_gauss(Z(i));
    %ci(i) = add confidence intervals

    RES{i} = struct('n', N, 'hits', est(i).*N, 'prop', est(i), 'SE', SDs(i), 'Z', Z(i), 'p_val', pval(i));

    varargout{i} = RES{i};
end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [result] = analyze_accuracy_robustness(experiment, trackers, sequences, varargin)
% analyze_accuracy_robustness Performs ranking analysis
%
% Performs ranking analysis for a given experiment on a set trackers and sequences.
%
% Input:
% - experiment (structure): A valid experiment structures.
% - trackers (cell): A cell array of valid tracker descriptor structures.
% - sequences (cell): A cell array of valid sequence descriptor structures.
% - varargin[Ranking] (boolean): Perform ranking analysis, used by default.
% - varargin[Tags] (cell): An array of tag names that should be used
% instead of sequences.
% - varargin[UsePractical] (boolean): Use practical difference for accuracy.
% - varargin[Alpha] (double): Statistical significance parameter.
% - varargin[Adaptation] (string): Type of rank adaptation. See
% adapter_ranks for more details.
%
% Output:
% - result (structure): A structure with the following fields
%     - accuracy
%          - values: average overlap matrix
%          - ranks: accuracy ranks matrix (if applicable)
%     - robustness
%          - values: number of failures matrix
%          - normalized: number of failures matrix (normalized)
%          - ranks: robustness ranks matrix (if applicable)
%     - lengths: number of frames for individual selectors
%     - tags: names of individual selectors

    usepractical = false;
    ranking = true;
    tags = {};
	adaptation = 'mean';
    alpha = 0.05;

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

    print_text('Ranking analysis for experiment %s ...', experiment.name);

    if ~any(strcmp(experiment.type, {'supervised', 'realtime'}))
        error('Ranking analysis can only be used in supervised experiment scenario.');
    end;

    if experiment.parameters.repetitions < 5 && ranking
        error('The experiment specifies less than 5 repetitions. Not valid for statistical consideration.');
    end

    if experiment.parameters.repetitions < 15 && ranking
        print_text('Warning: the experiment specifies less than 15 repetitions, the results may be statistically unstable.');
    end;

    print_indent(1);

    experiment_sequences = convert_sequences(sequences, experiment.converter);

    if ~isempty(tags)

        tags = unique(tags); % Remove any potential duplicates.

        selectors = sequence_tag_selectors(experiment, ...
            experiment_sequences, tags);

    else

        selectors = sequence_selectors(experiment, experiment_sequences);

    end;

    if ~ranking
        % Disable ranking
        alpha = -1;
    end;

    [accuracy, robustness, lengths] = trackers_ar(experiment, trackers, ...
        experiment_sequences, selectors, alpha, usepractical, adaptation);

    result = struct('accuracy', accuracy, 'robustness', robustness, 'lengths', lengths);
    result.tags = cellfun(@(x) x.name, selectors, 'UniformOutput', false);

    print_indent(-1);

end

function [accuracy, robustness, lengths] = trackers_ar(experiment, trackers, ...
    sequences, selectors, alpha, usepractical, adaptation)

    N_trackers = length(trackers) ;
    N_selectors = length(selectors) ;

    accuracy.values = zeros(N_selectors, N_trackers);
    robustness.values = zeros(N_selectors, N_trackers);
    robustness.normalized = zeros(N_selectors, N_trackers);

    if alpha > 0
        accuracy.ranks = zeros(N_selectors, N_trackers);
        robustness.ranks = zeros(N_selectors, N_trackers);
    end;

    lengths = zeros(N_selectors, 1);

    for a = 1:length(selectors)

	    print_indent(1);

	    print_text('Processing selector %s ...', selectors{a}.name);

        % rank trackers and calculate statistical significance of differences
        [average_overlap, average_failures, average_failurerate, HA, HR, available] = ...
            trackers_raw_scores_selector(experiment, trackers, sequences, selectors{a}, alpha, usepractical);

        accuracy.values(a, :) = average_overlap;
        robustness.values(a, :) = average_failures;
        robustness.normalized(a, :) = average_failurerate;

        if alpha > 0

            [~, order_by_accuracy] = sort(average_overlap(available), 'descend');
            accuracy_ranks = ones(size(available)) * length(available);
            [~, accuracy_ranks(available)] = sort(order_by_accuracy, 'ascend') ;

            [~, order_by_robustness] = sort(average_failures(available), 'ascend');
            robustness_ranks = ones(size(available)) * length(available);
            [~, robustness_ranks(available)] = sort(order_by_robustness,'ascend');

            % get adapted ranks
            adapted_accuracy_ranks = adapted_ranks(accuracy_ranks, HA, adaptation);
            adapted_robustness_ranks = adapted_ranks(robustness_ranks, HR, adaptation);

            % mask out results that are not available
            adapted_accuracy_ranks(~available) = nan;
            adapted_robustness_ranks(~available) = nan;
            accuracy.ranks(a, :) = adapted_accuracy_ranks;
            robustness.ranks(a, :) = adapted_robustness_ranks;

        end;

        lengths(a) = selectors{a}.length(sequences);

	    print_indent(-1);

    end

end

function [average_accuracy, average_failures, average_failurerate, HA, HR, available] ...
    = trackers_raw_scores_selector(experiment, trackers, sequences, selector, alpha, usepractical)

    cacheA = cell(length(trackers), 1);
    cacheR = cell(length(trackers), 1);

    HA = false(length(trackers)); % results of statistical testing
    HR = false(length(trackers)); % results of statistical testing

    average_accuracy = nan(length(trackers), 1);
    average_failures = nan(length(trackers), 1);
    average_failurerate = nan(length(trackers), 1);

    available = true(length(trackers), 1);

    if alpha > 0 && usepractical
        practical = selector.groundtruth_values(sequences, 'practical');
        practical = cat(1, practical{:});
    else
        practical = [];
    end

	print_indent(1);
    [~, lengths] = selector.length(sequences);

    for t1 = 1:length(trackers)

		print_text('Processing tracker %s ...', trackers{t1}.identifier);

        if isempty(cacheA{t1})
            [O1, F1] = calculate_accuracy_overlap(selector, experiment, trackers{t1}, sequences);
            cacheA{t1} = O1; cacheR{t1} = F1;
        else
            O1 = cacheA{t1}; F1 = cacheR{t1};
        end;

        if isempty(O1)
            available(t1) = false;
			HA(t1, :) = true; HA(:, t1) = true;
			HR(t1, :) = true; HR(:, t1) = true;
            HA(t1, t1) = false; HR(t1, t1) = false;
            continue;
        end

        valid_frames = ~isnan(O1) ;

        % O1 ... stacked per-frame overlaps (already averaged over
        % repeats).
        %
        % F1 ... fragments (rows) x repeats (columns) of raw failure count.

        % Average accuracy is average over valid frames (non NaN).
        if all(valid_frames == 0)
            average_accuracy(t1) = 0;
        else
            average_accuracy(t1) = mean(O1(valid_frames));
        end;

        % Average failures are sum of failures in fragments averaged over
        % repetitions
        average_failures(t1) = mean(sum(F1, 1));

        % Average failure rate is sum of failures in fragments divided by
        % total length of selector, averaged over repetitions
        average_failurerate(t1) = mean(sum(F1, 1) ./ sum(lengths));

        if alpha > 0

            for t2 = t1+1:length(trackers)

                t2

                if isempty(cacheA{t1})
                    [O1, F1] = calculate_accuracy_overlap(selector, experiment, trackers{t1}, sequences);
                    cacheA{t1} = O1; cacheR{t1} = F1;
                else
                    O1 = cacheA{t1}; F1 = cacheR{t1};
                end;

                if isempty(cacheA{t2})
                    [O2, F2] = calculate_accuracy_overlap(selector, experiment, trackers{t2}, sequences);
                    cacheA{t2} = O2; cacheR{t2} = F2;
                else
                    O2 = cacheA{t2}; F2 = cacheR{t2};
                end;

                if isempty(O2)
                    available(t2) = false;
                    continue;
                end

                % If alpha is 0 then we disable the equivalence testing
                if alpha == 0

                    ha = true; hr = true;

                else

                    [ha, hr] = test_significance(O1, F1, O2, F2, alpha, practical);

                end;

                HA(t1, t2) = ha; HA(t2, t1) = HA(t1, t2);
                HR(t1, t2) = hr; HR(t2, t1) = HR(t1, t2);
            end;
        end;
    end;

	print_indent(-1);

    average_accuracy(isnan(average_accuracy)) = 0;

end

function [aggregated_overlap, aggregated_failures] = calculate_accuracy_overlap(selector, experiment, tracker, sequences)

    aggregated_overlap = [];
    aggregated_failures = [];

    burnin = experiment.parameters.burnin;

    groundtruth = selector.groundtruth(sequences);
    trajectories = selector.results(experiment, tracker, sequences);

    repeat = experiment.parameters.repetitions;

    for s = 1:numel(groundtruth)

        accuracy = nan(repeat, length(groundtruth{s}));
        failures = nan(repeat, 1);

        for r = 1:size(trajectories, 2)

            if isempty(trajectories{s, r})
                continue;
            end;

            [~, frames] = estimate_accuracy(trajectories{s, r}, groundtruth{s}, 'burnin', burnin, 'BindWithin', [sequences{s}.width, sequences{s}.height]);

            accuracy(r, :) = frames;

            failures(r) = estimate_failures(trajectories{s, r}, groundtruth{s});

        end;

        frames = num2cell(accuracy, 1);
        sequence_overlaps = cellfun(@(frame) nanmean(frame), frames);

        failures(isnan(failures)) = nanmean(failures);
        sequence_failures = failures';

        if ~isempty(sequence_overlaps)
            aggregated_overlap = [aggregated_overlap, sequence_overlaps]; %#ok<AGROW>
        end;

        if ~isempty(sequence_failures)
            aggregated_failures = [aggregated_failures; sequence_failures]; %#ok<AGROW>
        end;

    end;

    aggregated_failures = aggregated_failures(~isnan(aggregated_failures(:, 1)), :);

end

function [ha, hr, hp] = test_significance(A1, R1, A2, R2, alpha, practical)
% test_significance Verify difference of A-R performance for two trackers
%
% Compare A-R performance of two trackers taking into account statistical
% and practical difference of results.
%
% Input:
% - A1 (double matrix): Per-frame accuracy for first tracker
% - R1 (double matrix): Per-segment robustness for first tracker
% - A2 (double matrix): Per-frame accuracy for second tracker
% - R2 (double matrix): Per-segment robustness for second tracker
% - alpha (double): Confidence parameter
% - practical (boolean): Take into account practical difference for the
% frames
%
% Output:
% - ha (boolean): Is accuracy different
% - hr (boolean): Is robustness different
% - hp (number): Number of frames for which the practical difference test
% was positive
%

    % Testing accuracy significance

    % Statistical test
    dif = A1 - A2;
    valid = ~isnan(dif);
    dif = dif(valid) ;
    if (length(dif) < 25)
        print_text('Warning: less than 25 samples when comparing trackers. Cannot reject hypothesis');
        ha = 1;
    else
        if (is_octave)
            try
            pa = wilcoxon_test(A1(valid), A2(valid));
            ha = (pa <= alpha);
            catch
                %A1(valid) - A2(valid)
                %pa = 0;
                ha = 1;
            end;
        else
            [~, ha, ~] = signrank(dif, [], 'alpha', alpha ) ;
        end;
    end;

    hp = 0;

    % Practical difference of accuracy
    if ~isempty(practical)
        hp = sum(dif' < practical(valid));
        if abs(mean(dif' ./ practical(valid))) < 1
            ha = 0;
        end;

    end;

    % Testing robustness significance
    R1 = R1(:);
    R2 = R2(:);
	%R1 = mean(R1,2);
    %R2 = mean(R2,2);

    if (is_octave)
       pr = u_test(R1, R2);
       hr = (pr <= alpha);
    else
       [~, hr] = ranksum(R1, R2, 'alpha', alpha) ;
    end;

end

author: votchallenge

repo: toolkit-legacy


text:
%kkpoisson 'Introduce Poisson Noise in Input Object'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kpoisson.pane file
%
% Parameters:
% InputFile: i 'Input File', required: 'input file'
% OutputFile: o 'Output File', required: 'output file'
% Toggle: replace 'Replace Data with Noise', default: 0: 'Replace existing values with noise'
% Toggle: add 'Add Noise to Data', default: 0: 'Add noise to the specified segment'
%
% Example: o = kkpoisson(i, {'i','';'o','';'replace',0;'add',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kpoisson - Introduce Poisson Noise in Input Object
%
%  DESCRIPTION
% "kpoisson" introduces Poisson noise of the specified parameters in the
% input data object (-i).  The noise can either replace the existing value
% data or it can be added to the existing value data by specifying
% the appropriate flag (-replace or -add).
% The time (-ptime) and the variance (-pvar) of the Poisson distribution
% must be specified.
%
%  "Data Type" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Map Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "Validity Mask" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Explicit Location and Time Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Failure Modes"
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::knoise, DATAMANIP::kexpon, DATAMANIP::kgauss,
% DATAMANIP::krayleigh, DATAMANIP::kuniform
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkpoisson(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkpoisson(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'replace', 0;'add', 0};
maxval={0,0,0,0};
minval={0,0,0,0};
istoggle=[0,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'knoise"  -poiss'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function e=newoenvel(f,u,j,i);
%
% Finds the upper envelope of a 'function' f(u) between u(i) and u(j)
%
if i>j,
  [j i]
  error('i>j');
end;
if i==j,
  e=[i]; 	    	%Terminate with a 'shock'
elseif i+1==j,
  e=[i j]; 		%Terminate with a 'rarefaction'
else
  du=u(i)-u(i+1:j);
  df=f(i)-f(i+1:j);
  s=df./du;
  [m k]=max(s);			%New point on envelope
  e=[i newoenvel(f,u,j,i+k)];	%Find the rest of the envelope
end;

repo: Aero-matlab


text:
 function ni = trl_curvature(yi, bi, ri, li, ctype)
%function ni = trl_curvature(yi, bi, ri, li, ctype)
%|
%| Compute surrogate parabola curvatures for Poisson transmission model
%| ctype:
%|	oc	erdogan's optimal curvatures
%|	pc	precomputed curvatures, ala trpl3
%|	nc	newton curvatures
%|
%| fix: align with C version of trpl2,3
%| The minimum returned curvature will be zero.
%| It is the user's responsibility to impose additional bounds
%| if desired for certain algorithms.
%|
%| Copyright 2002-1-28, Jeff Fessler, University of Michigan

if nargin == 1 && streq(yi, 'test'), trl_curvature_test, return, end
if nargin < 4, ir_usage, end

[h dh] = trl_h_dh;

if ~isvar('ctype') || isempty(ctype)
	ctype = 'oc'
end

% Compute optimal surrogate parabola curvatures
% for Poisson transmission model based on Erdogan's formula.
switch ctype
case 'oc'

	% compute curvature at l=0
	ni_max = zeros(size(yi));
	if numel(bi) == 1 % scalar bi (must be positive!)
		ni_max = bi .* (1 - yi .* ri ./ (bi + ri).^2);
	else
		i0 = bi > 0;
		if numel(ri) == 1
			rii = 1;
		else
			rii = ri(i0);
		end
		ni_max(i0) = bi(i0) .* (1 - yi(i0) .* rii ./ (bi(i0) + rii).^2);
	end
	ni_max = max(ni_max, 0);
	ni = ni_max;

	if 0
		il0 = li <= 0;
	else % trick in C program due to numerical precision issues
		il0 = li < 0.1;
	end

	tmp = h(yi,bi,ri,li) - h(yi,bi,ri,0) - li .* dh(yi,bi,ri,li);
	i = ~il0;
	ni(i) = 2 ./ li(i).^2 .* max(tmp(i),0);

	if any(ni > ni_max)
	%	plot([ni_max(:) ni(:) ni(:)>ni_max(:)])
		warning 'large ni'
	end

% Precomputed approximating parabola curvatures
% for Poisson transmission model.
% The minimum returned curvature will be zero.
% This is compatible with trpl/trp_init_der02_sino() in aspire.
case 'pc'

	% ni = (yi-ri)^2 / yi, if yi > ri >= 0 and bi > 0
	ii = (yi > ri) & (ri >= 0) & (bi > 0); % good rays
	ni = zeros(size(yi));
	ni(ii) = (yi(ii) - ri(ii)).^2 ./ yi(ii);

% newton curvatures (current 2nd derivative)
case 'nc'
	bel = bi .* exp(-li);
	yb = bel + ri;
	ni = (1 - ri.*yi./yb.^2) .* bel;

otherwise
	fail 'bug'
end

% trl_h_dh()
% transmission Poisson likelihood function
function [h, dh] = trl_h_dh
h = @(y,b,r,l) y .* log(b.*exp(-l)+r) - (b.*exp(-l)+r);
dh = @(y,b,r,l) (1 - y ./ (b.*exp(-l)+r)) .* b.*exp(-l);

% trl_curvature_test()
% demonstrate an example surrogate parabola!
function trl_curvature_test
[h dh] = trl_h_dh;
if 0
	l = linspace(0,2,101)';
	n = trl_curvature(2+0*l, 3, 3, l, 'oc');
	plot(l, n, '-o'), xlabel l, ylabel n
else
%	y = 4; b = 3; r = 1; ln = 3.8;
%	y = 2; b = 3; r = 3; ln = 0.05;
	y = 2; b = 3; r = 1; ln = 1;
	n.oc = trl_curvature(y, b, r, ln, 'oc');
	n.pc = trl_curvature(y, b, r, ln, 'pc');
	l = linspace(-0.2,5,101)';
	hn = h(y,b,r,ln);
	dhn = dh(y, b, r, ln);%
	hl = h(y,b,r,l);
	q.oc = hn + dhn * (l-ln) - 0.5 * n.oc * (l-ln).^2;
	q.pc = hn + dhn * (l-ln) - 0.5 * n.pc * (l-ln).^2;
	if im
		clf, plot(l, hl, '-', l, q.oc, '--', l, q.pc, '-.', ln, hn, 'o')
		grid
		axis([minmax(l)' minmax(hl)'])
		legend('h(l)',	sprintf('q_o(l) n_i=%g', n.oc), ...
				sprintf('q_p(l) n_i=%g', n.pc))
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [result_features]=back_find_low_node_para(result_dir,nSubj,k_times,nROI,w,cross_val,ttest_p,midw_lasso)

% This function aims to find features used in the classification, suitable
% for the SR,SLR,SGR,GSR,WSR,WSGR,SSGSR network construction method. And it is also
% suitable for the LOOCV or the 10-fold cross validation. The found
% features can be applied on the visualization software to present the
% important brain regions or links.
% Input:
%         result_dir: the directory you want to store all the result files;
%         nSubj: number of subjects;
%         k_times: times of 10-fold cross validation;
%         nROI: number of ROIs
%         w: weight of each selected features;
%         cross_val: loocv or 10-fold;
%         ttest_pm,midw_lasso: the feature selection indexes of ttest or
%         lasso;
%
% Output:
%        result_features: cell array, consisting of the averaged weight and occurrence of the features;

% Written by Zhen Zhou, zzstefan@email.unc.edu
% IDEA lab, https://www.med.unc.edu/bric/ideagroup
% Department of Radiology and BRIC, University of North Carolina at Chapel Hill

% clear all
% load SR_10_fold_middle.mat;
% load SR_loocv_middle.mat;
% cross_val='loocv';
% ttest_p=feature_index_ttest;
% midw_lasso=feature_index_lasso;
fprintf('Begin contributing feature identification\n');
if strcmpi(cross_val,'10-fold')
    for i=1:size(ttest_p,1)
        for j=1:size(ttest_p,2)
            tmp=find(ttest_p{i,j}<0.05);
            tmp_b{i,j}=find(midw_lasso{i,j});
            index=tmp(tmp_b{i,j});
            for k=1:length(index)
                [first{i,j,k}(1),first{i,j,k}(2)]=find_elements(nROI,index(k));
            end
        end
    end
elseif strcmpi(cross_val,'loocv')
    for i=1:length(ttest_p)
        tmp=find(ttest_p{i}<0.05);
        tmp_b{i}=find(midw_lasso{i});
        index=tmp(tmp_b{i});
        for k=1:length(index)
            [first{i,k}(1),first{i,k}(2)]=find_elements(nROI,index(k));
        end
    end
end

%first=cellfun(@(x) x',first,'UniformOutput',false);
temp=first(:);

[a1,b,c]=unique(cellfun(@char,temp,'un',0));
lo=histc(c,1:max(c));
loo=lo(:)>1;
out=[temp(b(loo)),num2cell(lo(loo))];
out(1,:)=[];

All_link=out(find(cell2mat(out(:,2))>0),1);%find all links which exists in any one  LOOCV fold
if strcmpi(cross_val,'10-fold')
    for j=1:length(All_link)
        weight=[];
        temp_index=cell2mat(cellfun(@(x)isequal(x,All_link{j}),first,'un',0));
        ind=find(temp_index);
        [A,B,C]=ind2sub(size(temp_index),ind);
        for k=1:length(A)
            weight=[weight,w{A(k),B(k)}(C(k))];
        end
        all_weight(j)=mean(weight);
    end
elseif strcmpi(cross_val,'loocv')
    for j=1:length(All_link)
        weight=[];
        temp_index=cell2mat(cellfun(@(x)isequal(x,All_link{j}),first,'un',0));
        [A,B]=find(temp_index);
        for k=1:length(A)
            weight=[weight,w{A(k)}(B(k))];
        end
        all_weight(j)=mean(weight);
    end
end

matrix_1=zeros(nROI,nROI);
matrix_2=zeros(nROI,nROI);
for m=1:size(out,1)
    matrix_1(out{m,1}(1),out{m,1}(2))=all_weight(m);
    matrix_1(out{m,1}(2),out{m,1}(1))=all_weight(m);
    matrix_2(out{m,1}(1),out{m,1}(2))=out{m,2};
    matrix_2(out{m,1}(2),out{m,1}(1))=out{m,2};
end
result_features{1}=matrix_1;

if strcmpi(cross_val,'loocv')
    matrix_2=matrix_2/nSubj;
elseif strcmpi(cross_val,'10-fold')
    matrix_2=matrix_2/(10*k_times);
end
result_features{2}=matrix_2;

figure('visible','off');
subplot(1,2,1);
imagesc(matrix_1);
colormap jet
colorbar
axis square
xlabel('ROI');
ylabel('ROI');
title('Averaged weight');

subplot(1,2,2);
imagesc(matrix_2);
colormap jet
colorbar
axis square
xlabel('ROI');
ylabel('ROI');
title('Normalized occurence');
print(gcf,'-r1000','-dtiff',char(strcat(result_dir,'/result_features_weigthAndOccurence.tiff')));

fprintf('End contributing feature identification\n');

author: zzstefan

repo: BrainNetClass


text:
function [ptbifu, ptroot] = points_init(bw)

% This function initializes the feature points from the binary image

[M, N]= size(bw);
imdim = M*N + 1;

Neighbor = [-M, 1, M, -1, -1-M, 1-M,  1+M, -1+M];
Len = prod(size(Neighbor));

seeds = find(bw == 1);
npix = prod(size(seeds));
countmap = zeros(size(bw));

for k =1:npix
    localidx = seeds(k);
    neighidx = localidx + Neighbor;
    for i=1:Len
        idx = neighidx(i);
        if (idx>0) & (idx<imdim) & (bw(idx) ==1)
            countmap(localidx) = countmap(localidx)+1;
        end
    end
end

% The bifurcation candidates may connect, but only the central is true
% bifurcation point
bw1 = (countmap>=3);
[labelmap, numlabel] = bwlabel(bw1, 8);
for i=1:numlabel
    candidates = find(labelmap==i);
    ptbifu(i) = candidates(fix((end+1)/2));
end

ptroot = find(countmap == 1);

ptbifu = ptbifu(:);
ptroot = ptroot(:);

author: Sable

repo: mcbench-benchmarks


text:
function report_sectors(XS,str)

ntickers = size(XS,1);

[Ctry,ia,ic] = unique(XS(:,3)); % ia (index of first appearance)
N_sectors = length(Ctry);
fprintf('\n');
fprintf('%d %s Sectors for %d stocks:\n',N_sectors,str,ntickers);
for i=1:N_sectors
  n_ic = sum(ic==i);
  fprintf(' %s (%d = %.1f%%)\n',Ctry{i},n_ic,100*n_ic/length(ic));
end
fprintf('\n');

%{
[Ctry,ia,ic] = unique(XS(:,4)); % ia (index of first appearance)
N_industries = length(Ctry);
fprintf('\n');
fprintf('%d %s Industries for %d stocks:\n',N_industries,str,ntickers);
for i=1:N_industries
  n_ic = sum(ic==i);
  fprintf(' %s (%d = %.1f%%)\n',Ctry{i},n_ic,100*n_ic/length(ic));
end
fprintf('\n');
%}

end

author: gudbrandtandberg


text:
function [ss,ssfun]=moesp(y,u,d)
% MOESP Multivariable Output Error State Space Approach of Subspace Identification
%   [SS,SSMAT] = MOESP(Y,U,D) identifies the observable subspace based on
%   the measured (N by ny) output matrix, Y with N  sampling points and ny
%   variables, and the corresponding N by nu input matrix, U. The third
%   parameter d is the embeded dimension, which should be larger than the
%   order of the system to be identified.
%
%   The function returns the scores of the subspace, SS and a function
%   handle SSMAT, for further identification of state space matrices.
%
%   The system order can determined from the returned score vector, SS so
%   that sum(SS(1:n)) ~ sum(SS).
%
%   Once the system order is determined, the underline dynamic system is
%   identified by calling the retruned function handle:
%
%   [A,B,C,D] = SSMAT(n)
%
%   to represent a state space model:
%
%       x(k+1) = Ax(k) + Bu(k)
%       y(k)   = Cx(k) + Du(k)
%
% See also: n4sid, subid

% Version 1.0 by Yi Cao at Cranfield University on 27th April 2008

% Reference
% Tohru Katayama, "Subspace Methods for System Identification", Springer
% 2005.

% Example
%{
%   Consider a multivariable fourth order system a,b,c,d
%   with two inputs and two outputs:
    a = [0.603 0.603 0 0;-0.603 0.603 0 0;0 0 -0.603 -0.603;0 0 0.603 -0.603];
    b = [1.1650,-0.6965;0.6268 1.6961;0.0751,0.0591;0.3516 1.7971];
    c = [0.2641,-1.4462,1.2460,0.5774;0.8717,-0.7012,-0.6390,-0.3600];
    d = [-0.1356,-1.2704;-1.3493,0.9846];
%
%   We take a white noise sequence of 1000 points as input u.
    N = 1000;
    u = randn(N,2);
%
%   With noise added, the state space system equations become:
%                  x_{k+1) = A x_k + B u_k + K e_k
%                    y_k   = C x_k + D u_k + e_k
%                 cov(e_k) = R
%
    k = [0.1242,-0.0895;-0.0828,-0.0128;0.0390,-0.0968;-0.0225,0.1459]*4;
    r = [0.0176,-0.0267;-0.0267,0.0497];
%
%   The noise input thus is equal to (the extra chol(r) makes cov(e) = r):
    e = randn(N,2)*chol(r);
%
%   And the simulated noisy output:
    y = dlsim(a,b,c,d,u) + dlsim(a,k,c,eye(2),e);
%
%   Using this output in subid returns a more realistic image of
%   the singular value plot:
    k = 10;
    [ss,ssfun] = moesp(y,u,k);
%
%   To determine the order, check the score vector
%
bar(ss)
%
%   Clearly, the order should be 4. Therefore, the state space model is
%   obtained by calling ssfun:
%
    [A,B,C,D]=ssfun(4);
%
% To compare with the actual system, we use the bode diagram:
    w = [0:0.005:0.5]*(2*pi); 		% Frequency vector
    m1 = dbode(a,b,c,d,1,1,w);
    m2 = dbode(a,b,c,d,1,2,w);
    M1 = dbode(A,B,C,D,1,1,w);
    M2 = dbode(A,B,C,D,1,2,w);
%
% Plot comparison
    figure(1)
    hold off;subplot;clg;
    subplot(221);plot(w/(2*pi),[m1(:,1),M1(:,1)]);title('Input 1 -> Output 1');
    subplot(222);plot(w/(2*pi),[m2(:,1),M2(:,1)]);title('Input 2 -> Output 1');
    subplot(223);plot(w/(2*pi),[m1(:,2),M1(:,2)]);title('Input 1 -> Output 2');
    subplot(224);plot(w/(2*pi),[m2(:,2),M2(:,2)]);title('Input 2 -> Output 2');
%}

% Input and output check
error(nargchk(1,3,nargin));
error(nargoutchk(0,4,nargout));

[ndat,ny]=size(y);
[mdat,nu]=size(u);
if ndat~=mdat
    error('Y and U have different length.')
end

% block Hankel matrix
N=ndat-d+1;
Y = zeros(d*ny,N);
U = zeros(d*nu,N);
sN=sqrt(N);
sy=y'/sN;
su=u'/sN;
for s=1:d
    Y((s-1)*ny+1:s*ny,:)=sy(:,s:s+N-1);
    U((s-1)*nu+1:s*nu,:)=su(:,s:s+N-1);
end

% LQ decomposition
R=triu(qr([U;Y]'))';
R=R(1:d*(ny+nu),:);

% SVD
R22 = R(d*nu+1:end,d*nu+1:end);
[U1,S1]=svd(R22);

% sigular value
ss = diag(S1);
% n=find(cumsum(ss)>0.85*sum(ss),1);

ssfun = @ssmat;

    function [A,B,C,D]=ssmat(n)
        % C and A
        Ok = U1(:,1:n)*diag(sqrt(ss(1:n)));
        C=Ok(1:ny,:);
        A=Ok(1:ny*(d-1),:)\Ok(ny+1:d*ny,:);

        % B and D
        L1 = U1(:,n+1:end)';
        R11 = R(1:d*nu,1:d*nu);
        R21 = R(d*nu+1:end,1:d*nu);
        M1 = L1*R21/R11;
        m = ny*d-n;
        M = zeros(m*d,nu);
        L = zeros(m*d,ny+n);
        for k=1:d
            M((k-1)*m+1:k*m,:)=M1(:,(k-1)*nu+1:k*nu);
            L((k-1)*m+1:k*m,:)=[L1(:,(k-1)*ny+1:k*ny) L1(:,k*ny+1:end)*Ok(1:end-k*ny,:)];
        end
        DB=L\M;
        D=DB(1:ny,:);
        B=DB(ny+1:end,:);
    end
end


text:
function xhfd=hfd(x,kmax)
%function xhfd=hfd(x,kmax)
%Input:
%x: (either column or row) vector of length N
%kmax: maximum value of k
%Output:
%xhfd: Higuchi fractal dimension of x

if ~exist('kmax','var')||isempty(kmax),
    kmax=5;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=x(:)';
N=length(x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lmk=zeros(kmax,kmax);
for k=1:kmax,
    for m=1:k,
        Lmki=0;
        for i=1:fix((N-m)/k),
            Lmki=Lmki+abs(x(m+i*k)-x(m+(i-1)*k));
        end;
        Ng=(N-1)/(fix((N-m)/k)*k);
        Lmk(m,k)=(Lmki*Ng)/k; % Here is the problem in the code by Mr. Tikkuhirvi & Mr. Aino
    end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lk=zeros(1,kmax);
for k=1:kmax,
    Lk(1,k)=sum(Lmk(1:k,k))/k;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lnLk=log(Lk);
lnk=log(1./[1:kmax]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b=polyfit(lnk,lnLk,1);
xhfd=b(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function [ dz ] = odesys( ~,z,param )
% Copyright 2011 The MathWorks, Inc.
% Uses natural frequency (wn) and damping ratio (zeta) to form state
% equations for a standard 2nd order system.  The time derivative of
% the state vector is:

dz = [0 1; -param.wn^2 -2*param.zeta*param.wn]*z;

end