text:
classdef Residuals < Exportable

    %--- * --. --- --. .--. ... * ---------------------------------------------
    %               ___ ___ ___
    %     __ _ ___ / __| _ | __|
    %    / _` / _ \ (_ |  _|__ \
    %    \__, \___/\___|_| |___/
    %    |___/                    v 1.0RC1
    %
    %--------------------------------------------------------------------------
    %  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
    %  Written by:        Andrea Gatti
    %  Contributors:      Andrea Gatti, ...
    %  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
    %--------------------------------------------------------------------------

    %% CONSTANTS
    properties (Constant)
        RES_TYPE = {'0: no residuals', '1: PREPRO', '2: U1 engine', '3: U2 engine'};
    end

    properties
        type        % 0,1,2,3 see RES_TYPE
        time        % time as GPS_Time                    GPS_Time [1 x 1] stores n_epoch
        value       % matrix of residuals

        prn         % prn of the satellite (1 x col of pr/ph)
        obs_code    % type of tracking of the column (e.g. GL1C, GL1CL2WI, ...)

        rec_coo     % <optional> Coordinates of the receiver
    end

    methods
        % Creator
        function this = Residuals()
            % Object creator
            this.reset();
        end

    end

    % =========================================================================
    %%  METHODS OBJ MANAGEMENT
    % =========================================================================

    methods % Public Access
        function is_empty = isEmpty(this)
            % Return the empty status
            %
            % SYNTAX
            %   is_empty = this.isEmpty()

            is_empty = this.type == 0 || this.time.isEmpty;
        end

        function reset(this)
            % Reset the stored residuals
            % Empty the object
            %
            % SYNTAX
            %   this.reset

            this.type = 0;
            this.time = GPS_Time();
            this.value = [];

            this.prn = [];
            this.obs_code = '';
            this.rec_coo = Coordinates();
        end

        function import(this, type, time, value, prn, obs_code, rec_coo)
            % Import new residuals (and delete the previous content)
            %
            % INPUT
            %   type        % 0,1,2,3 see RES_TYPE
            %   time        % time as GPS_Time                        GPS_Time [1 x 1] stores n_epoch
            %   value       % matrix of residuals
            %
            %   prn         % prn of the satellite (1 x col of pr/ph)
            %   obs_code    % type of tracking of the column (e.g. GL1C, GL1CL2WI, ...)
            %
            %   rec_coo     % <optional> Coordinates of the receiver
            %
            % SYNTAX
            %   this.import(type, time, value, prn, obs_code, rec_coo)
            this.init(type, time, value, prn, obs_code, rec_coo)
        end

        function append(this, type, time, value, prn, obs_code, rec_coo)
            % Append new residuals to the one already stored
            %
            % INPUT
            %   type        % 0,1,2,3 see RES_TYPE
            %   time        % time as GPS_Time                        GPS_Time [1 x 1] stores n_epoch
            %   value       % matrix of residuals
            %
            %   prn         % prn of the satellite (1 x col of pr/ph)
            %   obs_code    % type of tracking of the column (e.g. GL1C, GL1CL2WI, ...)
            %
            %   rec_coo     % <optional> Coordinates of the receiver
            %
            % SYNTAX
            %   this.import(type, time, value, prn, obs_code, rec_coo)

            res = Residuals();
            res.init(type, time, value, prn, obs_code, rec_coo);

            this.injest(res);
        end

        function injest(this, res)
            % import and append from a residual object file

            % 1) Remove the old overlapped residuals
            % Find the first epoch to remove from the existing res
            if ~res.time.isEmpty
                if isempty(this.time)
                    id_start = this.time.length + 1;
                    id_stop = 0;

                else
                    id_start = find(this.time.getNominalTime.getRefTime(res.time.getNominalTime.first.getMatlabTime) > 0, 1, 'first');
                    id_stop = find(this.time.getNominalTime.getRefTime(res.time.getNominalTime.last.getMatlabTime) <= 0, 1, 'last');

                    if isempty(id_start)
                        id_start = this.time.length + 1;
                    end
                    if isempty(id_stop)
                        id_stop = 0;
                    end
                end
                % Find the last epoch to remove from the existing res
                id_ko = id_start : id_stop;
                if ~isempty(id_ko)
                    this.remEpoch(id_ko);
                end

                % 2) Insert data

                % Insert time
                time = this.time.getEpoch(1 : id_start-1);
                time.append(res.time);
                time.append(this.time.getEpoch(id_start : this.time.length));

                if isempty(this.prn)
                    code_old = [];
                else
                    code_old = Constellation_Collector.obsCode2num(this.obs_code, this.prn);
                    this.remEntry(code_old == 0);
                end
                if isempty(res.prn)
                    code_new = [];
                else
                    % uniform the size of obs_code
                    if size(res.obs_code,2) < size(this.obs_code,2)
                        n_el = size(this.obs_code,2);
                        res.obs_code = [res.obs_code, ' ' * char(ones(size(res.obs_code,1), n_el - size(res.obs_code,2), 'uint8'))];
                    end
                    if size(res.obs_code,2) > size(this.obs_code,2)
                        n_el = size(res.obs_code,2);
                        this.obs_code = [this.obs_code, ' ' * char(ones(size(this.obs_code,1), n_el - size(this.obs_code,2), 'uint8'))];
                    end
                    code_new = Constellation_Collector.obsCode2num(res.obs_code, res.prn);
                    res.remEntry(code_new == 0);
                end

                % new satellites to add
                [code_add, id_add] = setdiff(code_new, code_old);
                [code_common, id_new, id_old] = intersect( code_new, code_old);

                n_obs_new = size(res.value, 1);

                % resize data to add new epochs
                this.value = [this.value(1 : id_start-1, :); nan(n_obs_new, size(this.value, 2)); this.value(id_start : this.time.length, :)];

                % resize data to add columns for the new observables
                this.value = [this.value nan(size(this.value, 1), numel(code_add))];

                % add new data
                this.value(id_start + (0 : n_obs_new - 1) , [id_old; ((end - numel(code_add) + 1) : end)']) = res.value(:, [id_new; id_add]);

                nch_old = size(this.obs_code,2);
                nch_new = size(res.obs_code,2);
                if nch_old == nch_new % length of obs code might not be compatible
                    this.obs_code = [this.obs_code; res.obs_code(id_add, :)];
                elseif nch_old > nch_new
                    this.obs_code = [this.obs_code ; [res.obs_code(id_add, :) char(32*ones(size(res.obs_code,1),nch_old-nch_new,'uint8'))]];
                else
                    this.obs_code = [[this.obs_code char(32*ones(size(this.obs_code,1),nch_new-nch_old,'uint8'))]; res.obs_code(id_add, :)];
                end
                this.prn = [this.prn; res.prn(id_add)];

                this.time = time;

                this.type = res.type;       % save the last type
                this.rec_coo = res.rec_coo.getCopy; % save the last coo
            end
        end

        function remEpoch(this, lid_ko)
            % Remove an epoch from the residuals
            %
            % INPUT
            %   lid_ko  logical array of ko epochs
            %
            % SYNTAX
            %   this.remEpoch(lid_ko);

            if any(lid_ko)
                if ~islogical(lid_ko)
                    id_ko = lid_ko;
                    lid_ko = false(1,this.time.length);
                    lid_ko(id_ko) = true;
                end
                id_ok = ~lid_ko;
                this.time = this.time.getEpoch(id_ok);
                this.value(lid_ko, :) = [];
            end
        end

        function remEntry(this, lid_ko)
            % Remove an entry from the residuals
            %
            % INPUT
            %   lid_ko  logical array of ko entry
            %
            % SYNTAX
            %   this.remEntry(lid_ko);

            if any(lid_ko)
                if ~islogical(lid_ko)
                    id_ko = lid_ko;
                    lid_ko = false(1, size(this.value,2));
                    lid_ko(id_ko) = true;
                end
                this.value(:, lid_ko) = [];
                this.prn(lid_ko) = [];
                this.obs_code(lid_ko, :) = [];
            end
        end

        function cutEpochs(this, new_lim, end_lim)
            % Get the residual only in the time span given
            %
            % SYNTAX
            %   this.cutEpochs(new_limits)
            %   this.cutEpochs(lim_start, lim_stop)

            if nargin == 3
                new_lim = new_lim.getCopy;
                new_lim.append(end_lim);
            end
            time_res = this.time.getNominalTime();
            sss_start = find(time_res.getMatlabTime >= round(new_lim.first.getMatlabTime * 86400 * time_res.getRate) / (86400 * time_res.getRate), 1, 'first');
            sss_stop = find(time_res.getMatlabTime > round(new_lim.last.getMatlabTime * 86400 * time_res.getRate) / (86400 * time_res.getRate), 1, 'first');

            lid = false(this.time.length(), 1);
            if ~isempty(sss_start)
                lid(1 : (sss_start - 1)) = true;
            end
            if ~isempty(sss_stop)
                lid(sss_stop : end) = true;
            end
            this.remEpoch(lid);
        end
    end

    % =========================================================================
    %%  GETTERS
    % =========================================================================
    methods
        function [res, obs_code, prn] = getU1(this, sys_c, freq_c)
            % Get residual matrix of the combined/single_freq processing
            %
            % SYNTAX
            %    [res, obs_code, prn] = this.getU1()
            if this.type < 3
                id_ok =  true(numel(this.prn), 1);
                if nargin > 1 && ~isempty(sys_c)
                    id_ok(this.obs_code(:,1) ~= sys_c) = false;
                end
                if nargin > 2 && ~isempty(freq_c)
                    id_ok(this.obs_code(:,3) ~= freq_c) = false;
                end

                prn = this.prn(id_ok);
                obs_code = this.obs_code(id_ok,:);
                res = this.value(:, id_ok);
            else
                prn = [];
                obs_code = '';
                res = [];
            end
        end

        function [res, obs_code, prn] = getPrU2(this, sys_c, freq_c)
            % Get residual matrix of the uncombined processing
            % Pseudo-codes residuals
            %
            % SYNTAX
            %    [res, obs_code, prn] =  = this.getPrU2()
            if this.type == 3
                id_ok =  this.obs_code(:,2) == 'C';
                if nargin > 1 && ~isempty(sys_c)
                    id_ok(this.obs_code(:,1) ~= sys_c) = false;
                end
                if nargin > 2 && ~isempty(freq_c)
                    id_ok(this.obs_code(:,3) ~= freq_c) = false;
                end

                prn = this.prn(id_ok);
                obs_code = this.obs_code(id_ok,:);
                res = this.value(:, id_ok);
            else
                prn = [];
                obs_code = '';
                res = [];
            end
        end

        function [is_ph] = isPhase(this)
            % get an index that tell which residual are phase
            %
            % SYNTAX:
            %    [is_ph] = this.isPhase()
            if isempty(this.obs_code)
                is_ph = false(size(this.obs_code,1),1);
            else
                is_ph = this.obs_code(:,2) == 'L';
            end
        end

        function [is_co] = isCombined(this)
            % get an index that tell which residual are combined phases
            %
            % SYNTAX:
            %    [is_co] = this.isCombined()
            if size(this.obs_code,2) > 4
                is_co = this.obs_code(:,5) ~= ' ';
            else
                is_co = false(size( this.obs_code,1),1);
            end
        end

        function [res, obs_code, prn] = getPhU2(this, sys_c, freq_c)
            % Get residual matrix of the combined processing
            % Carrier-phase residuals
            %
            % SYNTAX
            %    [res, obs_code, prn] = this.getPhU2()
            if this.type == 3
                id_ok =  this.obs_code(:,2) == 'L';
                if nargin > 1 && ~isempty(sys_c)
                    id_ok(this.obs_code(:,1) ~= sys_c) = false;
                end
                if nargin > 2 && ~isempty(freq_c)
                    id_ok(this.obs_code(:,3) ~= freq_c) = false;
                end

                prn = this.prn(id_ok);
                obs_code = this.obs_code(id_ok,:);
                res = this.value(:, id_ok);
            else
                prn = [];
                obs_code = '';
                res = [];
            end
        end

        function [res, obs_code, prn, time] = getRangeResiduals(this, sys_c)
            % Get range residuals
            %
            % SYNTAX
            %    [res, obs_code, prn] = this.getU1()
            if this.type < 3
                if nargin == 1
                    [res, obs_code, prn] = this.getU1();
                else
                    [res, obs_code, prn] = this.getU1(sys_c);
                end
                time = this.time;
            else
                % To be done!!! in case of uncombined residuals
                prn = [];
                obs_code = '';
                res = [];
                time = GPS_Time;
            end
        end

        function [res, obs_code, prn, type] = get(this, sys_c, freq_c)
            % Get residual matrix stored in residuals
            %
            % INPUT
            %   sys_c   single character describing the constellation e.g. 'G', 'R', 'E', ...
            %   freq_c  single character describing the frequency number e.g. '1', '2', ....
            %
            % SYNTAX
            %    [res, obs_code, prn, type] = this.get(sys_c, freq_c)
            type = this.type;
            switch type
                case 0
                    prn = [];
                    obs_code = '';
                    res = [];
                case {1, 2} % Prepro or unconbined residuals (both uses U1)
                    if nargin == 1
                        [res, obs_code, prn] = this.getU1();
                    elseif nargin == 2
                        [res, obs_code, prn] = this.getU1(sys_c);
                    elseif nargin == 3
                        [res, obs_code, prn] = this.getU1(sys_c, freq_c);
                    end
                case 3 % if I have uncombined residuals, return just phases
                    if nargin == 1
                        [res, obs_code, prn] = this.getPhU2();
                        if isempty(res)
                            [res, obs_code, prn] = this.getPrU2();
                        end
                    elseif nargin == 2
                        [res, obs_code, prn] = this.getPhU2(sys_c);
                        if isempty(res)
                            [res, obs_code, prn] = this.getPrU2(sys_c);
                        end
                    elseif nargin == 3
                        [res, obs_code, prn] = this.getPhU2(sys_c, freq_c);
                        if isempty(res)
                            [res, obs_code, prn] = this.getPrU2(sys_c, freq_c);
                        end
                    end
            end
        end

        function res = getCopy(this)
            % Get a copy of the object
            %
            % SYNTAX
            %   res = this.getCopy();

            res = Residuals;
            res.importFromStruct(this.toStruct);
            res.time = this.time.getCopy;
            res.rec_coo = this.rec_coo.getCopy;
        end

        function sigma = getStd(this)
            % Get std of all the stored residuals
            % WARNING this is still a very rough estimation,
            %         different frequencies have different noise
            %
            % SINTAX
            %   sigma = this.getStd()
            sigma = std(zero2nan(serialize(this.get())), 'omitnan');
        end

    end

    % =========================================================================
    %%  AUXILLIARY
    % =========================================================================
    methods
        function [az, el, sat_coo, sat_name, go_id] = getAzimuthElevation(this, id_ok)
            % Get azimuth and elevation of each satellite stored in residuals
            %
            %
            %   [az, el, sat_coo, sat_name, go_id] = this.getAzimuthElevation();

            core = Core.getCurrentCore;
            sky = core.sky;
            if isempty(core.state.eph_name)
                fw = File_Wizard();
                fw.conjureNavFiles(this.time.first, this.time.last);
            end
            if nargin == 2
                time = this.time.getEpoch(id_ok);
            else
                time = this.time;
            end
            lim = time.first.getCopy;
            lim.append(this.time.last);
            flag_no_clock = true;
            core.initSkySession(lim, flag_no_clock);
            cc = core.getConstellationCollector;
            go_id = unique(cc.getIndex(this.obs_code(:,1), this.prn));
            sat_name = cc.getSatName(go_id);
            [az, el, sat_coo] = sky.getAzimuthElevation(this.rec_coo, time, go_id);
        end
    end

    % =========================================================================
    %%  MULTIPATH
    % =========================================================================
    methods
        function ant_mp = computeMultiPath(this, marker_name, sys_grp, l_max, flag_reg, is_ph, mode, time_lim)
            % Get multi path maps in different modes
            %
            % z_map  Zernike
            % r_map  Zernike + (the methods specified on mode)
            % g_map  Simple Gridding of size [stk_grid_step]
            % c_map  Congruent cells gridding of size [stk_grid_step]
            % g1_map Simple Gridding of size [1x1]
            % c1_map Congruent cells gridding of size [1x1]
            %
            %
            % INPUT
            %   marker_name     name of the station
            %   sys_grp         constellation grouping for MP estimation
            %   l_max           maximum degree of the 3 steps of the zernike interpolation [l_max1, l_max2, l_max3]
            %                   to disable Zernike use l_max = 0
            %   flag_reg        add regularization points (pseudo obs at 0) in the empty areas of the sky
            %   is_ph           use phases instead of pseudo-ranged (default = true)
            %   mode            - 0 use Zernike + staking maps using congruent cells (variable azimuthal resolution)
            %                                   n_step_az = (360/max_n_step_az * cosd(el))
            %                   - [0, n, m] use Zernike + stacking maps using congruent cells with maximum size of [n x m] note that the output will always be 0.5 x 0.5 degrees
            %                   - 1 use Zernike + stacking map of [5 x 1]
            %                   - [1, n, m] use Zernike + stacking map of [n x m]
            % NOTE
            %   For mode 0 and 1 the output matrix will always be 0.5 x 0.5 degrees
            %
            % SYNTAX
            %   this.computeMultiPath(marker_name, sys_grp, <l_max=[43,43,43]>, <flag_reg=true>, <is_ph=true>, <mode=[0 5 1]>)

            state = Core.getCurrentSettings;
            flag_discard_co = false;

            if nargin < 7 || isempty(mode)
                mode = 0; % Z + stacking
            end
            if numel(mode) == 3
                r_grid_step = mode(2,3);
                mode = mode(1);
            else
                r_grid_step = state.mp_zcongruent_up_nxm;
            end

            n_min = state.mp_n_min; % minimum  number of points per cell
            n_min = 15; %  <== DEBUG

            % z_map  Zernike
            % r_map  Zernike + (the methods specified on mode)
            % g_map  Simple Gridding of size [stk_grid_step]
            % c_map  Congruent cells gridding of size [stk_grid_step]
            % g1_map Simple Gridding of size [1x1]
            % c1_map Congruent cells gridding of size [1x1]
            ltype_of_grids = [...
                sum(state.mp_l_max) > 0 ...
                state.mp_zcongruent_up_nxm(1) > 0 ...
                state.mp_regular_up_nxm(1) > 0 ...
                state.mp_congruent_up_nxm(1) > 0 ...
                state.mp_regular_nxm(1) > 0 ...
                state.mp_congruent_nxm(1) > 0];

            log = Core.getLogger();
            ant_mp = struct();
            if this.isEmpty
                log.addWarning('Residuals have not been computed');
            else
                flag_debug = false;
                if flag_debug
                    % Enable all the grid types
                    ltype_of_grids = logical([1 1 1 1 1 1]); % All enabled
                end

                if nargin < 4 || isempty(l_max)
                    l_max = state.mp_l_max;
                end
                % Legacy support
                if numel(l_max) == 3
                    l_max = [l_max(1) l_max];
                end
                if numel(l_max) == 1
                    l_max = [l_max 0 0 0];
                end

                % Depending on the maximum zernike degree change the map resolution
                if max(l_max) > 180
                    grid_step = 0.1;
                elseif max(l_max) > 90
                    grid_step = 0.25;
                else
                    grid_step = 0.5;
                end

                if nargin < 2
                    marker_name = 'UNKN';
                end

                if nargin < 5 || isempty(flag_reg)
                    flag_reg = true;
                end

                if nargin < 6 || isempty(is_ph)
                    % If there are phases use phases
                    is_ph = any((serialize(this.obs_code(:,2:3:end-1))) == 'L');
                end
                if is_ph
                    name = 'Carrier-phase residuals';
                    search_obs = 'L';
                else
                    name = 'Pseudo-ranges residuals';
                    search_obs = 'C';
                end

                deg2rad = pi/180;

                cc = Core.getConstellationCollector();
                log.addMarkedMessage(sprintf('Computing azimuth and elevation for "%s"', marker_name));
                if nargin == 8 && ~isempty(time_lim)
                    id_span = this.time.getNominalTime(1) >= time_lim.first & this.time.getNominalTime(1) <= time_lim.last;
                    [az, el, ~, ~, go_id] = this.getAzimuthElevation(id_span);
                else
                    id_span = true(this.time.length, 1);
                    [az, el, ~, ~, go_id] = this.getAzimuthElevation();
                end
                sys_c_list = cc.getAvailableSys;

                log = Core.getLogger;
                log.addMarkedMessage(sprintf('Computing multipath mitigation coefficients for "%s"', marker_name));

                obs_code = this.obs_code;
                if Core.getCurrentSettings.FLAG_MP_IGNORE_TRK
                    for i = 1 : size(obs_code, 1)
                        obs_code(i, 4:3:end) = '_';
                    end
                end

                for sys_c = sys_c_list(:)'
                    ids = find(ismember(obs_code(:,1), sys_c) & any((obs_code(:,2:3:end-1)) == search_obs, 2));
                    if ~any(ids)
                        log.addWarning(sprintf('No %s found in %s for constellation %s', name, marker_name, cc.getSysName(sys_c)));
                    else
                        obs_id_num = cc.obsCode2num([repmat('G', numel(ids),1) obs_code(ids, 2:end)], zeros(size(ids, 1), 1)); % get all the data of the same frequency - all the satellites
                        uobs_id = unique(obs_id_num);
                        for  t = 1 : numel(uobs_id)
                            id = ids(obs_id_num == uobs_id(t)); % tracking for the specific obs_code
                            trk_code = obs_code(id(1),2:end);
                            if (numel(trk_code) > 5) && (all((trk_code(5:end)) == '_'))
                                trk_code = trk_code(1:4);
                            end
                            % recompute id if grouping is present
                            if numel(sys_grp.(sys_c)) > 1
                                ext_ids = find(ismember(obs_code(:,1), sys_grp.(sys_c)) & any((obs_code(:,2:3:end-1)) == search_obs, 2));
                                obs_id_num = cc.obsCode2num([repmat('G', numel(ext_ids),1) obs_code(ext_ids, 2:end)], zeros(size(ext_ids, 1), 1)); % get all the data of the same frequency - all the satellites
                                id = ext_ids(obs_id_num == uobs_id(t)); % tracking for the specific obs_code
                            end

                            data_found = false;

                            res = zero2nan(this.value(id_span, id));

                            % time filtering, search for badly estimated residuals
                            tmp = zero2nan(movmedian(std(res, 0, 2, 'omitnan'), 21, 'omitnan'));
                            thr = 5 * median(tmp, 'omitnan');
                            id_ko = flagMerge(tmp > thr, 11);
                            res(id_ko, :) = nan;

                            res_go_id = cc.getIndex(obs_code(id, 1), this.prn(id));

                            % Get all the data to interpolate
                            [~, id_sat] = ismember(res_go_id,go_id);

                            az_all = [];
                            el_all = [];

                            % Propagate orbit nans
                            for s = 1 : numel(res_go_id)
                                id_ko = isnan(el(:,id_sat(s))) | isnan(az(:,id_sat(s)));
                                res(id_ko, s) = nan;
                            end

                            res_all = res(~isnan(res(:)));
                            res_smt = Receiver_Commons.smoothMat(res, 'spline', 120/this.time.getRate);
                            res_smt = res_smt(~isnan(res(:)));

                            go_id_list = [];
                            for s = 1 : numel(res_go_id)
                                id_ok = ~isnan(res(:,s));
                                if any(id_ok)
                                    data_found = true;
                                    % res_all = [res_all; serialize(res(id_ok, s))]; %#ok<AGROW>
                                    az_all = [az_all; az(id_ok, id_sat(s)) .* deg2rad]; %#ok<AGROW>
                                    el_all = [el_all; el(id_ok, id_sat(s)) .* deg2rad]; %#ok<AGROW>
                                    go_id_list = [go_id_list; res_go_id(s)]; %#ok<AGROW>
                                end
                            end

                            if data_found
                                m_max = l_max;
                                % Remove outliers
                                id_ok = Core_Utils.polarCleaner(az_all, el_all, res_all, [360, 1; 3 3]) & Core_Utils.polarCleaner(az_all, el_all, res_smt, [360, 1; 3 3]);
                                log.addMessage(log.indent(sprintf('1. Outlier rejection (%.3f%%)', (sum(~id_ok) / numel(id_ok)) * 100), 9));
                                if flag_debug
                                    figure; plot(el_all/pi*180, res_all*1e3, '.'); hold on; plot(el_all(~id_ok)/pi*180, res_all(~id_ok)*1e3, 'o');
                                    legend('residuals', 'outliers');
                                    title((sprintf('Residuals of  %s %s%s [mm]', marker_name, sys_c, trk_code))); drawnow
                                    grid on;
                                end
                                clear res_smt;
                                az_all = az_all(id_ok);
                                el_all = el_all(id_ok);
                                res_all = res_all(id_ok);
                                n_obs = numel(res_all);

                                if flag_reg
                                    log.addMessage(log.indent('2. Preparing regularization', 9));
                                    % Get regularization points based on empty sky areas
                                    [data_map, n_data_map, az_grid, el_grid] = Core_Utils.hemiGridder(az_all, el_all, res_all, [1 1]);
                                    [az_grid, el_grid] = meshgrid(az_grid, el_grid);
                                    az_reg = az_grid(n_data_map <= n_min);
                                    el_reg = el_grid(n_data_map <= n_min);

                                    % In the knots with few data add zero
                                    az_all = [az_all; az_reg];
                                    el_all = [el_all; el_reg];
                                    res_all = [res_all; zeros(size(el_reg))];

                                    if flag_discard_co
                                        % First approach
                                        % Do not consider observations under
                                        % cut-off => map the radius starting from cut-off
                                        id_ko = (el_all * 180/pi) < Core.getState.getCutOff;
                                        az_all(id_ko) = [];
                                        el_all(id_ko) = [];
                                        res_all(id_ko) = [];
                                    else
                                        % Second approach
                                        % Add additional points at the border close to radius 1
                                        % This regularization is needed if the mapping of the radius
                                        % have a cut-off
                                        for i = 0 : 0.5 : (Core.getState.getCutOff - 2.5)
                                            az_all = [az_all; (-pi : 0.05 : pi)'];
                                            el_all = [el_all; i/180*pi + (-pi : 0.05 : pi)'*0];
                                            res_all = [res_all; (-pi : 0.05 : pi)'*0];
                                        end
                                    end
                                end

                                % Ignore data under cut-off
                                if flag_discard_co
                                    Zernike.setCutOff(max(0, (Core.getState.getCutOff - 0.5)));
                                end
                                res_work = res_all;

                                if ~any(ltype_of_grids(1:2))
                                    % Zernike maps are not requested
                                    z_map = 0;
                                    r_map = 0;
                                else

                                    step = 2 + flag_reg*1;
                                    % Perform the first of 3 Zernike steps
                                    Zernike.setMode(0); % Set Zernike engine to recursive
                                    if l_max(1) > 0
                                        log.addMessage(log.indent(sprintf('%d. Zernike coef. estimation (l_max = %d) (1/4)', step, l_max(1)), 9));
                                        Zernike.setModeMF(0);
                                        el2radius = Zernike.getElFun;
                                        [z_par, l, m] = Zernike.analysisAllBlock(l_max(1), m_max(1), az_all, el2radius(el_all), res_work, 1e-5);
                                        [z_map1, az_grid, el_grid] = Zernike.synthesisGrid(l, m, z_par, grid_step);
                                        z_map1((el_grid * 180/pi) < Core.getState.getCutOff, :) = 0; % remove cutoff;
                                        res_work = res_work - Core_Utils.hgrid2scatter(az_all, el_all, z_map1, false, 'spline');
                                        step = step + 1;
                                    else
                                        z_map1 = 0;
                                    end

                                    Zernike.setMode(0); % Set Zernike engine to recursive
                                    if l_max(2) > 0
                                        log.addMessage(log.indent(sprintf('%d. Zernike coef. estimation (l_max = %d) (2/4)', step, l_max(1)), 9));
                                        Zernike.setModeMF(1);
                                        el2radius = Zernike.getElFun;
                                        [z_par, l, m] = Zernike.analysisAllBlock(l_max(2), m_max(2), az_all, el2radius(el_all), res_work, 1e-5);
                                        [z_map2, az_grid, el_grid] = Zernike.synthesisGrid(l, m, z_par, grid_step);
                                        z_map2((el_grid * 180/pi) < Core.getState.getCutOff, :) = 0; % remove cutoff;
                                        res_work = res_work - Core_Utils.hgrid2scatter(az_all, el_all, z_map2, false, 'spline');
                                        step = step + 1;
                                    else
                                        z_map2 = 0;
                                    end

                                    % Perform the second of 3 Zernike steps
                                    if l_max(3) > 0
                                        log.addMessage(log.indent(sprintf('%d. Zernike coef. estimation (l_max = %d) (3/4)', step, l_max(2)), 9));
                                        Zernike.setModeMF(2);
                                        Zernike.setCutOff(0); % This mapping function is more unstable at low elevations => use polar regularization

                                        el2radius = Zernike.getElFun;
                                        [z_par, l, m] = Zernike.analysisAllBlock(l_max(3), m_max(3), az_all, el2radius(el_all), res_work, 1e-5);
                                        [z_map3, az_grid, el_grid] = Zernike.synthesisGrid(l, m, z_par, grid_step);
                                        z_map3((el_grid * 180/pi) < Core.getState.getCutOff, :) = 0; % remove cutoff;
                                        res_work = res_work - Core_Utils.hgrid2scatter(az_all, el_all, z_map3, false, 'spline');
                                        step = step + 1;
                                    else
                                        z_map3 = 0;
                                    end

                                    % Perform the third of 3 Zernike steps
                                    if l_max(4) > 0
                                        log.addMessage(log.indent(sprintf('%d. Zernike coef. estimation (l_max = %d) (4/4)', step, l_max(3)), 9));
                                        Zernike.setModeMF(3);
                                        Zernike.setCutOff(0); % This mapping function is more unstable at low elevations => use polar regularization
                                        [z_par, l, m] = Zernike.analysisAllBlock(l_max(4), m_max(4), az_all, el2radius(el_all), res_work, 1e-5);
                                        [z_map4, az_grid, el_grid] = Zernike.synthesisGrid(l, m, z_par, grid_step);
                                        z_map4((el_grid * 180/pi) < Core.getState.getCutOff, :) = 0; % remove cutoff;
                                        res_work = res_work - Core_Utils.hgrid2scatter(az_all, el_all, z_map4, false, 'spline');
                                        step = step + 1;
                                    else
                                        z_map4 = 0;
                                    end

                                    % Generate maps
                                    log.addMessage(log.indent(sprintf('%d. Compute mitigation grids', step), 9));
                                    step = step + 1;
                                    z_map = z_map1 + z_map2 + z_map3 + z_map4; % z_map  Zernike only

                                    if ~ltype_of_grids(2) % r_map  Zernike + (the methods specified on mode)
                                        r_map = 0;
                                    else
                                        % In this mode a gridding on the residuals is performed
                                        % to retrieve the high frequency multipath a double step procedure is performed,
                                        % first low-res than high-res

                                         if mode == 1
                                            flag_congruent = false;
                                        else % if mode == 0
                                            flag_congruent = true;
                                        end

                                        % low-res step
                                        grid_step1 = [max(r_grid_step(1), 1.5) max(r_grid_step(end), 0.5)];
                                        out_step1 = [min([grid_step1(1), grid_step(1), 0.5]) min([grid_step1(end), grid_step(end) 0.5])];

                                        % high-res step
                                        grid_step2 = r_grid_step;
                                        out_step2 = [min([r_grid_step(1), grid_step(1)]) min([r_grid_step(end), grid_step(end)])];

                                        % Set the final size of the r_grid (it is generally [0.25 x 0.1] deg
                                        out_size = [90 360] ./ fliplr(out_step2);

                                        % Restore data with no regularization
                                        res_work = res_all(1 : n_obs);
                                        res_work = res_work - Core_Utils.hgrid2scatter(az_all(1 : n_obs), el_all(1 : n_obs), Core_Utils.resize2(z_map, out_size));

                                        % Compute low-res grid - STEP 1
                                        log.addMessage(log.indent(sprintf('  - Zernike + LoRes Congruent %g x %g', grid_step1(1), grid_step1(end)), 9));
                                        [r_map1, n_data_map, az_grid, el_grid] = Core_Utils.hemiGridder(az_all(1 : n_obs), el_all(1 : n_obs), res_work, grid_step1, out_step1, flag_congruent, n_min);

                                        % Resize it to the final grid size
                                        r_map = Core_Utils.resize2(r_map1, out_size);

                                        res_work = res_work - Core_Utils.hgrid2scatter(az_all(1 : n_obs), el_all(1 : n_obs), r_map);

                                        % Compute high-res grid - STEP 2 -- this is not CONGRUENT!
                                        % note that high latitude cells are usually under n_min thr
                                        log.addMessage(log.indent(sprintf('  - Zernike + HiRes %g x %g', grid_step2(1), grid_step2(end)), 9));
                                        [r_map2, n_data_map, az_grid, el_grid] = Core_Utils.hemiGridder(az_all(1 : n_obs), el_all(1 : n_obs), res_work, grid_step2, grid_step2, false, 7);
                                        r_map = Core_Utils.resize2(z_map, out_size) + r_map + Core_Utils.resize2(r_map2, out_size);
                                    end
                                end

                                % Compute normal and congruent maps as comparison (no regularization)
                                if ltype_of_grids(3) % g_map  Simple Gridding of size
                                    log.addMessage(log.indent(sprintf('  - Regular grid %g x %g', state.mp_regular_up_nxm(1), state.mp_regular_up_nxm(end)), 9));
                                    g_map = Core_Utils.hemiGridder(az_all, el_all, res_all, state.mp_regular_up_nxm , grid_step, false, n_min);
                                else
                                    g_map = 0;
                                end
                                if ltype_of_grids(4) % c_map  Congruent cells gridding of size [stk_grid_step]
                                    log.addMessage(log.indent(sprintf('  - Congruent grid %g x %g', state.mp_congruent_up_nxm(1), state.mp_congruent_up_nxm(end)), 9));
                                    c_map = Core_Utils.hemiGridder(az_all, el_all, res_all, state.mp_congruent_up_nxm, grid_step, true, n_min);
                                else
                                    c_map = 0;
                                end
                                if ltype_of_grids(5) % g1_map Simple Gridding of size [1x1]
                                    log.addMessage(log.indent(sprintf('  - Regular grid %g x %g', state.mp_regular_nxm(1), state.mp_regular_nxm(end)), 9));
                                    g1_map = Core_Utils.hemiGridder(az_all(res_all ~= 0), el_all(res_all ~= 0), res_all(res_all ~= 0), state.mp_regular_nxm, [], false, n_min);
                                else
                                    g1_map = 0;
                                end
                                if ltype_of_grids(6) % c1_map Congruent cells gridding of size [1x1]
                                    log.addMessage(log.indent(sprintf('  - Congruent grid %g x %g', state.mp_congruent_nxm(1), state.mp_congruent_nxm(end)), 9));
                                    c1_map = Core_Utils.hemiGridder(az_all(res_all ~= 0), el_all(res_all ~= 0), res_all(res_all ~= 0), state.mp_congruent_nxm, state.mp_congruent_nxm, true, n_min);
                                else
                                    c1_map = 0;
                                end

                                if flag_debug
                                    clim = [-1 1] * max(-perc(1e3*(r_map(:)), 0.003),perc(1e3*(r_map(:)), 0.997));
                                    mp_map = z_map1;
                                    [az_grid, el_grid] = Core_Utils.getPolarGrid(360 / size(mp_map, 2), 90 / size(mp_map, 1));
                                    az_grid = Core_Utils.deg2rad(az_grid)';
                                    el_grid = Core_Utils.deg2rad(el_grid);

                                    %figure; imagesc(1e3*(z_map)); colormap((Cmap.get('PuOr', 2^11))); caxis([-5 5]); colorbar;
                                    if numel(z_map1) > 1
                                        figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*(z_map1)); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    end
                                    title((sprintf('Zernike expansion (1) of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow

                                    if numel(z_map2) > 1
                                        figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*(z_map2)); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    end
                                    title((sprintf('Zernike expansion (2) of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow

                                    if numel(z_map3) > 1
                                        figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*(z_map3)); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    end
                                    title((sprintf('Zernike expansion (3) of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow

                                    if numel(z_map4) > 1
                                        figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*(z_map4)); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    end
                                    title((sprintf('Zernike expansion (4) of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow

                                    figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*(z_map)); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    title((sprintf('Zernike expansion of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow

                                    mp_map = r_map;
                                    [az_grid, el_grid] = Core_Utils.getPolarGrid(360 / size(mp_map, 2), 90 / size(mp_map, 1));
                                    az_grid = Core_Utils.deg2rad(az_grid)';
                                    el_grid = Core_Utils.deg2rad(el_grid);

                                    figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*r_map); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    title((sprintf('Final map of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow

                                    mp_map = c_map;
                                    [az_grid, el_grid] = Core_Utils.getPolarGrid(360 / size(mp_map, 2), 90 / size(mp_map, 1));
                                    az_grid = Core_Utils.deg2rad(az_grid)';
                                    el_grid = Core_Utils.deg2rad(el_grid);
                                    figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*(mp_map)); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    title((sprintf('Gridded map with congruent cells of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow

                                    mp_map = g_map;
                                    [az_grid, el_grid] = Core_Utils.getPolarGrid(360 / size(mp_map, 2), 90 / size(mp_map, 1));
                                    az_grid = Core_Utils.deg2rad(az_grid)';
                                    el_grid = Core_Utils.deg2rad(el_grid);
                                    figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*(mp_map)); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    title((sprintf('Gridded map of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow

                                    mp_map = c1_map;
                                    [az_grid, el_grid] = Core_Utils.getPolarGrid(360 / size(mp_map, 2), 90 / size(mp_map, 1));
                                    az_grid = Core_Utils.deg2rad(az_grid);
                                    el_grid = Core_Utils.deg2rad(el_grid);

                                    figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*(mp_map)); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    title((sprintf('Gridded map with congruent cells [1 x 1] of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow

                                    mp_map = g1_map;
                                    [az_grid, el_grid] = Core_Utils.getPolarGrid(360 / size(mp_map, 2), 90 / size(mp_map, 1));
                                    az_grid = Core_Utils.deg2rad(az_grid)';
                                    el_grid = Core_Utils.deg2rad(el_grid);
                                    figure; polarImagesc(az_grid, (pi/2 - el_grid), 1e3*(mp_map)); colormap((Cmap.get('PuOr', 2^11))); caxis(clim); colorbar;
                                    title((sprintf('Gridded map [1 x 1] of %s %s%s [mm]', marker_name, sys_c, trk_code)), 'interpreter', 'none'); drawnow
                                end

                                if ~isfield(ant_mp, sys_c)
                                    ant_mp.(sys_c) = struct;
                                end
                                if ~isfield(ant_mp.(sys_c), trk_code)
                                    if (numel(trk_code) == 8) && sum(trk_code(5:end) == '_') == 4
                                        trk_code = trk_code(1:4);
                                    end
                                    trk_code = strrep(trk_code, ' ', '_'); % structures do not support spaces
                                    ant_mp.(sys_c).(trk_code) = struct;
                                end
                                % Keep multiple solutions in the struct
                                % decide a-posteriori what it's better

                                % Save grids of multi-path
                                ant_mp.(sys_c).(trk_code).z_map = single(z_map);             % Zernike map
                                ant_mp.(sys_c).(trk_code).r_map = single(r_map);             % Zernike map + gridded residuals
                                ant_mp.(sys_c).(trk_code).g_map = single(g_map);             % Simple Gridding of size [stk_grid_step]
                                ant_mp.(sys_c).(trk_code).c_map = single(c_map);             % Congruent cells gridding of size [stk_grid_step]
                                ant_mp.(sys_c).(trk_code).g1_map = single(g1_map);           % Simple Gridding of size [1x1]
                                ant_mp.(sys_c).(trk_code).c1_map = single(c1_map);           % c1_map Congruent cells gridding of size [1x1]
                                % Save Zernike coefficients
                                %ant_mp.(sys_c).(trk_code).z_par = [z_par1 z_par2];
                                %ant_mp.(sys_c).(trk_code).l = l;
                                %ant_mp.(sys_c).(trk_code).m = m;
                            else
                                if ~data_found
                                    log = Core.getLogger;
                                    log.addWarning(sprintf('No %s %s found in %s for constellation %s', name, trk_code, marker_name, cc.getSysName(sys_c)));
                                end

                            end
                        end
                    end
                end
                % Get the time limit of the map solution
                ant_mp.time_lim = this.time.getEpoch(minMax(find(id_span)));
            end
        end
    end

    % =========================================================================
    %%  SHOW
    % =========================================================================
    methods
        function fh_list = showResSkyCartScatter(this, marker_name, sys_c_list, is_ph)
            % Plot residuals of the solution on cartesian axes
            %
            % SYNTAX
            %   this.showResSkyCartScatter(marker_name, sys_c_list, is_ph)
            log = Core.getLogger();
            fh_list = [];
            if this.isEmpty
                log.addWarning('Residuals have not been computed');
            else
                if nargin < 3 || isempty(sys_c_list)
                    sys_c_list = unique(this.obs_code(:,1));
                end
                if nargin < 4 || isempty(is_ph)
                    % If there are phases use phases
                    is_ph = any((serialize(this.obs_code(:,2:3:end-1))) == 'L');
                end
                if is_ph
                    name = 'Carrier-phase residuals';
                    search_obs = 'L';
                    scale = 1e3;
                else
                    name = 'Pseudo-ranges residuals';
                    search_obs = 'C';
                    scale = 1e2;
                end

                if nargin < 2
                    marker_name = '';
                end

                cc = Core.getConstellationCollector();
                [az, el, sat_coo, sat_name, go_id] = this.getAzimuthElevation();
                for sys_c = sys_c_list(:)'
                    ids = find(this.obs_code(:,1) == sys_c & any((this.obs_code(:,2:3:end-1)) == search_obs, 2));
                    if ~any(ids)
                        log.addWarning(sprintf('No %s found in %s for constellation %s', name, marker_name, cc.getSysName(sys_c)));
                    else
                        obs_id_num = cc.obsCode2num(this.obs_code(ids,:), zeros(size(ids, 1), 1));
                        uobs_id = unique(obs_id_num);
                        for  t = 1 : numel(uobs_id)
                            id = ids(obs_id_num == uobs_id(t)); % tracking for the specific obs_code
                            trk_code = this.obs_code(id(1),:);

                            res = this.value(:, id) * scale;
                            res_go_id = cc.getIndex(this.obs_code(id, 1), this.prn(id));

                            fh = figure('Visible', 'off'); fh.Name = sprintf('%03d: %s Res Cart %s', fh.Number, marker_name, strtrim(trk_code)); fh.NumberTitle = 'off';
                            Core_UI.beautifyFig(fh);

                            fh_list = [fh_list; fh]; %#ok<AGROW>
                            fig_name = sprintf('Res_cart_%s_%s_%s_%s', marker_name, cc.getSysName(sys_c), this.time.first.toString('yyyymmdd_HHMM'));
                            fh.UserData = struct('fig_name', fig_name);

                            [~, id_sat] = intersect(go_id, res_go_id);

                            figure(fh); % get focus;
                            hold on;
                            go_id_list = [];
                            for s = 1 : numel(res_go_id)
                                id_ok = find(res(:,s) ~= 0);
                                if any(id_ok)
                                    [~, id_sort] = sort(abs(res(id_ok,s)));
                                    id_ok = id_ok(id_sort);
                                    line = scatter(az(id_ok, id_sat(s)), el(id_ok, id_sat(s)), 45, serialize(res(id_ok, s)), 'filled');
                                    line.UserData = res_go_id(s);
                                    go_id_list = [go_id_list; res_go_id(s)];
                                end
                            end
                            if isempty(go_id_list)
                                delete(fh);
                            else
                                ylim([0 90]); xlim([-180 180]);
                                caxis([-1 1] * max(2, min(6*std(noZero(res),'omitnan'), max(abs(noZero(res(:)))))));
                                colormap((Cmap.get('PuOr', 2^11)));
                                fh.Color = [.95 .95 .95]; cb = colorbar(); cbt = title(cb, iif(scale == 1e2, '[cm]', '[mm]')); cbt.Parent.UserData = cbt; ax = gca; ax.Color = 'none';
                                h = title(sprintf('Satellites residuals - receiver %s - %s\\fontsize{5} \n', strrep(marker_name,'_','\_'), cc.getSysExtName(sys_c)));  h.FontWeight = 'bold';
                                hl = xlabel('Azimuth [deg]'); hl.FontWeight = 'bold';
                                hl = ylabel('Elevation [deg]'); hl.FontWeight = 'bold';

                                Core_UI.addSatMenu(fh, go_id_list);
                                Core_UI.beautifyFig(fh);
                                Core_UI.addExportMenu(fh);
                                Core_UI.addBeautifyMenu(fh);
                                fh.Visible = 'on'; drawnow;
                            end
                        end
                    end
                end
                if isempty(fh_list)
                    log.addWarning('Residuals have not been computed');
                end
            end
        end

        function fh_list = showResSkyPolarScatter(this, marker_name, sys_c_list, is_ph)
            % Plot residuals of the solution on polar axes
            %
            % SYNTAX
            %   this.showResSkyPolarScatter(marker_name, sys_c_list, is_ph)
            log = Core.getLogger();
            fh_list = [];
            if this.isEmpty
                log.addWarning('Residuals have not been computed');
            else
                if nargin < 3 || isempty(sys_c_list)
                    sys_c_list = unique(this.obs_code(:,1));
                end
                if nargin < 4 || isempty(is_ph)
                    % If there are phases use phases
                    is_ph = any((serialize(this.obs_code(:,2:3:end-1))) == 'L');
                end
                if is_ph
                    name = 'Carrier-phase residuals';
                    search_obs = 'L';
                    scale = 1e3;
                else
                    name = 'Pseudo-ranges residuals';
                    search_obs = 'C';
                    scale = 1e2;
                end

                if nargin < 2
                    marker_name = '';
                end

                cc = Core.getConstellationCollector();
                [az, el, ~, ~, go_id] = this.getAzimuthElevation();
                for sys_c = sys_c_list(:)'
                    ids = find(this.obs_code(:,1) == sys_c & any((this.obs_code(:,2:3:end-1)) == search_obs, 2));
                    if ~any(ids)
                        log.addWarning(sprintf('No %s found in %s for constellation %s', name, marker_name, cc.getSysName(sys_c)));
                    else
                        obs_id_num = cc.obsCode2num(this.obs_code(ids,:), zeros(size(ids, 1), 1));
                        uobs_id = unique(obs_id_num);
                        for  t = 1 : numel(uobs_id)
                            id = ids(obs_id_num == uobs_id(t)); % tracking for the specific obs_code
                            trk_code = this.obs_code(id(1),:);

                            res = this.value(:, id) * scale;
                            res_go_id = cc.getIndex(this.obs_code(id, 1), this.prn(id));

                            fh = figure('Visible', 'off'); fh.Name = sprintf('%03d: %s Res Polar %s', fh.Number, marker_name, strtrim(trk_code)); fh.NumberTitle = 'off';
                            Core_UI.beautifyFig(fh);

                            fh_list = [fh_list; fh]; %#ok<AGROW>
                            fig_name = sprintf('Res_polar_%s_%s_%s_%s', marker_name, cc.getSysName(sys_c), this.time.first.toString('yyyymmdd_HHMM'));
                            fh.UserData = struct('fig_name', fig_name);

                            [~, id_sat] = ismember(res_go_id,go_id);

                            figure(fh); % get focus;
                            go_id_list = [];
                            for s = 1 : numel(res_go_id)
                                id_ok = find(res(:,s) ~= 0);
                                if any(id_ok)
                                    [~, id_sort] = sort(abs(res(id_ok,s)));
                                    id_ok = id_ok(id_sort);
                                    line = polarScatter(az(id_ok, id_sat(s))/180*pi, (90 -el(id_ok, id_sat(s)))/180*pi, 45, serialize(res(id_ok, s)), 'filled');
                                    line.UserData = res_go_id(s);
                                    hold on;
                                    go_id_list = [go_id_list; res_go_id(s)]; %#ok<AGROW>
                                end
                            end
                            if isempty(go_id_list)
                                delete(fh);
                            else
                                caxis([-1 1] * max(2, min(6*std(noZero(res),'omitnan'), max(abs(noZero(res(:)))))));
                                colormap((Cmap.get('PuOr', 2^11)));
                                fh.Color = [.95 .95 .95]; cb = colorbar(); cbt = title(cb, iif(scale == 1e2, '[cm]', '[mm]')); cbt.Parent.UserData = cbt; ax = gca; ax.Color = 'none';
                                h = title(sprintf('Satellites residuals\nreceiver %s - %s %s\\fontsize{5} \n', strrep(marker_name,'_','\_'), cc.getSysExtName(sys_c), strtrim(trk_code(2:end))));  h.FontWeight = 'bold';

                                Core_UI.addSatMenu(fh, go_id_list);
                                Core_UI.beautifyFig(fh);
                                Core_UI.addExportMenu(fh);
                                Core_UI.addBeautifyMenu(fh);
                                fh.Visible = 'on'; drawnow;
                            end
                        end
                    end
                end
                if isempty(fh_list)
                    log.addWarning('Residuals have not been computed');
                end
            end
        end

        function fh_list = showRes(this, marker_name, sys_c_list, is_ph)
            % Plot residuals of the solution
            %
            % SYNTAX
            %   fh_list = this.showRes(marker_name, sys_c_list, is_ph)
            log = Core.getLogger();
            fh_list = [];
            if this.isEmpty
                log.addWarning('Residuals have not been computed');
            else
                if nargin < 3 || isempty(sys_c_list)
                    sys_c_list = unique(this.obs_code(:,1));
                end
                if nargin < 4 || isempty(is_ph)
                    % If there are phases use phases
                    is_ph = any((serialize(this.obs_code(:,2:3:end-1))) == 'L');
                end
                if is_ph
                    name = 'Carrier-phase residuals';
                    search_obs = 'L';
                    scale = 1e3;
                else
                    name = 'Pseudo-ranges residuals';
                    search_obs = 'C';
                    scale = 1e2;
                end

                if nargin < 2 || isempty(marker_name)
                    marker_name = '';
                end

                cc = Core.getConstellationCollector();
                for sys_c = sys_c_list(:)'
                    ids = find(this.obs_code(:,1) == sys_c & any((this.obs_code(:,2:3:end-1)) == search_obs, 2));
                    if ~any(ids)
                        log.addWarning(sprintf('No %s found in %s for constellation %s', name, marker_name, cc.getSysName(sys_c)));
                    else
                        obs_id_num = cc.obsCode2num(this.obs_code(ids,:), zeros(size(ids, 1), 1));
                        uobs_id = unique(obs_id_num);
                        for  t = 1 : numel(uobs_id)
                            id = ids(obs_id_num == uobs_id(t)); % tracking for the specific obs_code
                            trk_code = this.obs_code(id(1),:);
                            res = zero2nan(this.value(:, id) * scale);
                            %res = Receiver_Commons.smoothMat(res, 'spline', 10);
                            res_go_id = cc.getIndex(this.obs_code(id, 1), this.prn(id));

                            fh = figure('Visible', 'off'); fh.Name = sprintf('%03d: %s Res %s %s', fh.Number, marker_name, cc.getSysName(sys_c), strtrim(trk_code(2:end))); fh.NumberTitle = 'off';
                            Core_UI.beautifyFig(fh); Core_UI.beautifyFig(fh); drawnow;

                            fig_name = sprintf('Res_polar_%s_%s_%s_%s', marker_name, cc.getSysName(sys_c), strtrim(trk_code(2:end)), this.time.first.toString('yyyymmdd_HHMM'));
                            fh.UserData = struct('fig_name', fig_name);

                            time = this.time.getMatlabTime;

                            go_id_list = [];
                            sat_name_list = {};
                            figure(fh); % get focus;
                            for s = 1 : numel(res_go_id)
                                id_ok = ~isnan(res(:,s));
                                if any(id_ok)
                                    line = Core_Utils.plotSep(time(id_ok), serialize(res(id_ok, s)), '.-', 'Color', Core_UI.getColor(s, numel(res_go_id)));
                                    line.UserData = res_go_id(s);
                                    hold on;
                                    go_id_list = [go_id_list; res_go_id(s)];
                                    sat_name_list = [sat_name_list {cc.getSatName(res_go_id(s))}];
                                end
                            end
                            if isempty(go_id_list)
                                delete(fh);
                            else
                                xlim(minMax(time));
                                ylim([-1 1] * max(abs(ylim)));
                                setTimeTicks();
                                h = title(sprintf('Receiver %s - %s %s\\fontsize{5} \n', strrep(marker_name,'_','\_'), cc.getSysExtName(sys_c), strtrim(trk_code(2:end))));  h.FontWeight = 'bold';
                                [~, icons] = legend(sat_name_list, 'Location', 'NorthEastOutside');
                                icons = icons(numel(sat_name_list) + 2 : 2 : end);
                                for i = 1 : numel(icons)
                                    icons(i).MarkerSize = 18;
                                    icons(i).LineWidth = 2;
                                end
                                ylabel(sprintf('Satellite Residuals %s', iif(scale == 1e2, '[cm]', '[mm]')));
                                Core_UI.addSatMenu(fh, go_id_list);
                                Core_UI.beautifyFig(fh);
                                Core_UI.addExportMenu(fh);
                                Core_UI.addBeautifyMenu(fh);
                                fh_list = [fh_list; fh]; %#ok<AGROW>
                                fh.Visible = 'on'; drawnow;
                            end
                        end
                    end
                end
                if isempty(fh_list)
                    log.addWarning('Residuals have not been computed');
                end
            end
        end

        function fh_list = showResPerSat(this, marker_name, sys_c_list, is_ph)
            % Plot the residuals of phase per tracking
            %
            % INPUT
            %   res     is the matrix of residuals satellite by satellite and can be passed from e.g. NET
            %
            % SYNTAX
            %   fh_list = this.showResPerSat(marker_name, sys_c_list, is_ph)

            log = Core.getLogger();
            fh_list = [];
            if this.isEmpty
                log.addWarning('Residuals have not been computed');
            else
                if nargin < 3 || isempty(sys_c_list)
                    sys_c_list = unique(this.obs_code(:,1));
                end
                if nargin < 4 || isempty(is_ph)
                    % If there are phases use phases
                    is_ph = any((serialize(this.obs_code(:,2:3:end-1))) == 'L');
                end
                if is_ph
                    name = 'Carrier-phase residuals';
                    search_obs = 'L';
                    scale = 1e3;
                else
                    name = 'Pseudo-ranges residuals';
                    search_obs = 'C';
                    scale = 1e2;
                end

                if nargin < 2
                    marker_name = '';
                end

                cc = Core.getConstellationCollector();
                for sys_c = sys_c_list(:)'
                    ids = find(this.obs_code(:,1) == sys_c & any((this.obs_code(:,2:3:end-1)) == search_obs, 2));
                    if ~any(ids)
                        log.addWarning(sprintf('No %s found in %s for constellation %s', name, marker_name, cc.getSysName(sys_c)));
                    else
                        obs_id_num = cc.obsCode2num(this.obs_code(ids,:), zeros(size(ids, 1), 1));
                        uobs_id = unique(obs_id_num);
                        for  t = 1 : numel(uobs_id)
                            id = ids(obs_id_num == uobs_id(t)); % tracking for the specific obs_code
                            trk_code = this.obs_code(id(1),:);

                            res = this.value(:, id) * scale;

                            fh = figure('Visible', 'off'); fh.Name = sprintf('%03d: %s Res %s', fh.Number, marker_name, trk_code); fh.NumberTitle = 'off';
                            Core_UI.beautifyFig(fh); drawnow;

                            fh_list = [fh_list; fh]; %#ok<AGROW>
                            fig_name = sprintf('Res_Per_Sat_%s_%s_%s_%s_%s', trk_code(2:end), marker_name, cc.getSysName(sys_c), trk_code, this.time.first.toString('yyyymmdd_HHMM'));
                            fh.UserData = struct('fig_name', fig_name);

                            ax2 = subplot(1, 24, 19:24);
                            ax1 = subplot(1, 24, 1:16);

                            data_found = false;
                            figure(fh); % get focus;
                            for s = 1 : numel(id)
                                id_ok = find(~isnan(zero2nan(res(:,s))));
                                if any(id_ok)
                                    data_found = true;
                                    [~, id_sort] = sort(abs(res(id_ok, s)));
                                    scatter(ax1, id_ok(id_sort),  this.prn(id(s)) * ones(size(id_ok)), 80, (res(id_ok(id_sort), s)), 'filled');
                                    hold(ax1, 'on');
                                    err = std(zero2nan(res(:,s)), 'omitnan');
                                    if  verLessThan('matlab', '9.4')
                                        plot(ax2, mean(zero2nan(res(:,s)), 'omitnan') + [-err err], this.prn(id(s)) * [1 1], '.-', 'MarkerSize', 15, 'LineWidth', 3, 'Color', [0.6 0.6 0.6]);
                                        plot(ax2, mean(zero2nan(res(:,s)), 'omitnan'), this.prn(id(s)), '.', 'MarkerSize', 30, 'Color', [0.6 0.6 0.6]);
                                    else
                                        errorbar(ax2, mean(zero2nan(res(:,s)), 'omitnan'), this.prn(id(s)), err, '.', 'horizontal', 'MarkerSize', 30, 'LineWidth', 3, 'Color', [0.6 0.6 0.6]);
                                    end
                                    hold(ax2, 'on');
                                end
                            end

                            if ~data_found
                                close(fh)
                                log = Core.getLogger;
                                log.addWarning(sprintf('No %s %s found in %s for constellation %s', name, trk_code, marker_name, cc.getSysName(sys_c)));
                            else
                                cax = caxis(ax1); caxis(ax1, [-1 1] * max(abs(cax)));
                                colormap(Cmap.get('PuOr', 2^11));
                                if min(abs(cax)) > 5
                                    setColorMap('PuOr', caxis(), 0.90, [-5 5])
                                end
                                cb = colorbar(ax1); cb.UserData = title(cb, iif(scale == 1e2, '[cm]', '[mm]')); ax1.Color = [0.9 0.9 0.9];
                                prn_ss = unique(cc.prn(cc.system == sys_c));
                                xlim(ax1, [1 size(res,1)]);
                                ylim(ax1, [min(prn_ss) - 1 max(prn_ss) + 1]);
                                h = ylabel(ax1, 'PRN'); h.FontWeight = 'bold';
                                ax1.YTick = prn_ss;
                                grid(ax1, 'on');
                                h = xlabel(ax1, 'epoch'); h.FontWeight = 'bold';
                                h = title(ax1, sprintf('%s %s %s\\fontsize{5} \n', cc.getSysName(sys_c), strrep(marker_name, '_', '\_'), trk_code(2:end)), 'interpreter', 'tex'); h.FontWeight = 'bold';

                                ylim(ax2, [min(prn_ss) - 1 max(prn_ss) + 1]);
                                xlim(ax2, [-1 1] * (max(max(abs(mean(zero2nan(res(:,:)), 'omitnan'))), ...
                                    max(std(zero2nan(res(:,:)), 'omitnan'))) + 1));
                                ax2.YTick = prn_ss; ax2.Color = [1 1 1];
                                grid(ax2, 'on');
                                xlabel(ax2, sprintf('mean %s', iif(scale == 1e2, 'cm', 'mm')));
                                h = title(ax2, sprintf('mean\\fontsize{5} \n'), 'interpreter', 'tex'); h.FontWeight = 'bold';
                                linkaxes([ax1, ax2], 'y');

                                Core_UI.beautifyFig(fh, 'dark');
                                Core_UI.addBeautifyMenu(fh);
                                fh.Visible = 'on'; drawnow;
                            end
                        end
                    end
                end
            end
        end
    end

    % =========================================================================
    %%  PRIVATE
    % =========================================================================
    methods (Access = private)
        function init(this, type, time, value, prn, obs_code, rec_coo)
            % Init the residual object with new residuals (destroy the prevous content
            %
            % INPUT
            %   type        % 0,1,2,3 see RES_TYPE
            %   time        % time as GPS_Time                        GPS_Time [1 x 1] stores n_epoch
            %   pr          % matrix of pseudorange residuals
            %   ph          % matrix of carrier-phase residuals
            %
            %   value       % matrix of residuals
            %   obs_code    % type of tracking of the column (e.g. GL1C, GL1CL2WI, ...)
            %
            %   rec_coo     % <optional> Coordinates of the receiver
            %
            % SYNTAX
            %   this.import(type, time, pr, ph, prn, obs_code, rec_coo)

            this.type = type;
            this.time = time;
            this.value = value;
            this.prn = prn;
            this.obs_code = obs_code;
            this.rec_coo = rec_coo;

            % Remove entry with no obs_code
            code_ko = Constellation_Collector.obsCode2num(this.obs_code, this.prn);
            this.remEntry(code_ko == 0);
        end
    end

    % =========================================================================
    %%  STATIC
    % =========================================================================
    methods (Static)

    end
end

author: goGPS-Project


text:
function y=packr(x)
% PACKR
% PACKR(x) deletes all the rows of x that contain a MATLAB missing element (NaN)
t=isnan(x);
c=sum(t')'==0;
y=x(c==1,:);

author: Sable

repo: mcbench-benchmarks


text:
%GETWINDOWS Get pixel feature vectors around given pixels in image dataset
%
%		L = GETWINDOWS(A,INDEX,WSIZE,INCLUDE)
%		L = GETWINDOWS(A,[ROW,COL],WSIZE,INCLUDE)
%
% INPUT
%   A       Dataset containing feature images
%   INDEX   Index vector of target pixels in the images (Objects in A)
%   ROW     Column vector of row-coordinates of target pixels
%   COL     Column vector of column-coordinates of target pixels
%   WSIZE   Desired size of rectangular window around target pixels
%   INCLUDE Flag (0/1), indicating whether target pixels should be included
%           (1,default), or not (0) in result.
% OUTPUT
%   L       Index in A of window pixels
%
% DESCRIPTION
% This routine generates all objects in a dataset constructed by image
% features that are in a window of size WSIZE around the target pixels
% given by INDEX, or by [ROW,COL]. If WSIZE is omitted or empty ([],
% default), just the 4 4-conntected neighbors of the target pixels are
% returned. L points to the window pixels, such that A(L,:) is a dataset
% of the corresponding objects.
%
% SEE ALSO (PRTools Guide)
% DATASETS, IM2FEAT

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function L = getwindows(a,index,wsize,include)

if nargin < 4, includeflag = 1; end
if nargin < 3, wsize = []; end

isdataset(a);
isfeatim(a);

imsize = getobjsize(a);

if size(index,2) == 1
	[row,col] = ind2sub(imsize,index);
elseif size(index,2) == 2
	row = index(:,1);
	col = index(:,2);
end
n = length(row);

imsize = getobjsize(a);

if ~isempty(wsize) % rectangular neighborhoods
	if length(wsize) == 1 % square windows
		row1 = ceil(wsize/2)-1; col1 = row1; % # pixels above / left of target pixel
		row2 = floor(wsize/2);  col2 = row2; % # pixels below / right of target pixel
	elseif length(wsize) == 2 % possibly non-square windows
		row1 = ceil(wsize(1)/2)-1;
		row2 = floor(wsize(1)/2);
		col1 = ceil(wsize(2)/2)-1;
		col2 = floor(wsize(2)/2);
	else
		error('Window size should be 1D or 2D')
	end
	[R,C] = meshgrid([-row1:row2],[-col1:col2]);
	k = length(R(:));
	R = repmat(row(:),1,k)+repmat(R(:)',length(row),1);
	C = repmat(col(:),1,k)+repmat(C(:)',length(col),1);
else
	R = repmat(row(:),1,5)+repmat([0 -1 0 1 0],length(row),1);
	C = repmat(col(:),1,5)+repmat([-1 0 0 0 1],length(col),1);
end
JR = [find(R <= 0); find(R > imsize(1))];
R(JR) = []; C(JR) = [];
JC = [find(C <= 0); find(C > imsize(2))];
R(JC) = []; C(JC) = [];
L = sub2ind(imsize,R,C);
%b(L) = ones(length(L),1);

if include % we are done
	L = unique(L);
else       % remove given objects
	b = zeros(imsize);             % create an image of the right size
	b(L) = ones(length(L),1);      % flag the objects we found.
	Z = sub2ind(imsize,row,col);   % for all original objects,
	b(Z) = zeros(length(Z),1);     % remove flags
	L = find(b > 0);               % and see what is left
end
%b = a(L,:);

L = unique(L);
J = [find(L<1); find(L>prod(imsize))];
L(J) = [];
%b = a(L,:);

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% Mathematics Q454519
% https://math.stackexchange.com/questions/454519
% Least Squares with L2 Linear Norm Regularization (Not Squared)
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     22/08/2017
%   *   First release.

%% General Parameters

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = OFF;

%% Simulation Parameters

numRows     = 4;
paramBeta   = 0.5;

numIterations   = 25;

%% Generate Data

vU = randn([numRows, 1]);
vX = 1 + rand([numRows, 1]);

mX = diag(1 ./ vX);

%% Solution by CVX

cvx_begin('quiet')
    cvx_precision('best');
    variable vW(numRows)
    minimize( (0.5 * sum_square(vW - vU)) + (paramBeta * norm(vW ./ vX)) )
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vW.'), ' ]']);
disp([' ']);

%% Solution by Iterative Reweighted Least Squares (IRLS)

hObjFun = @(vW) (0.5 * sum((vW - vU) .^ 2)) + (paramBeta * norm(vW ./ vX));

vObjVal = zeros([numIterations, 1]);

vW = vU;
vObjVal(1) = hObjFun(vW);

for ii = 2:numIterations

    % Calculation of the inverse term. Since all diagonlas terms can be
    % done using vectors.
    fctrTerm    = paramBeta / norm(vW ./ vX);
    vXX         = (1 ./ (vX .^ 2));
    vInvXX      = 1 ./ ((fctrTerm * vXX) + ones([numRows, 1]));

    vW = vInvXX .* vU;

    vObjVal(ii) = hObjFun(vW);
end

disp([' ']);
disp(['Iterative Reweighted Least Squares (IRLS) Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(vObjVal(numIterations))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vW.'), ' ]']);
disp([' ']);

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
hLineSeries = plot(1:numIterations, [vObjVal, cvx_optval * ones([numIterations, 1])]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', ['Objective Function Value vs. Iteration'], ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'Iteration Number', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Objective Function Value', ...
    'FontSize', fontSizeAxis);
set(hAxes, 'XLim', [1, numIterations]);
hLegend = ClickableLegend({['IRLS'], ['Optimal Value (CVX)']});
set(hAxes, 'LooseInset', [0.07, 0.07, 0.07, 0.07]);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function obj = removeEvent(obj, name)
    % Remove event functions defined for the system
    %
    % Parameters:
    % name: the name of the constraint @type cellstr

    assert(ischar(name) || iscellstr(name), ...
        'The name must be a character vector or cellstr.');
    if ischar(name), name = cellstr(name); end

    for i=1:length(name)
        constr = name{i};

        if isfield(obj.EventFuncs, constr)
            obj.EventFuncs = rmfield(obj.EventFuncs,constr);
        else
            error('The event (%s) does not exist.\n',constr);
        end
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [sigma,shrinkage]=cov1para(x,shrink)

% function sigma=cov1para(x)
% x (t*n): t iid observations on n random variables
% sigma (n*n): invertible covariance matrix estimator
%
% Shrinks towards one-parameter matrix:
%    all variances are the same
%    all covariances are zero
% if shrink is specified, then this value is used for shrinkage

% This version: 04/2014

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is released under the BSD 2-clause license.

% Copyright (c) 2014, Olivier Ledoit and Michael Wolf
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in the
% documentation and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% de-mean returns
[t,n]=size(x);
meanx=mean(x);
x=x-meanx(ones(t,1),:);
%x = bsxfun(@minus,x,meanx);

% compute sample covariance matrix
sample=(1/t).*(x'*x);

% compute prior
%meanvar=mean(diag(sample));
meanvar = trace(sample)/n;
prior=meanvar*eye(n);

%if (nargin < 2 | shrink == -1) % compute shrinkage parameters

  % what we call p
  y=x.^2;
  phiMat=y'*y/t-sample.^2;
  phi=sum(sum(phiMat));

  % what we call r is not needed for this shrinkage target

  % what we call c
  %gamma=norm(sample-prior,'fro')^2;
  gamma = sample-prior;
  gamma = sum(gamma(:).^2);

  % compute shrinkage constant
  kappa=phi/gamma;
  shrinkage=max(0,min(1,kappa/t));

%else % use specified number
%  shrinkage=shrink;
%end

% compute shrinkage estimator
sigma=shrinkage*prior+(1-shrinkage)*sample;

author: brainstorm-tools


text:
function h = scattercloud(x,y,n,l,clm,cmap)
%SCATTERCLOUD display density of scatter data
%   SCATTERCLOUD(X,Y) creates a scatterplot of X and Y, displayed over a
%   surface representing the smoothed density of the points.  The density is
%   determined with a 2D histogram, using 25 equally spaced bins in both
%   directions.
%   SCATTERCLOUD(X,Y,N) uses N equally spaced bins.
%   SCATTERCLOUD(X,Y,N,L) uses L as a parameter to the smoothing algorithm.
%    Defaults to 1.  Larger values of L lead to a smoother density, but a
%    worse fit to the original data.
%   SCATTERCLOUD(X,Y,N,L,CLM) uses CLM as the color/linestyle/marker for
%    the scatter plot.  Defaults to 'k+'.
%   SCATTERCLOUD(X,Y,N,L,CLM,CMAP) uses CMAP as the figure's colormap.  The
%    default is 'flipud(gray(256))'.
%   H = SCATTERCLOUD(...) returns the handles for the surface and line
%    objects created.
%
%   Example:
%
%     scattercloud(1:100 + randn(1,100), sin(1:100) + randn(1,100),...
%                  50,.5,'rx',jet(256))
%
%   References:
%     Eilers, Paul H. C. & Goeman, Jelle J. (2004). Enhancing scatterplots
%   with smoothed densities. Bioinformatics 20(5), 623-628.

error(nargchk(2,6,nargin),'struct');

x = x(:);
y = y(:);

if length(x) ~= length(y)
    error('SCATTERCLOUDDataVectorSizesDoNotMatch','The number of elements in x and y do not match')
end

if nargin < 6
    cmap = flipud(gray(256));
end

if nargin < 5
    clm = 'k+';
end

if nargin < 4
    l = 1;
end

if nargin < 3
    n = 25;
end

% min/max of x and y
minX = min(x);
maxX = max(x);
minY = min(y);
maxY = max(y);

% edge locations
xEdges = linspace(minX,maxX,n);
yEdges = linspace(minY,maxY,n);

% shift edges
xDiff = xEdges(2) - xEdges(1);
yDiff = yEdges(2) - yEdges(1);
xEdges = [-Inf, xEdges(2:end) - xDiff/2, Inf];
yEdges = [-Inf, yEdges(2:end) - yDiff/2, Inf];

% number of edges
numX = numel(xEdges);
numY = numel(yEdges);

% hold counts
C = zeros(numY,numX);

% do counts
for i = 1:numY-1
    for j = 1:numX-1
        C(i,j) = length(find(x >= xEdges(j) & x < xEdges(j+1) &...
                             y >= yEdges(i) & y < yEdges(i+1)));
    end
end

% get rid of Infs from the edges
xEdges = [xEdges(2) - xDiff,xEdges(2:end-1), xEdges(end-1) + xDiff];
yEdges = [yEdges(2) - yDiff,yEdges(2:end-1), yEdges(end-1) + yDiff];

% smooth the density data, in both directions.
C = localSmooth(localSmooth(C,l)',l)';

% create the graphics
ax = newplot;
s = surf(xEdges,yEdges,zeros(numY,numX),C,...
         'EdgeColor','none',...
         'FaceColor','interp');
view(ax,2);
colormap(ax,cmap);
grid(ax,'off');
holdstate = get(ax,'NextPlot');
set(ax,'NextPlot','add');
p = plot(x,y,clm);
axis(ax,'tight');
set(ax,'NextPlot',holdstate)

% outputs
if nargout
    h = [s;p];
end

function B = localSmooth(A,L)
r = size(A,1);
I = eye(r);
D1 = diff(I);
D2 = diff(I,2);
B = (I + L ^ 2 * D2' * D2 + 2 * L * D1' * D1) \ A;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode