text:
%   CLASS File_Rinex
% =========================================================================
%
% DESCRIPTION
%   Class to store file_paths for RINEX files
%
% EXAMPLE
%   fr = File_Rinex(file_name, verbosity_lev);
%
% FOR A LIST OF CONSTANTS and METHODS use doc File_Rinex

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    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_Rinex < Exportable

    properties (SetAccess = protected, GetAccess = protected)
        log = Logger.getInstance(); % Handler to the log object
    end

    properties (SetAccess = protected, GetAccess = public)
        is_valid = false;                            % flag, if true it means that the object contains at least one valid rinex file
        base_dir = {'../data/project/default_DD/RINEX/'};           % directory containing all the files
        file_name_list = {'yamatogawa_rover', 'yamatogawa_master'}; % file names (they can be multiple files for different days)
        ext = {'.obs', '.obs'};                                     % file names extension (they can be multiple files for different days)
        is_valid_list = false(1, 2);                 % for each element of file_name_list check the validity of the file

        is_composed = false;                         % when this flag is set, it means that the file_name depends on variables such as DOY DOW YYYY SSSS MM ecc...

        trck_availability                            % boolean to get tracking availabilty the tracking are in Core_sky.group_delays_flags

        marker_name = {};                            % marker name of the files
        coo = Coordinates.fromXYZ([0 0 0])           % receiver coordinates
        first_epoch = GPS_Time();                    % first epoch stored in the RINEX (updated after checkValidity)
        last_epoch = GPS_Time();                     % last epoch stored in the RINEX (updated after checkValidity)

        verbosity_lev = 0;                           % Level of verbosity  0 = all messages, see Logger for more information (if negative) no messages at all

        eoh = 0;                                     % end of header line - store the last line of the header
    end

    properties (SetAccess = private, GetAccess = private)
        id_date = 2:28;                              % Last character containing the 6 fields of the date (in case of 4digits year), it the pends on the type of rinex files (28 -> OBS RINEX 3)
    end

    methods
        function this = File_Rinex(file_name, verbosity_lev, flag_header_only)
            % Creator of file_rinex simple parser
            %
            % INPUT
            %   file_name        file path (can be a cell array)
            %   verbosity_lev    verbosity level for logging purposes
            %                    (default = 100)
            %   flag_header_only if true and first/last epochs are not found
            %                    in header do not search them in the file
            %                    (default = false)
            %                    this option may leave the obj corrupted (use with care)
            %
            % SYNTAX
            %   File_Rinex (file_name , verbosity_level, flag_header_only)
            %
            if nargin < 3 || isempty(flag_header_only)
                flag_header_only = false;
            end

            if nargin == 0
                % Empty File_Rinex;
            else
                % fill the path with the imported file names
                if ~iscellstr(file_name)
                    file_name = {file_name};
                end
                this.file_name_list = {};
                this.ext = {};
                if nargin == 1 || isempty(verbosity_lev)
                    verbosity_lev = 100;
                end

                if verbosity_lev < 50 && numel(file_name) > 20
                    w_bar = Go_Wait_Bar.getInstance(numel(file_name), 'Checking rinex files');
                    w_bar.createNewBar
                    for f = 1 : numel(file_name)
                        [this.base_dir{f}, this.file_name_list{f}, this.ext{f}] = fileparts(checkPath(file_name{f}));
                        w_bar.go();
                    end
                    w_bar.close();
                else
                    for f = 1 : numel(file_name)
                        [this.base_dir{f}, this.file_name_list{f}, this.ext{f}] = fileparts(checkPath(file_name{f}));
                    end
                end

                if nargin >= 2
                    this.verbosity_lev = verbosity_lev;
                end

                this.checkValidity(flag_header_only);
            end
        end

        function copy = getCopy(this)
            for r = numel(this):-1:1
                copy(r) = File_Rinex();
                copy(r).copyFrom(this(r));
            end
        end

        function copyFrom(this, file_rinex, id)
            % Copy from an object of the same type
            %
            % SYNTAX
            %   this.copyFrom(time)
            if nargin < 3 || isempty(id)
                id = 1 : numel(file_rinex.file_name_list);
            end
            % Get id within coo and time limits objects
            valid_lid = false(size(file_rinex.is_valid_list));
            valid_lid(id) = true;
            valid_lid = valid_lid & file_rinex.is_valid_list;
            valid_id = cumsum(file_rinex.is_valid_list);
            valid_id = valid_id(valid_lid);
            % Extract fields to copy
            this.is_valid       = all(file_rinex.is_valid_list(id));
            this.base_dir       = file_rinex.base_dir(id);
            this.file_name_list = file_rinex.file_name_list(id);
            this.ext            = file_rinex.ext(id);
            this.is_valid_list  = file_rinex.is_valid_list(id);
            this.is_composed    = file_rinex.is_composed;
            this.first_epoch    = file_rinex.first_epoch.getEpoch(valid_id);
            this.last_epoch     = file_rinex.last_epoch.getEpoch(valid_id);
            this.verbosity_lev  = file_rinex.verbosity_lev;
            this.eoh            = file_rinex.eoh(id);
            this.coo            = file_rinex.coo.getCopy();
        end
    end

    methods
        function checkValidity(this, flag_header_only)
            % Update the status of validity of the files here pointed
            % Buffered read
            %
            % SYNTAX
            %   this.checkValidity();

            % pre allocate
            this.is_valid_list = false(1, numel(this.file_name_list));
            this.eoh = zeros(1, numel(this.file_name_list));
            this.first_epoch = GPS_Time();
            this.last_epoch = GPS_Time();
            this.coo = Coordinates;           % receiver coordinates
            log = Core.getLogger;
            % for each file present in the list
            for f = 1 : numel(this.file_name_list)
                marker_name = '';
                full_path = fullfile(this.base_dir{f}, [this.file_name_list{f} this.ext{f}]);

                % try to find the first and the last epoch stored in the file
                try
                    %%
                    this.is_valid_list(f) = true;
                    fid = fopen(full_path, 'rt');
                    if fid < 0
                        log.addError(['"' this.file_name_list{f} this.ext{f} '" appears to be missing'], this.verbosity_lev);
                        this.is_valid_list(f) = false;
                    else
                        buf = fread(fid, 1e4, '*char')';
                        if length(buf) > 65 && strcmp(buf(61:64), 'CRIN')
                            log.addError(sprintf('Check the following file, it seems to be hatanaka compressed\nDecompress "%s"', full_path));
                            this.is_valid_list(f) = false;
                            fclose(fid);
                        else
                            l = 1;
                            % read by buffer 10k char at a time
                            % fseek(fid, 0, 'bof'); buf = fread(fid, 1e4, '*char')';

                            % detect windows carriage return
                            if ~isempty(find(buf(1:min(1000,numel(buf))) == 13, 1, 'first'))
                                has_cr = true;  % The file has carriage return - I hate you Bill!
                            else
                                has_cr = false;  % The file is UNIX standard
                            end

                            % get new line separators
                            nl = regexp(buf, '\n')';
                            if nl(end) <  (numel(buf) - double(has_cr))
                                nl = [nl; numel(buf)];
                            end
                            lim = [[1; nl(1 : end - 1) + 1] (nl - 1 - double(has_cr))];
                            lim = [lim lim(:,2) - lim(:,1)];
                            while lim(end,3) < 3
                                lim(end,:) = [];
                            end

                            % Check for meteo rinex
                            flag_met = buf(21) == 'M';

                            line = buf(lim(l,1) : lim(l,2));
                            date_start = '';
                            trck_name = Core_Sky.group_delays_flags;
                            trck_availability = false(size(trck_name,1),1);
                            cur_trck_sys = '';
                            date_stop = '';
                            coo = '';
                            eof = false;
                            par_to_find = 4;
                            %while (par_to_find > 0) && isempty(strfind(line,'END OF HEADER')) && not(eof)
                            while ~((length(line) > 61) && (line(61) == 'E')) && not(eof)
                                l = l + 1;
                                if l >= size(lim, 1)
                                    % out of buffer
                                    % read block:
                                    buf = [buf char(fread(fid, 1e4))'];

                                    % get new line separators
                                    nl = regexp(buf, '\n')';
                                    if nl(end) <  (numel(buf) - double(has_cr))
                                        nl = [nl; numel(buf)];
                                    end
                                    lim = [[1; nl(1 : end - 1) + 1] (nl - 1 - double(has_cr))];
                                    lim = [lim lim(:,2) - lim(:,1)];
                                    while lim(end,3) < 3
                                        lim(end,:) = [];
                                    end
                                end
                                if l > size(lim, 1)
                                    eof = true;
                                else
                                    line = buf(lim(l,1) : lim(l,2));
                                    if par_to_find > 0
                                        if numel(line) >= 70 && isempty(marker_name) && line(64) == 'K' % read marker_name
                                            marker_name = strtrim(regexp(line, '.*(?=MARKER NAME)', 'match', 'once'));
                                            par_to_find = par_to_find - 1;
                                        end
                                        if numel(line) >= 76
                                            if line(61) == 'T'
                                                %tmp = regexp(line, '.*(?=GPS         TIME OF FIRST OBS)', 'match', 'once');
                                                if line(69) == 'F'
                                                    date_start = strtrim(line(1:48));
                                                    par_to_find = par_to_find - 1;
                                                elseif line(69) == 'L'
                                                    %tmp = regexp(line, '.*(?=GPS         TIME OF LAST OBS)', 'match', 'once');
                                                    date_stop = strtrim(line(1:48));
                                                    par_to_find = par_to_find - 1;
                                                end
                                            elseif line(76) == 'Y' || (numel(line) >= 79 && line(79) == 'V') %  SYS / # / OBS TYPES || # / TYPES OF OBSERV
                                                if line(1) ~= ' '
                                                    cur_trck_sys = line(1);
                                                else
                                                    cur_trck_sys = ' ';
                                                end
                                                trcks = strsplit(line(8:60),' ');
                                                for t = trcks
                                                    if ~isempty(t{1})
                                                        t = t{1};
                                                        t = [t repmat(' ',1,3-length(t))];
                                                        idx = trck_name(:,1) == cur_trck_sys & trck_name(:,2) == t(1) & trck_name(:,3) == t(2) & trck_name(:,4) == t(3);
                                                        trck_availability(idx) = true;
                                                        if strcmp(t,'P1 ')% rinex 2 codes does not specify contellation, assuming GPS
                                                            trck_availability(strLineMatch(trck_name,'GC1W')) = true;
                                                        elseif strcmp(t,'P2 ')
                                                            trck_availability(strLineMatch(trck_name,'GC2W')) = true;
                                                        elseif strcmp(t,'C1 ')
                                                            trck_availability(strLineMatch(trck_name,'GC1C')) = true;
                                                        elseif strcmp(t,'C2 ')
                                                            trck_availability(strLineMatch(trck_name,'GC2C')) = true;
                                                        end
                                                    end
                                                end
                                            elseif line(66) == 'Y' % "# / TYPES OF OBSERV"
                                                % temporary solution for RINEX 2
                                                % to be improved
                                                trcks = strsplit(line(8:60),' ');
                                                for cur_trck_sys = Core.getConstellationCollector.getActiveSysChar
                                                    for t = trcks
                                                        if ~isempty(t{1})
                                                            t = t{1};
                                                            if numel(t) == 2 % correct missing character for RINEX 2
                                                                % It should be done multi constellation
                                                                % GPS C1 -> C1C
                                                                % GPS C2 -> C2C
                                                                t = [t 'C'];
                                                            end
                                                            idx = trck_name(:,1) == cur_trck_sys & trck_name(:,2) == t(1) & trck_name(:,3) == t(2) & trck_name(:,4) == t(3);
                                                            trck_availability(idx) = true;
                                                        end
                                                    end
                                                end
                                            else
                                                %tmp = regexp(line, '.*(?=APPROX POSITION XYZ)', 'match', 'once');
                                                % character to recognize approximate position for met file: 'E' => sensor pos
                                                % character to recognize approximate position for met file: 'P' => approx position
                                                pos_char = iif(flag_met, 'E', 'P');
                                                if line(62) == pos_char % APPROXIMATE POSITION
                                                    coo = strtrim(line(1:43));
                                                    par_to_find = par_to_find - 1;
                                                end
                                            end
                                        end
                                    end
                                end
                            end

                            eoh = l; %end of header
                            if ~isempty(trck_availability)
                                this.trck_availability = trck_availability;
                            end

                            % If I did not found date of start and date of end in the file header
                            % Try to read them
                            if ~isempty(date_start)
                                this.first_epoch.addEpoch(date_start, [], true);
                            else
                                l = l + 1;
                                if l > size(lim, 1)
                                    % out of buffer
                                    % read block:
                                    buf = [buf char(fread(fid, 1e4))'];

                                    % get new line separators
                                    nl = regexp(buf, '\n')';
                                    if nl(end) <  (numel(buf) - double(has_cr))
                                        nl = [nl; numel(buf)];
                                    end
                                    lim = [[1; nl(1 : end - 1) + 1] (nl - 1 - double(has_cr))];
                                    lim = [lim lim(:,2) - lim(:,1)];
                                    while lim(end,3) < 3
                                        lim(end,:) = [];
                                    end
                                end
                                if l <= size(lim,1)
                                    epoch_line = buf(lim(l,1) : lim(l,2));
                                    % try to guess the time format
                                    [id_start, id_stop] = regexp(epoch_line, '[.0-9]*');
                                    this.id_date = id_start(1) : id_stop(6); % save first and last char limits of the date in the line -> suppose it composed by 6 fields
                                    this.first_epoch.addEpoch(epoch_line(this.id_date), [], true);
                                else
                                    error(sprintf('"%s" seems corrupted', full_path));
                                end
                            end
                            log.addStatusOk(['"' this.file_name_list{f} this.ext{f} '" appears to be a valid RINEX'], this.verbosity_lev);
                            log.addMessage(sprintf('        first epoch found at: %s', this.first_epoch.last.toString()), this.verbosity_lev);

                            if ~isempty(date_stop)
                                this.last_epoch.addEpoch(date_stop, [], true);
                                log.addMessage(sprintf('        last  epoch found at: %s', this.last_epoch.last.toString()), this.verbosity_lev);
                            else
                                if flag_header_only
                                    log.addWarning('Last epoch not found in header, search in file is not enabled\nThe last epoch has not been saved within the object', this.verbosity_lev);
                                    if ~isempty(this.first_epoch.last.getMatlabTime)
                                        this.last_epoch.addEpoch(this.first_epoch.last.getMatlabTime, [], true);
                                    end
                                else
                                    % go to the end of the file to search for the last epoch
                                    % to be sure to find at least one line containing a valid epoch, go to the end of the file minus 10000 characters
                                    fseek(fid, 0, 'eof');
                                    f_size = ftell(fid);
                                    fseek(fid, - min(1e4, f_size - 1), 'eof');
                                    buf = fread(fid, min(1e4, f_size - 1),'*char')';

                                    % Start searching for a valid epoch
                                    time = [];
                                    loop_n = 1;
                                    file_ko = false;
                                    while isempty(time) && ~isempty(buf) && ~file_ko
                                        if ~isempty(find(buf(1:min(1000,numel(buf))) == 13, 1, 'first'))
                                            has_cr = true;  % The file has carriage return - I hate you Bill!
                                        else
                                            has_cr = false;  % The file is UNIX standard
                                        end

                                        % get new line separators
                                        nl = regexp(buf, '\n')';
                                        if isempty(nl)
                                            file_ko = true;
                                        else
                                            if nl(end) <  (numel(buf) - double(has_cr))
                                                nl = [nl; numel(buf)];
                                            end
                                            lim = [[1; nl(1 : end - 1) + 1] (nl - 1 - double(has_cr))];
                                            lim = [lim lim(:,2) - lim(:,1)];
                                            lim(lim(:,3) < this.id_date(end), :) = [];

                                            l = size(lim, 1);
                                            while (l >= 1) && isempty(time)
                                                line = buf(lim(l,1) : lim(l,2));
                                                if ~isempty(regexp(line(1:15),'( [0-9]{1,4} [ 0-9]{1}[0-9]{1} [ 0-9]{1}[0-9]{1})', 'once'))
                                                    % Check that the read epoch is within 30 days from the first epoch
                                                    % (if it's further in time it's probably a misleading false epoch line)
                                                    time = GPS_Time(line(this.id_date), [], true);
                                                end
                                                l = l - 1;
                                            end
                                            if ~file_ko
                                                if isempty(time) || time.isEmpty()
                                                    loop_n = loop_n + 1;
                                                    fseek(fid, loop_n * -10000, 'eof'); % If no valid time have been found try to go back more...
                                                    buf = fread(fid, 10000, '*char')';
                                                end
                                            end
                                        end
                                    end
                                    if file_ko
                                        log.addError(sprintf('Check the following file, it seems to be corrupted\nPath: "%s"', full_path));
                                        this.first_epoch = this.first_epoch.getEpoch(1 : this.last_epoch.length);
                                        this.is_valid_list(f) = false;
                                    elseif isempty(time)
                                        log.addError(sprintf('Check the following file, it seems to be empty\nPath: "%s"', full_path));
                                        this.first_epoch = this.first_epoch.getEpoch(1 : this.last_epoch.length);
                                        this.is_valid_list(f) = false;
                                    else
                                        this.last_epoch.append(time);
                                        log.addMessage(sprintf('        last  epoch found at: %s', this.last_epoch.last.toString()), this.verbosity_lev);
                                    end
                                end
                            end

                            if this.is_valid_list(f)
                                if ~isempty(coo)
                                    this.coo.append(Coordinates.fromStringXYZ(coo));
                                else
                                    this.coo.append(Coordinates.fromXYZ([0 0 0]));
                                end
                                this.eoh(f) = eoh;
                            end

                            fclose(fid);
                        end
                    end
                catch ex
                    log.addWarning(['"' this.file_name_list{f} this.ext{f} '" appears to be a corrupted RINEX file or missing'], this.verbosity_lev);
                    %Core_Utils.printEx(ex);
                    this.is_valid_list(f) = false;
                    fclose(fid);
                end
                % store marker_name
                if isempty(marker_name)
                    this.marker_name{f} = upper(this.file_name_list{f}(1:min(4, numel(this.file_name_list{f}))));
                else
                    this.marker_name{f} = marker_name;
                end
            end
            this.is_valid = all(this.is_valid_list);
            if (~this.is_valid)
                log.addWarning('Some or all the RINEX files are corrupted or missing!!!', this.verbosity_lev);
            end
        end

        function checkValidityUnbuffered(this, flag_header_only)
            % Update the status of validity of the files here pointed
            % Unbuffered read => read line by line (old approach)
            %
            % SYNTAX
            %   this.checkValidity();

            % pre allocate
            this.is_valid_list = false(1, numel(this.file_name_list));
            this.eoh = zeros(1, numel(this.file_name_list));
            this.first_epoch = GPS_Time();
            this.last_epoch = GPS_Time();
            this.coo = Coordinates;           % receiver coordinates
            % for each file present in the list
            for f = 1 : numel(this.file_name_list)
                marker_name = '';

                % try to find the first and the last epoch stored in the file
                try
                    %%
                    fid = fopen(fullfile(this.base_dir{f}, [this.file_name_list{f} this.ext{f}]));
                    if fid < 0
                        this.log.addError(['"' this.file_name_list{f} this.ext{f} '" appears to be missing'], this.verbosity_lev);
                        this.is_valid_list(f) = false;
                    else
                        l = 1;
                        line = fgetl(fid);
                        date_start = '';
                        date_stop = '';
                        coo = '';
                        par_to_find = 4;
                        while par_to_find > 0 && isempty(strfind(line,'END OF HEADER')) && ischar(line) %#ok<*STREMP>
                            l = l + 1;
                            line = fgetl(fid);
                            if numel(line) > 70 && isempty(marker_name) && line(64) == 'K' % read marker_name
                                marker_name = strtrim(regexp(line, '.*(?=MARKER NAME)', 'match', 'once'));
                                par_to_find = par_to_find - 1;
                            end
                            if numel(line) > 76
                                if line(61) == 'T'
                                    %tmp = regexp(line, '.*(?=GPS         TIME OF FIRST OBS)', 'match', 'once');
                                    if line(69) == 'F'
                                        date_start = strtrim(line(1:48));
                                        par_to_find = par_to_find - 1;
                                    elseif line(69) == 'L'
                                        %tmp = regexp(line, '.*(?=GPS         TIME OF LAST OBS)', 'match', 'once');
                                        date_stop = strtrim(line(1:48));
                                        par_to_find = par_to_find - 1;
                                    end
                                else
                                    %tmp = regexp(line, '.*(?=APPROX POSITION XYZ)', 'match', 'once');
                                    if line(62) == 'P' % APPROXIMATE POSITION
                                        coo = strtrim(line(1:60));
                                        par_to_find = par_to_find - 1;
                                    end
                                end
                            end
                        end

                        if ~isempty(coo)
                            this.coo.append(Coordinates.fromStringXYZ(coo));
                        else
                            this.coo.append(Coordinates.fromXYZ([0 0 0]));
                        end
                        this.eoh(f) = l;

                        % If I did not found date of start and date of end in the file header
                        % Try to read them
                        if ~isempty(date_start)
                            this.first_epoch.addEpoch(date_start, [], true);
                        else
                            epoch_line = fgetl(fid);
                            if epoch_line == -1
                                error('Empty file');
                            end
                            % try to guess the time format
                            [id_start, id_stop] = regexp(epoch_line, '[.0-9]*');
                            this.id_date = id_start(1) : id_stop(6); % save first and last char limits of the date in the line -> suppose it composed by 6 fields

                            this.first_epoch.addEpoch(epoch_line(this.id_date), [], true);
                        end
                        this.log.addStatusOk(['"' this.file_name_list{f} this.ext{f} '" appears to be a valid RINEX'], this.verbosity_lev);
                        this.log.addMessage(sprintf('        first epoch found at: %s', this.first_epoch.last.toString()), this.verbosity_lev);

                        if ~isempty(date_stop)
                            this.last_epoch.addEpoch(date_stop, [], true);
                            this.log.addMessage(sprintf('        last  epoch found at: %s', this.last_epoch.last.toString()), this.verbosity_lev);
                        else
                            if flag_header_only
                                this.log.addWarning('Last epoch not found in header, search in file is not enabled\nThe last epoch has not been saved within the  object', this.verbosity_lev);
                            else
                                % go to the end of the file to search for the last epoch
                                % to be sure to find at least one line containing a valid epoch, go to the end of the file minus 5000 characters
                                fseek(fid, -10000, 'eof');
                                txt = fread(fid, 10000,'*char')';

                                % Start searching for a valid epoch
                                time = [];
                                loop_n = 1;
                                while isempty(time) && ~isempty(txt)
                                    if ~isempty(find(txt(1:min(1000,numel(txt))) == 13, 1, 'first'))
                                        has_cr = true;  % The file has carriage return - I hate you Bill!
                                    else
                                        has_cr = false;  % The file is UNIX standard
                                    end

                                    % get new line separators
                                    nl = regexp(txt, '\n')';
                                    if nl(end) <  (numel(txt) - double(has_cr))
                                        nl = [nl; numel(txt)];
                                    end
                                    lim = [[1; nl(1 : end - 1) + 1] (nl - 1 - double(has_cr))];
                                    lim = [lim lim(:,2) - lim(:,1)];
                                    lim(lim(:,3) < 20, :) = [];

                                    l = size(lim, 1);
                                    while (l >= 1) && isempty(time)
                                        line = txt(lim(l,1) : lim(l,2));
                                        if ~isempty(regexp(line(1:15),'( [0-9]{1,4} [ 0-9]{1}[0-9]{1} [ 0-9]{1}[0-9]{1})', 'once'))
                                            % Check that the read epoch is within 30 days from the first epoch
                                            % (if it's further in time it's probably a misleading false epoch line)
                                            time = GPS_Time(line(this.id_date), [], true);
                                        end
                                        l = l - 1;
                                    end
                                    if isempty(time) || time.isEmpty()
                                        loop_n = loop_n + 1;
                                        fseek(fid, loop_n * -10000, 'eof'); % If no valid time have been found try to go back more...
                                        txt = fread(fid, 10000,'*char')';
                                    end
                                end
                                this.last_epoch.append(time);
                                this.log.addMessage(sprintf('        last  epoch found at: %s', this.last_epoch.last.toString()), this.verbosity_lev);
                            end
                        end
                        fclose(fid);
                        this.is_valid_list(f) = true;
                    end
                catch ex
                    this.log.addWarning(['"' this.file_name_list{f} this.ext{f} '" appears to be a corrupted RINEX file'], this.verbosity_lev);
                    this.is_valid_list(f) = false;
                end

                % store marker_name
                if isempty(marker_name)
                    this.marker_name{f} = upper(this.file_name_list{f}(1:min(4, numel(this.file_name_list{f}))));
                else
                    this.marker_name{f} = marker_name;
                end
            end
            this.is_valid = all(this.is_valid_list);
            if (~this.is_valid)
                this.log.addWarning('Some or all the RINEX files are corrupted or missing!!!', this.verbosity_lev);
            end
        end

        function checkCoordinates(rin_list, force_delete)
            flag_show = false;
            dist_lim = 1; % [m]
            for r = 1 : numel(rin_list)
                xyz = rin_list(r).coo.xyz;
                id_ko = (all(xyz == 0,2));
                xyz(id_ko,:) = [];
                if not(isempty(xyz))
                    xyz = bsxfun(@minus, xyz, median(xyz, 1, 'omitnan'));
                    dist = 1e-3 * (sqrt(xyz(:,1).^2 + xyz(:,2).^2 + xyz(:,3).^2)); % [km]
                    if any(dist > dist_lim)
                        % This set of file RINEX are from a receiver that have been moved
                        % for more than 100 meters

                        marker_name = rin_list(r).marker_name{find(rin_list(r).is_valid_list, 1, 'first')};
                        if flag_show
                            rin_list(r).coo.showPositionENU
                            title(marker_name);
                        end

                        log = Core.getLogger;
                        log.addWarning(sprintf('The RINEX of the station %s are anomalous. The antenna could have been moved or the files are from different omonimous receivers\n > MAX detected distance: %.3g [km]\n\nThe following files will be ignored:', marker_name, max(dist)));
                        id_coo = find(not(id_ko));
                        id_valid = find(rin_list(r).is_valid_list);
                        id_ko = find(dist > dist_lim)';
                        for i = id_ko
                            file_name = fullfile(rin_list(r).base_dir{id_valid(id_coo(i))}, [rin_list(r).file_name_list{id_valid(id_coo(i))}, rin_list(r).ext{id_valid(id_coo(i))}]);
                            if nargin == 2 && force_delete
                                delete(file_name);
                            end
                            log.addMonoMessage(log.indent(sprintf('- %s', file_name)));
                        end
                        rin_list(r).first_epoch.remEpoch(id_coo(id_ko));
                        rin_list(r).last_epoch.remEpoch(id_coo(id_ko));
                        rin_list(r).is_valid_list(id_valid(id_coo(id_ko))) = false;
                        rin_list(r).coo.rem(id_coo(id_ko));
                    end
                end
            end
        end

        function file_name = getFileName(this, file_number)
            % Get the full path of a file (if the object contains a list of files, the id can be specified)
            % SYNTAX: file_name = this.getFileName(<file_number = 1>)
            if nargin == 1
                file_number = 1;
            end
            file_name = fullfile(this.base_dir{file_number}, [this.file_name_list{file_number} this.ext{file_number}]);
        end

        function first_epoch = getFirstEpoch(this, session)
            % Get the first epoch of a session
            %
            % SYNTAX
            %   first_epoch = this.getFirstEpoch(session)
            if nargin == 1
                if any(this.is_valid_list)
                    first_epoch = this.first_epoch.getCopy;
                else
                    first_epoch = GPS_Time();
                end
            else
                if Core.getState.isRinexSession
                    if this.is_valid_list(session)
                        id_ss = cumsum(this.is_valid_list);
                        first_epoch = this.first_epoch.getEpoch(id_ss(session));
                    else
                        first_epoch = GPS_Time();
                    end
                else
                    first_epoch = this.state.getSessionLimits(session).first;
                end
            end
        end

        function last_epoch = getLastEpoch(this, session)
            % Get the last epoch of a session
            %
            % SYNTAX
            %   first_epoch = this.getFirstEpoch(session)
            if nargin == 1
                if any(this.is_valid_list)
                    last_epoch = this.last_epoch.getCopy;
                else
                    last_epoch = GPS_Time();
                end
            else
                if Core.getState.isRinexSession
                    if this.is_valid_list(session)
                        id_ss = cumsum(this.is_valid_list);
                        last_epoch = this.last_epoch.getEpoch(id_ss(session));
                    else
                        last_epoch = GPS_Time();
                    end
                else
                    last_epoch = this.state.getSessionLimits(session).last;
                end
            end
        end

        function line_num = getEOH(this, file_number)
            % Get the end of header line of a RINEX file (if the object contains a list of files, the id can be specified)
            % SYNTAX: line_num = this.getEOH(<file_number = 1>)
            if nargin == 1
                file_number = 1;
            end
            line_num = this.eoh(file_number);
        end

        function validity = isValid(this, file_number)
            % Get the validity of a RINEX file or the object (if the object contains a list of files, the id can be specified)
            % SYNTAX: validity = isValid(<file_number>)
            validity = false;
            for r = 1 : numel(this)

                if (nargin == 1)
                    validity = validity || any([this(r).is_valid_list]);
                else
                    if file_number <= numel(this(r).is_valid_list)
                        validity = validity || this(r).is_valid_list(file_number);
                    end
                end
            end
        end

        function keepFiles(this, date_start, date_stop)
            % keep only files contained in the two dates
            %
            % SYNTAX
            % this.keepFiles(date_start, date_stop)
            for r = 1 : numel(this)
                [in_bound] = Core_Utils.timeIntersect(this(r).first_epoch, this(r).last_epoch, date_start, date_stop);
                to_keep = this(r).is_valid_list;
                to_keep(to_keep) = in_bound;
                this(r).keep(to_keep);
            end
        end

        function has_obs = hasObsInSession(this, date_start, date_stop)
            % has observation in the two dates
            %
            % SYNTAX
            % this.keepFiles(date_start, date_stop)
            has_obs = false(numel(this),1);
            for r = 1 : numel(this)
                if this(r).isValid
                    in_bound = Core_Utils.timeIntersect(this(r).first_epoch, this(r).last_epoch, date_start, date_stop);
                    has_obs(r) = any(in_bound);
                end
            end
        end

        function keep(this,id)
            % keep only element specified by id
            %
            % SYNTAX
            % this.keep(id)
            this.base_dir = this.base_dir(id);
            this.file_name_list = this.file_name_list(id);
            this.ext = this.ext(id);
            this.eoh = this.eoh(id);
            this.first_epoch = this.first_epoch.getEpoch(id(this.is_valid_list));
            this.last_epoch  = this.last_epoch.getEpoch(id(this.is_valid_list));
            this.is_valid_list = this.is_valid_list(id);
        end

        function printMissingFiles(rin_list, flag_download)
            % Show all the files that seems to be missing for the
            % processing
            if nargin == 1
                flag_download = false;
            end
            clc
            fprintf('These files are missing or currupted\n')
            for r = 1: numel(rin_list)
                id_ko = find(not(rin_list(r).is_valid_list));
                for i = id_ko
                    file_name = [rin_list(r).file_name_list{i}, rin_list(r).ext{i}];
                    if not(flag_download)
                        fprintf('%s\n', file_name);
                    end
                    if numel(file_name) == 12 % it is a RINEX2 file name
                        % Extracting year and doy
                        year = 2000 + str2double(file_name(end-2:end-1));
                        doy = str2double(file_name(5:7));
                        if flag_download
                            str_cmd = sprintf('/usr/bin/python3 /home/gred/Repositories/goget/goGet.py download -m %s --year-doy %04d %03d -o %s --no-wait', rin_list(r).marker_name{i}, year, doy, [rin_list(r).base_dir{i} filesep]);
                            fprintf('%s\n', str_cmd);
                        end
                    end
                end

            end

        end

    end
end

author: goGPS-Project


text:
function [opts, imdb] = os_setup(varargin)

global is_iasonas
if (exist('is_iasonas', 'var') ||isempty(is_iasonas) || ~is_iasonas)
    setup ;
end

opts.transf = {};
opts.seed = 1 ;
opts.batchSize = 128 ;
opts.useGpu = true ;
opts.regionBorder = 0.05 ;
opts.numDCNNWords = 64 ;
opts.numDSIFTWords = 256 ;
opts.numSamplesPerWord = 1000 ;
opts.printDatasetInfo = false ;
opts.excludeDifficult = true ;
opts.encoders = {struct('type', 'rcnn', 'opts', {})} ;
opts.dataset = 'os' ;
opts.osDir = 'data/os' ;
opts.fmdDir = 'data/fmd' ;
opts.dtdDir = 'data/dtd';
opts.kthDir = 'data/kth';
opts.alotDir = 'data/alot';
opts.mitDir = 'data/mit_indoor';
opts.msrcDir = 'data/msrc_c';
opts.cubDir = 'data/cub';
opts.vocDir = 'data/VOC2007';
opts.vocDir12 = 'data/VOC2012';
opts.writeResults = false;
opts.compid = 'comp2';
opts.suffix = 'baseline' ;
opts.prefix = 'v22' ;
opts.model = 'imagenet-vgg-m.mat';
opts.layer = 13 ; % for D-CNN (R-CNN is the penultimate)
opts.segProposalType = 'crisp' ;
opts.gpuId = 1;
opts.crf = [];
[opts, varargin] = vl_argparse(opts,varargin) ;

opts.expDir = sprintf('data/%s/%s-seed-%02d', opts.prefix, opts.dataset, opts.seed) ;
opts.imdbDir = fullfile(opts.expDir, 'imdb') ;
if ~iscell(opts.suffix)
  opts.resultPath = fullfile(opts.expDir, sprintf('result-%s.mat', opts.suffix)) ;
  opts.segResultPath = fullfile(opts.expDir, sprintf('%s/result-%s-seg.mat', opts.segProposalType, opts.suffix)) ;
  opts.segPublishDir = fullfile(opts.expDir, sprintf('%s/result-%s-seg-figures', opts.segProposalType, opts.suffix)) ;
  vl_xmkdir(opts.segPublishDir);
  if opts.writeResults
      opts.vocResultPath = fullfile(opts.expDir, [sprintf('result-%s/Main/%s_cls_', opts.suffix, opts.compid) '%s_%s.txt']);
  end
else
  for s = 1:numel(opts.suffix)
    opts.resultPath{s} = fullfile(opts.expDir, sprintf('result-%s.mat', opts.suffix{s})) ;
    opts.segResultPath{s} = fullfile(opts.expDir, sprintf('%s/result-%s-seg.mat', opts.segProposalType, opts.suffix{s})) ;
    opts.segPublishDir{s} = fullfile(opts.expDir, sprintf('%s/result-%s-seg-figures', opts.segProposalType, opts.suffix{s})) ;
    vl_xmkdir(opts.segPublishDir{s});
  end
end
opts = vl_argparse(opts,varargin) ;

if nargout <= 1, return ; end

% Setup GPU if needed
if opts.useGpu
  gpuDevice(opts.gpuId) ;
end

% -------------------------------------------------------------------------
%                                                            Setup encoders
% -------------------------------------------------------------------------

models = {} ;
for i = 1:numel(opts.encoders)
  if isstruct(opts.encoders{i})
    name = opts.encoders{i}.name ;
    opts.encoders{i}.path = fullfile(opts.expDir, [name '-encoder.mat']) ;
    opts.encoders{i}.codePath = fullfile(opts.expDir, [name '-codes.mat']) ;
    models = horzcat(models, get_cnn_model_from_encoder_opts(opts.encoders{i})) ;
  else
    for j = 1:numel(opts.encoders{i})
      name = opts.encoders{i}{j}.name ;
      opts.encoders{i}{j}.path = fullfile(opts.expDir, [name '-encoder.mat']) ;
      opts.encoders{i}{j}.codePath = fullfile(opts.expDir, [name '-codes.mat']) ;
      models = horzcat(models, get_cnn_model_from_encoder_opts(opts.encoders{i}{j})) ;
    end
  end
end

% -------------------------------------------------------------------------
%                                                       Download CNN models
% -------------------------------------------------------------------------

for i = 1:numel(models)
    crt_model_path = fullfile('data/models', models{i});
    vl_xmkdir('data/models');
    if ~exist(crt_model_path, 'file')
        fprintf('downloading model %s\n', models{i});
        urlwrite(fullfile('http://www.vlfeat.org/matconvnet/models', ...
            models{i}), crt_model_path);
        % backwards compatibility for the model:
        fix_load_nn(crt_model_path);
    end
end

% -------------------------------------------------------------------------
%                                                              Load dataset
% -------------------------------------------------------------------------

vl_xmkdir(opts.expDir) ;
vl_xmkdir(opts.imdbDir) ;

imdbPath = fullfile(opts.imdbDir, sprintf('imdb-seed-%d.mat', opts.seed)) ;
if exist(imdbPath)
  imdb = load(imdbPath) ;
  return ;
end

switch opts.dataset
  case 'os'
    imdb = os_get_database(opts.osDir) ;
  case 'os-a'
    imdb = os_attr_get_database(opts.osDir);
  case 'fmd'
    imdb = fmd_get_database(opts.fmdDir, 'seed', opts.seed) ;
  case 'dtd'
    imdb = dtd_get_database(opts.dtdDir, 'seed', opts.seed);
  case 'kth'
    imdb = kth_get_database(opts.kthDir, 'seed', opts.seed);
  case 'voc07'
    imdb = voc_get_database(opts.vocDir, 'seed', opts.seed);
  case 'voc12'
    imdb = voc_get_database(opts.vocDir12, 'seed', opts.seed);
  case 'voc12s'
    imdb = voc_get_seg_database(opts.vocDir12, 'seed', opts.seed);
  case 'mit'
    imdb = mit_indoor_get_database(opts.mitDir);
  case 'msrc'
    imdb = msrc_c_get_database(opts.msrcDir);
  case 'cubcrop'
    imdb = cub_get_database(opts.cubDir, true);
  case 'cub'
    imdb = cub_get_database(opts.cubDir, false);
  case 'alot'
    imdb = alot_get_database(opts.alotDir, 'seed', opts.seed);
  otherwise
    serror('Unknown dataset %s', opts.dataset) ;
end

save(imdbPath, '-struct', 'imdb') ;

if opts.printDatasetInfo
  print_dataset_info(imdb) ;
end

% -------------------------------------------------------------------------
function model = get_cnn_model_from_encoder_opts(encoder)
% -------------------------------------------------------------------------
p = find(strcmp('model', encoder.opts)) ;
if ~isempty(p)
  [~,m,e] = fileparts(encoder.opts{p+1}) ;
  model = {[m e]} ;
else
  model = {} ;
end

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
function dispScanStats(scanBinsV, volHistV, name, nameVol, planC, indexS, opt)
%Command line display of basic scan statistics
%scanBinsV is a vector of the midpoint scanBin values.
%volHistV is either a histogram of volumes or surface areas.
%LM: 14 Oct 02, JOD.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

% ESPEZI OCT 2014 added nameVol, changed order of printed items and added dose name

ud = get(findobj('Tag', 'IVHGui'),'userdata');
scanNum = get(ud.af.handles.scan,'value');
imageType = planC{indexS.scan}(scanNum).scanInfo(1).imageType;

if strcmpi(imageType,'CT')
    units = 'HU';
elseif strcmpi(imageType,'PET')
    units = 'SUV';
else
    units = '';
end

switch lower(opt)

  case 'standardscan'
    disp('-----------------------')
    disp('')
    disp(['Structure is:  ' name])

    totalVol = sum(volHistV);
    disp(['Total volume is:  ' num2str(totalVol) ' cubic cm.'])

    disp(['Scan name is:  ' nameVol])

    meanD = sum(scanBinsV.*volHistV)/sum(volHistV);
    disp(['Mean ' imageType ' ' units ' is:  ' num2str(meanD)])

    ind = max(find([volHistV~=0]));
    maxD = scanBinsV(ind);
    disp(['Max' imageType ' ' units ' is:  ' num2str(maxD)])

    ind = min(find([volHistV~=0]));
    minD = scanBinsV(ind);
    disp(['Min ' imageType ' ' units ' is:  ' num2str(minD)])
    disp('')
    disp('-----------------------')

  case 'dshscan'

    disp('-----------------------')
    disp('')
    disp(['Structure is:  ' name])

    areaV = volHistV;  %actually areas, not volumes.
    scansV = scanBinsV;

    totalArea = sum(areaV);
    disp(['Total surface area is:  ' num2str(totalArea) ' square cm.'])

    disp(['Scan name is:  ' nameVol])

    meanD = sum(scansV.*areaV)/sum(areaV);
    disp(['Mean surface ' imageType ' ' units ' is:  ' num2str(meanD)])

    maxScan = max(scansV);
    disp(['Max ' imageType ' ' units ' is:  ' num2str(maxScan)])

    minScan = min(scansV);
    disp(['Min ' imageType ' ' units ' is:  ' num2str(minScan)])
    disp('')
    disp('-----------------------')

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% DATES
% Converts between the following date formats:
%  - Year, Month, Day
%  - Year, Day-of-Year
%  - Year, Decimal Year
%  - GPS week, seconds of week
%
% Required M-files: cal2jd, doy2jd, gps2jd, jd2cal, jd2dow, jd2doy, jd2gpsa

% Created 24 May 96, M. Craymer
% Modified 03 Aug 96, Converted to MATLAB
% Modified 29 Sep 98, Convert each option to all other formats
% Modified 25 Mar 99, Converted to use isempty
% Modified 26 Apr 99, Rewrote to use new date conversion routines
%                     Added conversion of GPS week and day of week
% Modified 19 Feb 11, Replaced clear all with clear for specific variables used

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

while 1
clear ver day month buf fmt date ok iyr mn dy jd doy yr dow gpsweek sow rollover
warning on
ver = '1.0 (99.04.26)';
day = ['Sun';'Mon';'Tue';'Wed';'Thu';'Fri';'Sat'];
month = ['Jan';'Feb';'Mar';'Apr';'May';'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec'];

fprintf(1,'\n\n');
fprintf(1,'-------------------------------------------------\n');
fprintf(1,' DATES: Converts between different date formats.\n');
fprintf(1,'       M.Craymer, v%s\n',ver);
fprintf(1,'-------------------------------------------------\n');

fprintf(1,'\n');
fprintf(1,'Available input date formats:\n');
fprintf(1,' 1) Year, Month, Day\n');
fprintf(1,' 2) Year, Day of Year\n');
fprintf(1,' 3) Year (including decimal year)\n');
fprintf(1,' 4) GPS Week, Sec of Week, Number Rollovers\n');
fprintf(1,' 5) Julian Date\n');
fprintf(1,'\n');

buf = input('Enter date format and date [quit] > ','s');
if isempty(buf)
  fmt = 0;
else
  date = sscanf(buf,'%f');
  fmt = date(1);
end
ok = 1;

%----- Exit program
if fmt == 0
  clear ver day month buf fmt date ok iyr mn dy jd doy yr dow gpsweek sow rollover
  return;

%----- Day, month, year
elseif fmt == 1
  if length(date) < 4
    warning('Too few date arguments entered');
    ok = 0;
  else
    iyr = date(2);
    mn = date(3);
    dy = date(4);
    if mn<1 | mn>12
      warning('Invalid month');
      ok = 0;
    end
    if dy<1 | dy>=32
      warning('Invalid day');
      ok = 0;
    end
    if ok
      jd = cal2jd(iyr,mn,dy);
      doy = jd2doy(jd);
      yr = jd2yr(jd);
      dow = jd2dow(jd);
      [gpsweek,sow,rollover] = jd2gps(jd);
    end
  end

%----- Year, day of year
elseif fmt == 2
  if length(date) < 3
    warning('Too few date arguments entered')
    ok = 0;
  else
    iyr = date(2);
    doy = date(3);
    if doy<1 | doy>=367
      warning('Invliad day of year');
      ok = 0;
    end
    if ok
      jd=doy2jd(iyr,doy);
      [iyr,mn,dy]=jd2cal(jd);
      yr = jd2yr(jd);
      dow = jd2dow(jd);
      [gpsweek,sow,rollover] = jd2gps(jd);
    end
  end

%----- Year & decimal of year
elseif fmt == 3
  if length(date) < 2
    warning('No date arguments entered');
    ok = 0;
  else
    yr = date(2);
    jd=yr2jd(yr);
    [iyr,mn,dy]=jd2cal(jd);
    doy = jd2doy(jd);
    dow = jd2dow(jd);
    [gpsweek,sow,rollover] = jd2gps(jd);
  end

%----- GPS week, sec of week
elseif fmt == 4
  if length(date) < 4
    rollover = 0;
  else
    rollover = date(4);
  end
  if length(date) < 3
    sow = 0;
  else
    sow = date(3);
  end
  if length(date) < 2
    warning('No date arguments entered');
    ok = 0;
  else
    gpsweek = date(2);
    if gpsweek<1 | gpsweek>=1025
      warning('Invalid GPS week');
      ok = 0;
    end
    if sow<0 | sow>604800
      warning('Invliad sec of week');
      ok = 0;
    end
    if ok
      jd = gps2jd(gpsweek,sow,rollover);
      [iyr,mn,dy]=jd2cal(jd);
      doy = jd2doy(jd);
      yr = jd2yr(jd);
      dow = jd2dow(jd);
    end
  end

%----- Julian date
elseif fmt == 5
  if length(date) < 2
    warning('No date arguments entered');
    ok = 0;
  else
    jd = date(2);
    if jd<1
      warning('Invalid Julian date');
      ok = 0;
    end
    if ok
      [iyr,mn,dy]=jd2cal(jd);
      doy = jd2doy(jd);
      yr = jd2yr(jd);
      dow = jd2dow(jd);
      [gpsweek,sow,rollover] = jd2gps(jd);
    end
  end

%----- Invalid date format
else
  warning('Undefined date format entered');
  ok = 0;
end

%----- List results
if ok
  fprintf(1,'\n');
  fprintf(1,'Date    %4d-%2d-%4.1f\n',[iyr,mn,dy]);
  fprintf(1,'Year    %12.4f\n',yr);
  fprintf(1,'Day of Year   %6.1f\n',doy);
  fprintf(1,'Day of Week      %s\n',day(dow,:));
  fprintf(1,'GPS Week      %6d\n',gpsweek);
  fprintf(1,'Sec of Week   %6.0f\n',sow);
  fprintf(1,'Julian Date   %9.1f\n',jd);
end

end %while

author: Sable

repo: mcbench-benchmarks


text:
function [result]=sv_NodeCalc(mCatalog, fSplitTime, bTimePeriod, fTimePeriod, nCalculateMC, vCluster, bDecluster, fMcOrg, fBValueOrg, bSeisModel)
% function [result]=sv_NodeCalc(mCatalog, fSplitTime, bTimeperiod, fTimePeriod, nCalculateMC, vCluster, bDecluster, fMcOrg, fBValueOrg, bSeisModel)
% -------------------------------------------------------------------------------------------------------------------------------------------------
% Function to calculate absolute difference of number of events between
% to time periods of an earthquake catalog normalized to a year
% Incoming variables:
% mCatalog     : current earthquake catalog
% fSplitTime   : Splitting time
% bTimePeriod  : Use specific time period to compare (1) or from beginning till end of catalog (0)
% fTimePeriod  : Time period in days to be compared before and after fSplitTime
% nCalculateMC : Number for Method of Mc determination
% vCluster     : Vector of cluster numbers
% bDecluster   : Decluster catalog or not
% fMcOrg       : Mc of entire catalog
% fBValueOrg   : b-value of entire catalog
% bSeisModel   : Method to model seismicity varation: 0 - b-value fittng, 1 - grid search
%
% Outgoing variable:
% result.dNdiffsum      : total difference of number of events in the two time periods
% result.fMc            : magnitude of completeness for entire catalog
% result.fMcFirstPeriod : magnitude of completeness for first period defined
% result.fMcSecondPeriod : magnitude of completeness for second period defined
% result.fdMc            : Difference in Mc: Mc(Per1)-Mc(Per2)
% result.fMshift         : Simple magnitude shift
% result.fMshiftFit      : Goodness of fit in percent to modeling second period with simple magnitude shift of first period
% result.fMshiftSig      : Simple magnitude shift at 99% significance level of z-statistic
% result.fFactorHi       : Mean rate factor for EQ >= Mc(Per1) between two periods
% result.fPerHi          : Percental seismicity change
% result.fFactorLow      : Mean rate factor for EQ < Mc(Per1) between two periods
% result.fPerLow         : Percental seismicity change
% result.fMRateFit       : Goodness of fit by modeling second period with rate factor multiplication of first period
% result.fMshiftCom      : Magnitude shift combined with stretch
% result.fStretch        : Magnitude stretch
% result.fMTransFit      : Goodness of fit by modeling second period manipulating the first period
% result.fdMag           : Magnitude shift for a M=1 EQ using result.fMshiftCom and result.fStretch
% result.mModelCat       : Magnitude alternation for a magnitude M=1 EQ
%
% Author: J. Woessner
% woessner@seismo.ifg.ethz.ch
% last update: 01.11.02
% Changes:
% 19.08.02: Replaced  fcumulsum.m with calc_cumulsum.m
% 08.10.02: Deleted return at the end.

% Changes:
% 19.08.02: Replaced  fcumulsum.m with calc_cumulsum.m
% 08.10.02: Deleted return at the end.

% Init variable
result=[];
fTimePeriod =fTimePeriod/365;

if bDecluster == 1
    % Determine degree of Poissonian distribution using Chi^2-Test
    [result.fChi2 result.fChi2_sig90 result.fChi2_sig95 result.Chi2_sig99 result.nPoissDeg]= calc_Chi2(mCatalog);
    vSel = (vCluster(:,1) > 0);
    mCatalogDecl = mCatalog(~vSel,:);
    [result.fChi2_Dec result.fChi2_sig90_Dec result.fChi2_sig95_Dec result.Chi2_sig99_Dec...
            result.nPoissDeg_Dec]= calc_Chi2(mCatalogDecl);
    % Determine degree of Clustering
    [result.fClusterDeg] = calc_ClusterDeg(mCatalog, vCluster);
elseif (bSeisModel == 1)
    % Create the catalogs for the two time peiods
    [result.mFirstCatalog, result.mSecondCatalog, result.fFirstPeriodExact, result.fSecondPeriodExact, result.fFirstPeriod,...
            result.fSecondPeriod] = ex_SplitCatalog(mCatalog, fSplitTime, bTimePeriod, fTimePeriod, bTimePeriod, fTimePeriod);

    %% Determination of difference between two time periods normalized to year
    [dEv_val dMags dEv_valsum dEv_valsum_rev,  dMags_rev] = calc_cumulsum(result.mFirstCatalog);
    [dEv_val2 dMags2 dEv_valsum2 dEv_valsum_rev2,  dMags_rev2] = calc_cumulsum(result.mSecondCatalog);
    result.dNdiff = max(dEv_val2)/result.fSecondPeriodExact-max(dEv_val)/result.fFirstPeriodExact; % Normalization
    result.dNdiffsum = cumsum(result.dNdiff');
    %result.dNdiffsumVal = result.dNdiffsum(length(result.dNdiffsum));

    % Determine Mc for entire catalog, two time periods and difference in Mc
    result.fMc = calc_Mc(mCatalog, nCalculateMC);
    if isempty(result.fMc)
        result.fMc = NaN;
    end
    result.fMcFirstPeriod = calc_Mc(result.mFirstCatalog, nCalculateMC);
    if isempty(result.fMcFirstPeriod)
        result.fMcFirstPeriod = NaN;
    end
    result.fMcSecondPeriod = calc_Mc(result.mSecondCatalog, nCalculateMC);
    if isempty(result.fMcSecondPeriod)
        result.fMcSecondPeriod = NaN;
    end
    result.fdMc = result.fMcSecondPeriod - result.fMcFirstPeriod;

    % Determine simple magnitude shift and Goodness of Fit
    [result.fMshift result.fMshiftFit] = calc_Magshift(mCatalog, fSplitTime, bTimePeriod, fTimePeriod, nCalculateMC);
    if isempty(result.fMshift)
        result.fMshift = NaN;
        result.fMshiftFit = NaN;
    end
    % Determine significant Mshift at 99% level
    [mLMagsig mHMagsig fLZmax fLZmean fLZmin fHZmax fHZmean,  fHZmin] =...
        calc_Magsig(result.mFirstCatalog, result.mSecondCatalog , result.fFirstPeriodExact, result.fSecondPeriodExact, 0.1);
    if ((fLZmax >= 2.57 & fHZmin <= -2.57) | fLZmin <= -2.57 & fHZmax >= 2.57)
        result.fMshiftSig = result.fMshift;
    else
        result.fMshiftSig = NaN;
    end
    % Determine rate factors and percentage change of EQ for M<Mc and M>=Mc
    [result.fFactorHi fStdHi fResHi result.fPerHi result.fFactorLow fStdLow fResLow result.fPerLow result.fMRateFit]...
        = calc_ratefac(result.mFirstCatalog, result.mSecondCatalog, fTimePeriod, fTimePeriod,...
        result.fMcFirstPeriod, result.fMcSecondPeriod);
    [result.fMshiftCom result.fStretch result.fMTransFit result.fARate] = calc_Shiftstretch(mCatalog, fSplitTime, bTimePeriod,...
        fTimePeriod, nCalculateMC, fMcOrg, fBValueOrg);
    %% Magnitude fit for a magnitude M=1 EQ
    result.fdMag = 1+((result.fStretch-1)+result.fMshiftCom); % Theoretical map
    result.mModelCat = result.mFirstCatalog;
    result.mModelCat(:,6) = result.fStretch.*result.mModelCat(:,6)+result.fMshiftCom;
else
    % Create the catalogs for the two time peiods
    [result.mFirstCatalog, result.mSecondCatalog, result.fFirstPeriodExact, result.fSecondPeriodExact, result.fFirstPeriod,...
            result.fSecondPeriod] = ex_SplitCatalog(mCatalog, fSplitTime, bTimePeriod, fTimePeriod, bTimePeriod, fTimePeriod);

    %% Determination of difference between two time periods normalized to year
    [dEv_val dMags dEv_valsum dEv_valsum_rev,  dMags_rev] = calc_cumulsum(result.mFirstCatalog);
    [dEv_val2 dMags2 dEv_valsum2 dEv_valsum_rev2,  dMags_rev2] = calc_cumulsum(result.mSecondCatalog);
%     result.dNdiff = max(dEv_val2)/result.fSecondPeriodExact-max(dEv_val)/result.fFirstPeriodExact; % Normalization
    result.dNdiff = max(dEv_valsum2)/result.fSecondPeriodExact-max(dEv_valsum)/result.fFirstPeriodExact;
    result.dNdiffsum = max(cumsum(result.dNdiff'));
    %result.dNdiffsumVal = result.dNdiffsum(length(result.dNdiffsum));
    if isempty(result.dNdiff)
        result.dNdiff = NaN;
    end
    if isempty(result.dNdiffsum)
        result.dNdiffsum = NaN;
    end
    % Determine Mc for entire catalog, two time periods and difference in Mc
    result.fMc = calc_Mc(mCatalog, nCalculateMC);
    if isempty(result.fMc)
        result.fMc = NaN;
    end
    result.fMcFirstPeriod = calc_Mc(result.mFirstCatalog, nCalculateMC);
    if isempty(result.fMcFirstPeriod)
        result.fMcFirstPeriod = NaN;
    end
    result.fMcSecondPeriod = calc_Mc(result.mSecondCatalog, nCalculateMC);
    if isempty(result.fMcSecondPeriod)
        result.fMcSecondPeriod = NaN;
    end
    result.fdMc = result.fMcSecondPeriod - result.fMcFirstPeriod;

    % Mc determination by grid search
    %[result.fProbMcGrid result.fMcGrid] = calc_McGridN(mCatalog,0.1);
    [result.fProbMcGrid result.fMcGrid] = calc_McCdf2(mCatalog,0.1)
    if (isempty(result.fProbMcGrid) | isempty(result.fMcGrid))
        result.fProbMcGrid = NaN;
        result.fMcGrid = NaN;
    end
    %% Use data only above completeness
    %fMinMc = min([result.fMcFirstPeriod result.fMcSecondPeriod]);
    %     vSel1 = (result.mFirstCatalog(:,6) >= result.fMcFirstPeriod);
    %     vSel2 = (result.mSecondCatalog(:,6) >= result.fMcFirstPeriod);
    %     result.mFirstCatalog = result.mFirstCatalog(vSel1,:);
    %     result.mSecondCatalog = result.mSecondCatalog(vSel2,:);

    %% Seismicity variation modelling
    %!!!!! NO MODELS WITH STRETCH ALLOWED !!!!!!!!!!!!!!!!!!!!!!!!!
    %%% Maximum likelihood scores and BICs
    [fProb_nochange, fBic_nochange] = calc_loglikelihood_nochange2(result.mFirstCatalog, result.mSecondCatalog);
    [fdM, fProb_dM, fBic_dM, mLikeli_dM] = calc_loglikelihood_dM2(result.mFirstCatalog, result.mSecondCatalog);
    %[fS, fProb_stretch, fBic_stretch, mLikeli_dS] = calc_loglikelihood_stretch(result.mFirstCatalog, result.mSecondCatalog);
    [fFac, fProb_Rate, fBic_Rate, mLikeli_Rate] = calc_loglikelihood_rate2(result.mFirstCatalog, result.mSecondCatalog);
    [fdM_rate, fdM_Fac, fProb_dMrate, fBic_dMrate, mLikeli_dMrate] = calc_loglikelihood_dM_rate2(result.mFirstCatalog, result.mSecondCatalog);
    %[fdS_rate, fdS_Fac, fProb_dSrate, fBic_dSrate mLikeli_dSrate] = calc_loglikelihood_stretch_rate(result.mFirstCatalog, result.mSecondCatalog);
    %[fdM_st, fStretch, fProb_Trans, fBic_Trans, mLikeli_Trans] = calc_loglikelihood_Trans(result.mFirstCatalog, result.mSecondCatalog);
    %[fdM_all, fdS_all, fFac_all, fProb_all, fBic_all, mLikeli_all] = calc_loglikelihood_dMdSrate(result.mFirstCatalog, result.mSecondCatalog);

    %vBic = [fBic_nochange; fBic_dM; fBic_stretch; fBic_Rate; fBic_dMrate; fBic_Trans; fBic_dSrate; fBic_all];
    vBic = [fBic_nochange; fBic_dM;  fBic_Rate; fBic_dMrate];
    [result.nType] = find(vBic == min(vBic));
    %% Initialize result
    result.fdM = nan; % Shift
    result.fdS = nan; % Stretch
    result.fRf = nan; % rate factor
    if length(result.nType) > 1
        vBic;
        result.nModelChoice = NaN;
    elseif length(result.nType) == 0
        result.nModelChoice = nan;
    else
        switch result.nType
        case 1
            result.nModelChoice = 1;
            %     result.fdM = 0; % Shift
            %     result.fdS = 1; % Stretch
            %     result.fRf = 1; % rate factor
        case 2
            result.nModelChoice = 2;
            result.fdM = fdM;
%         case 3
%             result.nModelChoice = 3;
%             result.fdS = fS;
        case 3
            result.nModelChoice = 3;
            result.fRf = fFac;
        case 4
            result.nModelChoice = 4;
            result.fdM = fdM_rate;
            result.fRf = fdM_Fac;
%         case 6
%             result.nModelChoice = 6;
%             result.fdM = fdM_st;
%             result.fdS = fStretch;
%         case 7
%             result.nModelChoice = 7;
%             result.fdS = fdS_rate;
%             result.fRf = fdS_Fac;
%         case 8
%             result.nModelChoice = 8;
%             result.fdM = fdM_all;
%             result.fdS = fdS_all;
%             result.fRf = fFac_all;
        otherwise
            disp('Something is equal');
        end; % END of Switch result.Type
    end;% END of IF result.Type
end; % END of IF bSeismodel
end % End of IF bDecluster

author: CelsoReyes


text:
function [net, state, stats, db] = train_network(net, getBatch, opts)
% TRAIN_NETWORK  Adapted version of CNN_TRAIN_DAG from MatConvNet package.
% It performs the training of a CNN using the DagNN wrapper.

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

opts.plotStatistics = true;
opts.prefetch = false ;
opts.postEpochFn = [] ;  % postEpochFn(net,params,state) called after each epoch; can return a new learning rate, 0 to stop, [] for no change
opts.extractStatsFn = @extractStats ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end

% -------------------------------------------------------------------------
%                                                         		Prepare data
% -------------------------------------------------------------------------
%
% Load db.mat for the training data and initialize negatives
if exist(opts.dbPath, 'file')
  fprintf('>> Loading DB file...\n'); t=tic;
  db = load(opts.dbPath) ;
  db.train.nidxs = zeros(opts.numNegative, numel(db.train.qidxs));
  db.val.nidxs = zeros(opts.numNegative, numel(db.val.qidxs));
  fprintf('>>>> done in %s\n', htime(toc(t)));
  if ~isfield(db.val, 'data')
    fprintf('>> Validation images not provided, loading...\n');
    fprintf('>>>> imageDir: %s\n', opts.augment.imageDir);
    % we dont use augmentation for validation data
    optsValData.imageDir = opts.augment.imageDir;
    optsValData.crop = opts.augment.crop;
    optsValData.gpus = opts.augment.gpus;
    optsValData.imageSize = opts.augment.imageSize;
    optsValData.jitterFlip = false;
    optsValData.jitterLocation = false;
    optsValData.jitterScale = 1;
    optsValData.jitterBrightness = 0;
    optsValData.jitterContrast = 0;
    optsValData.jitterSaturation = 0;
    images = cellfun(@(x) cid2filename(x, optsValData.imageDir), db.val.cids, 'UniformOutput', false);
    db.val.data = get_images(images, optsValData);
  end
else
  disp('!!! DB file doesnt exist. KEYBOARD');
  keyboard
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

if isempty(opts.derOutputs)
  error('!!! DEROUTPUTS must be specified when training.\n') ;
end

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('>> Resuming by loading epoch %d\n', start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
  saveState(modelPath(0), net, state) ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.
  rng(epoch + opts.randomSeed) ;
  prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.getBatch = getBatch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;

  % Augment data if necessary: epochSize, imread with jitter, flip pairs, etc
  params.db = augment_database(db, opts.augment);

  if numel(opts.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train') ;
    [net, state] = processEpoch(net, state, params, 'val') ;
    saveState(modelPath(epoch), net, state) ;
    lastStats = state.stats ;
  else
    spmd
      [net, state] = processEpoch(net, state, params, 'train') ;
      [net, state] = processEpoch(net, state, params, 'val') ;
      if labindex == 1
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  saveStats(modelPath(epoch), stats) ;

  if opts.plotStatistics
    switchFigure(1) ; clf ;
    plots = setdiff(cat(2,fieldnames(stats.train)',fieldnames(stats.val)'), {'num', 'time'}) ;
    for p = plots
      p = char(p) ; values = zeros(0, epoch) ; leg = {} ;
      for f = {'train', 'val'}
        f = char(f) ;
        if isfield(stats.(f), p)
          tmp = [stats.(f).(p)] ; values(end+1,:) = tmp(1,:)' ; leg{end+1} = f ;
        end
      end
      subplot(1,numel(plots),find(strcmp(p,plots))) ; plot(1:epoch, values','o-') ;
      xlabel('epoch') ; title(p) ; legend(leg{:}) ; grid on ;
    end
    drawnow ; print(1, modelFigPath, '-dpdf') ;
  end

  if ~isempty(opts.postEpochFn)
    if nargout(opts.postEpochFn) == 0
      opts.postEpochFn(net, params, state) ;
    else
      lr = opts.postEpochFn(net, params, state) ;
      if ~isempty(lr), opts.learningRate = lr; end
      if opts.learningRate == 0, break; end
    end
  end

  % With multiple GPUs, return one copy
  if isa(net, 'Composite'), net = net{1} ; end

end

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% create opts needed for getBatch (so that we dont send all)
getbatchopts.batchSize = params.batchSize;
getbatchopts.numSubBatches = params.numSubBatches;
getbatchopts.numRemine = params.numRemine;
getbatchopts.memoryMapRemine = params.memoryMapRemine;
getbatchopts.gpus = params.gpus;

% initialize with momentum 0
if isempty(state) || isempty(state.solverState)
  state.solverState = cell(1, numel(net.params)) ;
  state.solverState(:) = {0} ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net.move('gpu') ;
  for i = 1:numel(state.solverState)
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gpuArray(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gpuArray, s, 'UniformOutput', false) ;
    end
  end
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  net.setParameterServer(parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

num = 0 ;
epoch = params.epoch ;
adjustTime = 0 ;

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

start = tic ;
time_remine = 0;

nbatch = floor(numel(params.db.(mode).qidxs)/params.batchSize); % floor so only complete batches
for batch=1:nbatch

  % #queries done in this lab and this batch
  bqnum = numel(labindex:numlabs:params.batchSize);

  for subbatch=1:params.numSubBatches
    % #queries done in this lab and this subbatch
    sbqnum = numel(subbatch:params.numSubBatches:bqnum);
    num = num + sbqnum;

    % get inputs and latest negatives
    [inputs, params.db.(mode).nidxs, time_r] = params.getBatch(getbatchopts, params.db, net, batch, subbatch, mode, epoch) ;
    time_remine = time_remine + time_r;

    if params.prefetch
      disp('Prefetch not implemented!!! Keyboard...\n');
      keyboard
    end

    if strcmp(mode, 'train')
      net.mode = 'normal' ;
      net.accumulateParamDers = (subbatch ~= 1) ;
      net.eval(inputs, params.derOutputs, 'holdOn', subbatch < params.numSubBatches) ;
    else
      net.mode = 'test' ;
      net.eval(inputs) ;
    end
  end

  % Accumulate gradient.
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    posnegNum = size(params.db.(mode).pidxs,1) + size(params.db.(mode).nidxs,1);
    pairNum = params.batchSize * posnegNum; % number of pairs
    imNum   = params.batchSize * (1 + posnegNum); % number of images (for batchNorm)
    state = accumulateGradients(net, state, params, pairNum, imNum, parserv) ;
  end

  % Get statistics.
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats.num = num ;
  stats.time = time ;
  stats = params.extractStatsFn(stats,net) ;
  batchSpeed = batch / (time - time_remine);
  % if t == 3*params.batchSize + 1
  %   % compensate for the first three iterations, which are outliers
  %   adjustTime = 4*batchTime - time ;
  %   stats.time = time + adjustTime ;
  % end

  fprintf('%s: ep %02d: bch %4d/%4d:', mode, epoch, batch, nbatch) ;

  fprintf(' %.1f bch/s:', batchSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s %.3f:', f(1:3), stats.(f)) ;
  end
  fprintf('\n') ;
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveSolverState
  state.solverState = [] ;
else
  for i = 1:numel(state.solverState)
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gather(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gather, s, 'UniformOutput', false) ;
    end
  end
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, pairNum, imNum, parserv)
% pairNum : number of pairs per batch #query x (#pos + #neg)
% imNum   : number of images that pass per batch #query x (1 + #pos + #neg)
% -------------------------------------------------------------------------------

numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for p=1:numel(net.params)

  if ~isempty(parserv)
    parDer = parserv.pullWithIndex(p) ;
  else
    parDer = net.params(p).der ;
  end

  if isempty(parDer) % for GOOGLENET part that is not used
    continue;
  end

  switch net.params(p).trainMethod

    case 'average' % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      net.params(p).value = vl_taccum(...
          1 - thisLR, net.params(p).value, ...
          (thisLR/imNum/net.params(p).fanout),  parDer) ;

    case 'gradient'
      thisDecay = params.weightDecay * net.params(p).weightDecay ;
      thisLR = params.learningRate * net.params(p).learningRate ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/pairNum, parDer, ...
                           thisDecay, net.params(p).value) ;

        if isempty(params.solver)
          % Default solver is the optimised SGD.
          % Update momentum.
          state.solverState{p} = vl_taccum(...
            params.momentum, state.solverState{p}, ...
            -1, parDer) ;

          % Nesterov update (aka one step ahead).
          if params.nesterovUpdate
            delta = params.momentum * state.solverState{p} - parDer ;
          else
            delta = state.solverState{p} ;
          end

          % Update parameters.
          net.params(p).value = vl_taccum(...
            1,  net.params(p).value, thisLR, delta) ;

        else
          % call solver function to update weights
          [net.params(p).value, state.solverState{p}] = ...
            params.solver(net.params(p).value, state.solverState{p}, ...
            parDer, params.solverOpts, thisLR) ;
        end
      end
    otherwise
      error('Unknown training method ''%s'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) strcmp(x,'objective') || strcmp(x,'error'), {net.layers.name}));
for i = 1:numel(sel)
  % if net.layers(sel(i)).block.ignoreAverage, continue; end;
  stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end

end
if numGpus >= 1 && cold
  clearMex() ;
  if numGpus == 1
    gpuinfo = gpuDevice(opts.gpus);
    fprintf('>> Running on GPU %s with Index %d\n', gpuinfo.Name, gpuinfo.Index);
  else
    spmd
      clearMex() ;
      gpuinfo = gpuDevice(opts.gpus(labindex));
      fprintf('>> Running on GPU %s with Index %d\n', gpuinfo.Name, gpuinfo.Index);
    end
  end
end

author: filipradenovic

repo: cnnimageretrieval


text:
function z = quad_over_lin( varargin )

%QUAD_OVER_LIN Sum of squares over linear.
%   Z=QUAD_OVER_LIN(X,Y), where X is a vector and Y is a scalar, is equal to
%   SUM(ABS(X).^2)./Y if Y is positive, and +Inf otherwise. Y must be real.
%
%   If X is a matrix, QUAD_OVER_LIN(X,Y) is a row vector containing the values
%   of QUAD_OVER_LIN applied to each column. If X is an N-D array, the operation
%   is applied to the first non-singleton dimension of X.
%
%   QUAD_OVER_LIN(X,Y,DIM) takes the sum along the dimension DIM of X.
%   A special value of DIM == 0 is accepted here, which is automatically
%   replaced with DIM == NDIMS(X) + 1. This has the effect of eliminating
%   the sum; thus QUAD_OVER_LIN( X, Y, NDIMS(X) + 1 ) = ABS( X ).^2 ./ Y.
%
%   In all cases, Y must be compatible in the same sense as ./ with the squared
%   sum; that is, Y must be a scalar or the same size as SUM(ABS(X).^2,DIM).
%
%   Disciplined convex programming information:
%       QUAD_OVER_LIN is convex, nonmontonic in X, and nonincreasing in Y.
%       Thus when used with CVX expressions, X must be convex (or affine)
%       and Y must be concave (or affine).

%
% Check arguments
%

persistent P
if isempty( P ),
    P.map = cvx_remap( ...
        { { 'any' }, { 'nonpositive' } }, ...
        { { 'constant' }, { 'positive' } }, ...
        { { 'l_convex' },  { 'l_concave' } }, ...
        { { 'l_concave' }, { 'l_convex' } }, ...
        { { 'r_affine', 'p_convex', 'n_concave' }, { 'positive' } }, ...
        { { 'r_affine', 'p_convex', 'n_concave' }, { 'concave' } }, ...
        { { 'affine', 'p_convex', 'n_concave' }, { 'positive' } }, ...
        { { 'affine', 'p_convex', 'n_concave' }, { 'concave' } }, ...
        [ 0, 1, 2, 2, 3, 4, 5, 6 ] );
    P.funcs = { @qol_cnst, @qol_log, @qol_sqr, @qol_lin, @qol_sqa, @qol_lin, @qol_cpx };
    P.constant = 1;
    P.name = 'quad_over_lin';
end
[ sx, x, y, dim ] = cvx_get_dimension( varargin, 3, 'zero', true );
if sx(dim) > 1, x = norms( x, 2, dim ); end
z = cvx_binary_op( P, x, y );

function z = qol_cnst( x, y )
z = x .^ 2 / y;

function z = qol_log( x, y )
z = exp( 2 * log( x ) - log( y ) );

function z = qol_sqr( x, y )
z = square( x ) ./ y;

function z = qol_lin( x, y ) %#ok
sz = max( size(x), size(y) );
cvx_begin
    epigraph variable z( sz ) nonnegative_
    { linearize(x), linearize(y), 0.5 * z } == rotated_lorentz( sz, 0 ); %#ok
cvx_end

function z = qol_sqa( x, y )
z = square( abs( x ) ) ./ y;

function z = qol_cpx( x, y ) %#ok
sz = max( size(x), size(y) );
cvx_begin
    epigraph variable z( sz ) nonnegative_
    { linearize(x), linearize(y), 0.5 * z } == rotated_complex_lorentz( sz, 0 ); %#ok
cvx_end

% 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:
function cnnnumgradcheck(net, x, y)
    epsilon = 1e-4;
    er      = 1e-8;
    n = numel(net.layers);
    for j = 1 : numel(net.ffb)
        net_m = net; net_p = net;
        net_p.ffb(j) = net_m.ffb(j) + epsilon;
        net_m.ffb(j) = net_m.ffb(j) - epsilon;
        net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
        net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
        d = (net_p.L - net_m.L) / (2 * epsilon);
        e = abs(d - net.dffb(j));
        if e > er
            error('numerical gradient checking failed');
        end
    end

    for i = 1 : size(net.ffW, 1)
        for u = 1 : size(net.ffW, 2)
            net_m = net; net_p = net;
            net_p.ffW(i, u) = net_m.ffW(i, u) + epsilon;
            net_m.ffW(i, u) = net_m.ffW(i, u) - epsilon;
            net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
            net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
            d = (net_p.L - net_m.L) / (2 * epsilon);
            e = abs(d - net.dffW(i, u));
            if e > er
                error('numerical gradient checking failed');
            end
        end
    end

    for l = n : -1 : 2
        if strcmp(net.layers{l}.type, 'c')
            for j = 1 : numel(net.layers{l}.a)
                net_m = net; net_p = net;
                net_p.layers{l}.b{j} = net_m.layers{l}.b{j} + epsilon;
                net_m.layers{l}.b{j} = net_m.layers{l}.b{j} - epsilon;
                net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
                net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
                d = (net_p.L - net_m.L) / (2 * epsilon);
                e = abs(d - net.layers{l}.db{j});
                if e > er
                    error('numerical gradient checking failed');
                end
                for i = 1 : numel(net.layers{l - 1}.a)
                    for u = 1 : size(net.layers{l}.k{i}{j}, 1)
                        for v = 1 : size(net.layers{l}.k{i}{j}, 2)
                            net_m = net; net_p = net;
                            net_p.layers{l}.k{i}{j}(u, v) = net_p.layers{l}.k{i}{j}(u, v) + epsilon;
                            net_m.layers{l}.k{i}{j}(u, v) = net_m.layers{l}.k{i}{j}(u, v) - epsilon;
                            net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
                            net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
                            d = (net_p.L - net_m.L) / (2 * epsilon);
                            e = abs(d - net.layers{l}.dk{i}{j}(u, v));
                            if e > er
                                error('numerical gradient checking failed');
                            end
                        end
                    end
                end
            end
        elseif strcmp(net.layers{l}.type, 's')
%            for j = 1 : numel(net.layers{l}.a)
%                net_m = net; net_p = net;
%                net_p.layers{l}.b{j} = net_m.layers{l}.b{j} + epsilon;
%                net_m.layers{l}.b{j} = net_m.layers{l}.b{j} - epsilon;
%                net_m = cnnff(net_m, x); net_m = cnnbp(net_m, y);
%                net_p = cnnff(net_p, x); net_p = cnnbp(net_p, y);
%                d = (net_p.L - net_m.L) / (2 * epsilon);
%                e = abs(d - net.layers{l}.db{j});
%                if e > er
%                    error('numerical gradient checking failed');
%                end
%            end
        end
    end
%    keyboard
end

author: rasmusbergpalm

repo: DeepLearnToolbox