text:
%   CLASS Core Block
% =========================================================================
%
% DESCRIPTION
%   Class to manage pre processing utilities
%
% EXAMPLE
%   go_block = Core_Block();
%
% FOR A LIST OF CONSTANTs and METHODS use doc Core_PP
%
% Note for the future: the class uses the current obs storage of goGPS
% -> switch to objects for rover and master observations is suggested

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

classdef Core_PP < handle

    properties (Constant, Access = private)
    end

    properties (Access = public)% Public Access
        log
        state
    end

    % ==================================================================================================================================================
    %  CREATOR
    % ==================================================================================================================================================

    methods (Static)
        function this = Core_PP()
            this.init()
        end
    end

    % ==================================================================================================================================================
    %  PRE PROCESSING MAIN FUNCTIONS
    % ==================================================================================================================================================

    methods % Public Access
        function init(this)
            % Init singletons
            %
            % SYNTAX
            %   this.init()
            this.log = Core.getLogger();
            this.state = Core.getState();
        end

        function [ph] = jmpFix(this, ph, lambda)
            % after a missing data try to fix the jmp in the observations with an N * lambda correction
            cs_factor = this.state.cs_thr;
            ph = ph';
            for a = 1 : size(ph,2)
                arc = zero2nan(ph(:,a));
                x_arc = (1 : numel(arc))';
                arc_fit = zeros(size(arc));

                int_factor = zeros(size(arc));
                int_factor(~isnan(arc)) = lambda(a,1) * cs_factor;
                arc = arc ./ int_factor;
                if any(arc)
                    %PLOT:figure;
                    %PLOT:plot(x_arc, arc, 'LineWidth', 1); hold on;
                    lim = getFlagsLimits(~isnan(arc));
                    for l = 1 : size(lim, 1) - 1
                        id = lim(1, 1) : lim(l, 2);
                        id = id(~isnan(arc(id)));
                        id = id((max(1, numel(id) - this.state.getMinArc(this.time.getRate)+1) : end));

                        p = polyfit(x_arc(id), arc(id), 3);
                        if l < size(lim, 1)
                            arc_fit1 = polyval(p, x_arc([lim(l, 2) lim(l + 1, 1)]));

                            id2 = lim(l + 1, 1) : lim(l + 1, 2);
                            id2 = id2(1 : min(numel(id), this.state.getMinArc (this.time.getRate)));
                            p2 = polyfit(x_arc(id2), arc(id2), 3);
                            arc_fit2 = polyval(p2, x_arc([lim(l, 2) lim(l + 1, 1)]));

                            arc(lim(l + 1, 1) : end) = arc(lim(l + 1, 1) : end) + round(mean(arc_fit1 - arc_fit2));
                        end
                        %PLOT: plot(x_arc, arc);
                    end
                    %PLOT: plot(x_arc, arc, 'LineWidth', 2); hold on;
                end
                ph(:,a) = nan2zero(arc) .* int_factor;
            end
            ph = ph';

            %PLOT: figure; plot(zero2nan(ph));
            %PLOT: d_ph = diff(zero2nan([ph(1,:); ph])); figure; plot(d_ph - movmedian(d_ph, 3, 'omitnan'));
        end
    end

    % ==================================================================================================================================================
    %  STATIC LAUNCHERS
    % ==================================================================================================================================================

    methods (Static) % Public Access
        function [data, flagged_out] = flagRawObs(data, time_set, time_ref, thr_factor, thr_min, win_size)
            % Flag the data above max(thr_min, thr_factor) * mean(movstd))
            % SYNTAX: data = flagRawObs(data,  time_set, time_ref, thr_factor, <thr_min = 0 >)
            if nargin < 5
                thr_min = 0;
            end
            if nargin < 6
                win_size = 5;
            end

            flag = 1;
            i = 0;
            % Interpolate to reference time to find outliers in the 4th derivative
            if ~isempty(time_set)
                data_interp = zeros(size(data));
                for s = 1 : size(data,2)
                    tmp = nan2zero(data(:,s));
                    if any(tmp)
                        id_interp = (conv(tmp ~= 0, ones(9,1), 'same') > 0);
                        data_interp(id_interp,s) = interp1(time_set, zero2nan(data(:,s)), time_ref(id_interp), 'spline', 'extrap');
                    end
                end
            else
                data_interp = data;
            end

            % loop for outliers rejection
            flagged_out = false(size(data));
            while (sum(flag(:)) > 0) && (i < 5)
                sensor = zero2nan(data_interp);
                sensor = bsxfun(@minus, sensor, median(sensor, 2, 'omitnan'));
                i = i + 1;
                % expand the set of observations to allow flogging at the border of the valid intervals
                sensor = movstd(sensor, win_size, 'omitnan');
                flag = sensor > max(thr_min, thr_factor * mean(serialize(zero2nan(sensor .* ~isnan(zero2nan(data)))), 'omitnan'));
                data(flag) = NaN;
                data_interp(flag) = NaN;
                flagged_out  = flagged_out | flag;
            end
        end

        function [data, flagged_out] = flagRawObsD4(data, time_set, time_ref, thr_factor, thr_min, win_size)
            % Flag the data above max(thr_min, thr_factor) * mean(movstd))
            %
            % SYNTAX
            %   data = flagRawObsD4(data,  time_set, time_ref, thr_factor, <thr_min = 0 >)
            if nargin < 5
                thr_min = 0;
            end
            if nargin < 6
                win_size = 5;
            end

            flag = 1;
            i = 0;
            % Interpolate to reference time to find outliers in the 4th derivative
            if ~isempty(time_set)
                data_interp = zeros(size(data));
                for s = 1 : size(data,2)
                    tmp = zero2nan(data(:,s));
                    id_interp = ~isnan(tmp);
                    if sum(~isnan(tmp) > 2)
                        data_interp(id_interp,s) = interp1(time_set(id_interp), tmp(id_interp), time_ref(id_interp), 'spline', 'extrap');
                    end
                end
            else
                data_interp = data;
            end

            % loop for outliers rejection
            flagged_out = false(size(data));
            while (sum(flag(:)) > 0) && (i < 5)
                sensor = Core_Utils.diffAndPred(zero2nan(data_interp), 4);
                sensor = bsxfun(@minus, sensor, median(sensor, 2, 'omitnan'));
                i = i + 1;
                % expand the set of observations to allow flogging at the border of the valid intervals
                sensor = movstd(sensor, win_size, 'omitnan');
                flag = sensor > max(thr_min, thr_factor * mean(serialize(zero2nan(sensor .* ~isnan(zero2nan(data)))), 'omitnan'));
                data(flag) = NaN;
                data_interp(flag) = NaN;
                flagged_out  = flagged_out | flag;
            end
        end

        function [data, flagged_out] = flagBorders(data, win_size)
            % Flag the data above max(thr_min, thr_factor) * mean(movstd))
            % SYNTAX: [data, flagged_out] = flagBorders(data, win_size)
            if nargin < 2
                win_size = 20;
            end

            % Interpolate to reference time
            time_set = 1 : size(data,1);
            id_interp = flagExpand(~isnan(data), 5) & isnan(data);
            data_interp = data;
            for s = 1 : size(data,2)
                tmp = data(~isnan(data(:,s)),s);
                if numel(tmp) > 5
                    data_interp(id_interp(:,s),s) = interp1(time_set(~isnan(data(:,s))), tmp, time_set(id_interp(:,s)), 'spline', 'extrap');
                end
            end

            % loop for outliers rejection
            sensor = Core_Utils.diffAndPred(zero2nan(data_interp), 4);
            sensor = bsxfun(@minus, sensor, median(sensor, 2, 'omitnan'));
            % expand the set of observations to allow flagging at the border of the valid intervals
            sensor(isnan(data)) = 0;
            sensor = maxfilt_mat(movstd(sensor, ceil(win_size/2), 'omitnan'),5);
            sensor(isnan(data)) = 0;
            std_max = sensor; std_max(flagExpand(sensor == 0, win_size)) = 0;
            sensor_max = maxfilt_mat(max(std_max, [], 2), ceil(win_size/2));
            sensor_max = simpleFill1D(sensor_max, sensor_max == 0);
            % remove possible outliers over thr and small arcs < 5 epochs
            flagged_out = ~Core_PP.remShortArcs(~(sensor > sensor_max | isnan(data))', 5)' & ~isnan(data);
            data(flagged_out) = NaN;
        end
    end

    % ==================================================================================================================================================
    %  STATIC public misc utilities
    % ==================================================================================================================================================

    methods (Static) % Public Access
        function [obs, dt_dj, is_jumping] = remDtJumps(obs)
            % remove a jumps of the clock from ph/pr
            % this piece of code is very very criptic, but it seems to work
            % review this whenever possible
            %
            % SYNTAX:
            %   [obs, dt_dj, is_jumping] = remDtJumps(obs)
            obs = zero2nan(obs);
            dt_dj = 0;
            is_jumping = false;

            % Iterate for very particular instable cases
            for i = 1 : 2
                % Try to find the biggest jumps (experimental)
                ddt = median(Core_Utils.diffAndPred(zero2nan(obs),1), 2, 'omitnan');
                pos_jmp = abs(ddt) > 1e5;
                d3dt = median(Core_Utils.diffAndPred(zero2nan(obs),3), 2, 'omitnan');
                ddt = cumsum(cumsum(nan2zero(d3dt)));
                dt0 = 0;
                if sum(pos_jmp) > 0
                    is_jumping = true;
                    ddt = ddt - simpleFill1D(ddt, pos_jmp);
                    dt0 = cumsum(nan2zero(ddt));
                    ph_bk = obs;
                    obs = bsxfun(@minus, obs, dt0);

                    % approach on 3rd and 4th derivate (classic)
                    ddt_sensor = cumsum(cumsum(nan2zero(d3dt)) - medfilt_mat(cumsum(nan2zero(d3dt)), 3));
                    % check if there is any discontinuity in the clock drift
                    clock_thresh = 1e3;
                    pos_jmp = abs(ddt_sensor - medfilt_mat(ddt_sensor, 3)) > clock_thresh;
                    pos_jmp(1:2) = 0;
                    if sum(pos_jmp) > 0
                        ddt = ddt - simpleFill1D(ddt, pos_jmp);
                        dt = cumsum(ddt);
                        ph_bk = obs;
                        obs = bsxfun(@minus, obs, dt);
                        d4dt = median(Core_Utils.diffAndPred(zero2nan(obs),4), 2, 'omitnan');
                        dobs = Core_Utils.diffAndPred(obs);
                        % find jmps on the median 4th derivate
                        jmp_candidate = flagExpand(abs(nan2zero(d4dt)) > clock_thresh, 2);
                        % detect the real jmp index
                        lim = getFlagsLimits(jmp_candidate);
                        jmp = false(size(jmp_candidate));
                        for l = 1 : size(lim,1)
                            dtmp = dobs(lim(l,1) : lim(l,2), :);
                            [~, id_max] = max(abs(nan2zero(median(dtmp, 2 ,'omitnan')))');
                            jmp(id_max + lim(l,1) - 1) = true;

                            % find the magnitude of the jump
                            dtmp_f = simpleFill1D(dtmp, isnan(dtmp)); % fill temp data
                            tmp_diff = diff(dtmp_f);
                            tmp_diff(:, isnan(dtmp(id_max,:))) = nan; % do not use filled data for estimating the magnitude of jump

                            d4dt(id_max + lim(l,1) - 1) = median(sign(dobs(id_max + lim(l,1) - 1,:)),'omitnan') * ...
                                median(abs(max(tmp_diff)),'omitnan');
                        end
                        % velocity offsets beteween sat
                        d4dt(~jmp) = 0;
                        obs = bsxfun(@minus, obs, cumsum(nan2zero(d4dt)));
                        dt = dt + cumsum(d4dt);
                    else
                        dt = zeros(size(ddt));
                    end
                    dt_dj = dt_dj + dt0 + dt;
                end
            end
            dt_dj = dt_dj / 299792458;
        end

        function [ph_out, dt] = remDtHF(ph)
            % [ !! ] EXPERIMENTAL FUNCTION
            % Remove high frequency from phases estimated by the 4th derivative
            % WARNING: this code works but it has a lot of border effects, especially if the dataset it's not continuous
            %
            % SYNTAX:
            %   [ph_out, dt] = remDtHF(ph)
            d4dt = median(Core_Utils.diffAndPred(zero2nan(ph),4), 2, 'omitnan');

            % Filter low frequencies:
%             dt_hf = cumsum(cumsum(cumsum(cumsum(nan2zero(d4dt)))));
%             x = (1 : numel(dt_hf))';
%             ws = 5;
%             margin = 2 * round(ws/2);
%             xi = (1 - margin : numel(dt_hf) + margin)';
%             dt_lf = interp1(x(~isnan(dt_hf)), dt_hf(~isnan(dt_hf)), xi, 'spline', 'extrap');
%             dt_lf = splinerMat([],medfilt_mat(dt_lf, 5), 5);
%             dt_hf = dt_hf - dt_lf(margin + 1 : end - margin);
            dt = cumsum(Core_PP.highPass(nan2zero(d4dt), 313, 0.01, 101));
            dt = cumsum(Core_PP.highPass(nan2zero(dt), 177, 0.01, 101));
            dt = cumsum(Core_PP.highPass(nan2zero(dt), 143, 0.01, 101));
            dt = cumsum(Core_PP.highPass(nan2zero(dt), 107, 0.01, 101));
            ph_out = nan2zero(bsxfun(@minus, zero2nan(ph), dt));
            dt = dt / 299792458;
        end

        function [pr, ph, dt_pr, dt_ph] = correctTimeDesyncHF(time_ref, time, pr, ph)
            %   Correction of jumps in code and phase due to dtR and time de-sync
            % SYNTAX:
            %   [pr, ph, dt_pr, dt_ph] = correctTimeDesync(time_ref, time, pr, ph)

            log = Logger.getInstance();

            time_desync  = round((time_ref - time) * 1e7) / 1e7;
            %figure(1); clf; plot(diff(zero2nan(ph))); hold on;
            %figure(2); clf; plot(diff(zero2nan(pr))); hold on;

            ph_ds = bsxfun(@minus, ph, time_desync .* 299792458);
            pr_ds = bsxfun(@minus, pr, time_desync .* 299792458);
            % if adding the desync will improve the std it means that the receiver does not compensate for it
            [ph, flag] = Core_PP.testDesyncCorrection(ph, ph_ds);
            if flag
                time_desync_ph = time_desync;
                log.addMessage('Correcting phase for time desync', 100);
            else
                time_desync_ph = 0;
            end

            [pr, flag] = Core_PP.testDesyncCorrection(pr, pr_ds);
            if flag
                time_desync_pr = time_desync;
                log.addMessage('Correcting pseudo-ranges for time desync', 100);
            else
                time_desync_pr = 0;
            end
            clear pr_ds ph_ds;
            [ph, dt_ph_jumps] = Core_PP.remDtJumps(ph);
            [ph, dt] = Core_PP.remDtHF(ph);
            dt_ph = dt_ph_jumps + dt;
            %figure(1); plot(diff(zero2nan(ph)),'.k');

            % correct the pseudo-ranges for HF
            pr_lf = bsxfun(@minus, pr, dt .* 299792458);
            [pr, flag] = Core_PP.testDesyncCorrection(pr, pr_lf);
            if flag
                dt_pr = dt_ph;
                log.addMessage('Correcting pseudo-ranges for dt HF as estimated from phases observations', 100);
            else
                dt_pr = zeros(size(dt_ph));
            end

            % correct the pseudo-ranges for jumps
            pr_dj = bsxfun(@minus, pr, dt_ph_jumps .* 299792458);
            [pr, flag] = Core_PP.testDesyncCorrection(pr, pr_dj);
            if flag
                dt_pr = dt_pr + dt_ph_jumps;
                log.addMessage('Correcting pseudo-ranges for dt as estimated from phases observations', 100);
            end

            % in some receivers the pseudo-range is not continuous while the phase are ok
            [pr_ds, dt_pr_jumps] = Core_PP.remDtJumps(pr);
            [pr, flag] = Core_PP.testDesyncCorrection(pr, pr_ds);
            if flag
                dt_pr = dt_pr_jumps + dt_pr;
                log.addMessage('Correcting pseudo-ranges for dt as estimated from their observations', 100);
            end

            dt_ph = dt_ph + time_desync_ph;
            dt_pr = dt_pr + time_desync_pr;
            %figure(2); plot(diff(zero2nan(pr)),'.k');
        end

        function [pr, ph, dt_pr, dt_ph] = correctTimeDesync(time_ref, time, pr, ph)
            %   Correction of jumps in code and phase due to dtR and time de-sync
            % SYNTAX:
            %   [pr, ph, dt_pr, dt_ph] = correctTimeDesync(time_ref, time, pr, ph)

            log = Core.getLogger();

            id_not_empty = time ~= 0;

            time_desync = nan(size(time));
            time_desync(id_not_empty)  = round((time(id_not_empty) - time_ref(id_not_empty)) * 1e7) / 1e7; % the rinex time has a maximum of 7 significant decimal digits
            time_desync = simpleFill1D(time_desync,isnan(time_desync),'nearest');

            if any(time_desync)
                [ph_dj, dt_ph_dj] = Core_PP.remDtJumps(ph);
                [pr_dj, dt_pr_dj] = Core_PP.remDtJumps(pr);
                ddt_pr = Core_Utils.diffAndPred(dt_pr_dj);

                %% time_desync is a introduced by the receiver to maintain the drift of the clock into a certain range
                ddt = [0; diff(time_desync)];
                ddrifting = ddt - ddt_pr;
                drifting = cumsum(ddt - ddt_pr);

                % Linear interpolation of ddrifting
                jmp_reset = find(abs(ddt_pr) > 1e-7); % points where the clock is reset
                jmp_fit = setdiff(find(abs(ddrifting) > 1e-7), jmp_reset); % points where desync interpolate the clock
                d_points = [drifting(jmp_reset); drifting(jmp_fit) - ddrifting(jmp_fit)/2];
                jmp = [jmp_reset; jmp_fit];
                drifting = interp1(jmp, d_points, (1 : numel(drifting))', 'spline');

                dt_ph = drifting + dt_ph_dj;
                dt_pr = drifting + dt_pr_dj;

                t_offset = round(mean(dt_pr(jmp) - time_desync(jmp) + ddrifting(jmp)/2) * 1e7) * 1e-7;
                dt_ph = dt_ph - t_offset;
                dt_pr = dt_pr - t_offset;

                ph = bsxfun(@minus, ph, dt_ph .* 299792458);
                pr = bsxfun(@minus, pr, dt_pr .* 299792458);
            else
                [ph_dj, dt_ph_dj] = Core_PP.remDtJumps(ph);
                [pr_dj, dt_pr_dj] = Core_PP.remDtJumps(pr);
                ddt_pr = Core_Utils.diffAndPred(dt_pr_dj);
                jmp_reset = find(abs(ddt_pr) > 1e-7); % points where the clock is reset
                if numel(jmp_reset) > 2
                    drifting = interp1(jmp_reset, dt_pr_dj(jmp_reset), (1 : numel(ddt_pr))', 'spline');
                else
                    drifting = 0;
                end
                dt_pr = dt_pr_dj - drifting;
                t_offset = mean(dt_pr);
                dt_ph = dt_ph_dj - drifting - t_offset;
                dt_pr = dt_pr - t_offset;

                ph = bsxfun(@minus, ph, dt_ph .* 299792458);
                pr = bsxfun(@minus, pr, dt_pr .* 299792458);
            end

            if any(dt_ph_dj)
                log.addMessage(log.indent('Correcting carrier phases jumps', 6),100);
            else
                log.addMessage(log.indent('Correcting carrier phases for a dt drift estimated from desync interpolation', 6),100);
            end
            if any(dt_pr_dj)
                log.addMessage(log.indent('Correcting pseudo-ranges jumps', 6),100);
            else
                log.addMessage(log.indent('Correcting pseudo-ranges for a dt drift estimated from desync interpolation', 6),100);
            end
        end

        function [obs_out] = remShortArcs(obs_in, min_len)
            % Removal of observation arcs shorter than given threshold.
            % SYNTAX:
            %   [obs_out] = Core_PP.remShortArcs(obs_in, threshold_gap);
            %
            % INPUT:
            %   obs_in = input observation matrix (num_sat x epochs)
            %   min_len = threshold on the number of epochs for an arc to be
            %                   considered "short"
            %
            % OUTPUT:
            %   obs_out = input observation matrix (num_sat x epochs)

            if (nargin < 2)
                min_len = 10;
            end

            obs_out = obs_in';

            for s = 1 : size(obs_in,1)

                % find intervals of zeros
                lim = getFlagsLimits(obs_in(s,:) == 0 | isnan(obs_in(s, :)));

                % find the intervals of good obs_in
                lim = [[1; lim(:,2)+1] [lim(:,1)-1; size(obs_in,2)]];
                lim = [lim (lim(:,2) - lim(:,1) +1)];

                lim = lim(lim(:,3) <= min_len, :);

                for i = 1 : size(lim,1)
                    obs_out(lim(i,1) : lim(i,2), s) = 0;
                end
            end

            obs_out = obs_out';
        end
    end

    % ==================================================================================================================================================
    %  PRIVATE FUNCTIONS called by public calls
    % ==================================================================================================================================================

    methods (Access = public)
    end

    % ==================================================================================================================================================
    %  STATIC FUNCTIONS used as utilities
    % ==================================================================================================================================================
    methods (Static, Access = public)
        function [data_out] = highPass(data, spline_base, cut_off, filt_border)
            % Compute an high pass filter on a dataset
            % Additionally apply a spline interpolation before thee fft transformation
            %
            % SYNTAX:
            %   [filtered_data] = highPass(data, spline_base, cut_off, filt_border);
            %
            % INPUT:
            %   data            dataset to filter
            %   spline_base     spline_base for spline filtering (if 0 do not do it)
            %   cut_off         cut-off frequency for the filter
            %   filt_border     filter dimension (larger -> smoother)
            %
            % OUTPUT:
            %   filtered_data   data reduced by splines and spectral filter

            if nargin < 4
                filt_border = 1e100;
            end

            spline_base = min(round(numel(data)/2), spline_base);

            %%% DEBUG figure(133); clf; plot(data); hold on;

            % if the number of data is even interpolate one observation
            is_even = false;
            if (cut_off > 0) && (mod(numel(data),2) == 0) % is even
                is_even = true;
                x = (1 : numel(data))';
                data = [interp1(x(~isnan(data)), data(~isnan(data)), 0, 'pchip', 'extrap'); data];
            end
            data = data - (0 : numel(data)-1)' * (median(data(max(1,length(data)-4):end)) - median(data(1:min(5,length(data))))) / numel(data);

            n_obs = length(data);

            % filter low frequencies with splines to limit border effects:
            if (spline_base > 0)
                x = (1 : n_obs)';
                ws = 5;
                margin = 2 * round(ws/2);
                xi = (1 - margin : n_obs + margin)';

                dt_lf = [median(data(1:min(5,length(data)))) * ones(margin,1); data(~isnan(data)); ones(margin,1) * median(data(max(1,length(data)-4):end))];
                %%% DEBUG figure(133); plot(xi - 1, dt_lf);
                dt_lf = splinerMat([],medfilt_mat(dt_lf, ws), spline_base);

                data = data - dt_lf(margin + 1 : end - margin);
            end
            %%% DEBUG figure(133); plot(data((1 + is_even) : end));

            if (cut_off > 0)
                % replicate signal to limit the border effects
                data = [flipud(data); data; flipud(data)];

                %%% DEBUG figure(110); clf; plot(zero2nan(data));

                % Take fft
                fftx = fftshift(fft(data));

                % Calculate the number of unique points (find the zero frequency - mean of the signal)
                f0_id = ceil((length(data) + 1) / 2);

                % This is an evenly spaced frequency vector with f0_id points.
                f = (0:f0_id-1) ./ length(data);

                % FFT is symmetric, throw away second half
                fftx_half = fftx(f0_id : end);

                % design the filter
                w = zeros(size(fftx_half));
                f_cut = cut_off;
                id_cut = find(f > f_cut, 1, 'first');
                filt_border = min(2*id_cut - 1, filt_border);
                w(1 : (id_cut - ((filt_border - 1) / 2))) = 0;
                w((id_cut + ((filt_border - 1) / 2)) : end) = 1;
                w((id_cut - ((filt_border - 1) / 2)) : (id_cut + ((filt_border - 1) / 2))) = 1 - (1 ./ (1+exp(-10:20/(filt_border-1):10)))'; % use a sigmoid filter for the design

                % compute power spectra - for debug
                %%% DEBUG psd = abs(fftx_half)/length(data)*sqrt(length(data)); % Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
                %%% DEBUG figure(125); clf; loglog(f, psd)

                % apply the filter
                fftx_half = fftx_half .* w;

                % compute power spectra of the filtered data - for debug
                %%% DEBUG psd_f = abs(fftx_half) / length(data) * sqrt(length(data)); % Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x

                fftx = ifftshift([flipud(conj(fftx_half)); fftx_half(2:end)]);

                %%% DEBUG figure(125); hold on; loglog(f, psd_f)
                data_out = real(ifft(fftx));
                data_out = data_out(n_obs + ((1 + is_even) : n_obs));
                %%% DEBUG figure(133); plot(data_out); title('data');
            else
                data_out = data;
            end
        end
    end

end

author: goGPS-Project


text:
function show_sample(samples)
% a simple visualization tool using isosurface to show each 3D sample.

n = size(samples,1);
for i = 1 : n
    the_sample = squeeze(samples(i,:,:,:,:));

    figure;
    p = patch(isosurface(the_sample,0.05));
    set(p,'FaceColor','red','EdgeColor','none');
    daspect([1,1,1])
    view(3); axis tight
    camlight
    lighting gouraud;
    axis off;
    set(gcf,'Color','white');
    set(gca,'position',[0,0,1,1],'units','normalized');
    axis tight;
    %title(i);
    pause;
    close(gcf);
end

author: zhirongw


text:
function [W] = spm_voice_iQ(Q)
% Discrete cosine transform of formant coefficients
% FORMAT [W] = spm_voice_iQ(Q)
%
% Q     - log formant frequencies
% G(1)  - log formant (pitch) Tu
% G(2)  - log timing  (pitch) Tv

%
% W     - (Nu x Nv) log formant coeficients (weights)
%
%   Nu  - number of formant coefficients
%   Nv  - number of timing  coefficients
%   Tu  - log formant (pitch)
%   Tv  - log timing  (pitch)
%
% This  auxiliary routine scales and transforms log formant coefficients
% using a pair of discrete cosine transforms with logarithmic scaling.
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_voice_iQ.m 7750 2019-12-05 17:54:29Z spm $

% defaults and (logarithmic) scaling
%--------------------------------------------------------------------------
try, Nu = VOX.Nu; catch, Nu  = 32;         end    % DCT order   (formants)
try, Nv = VOX.Nv; catch, Nv  = 8;          end    % DCT order   (interval)
try, Tu = VOX.Tu; catch, Tu  = 4;          end    % log scaling (formants)
try, Tv = VOX.Tv; catch, Tv  = 1;          end    % log scaling (interval)

% sizes
%--------------------------------------------------------------------------
[Ni,ni] = size(Q);

%  inverse transform
%--------------------------------------------------------------------------
U  = spm_voice_dct(Ni,Nu,Tu);                % DCT over formants
V  = spm_voice_dct(ni,Nv,Tv);                % DCT over intervals
W  = U\Q/V';                                 % coeficients
W  = W/std(W(:));                            % normalise

author: spm


text:
function handle_wrapcallbacktrycatch = ncorr_util_wrapcallbacktrycatch(handle_callback,handle_figure)
% This function is a wrapper that wraps the callback in a try-catch
% statement. If an error is thrown after the figure handle is closed, then
% disregard the error; if not, then rethrow it. The reason for this function
% is that closing a figure always interrupts a function, which can cause it
% to throw an error. Make sure to only throw the error if the figure which
% the function is attached to is still open.
%
% Inputs -----------------------------------------------------------------%
%   handle_callback - function handle;
%   handle_figure - figure handle;
%
% Outputs ----------------------------------------------------------------%
%   handle_wrapcallbacktrycatch - function handle;

    handle_wrapcallbacktrycatch = @wrapcallbacktrycatch;

    function wrapcallbacktrycatch(varargin)
        try
            handle_callback(varargin{:});
        catch err
            if (ishandle(handle_figure))
                rethrow(err);
            end
        end
    end
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
% CVX_CLEAR   Clears all active CVX data.
%    CVX_CLEAR clears the current CVX model in progress. This is useful if, for
%    example, you have made an error typing in your model and wish to start
%    over. Typing this before entering another CVX_BEGIN again avoids the
%    warning message that occurs if CVX_BEGIN detects a model in progress.

prob = evalin( 'caller', 'cvx_problem', '[]' );
if isa( prob, 'cvxprob' ),
	evalin( 'caller', 'pop( cvx_problem, ''clear'' )' );
end
% cvx_clearpath( 1 );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function xyz_io_test02 ( )

%*****************************************************************************80
%
%% XYZ_IO_TEST02 tests XYZ_HEADER_READ, XYZ_DATA_READ.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2009
%
%  Author:
%
%    John Burkardt
%
  file_name = 'xyz_io_prb_02.xyz';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'XYZ_IO_TEST02\n' );
  fprintf ( 1, '  XYZ_HEADER_READ reads the header of an XYZ file.\n' );
  fprintf ( 1, '  XYZ_DATA_READ reads the data of an XYZ file.\n' );

  xyz_write_test ( file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  XYZ_WRITE_TEST created some data.\n' );

  point_num = xyz_header_read ( file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  XYZ_HEADER_READ has read the header.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points = %d\n', point_num );

  xyz = xyz_data_read ( file_name, point_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  XYZ_DATA_READ has read the data.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample data:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : 11
    i = floor ( ( ( 11 - k     ) * 1 ...
                + (      k - 1 ) * point_num ) ...
                / ( 11     - 1 ) );
    fprintf ( 1, '  %4d  %10f  %10f  %10f\n', i, xyz(1:3,i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% MHSWTRANS
%
%  MCMC Metropolis-Hastings transition function that
%  utilizes the Swendsen-Wang proposal distribution.
%  A - The current joint assignment.  This should be
%      updated to be the next assignment
%  G - The network
%  F - List of all factors
%  variant - a number (1 or 2) indicating the variant of Swendsen-Wang to use.  In variant 1,
%            all the q_{i,j}'s are equal
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function A = MHSWTrans(A, G, F, variant)

%%%%%%%%%%%%%% Get Proposal %%%%%%%%%%%%%%
% Prune edges from q_list if the nodes don't have the same current value
q_list = G.q_list;
q_keep_indx = find(A(q_list(:, 1)) == A(q_list(:, 2)));
q_list = q_list(q_keep_indx, :);
% Select edges at random based on q_list
selected_edges_q_list_indx = find(q_list(:, 3) > rand(size(q_list,1), 1));
selected_edges = q_list(selected_edges_q_list_indx, 1:2);
% Compute connected components over selected edges
SelEdgeMat = sparse([selected_edges(:,1)'; selected_edges(:,2)'],...
                    [selected_edges(:,2)'; selected_edges(:,1)'],...
                    1, length(G.names), length(G.names));

[var2comp, cc_sizes] = scomponents(SelEdgeMat);
num_cc = length(cc_sizes);

% Select a connected component (the book calls this Y)
selected_cc = ceil(rand() * num_cc);
selected_vars = find(var2comp == selected_cc);
% Check that the dimensions are all the same and they have the same current assignment
assert(length(unique(G.card(selected_vars))) == 1);
assert(length(unique(A(selected_vars))) == 1);

% Pick a new label via sampling
old_value = A(selected_vars(1));
d = G.card(selected_vars(1));
LogR = zeros(1, d);
if variant == 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % YOUR CODE HERE
    % Specify the log of the distribution (LogR) from
    % which a new label for Y is selected for variant 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    LogR = ones(1,d).*log(1/d);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif variant == 2
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % YOUR CODE HERE
    % Specify the log of the distribution (LogR) from
    % which a new label for Y is selected for variant 2
    %
    % We suggest you read through the preceding code
    % before implementing this, one of the generated
    % data structures may be useful in implementing this section
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    LogR = BlockLogDistribution(selected_vars,G,F,A);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
    disp('WARNING: Unrecognized Swendsen-Wang Variant');
end

% Sample the new value from the distribution R
new_value = randsample(d, 1, true, exp(LogR));
A_prop = A;
A_prop(selected_vars) = new_value;

% Get the log-ratio of the probability of picking the connected component Y given A_prop over A
log_QY_ratio = 0.0;
for i = 1:size(G.q_list, 1)  % Iterate through *all* edges, not just the ones we selected earlier
    u = G.q_list(i, 1);
    v = G.q_list(i, 2);
    if length(intersect([u, v], selected_vars)) == 1  % the edge is from Y to outside-Y
        if A(u) == old_value && A(v) == old_value
            log_QY_ratio = log_QY_ratio - log(1 - G.q_list(i, 3));
        end
        if A_prop(u) == new_value && A_prop(v) == new_value
            log_QY_ratio = log_QY_ratio + log(1 - G.q_list(i, 3));
        end
    end
end

p_acceptance = 0.0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
% Compute acceptance probability
%
% Read through the preceding code to understand
% how to find the previous and proposed assignments
% of variables, as well as some ratios used in computing
% the acceptance probabilitiy.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lnum = LogProbOfJointAssignment(F,A_prop);
lden = LogProbOfJointAssignment(F,A);

p_acceptance = min(1,exp(lnum-lden+log_QY_ratio+LogR(old_value)-LogR(new_value)));
%p_acceptance = min(1,exp(log_QY_ratio));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Accept or reject proposal
if rand() < p_acceptance
    %disp('Accepted');
    A = A_prop;
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function f_x = ParFor5(in1)
%PARFOR5
%    F_X = PARFOR5(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    20-Sep-2019 09:35:35

u = in1(:,1);
ux = in1(:,4);
uxx = in1(:,5);
uxxx = in1(:,6);
f_x = ((u.*7.604652919865344e15+ux.*1.439296854415953e19+uxx.*4.056295991410688e16-uxxx.*3.942428646215516e17-u.*ux.*1.138844935991067e19-u.*uxxx.*1.804892689519346e18+2.199453572808704e15).*-1.0)./(u.*4.656689920376832e15-1.151078656075366e17);

author: dynamicslab

repo: SINDy-PI


text:
% Hoer's code (P. O. Hoyer. Non-negative Matrix Factorization with sparseness constraints.
% Journal of Machine Learning Research  5:1457-1469, 2004.) is modified

function [W, H, objhistory, objhistory_v] = nmf( V, rdim, showflag, ttt )
%% INPUT
% V: data matrix
% rdim : matrix factorization rank
% showflag: if it is 1. than it plots the change of reconstruction error
% for each iteration
% ttt: maximum number of iterations allowed

%% OUTPUT
% W: Mixing matrix (V=WH)
% H: Encoding matrix (V=WH)
% objhistory: change of total reconstruction error
% objhistory_v: reconstruction error for each individual sample

% Check that we have non-negative data
if min(V(:))<0, error('Negative values in data!'); end

% Dimensions
vdim = size(V,1);
samples = size(V,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% initialization %%%
fname2=['initials/Initials' num2str(vdim) 'x' num2str(rdim) '.txt'];
fid=fopen(fname2,'r');
if (fid==-1)
    ssbinitial(vdim,rdim);
end
fidW2 = fopen(fname2,'r');
W = fscanf(fidW2,'%f', [vdim inf]);
fclose(fidW2);

clear fname2

fname2=['initials/Initials' num2str(rdim) 'x' num2str(samples) '.txt'];
fid=fopen(fname2,'r');
if (fid==-1)
    ssbinitial(rdim,samples);
end
fidW2 = fopen(fname2,'r');
H = fscanf(fidW2,'%f', [rdim inf]);
fclose(fidW2);

clear fname2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (showflag==1)
    objhistory = ((sum(sum((V-W*H).^2))))/(vdim*samples);
    figure(2); clf;
    drawnow;
end

% Start iteration
iter = 0;
while 1,

    if (iter==ttt)
        break
    end
    iter = iter+1;
    % Compute new W and H (Lee and Seung; NIPS*2000)
    H = H.*(W'*V)./(W'*W*H + 1e-9);
    W = W.*(V*H')./(W*H*H' + 1e-9);

%     norms = sqrt(sum(W.^2));
%     H = H.*(norms'*ones(1,samples));
%     W = W./(ones(vdim,1)*norms);

    newobj = ((sum(sum((V-W*H).^2))))/(vdim*samples);
    newobj_v = ((sum(sum((V(:,end)-W*H(:,end)).^2))))/(vdim);
    if iter==1
        objhistory = [newobj];
        objhistory_v = [newobj_v];
    else
        objhistory = [objhistory newobj];
        objhistory_v = [objhistory_v newobj_v];
    end

    if(showflag==1)
        if (iter>1)
            plot(objhistory(2:end));
        end
        drawnow;
    end
end

author: andrewssobral

repo: lrslibrary


text:
function varargout = process_zscore( varargin )
% PROCESS_ZSCORE: Compute Z-Score for a matrix A (normalization respect to a baseline).
%
% DESCRIPTION:  For each channel:
%     1) Compute mean m and variance v for baseline
%     2) For each time sample, subtract m and divide by v

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2010-2015

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Z-score normalization [DEPRECATED]';
    sProcess.FileTag     = 'zscore';
    sProcess.Category    = 'Filter';
    sProcess.SubGroup    = 'Standardize';
    sProcess.Index       = 0;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Z-score';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Default values for some options
    sProcess.isSourceAbsolute = 0;
    sProcess.processDim       = 1;    % Process channel by channel

    % Definition of the options
    sProcess.options.description.Comment = ['For each signal in input:
' ...
                                            '1) Compute mean m and variance v for the baseline
' ...
                                            '2) For each time sample, subtract m and divide by v
' ...
                                            'Z = (Data - m) / v

'];
    sProcess.options.description.Type    = 'label';
    % === Baseline time window
    sProcess.options.baseline.Comment = 'Baseline:';
    sProcess.options.baseline.Type    = 'baseline';
    sProcess.options.baseline.Value   = [];
    % === Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    sProcess.options.sensortypes.InputTypes = {'data'};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    % Get time window
    if isfield(sProcess.options, 'baseline') && isfield(sProcess.options.baseline, 'Value') && iscell(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value{1})
        Time = sProcess.options.baseline.Value{1};
    else
        Time = [];
    end
    % Add time window to the comment
    if isempty(Time)
        Comment = 'Z-score normalization: [All file]';
    elseif any(abs(Time) > 2)
        Comment = sprintf('Z-score normalization: [%1.3fs,%1.3fs]', Time(1), Time(2));
    else
        Comment = sprintf('Z-score normalization: [%dms,%dms]', round(Time(1)*1000), round(Time(2)*1000));
    end
end

%% ===== RUN =====
function sInput = Run(sProcess, sInput) %#ok<DEFNU>
    % Get options
    if isfield(sProcess.options, 'baseline') && isfield(sProcess.options.baseline, 'Value') && iscell(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value{1})
        BaselineBounds = sProcess.options.baseline.Value{1};
    else
        BaselineBounds = [];
    end
    % Get baseline indices
    if ~isempty(BaselineBounds)
        iBaseline = panel_time('GetTimeIndices', sInput.TimeVector, sProcess.options.baseline.Value{1});
        if isempty(iBaseline)
            bst_report('Error', sProcess, [], 'Invalid baseline definition.');
            sInput = [];
            return;
        end
    % Get all file
    else
        iBaseline = 1:size(sInput.A,2);
    end
    % Compute zscore
    sInput.A = Compute(sInput.A, iBaseline);
    % Change DataType
    if ~strcmpi(sInput.FileType, 'timefreq')
        sInput.DataType = 'zscore';
    end
    % Default colormap
    if strcmpi(sInput.FileType, 'results')
        sInput.ColormapType = 'stat1';
        sInput.Function     = 'zscore';
    else
        sInput.ColormapType = 'stat2';
    end
    % Do not keep the Std field in the output
    if isfield(sInput, 'Std') && ~isempty(sInput.Std)
        sInput.Std = [];
    end
end

%% ===== COMPUTE =====
function A = Compute(A, iBaseline)
    disp('BST> process_zscore.m is deprecated, use "Standardize > Baseline normalization" instead.');
    % Calculate mean and standard deviation
    [meanBaseline, stdBaseline] = ComputeStat(A(:, iBaseline,:));
    % Compute zscore
    A = bst_bsxfun(@minus, A, meanBaseline);
    A = bst_bsxfun(@rdivide, A, stdBaseline);
end

function [meanBaseline, stdBaseline] = ComputeStat(A)
    % Compute baseline statistics
    stdBaseline  = std(A, 0, 2);
    meanBaseline = mean(A, 2);
    % Remove null variance values
    stdBaseline(stdBaseline == 0) = 1e-12;
end

author: brainstorm-tools


text:
% This is used to set the VIF algorithms that to be run
%
% Author:Xingchen Zhang, Ping Ye, Gang Xiao
% Contact: xingchen.zhang@imperial.ac.uk

function methods=configMethods

 methodVIFB={struct('name','ADF'),...
        struct('name','CBF'),...
        struct('name','CNN'),...
        struct('name','DLF'),...
        struct('name','FPDE'),...
        struct('name','GFCE'),...
        struct('name','GFF'),...
        struct('name','GTF'),...
        struct('name','HMSD_GF'),...
        struct('name','Hybrid_MSD'),...
        struct('name','IFEVIP'),...
        struct('name','LatLRR'),...
        struct('name','MGFF'),...
        struct('name','MST_SR'),...
        struct('name','MSVD'),...
        struct('name','NSCT_SR'),...
        struct('name','ResNet'),...
        struct('name','RP_SR'),...
        struct('name','TIF'),...
        struct('name','VSMWLS'),...
     };

methods = [methodVIFB];

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [mesh labels labelsmap] = loadMesh(filestr)
% loads a mesh from filestr, supports obj or off

file = fopen( strtrim( filestr ), 'rt');
if file == -1
    warning(['Could not open mesh file: ' filestr]);
    mesh = [];
    return;
end
mesh.filename = strtrim( filestr );

if strcmp( filestr(end-3:end), '.off')
    line = strtrim(fgetl(file));
    skipline = 0;
    if strcmp(line, 'OFF')
        line = strtrim(fgetl(file));
        skipline = 2;
    else
        line = line(4:end);
        skipline = 1;
    end
    [token,line] = strtok(line);
    numverts = eval(token);
    [token,line] = strtok(line);
    numfaces = eval(token);
    mesh.V = zeros( 3, numverts, 'single' );
    mesh.F = zeros( 3, numfaces, 'single' );

    DATA = dlmread(filestr, ' ', skipline, 0);
    DATA = DATA(1:numverts+numfaces, :);
    mesh.V(1:3, 1:numverts) = DATA(1:numverts, 1:3)';
    mesh.F(1:3, 1:numfaces) = DATA(numverts+1:numverts+numfaces, 2:4)' + 1;
elseif strcmp( filestr(end-3:end), '.obj')
    mesh.V = zeros(3, 10^6, 'single');
    mesh.Nv = zeros(3, 10^6, 'single');
    mesh.F = zeros(3, 5*10^6, 'uint32');
    v = 0;
    f = 0;
    vn = 0;

    while(~feof(file))
        line_type = fscanf(file,'%c',2);
        switch line_type(1)
            case {'#', 'g'}
                fgets(file);
            case 'v'
                if line_type(2) == 'n'
                    vn = vn + 1;
                    normal  = fscanf(file, '%f%f%f');
                    mesh.Nv(:, vn) = normal;
                elseif isspace( line_type(2) )
                    v = v + 1;
                    point = fscanf(file, '%f%f%f');
                    mesh.V(:, v) = point;
                else
                    fgets(file);
                end
            case 'f'
                f = f + 1;
                face = fscanf(file, '%u%u%u');
                mesh.F(:, f) = face;
            otherwise
                if feof(file)
                    break;
                end
                if isspace(line_type(1))
                    fseek(file, -1, 'cof');
                    continue;
                end
                fprintf('last string read: %c%c %s', line_type(1), line_type(2), fgets(file));
                fclose(file);
                error('only triangular obj meshes are supported with vertices, normals and faces.');
        end
    end
    mesh.V = mesh.V(:, 1:v);
    mesh.F = mesh.F(:, 1:f);
    mesh.Nv = mesh.Nv(:, 1:v);
end

fclose(file);

end

author: suhangpro

repo: mvcnn

save_path: github-repos/MATLAB/suhangpro-mvcnn


text:
function erg = rsqu(q, r)
% Description:
%   rsqu(r, q) computes the r2-value for two one-dimensional distributions
%   given by the vectors q and r
%
% Example code:
%  [erg] = rsqu(q , r)
%
% Input:
%   q: Data structrue (ex) Epoched data structure
%   r: Data structrue (ex) Epoched data structure
%
% Options:
%
% Return:
%    erg:  r2-values
%
% See also:
%
% Reference:
%           G. Schalk, D.J. McFarland, T. Hinterberger, N. Birbaumer, and
%           J. R. Wolpaw,"BCI2000: A General-Purpose Brain-Computer
%           Interface (BCI) System, IEEE Transactions on Biomedical
%           Engineering, Vol. 51, No. 6, 2004, pp.1034-1043.

%         We used BCI2000 open source toolbox code related in r-square value (rsqu.m)
%
%
% Ji Hoon, Jeong
% jh_jeong@korea.ac.kr

%%

q=double(q);
r=double(r);

sum1=sum(q);
sum2=sum(r);
n1=length(q);
n2=length(r);
sumsqu1=sum(q.*q);
sumsqu2=sum(r.*r);

G=((sum1+sum2)^2)/(n1+n2);

erg=(sum1^2/n1+sum2^2/n2-G)/(sumsqu1+sumsqu2-G);
end

author: PatternRecognition

repo: OpenBMI


text:
%  Computes the hessian numerically
%
%  ::
%
%    [H,issue]=numerical(fh,xbest,hessian_type)
%
%  Args:
%
%     - **fh** [char|function handle]: one-dimensional objective function
%     - **xbest** [vector]: point at which the hessian has to be computed
%     - **hessian_type** [{'fd'}|'opg']: type of hessian computed : finite
%       differences or outer-product-gradient
%
%  Returns:
%     :
%
%     - **H** [d x d matrix]: hessian
%     - **issue** [''|char]: description of any problem encountered during the
%       calculation of the hessian.
%
%

author: jmaih


text:
function beta = lars(X, y, method, stop, useGram, Gram, Cardi, bSparse, trace)
% This function is provided at
% http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3897
% I have made some small modifications  -- Deng Cai, Feb/2008

% LARS  The LARS algorithm for performing LAR or LASSO.
%    BETA = LARS(X, Y) performs least angle regression on the variables in
%    X to approximate the response Y. Variables X are assumed to be
%    normalized (zero mean, unit length), the response Y is assumed to be
%    centered.
%    BETA = LARS(X, Y, METHOD), where METHOD is either 'LARS' or 'LASSO'
%    determines whether least angle regression or lasso regression should
%    be performed.
%    BETA = LARS(X, Y, METHOD, STOP) with nonzero STOP will perform least
%    angle or lasso regression with early stopping. If STOP is negative,
%    STOP is an integer that determines the desired number of variables. If
%    STOP is positive, it corresponds to an upper bound on the L1-norm of
%    the BETA coefficients.
%    BETA = LARS(X, Y, METHOD, STOP, USEGRAM) specifies whether the Gram
%    matrix X'X should be calculated (USEGRAM = 1) or not (USEGRAM = 0).
%    Calculation of the Gram matrix is suitable for low-dimensional
%    problems. By default, the Gram matrix is calculated.
%    BETA = LARS(X, Y, METHOD, STOP, USEGRAM, GRAM) makes it possible to
%    supply a pre-computed Gram matrix. Set USEGRAM to 1 to enable. If no
%    Gram matrix is available, exclude argument or set GRAM = [].
%    BETA = LARS(X, Y, METHOD, STOP, USEGRAM, GRAM, TRACE) with nonzero
%    TRACE will print the adding and subtracting of variables as all
%    LARS/lasso solutions are found.
%    Returns BETA where each row contains the predictor coefficients of
%    one iteration. A suitable row is chosen using e.g. cross-validation,
%    possibly including interpolation to achieve sub-iteration accuracy.
%
% Author: Karl Skoglund, IMM, DTU, kas@imm.dtu.dk
% Reference: 'Least Angle Regression' by Bradley Efron et al, 2003.

%% Input checking
% Set default values.
if nargin < 9
    trace = 0;
end
if nargin < 8
    bSparse = 1;
end
if nargin < 7
    Cardi = [];
end
if nargin < 6
    Gram = [];
end
if nargin < 5
    useGram = 0;
end
if nargin < 4
    stop = 0;
end
if nargin < 3
    method = 'lasso';
end
if strcmpi(method, 'lasso')
    lasso = 1;
else
    lasso = 0;
end

if isempty(X)
    error('The code has been updated. Please input the X');
end

%% LARS variable setup
[n p] = size(X);
% nvars = min(n-1,p); %
nvars = p; %

maxk = 512*nvars; % Maximum number of iterations

if isempty(Cardi)
    if stop == 0
        if bSparse
            beta = sparse(p,2*nvars);
        else
            beta = zeros(p,2*nvars);
        end
    elseif stop < 0
        if bSparse
            beta = sparse(p,2*round(-stop));
        else
            beta = zeros(p,2*round(-stop));
        end
    else
        if bSparse
            beta = sparse(p,100);
        else
            beta = zeros(p,100);
        end
    end
else
    Cardi = unique(Cardi);
    Cardi(Cardi>nvars) = [];
    stop = -max(Cardi);
    if bSparse
        beta = sparse(p,length(Cardi));
    else
        beta = zeros(p,length(Cardi));
    end
    betak = zeros(p,1);
end

mu = zeros(n, 1); % current "position" as LARS travels towards lsq solution
I = 1:p; % inactive set
A = []; % active set

% Calculate Gram matrix if necessary
if isempty(Gram) && useGram
    error('The code has been updated. Please input the Gram');
%     clear Gram;
%     global Gram;
    %   Gram = X'*X; % Precomputation of the Gram matrix. Fast but memory consuming.
end

if ~useGram
    R = []; % Cholesky factorization R'R = X'X where R is upper triangular
end

lassocond = 0; % LASSO condition boolean
stopcond = 0; % Early stopping condition boolean
k = 0; % Iteration count
vars = 0; % Current number of variables

if trace
    disp(sprintf('Step\tAdded\tDropped\t\tActive set size'));
end

% TimeLoop = zeros(2*nvars,1);
tmpT = cputime;

%% LARS main loop
while vars < nvars && ~stopcond && k < maxk
    k = k + 1;
    c = X'*(y - mu);
    [C j] = max(abs(c(I)));
    j = I(j);

    if ~lassocond % if a variable has been dropped, do one iteration with this configuration (don't add new one right away)
        if ~useGram
            diag_k = X(:,j)'*X(:,j); % diagonal element k in X'X matrix
            if isempty(R)
                R = sqrt(diag_k);
            else
                col_k = X(:,j)'*X(:,A); % elements of column k in X'X matrix
                R_k = R'\col_k'; % R'R_k = (X'X)_k, solve for R_k
                R_kk = sqrt(diag_k - R_k'*R_k); % norm(x'x) = norm(R'*R), find last element by exclusion
                R = [R R_k; [zeros(1,size(R,2)) R_kk]]; % update R
            end
        end
        A = [A j];
        I(I == j) = [];
        vars = vars + 1;
        if trace
            disp(sprintf('%d\t\t%d\t\t\t\t\t%d', k, j, vars));
        end
    end

    s = sign(c(A)); % get the signs of the correlations

    if useGram
        if vars <= 200
            R = chol(Gram(A,A));
        elseif lassocond
            if (rJ <= 200) & vars <= 1000
                R = chol(Gram(A,A));
            else
                R(:,rJ) = []; % remove column j
                tmpn = size(R,2);
                for tmpk = rJ:tmpn
                    tmpp = tmpk:tmpk+1;
                    [G,R(tmpp,tmpk)] = planerot(R(tmpp,tmpk)); % remove extra element in column
                    if tmpk < tmpn
                        R(tmpp,tmpk+1:tmpn) = G*R(tmpp,tmpk+1:tmpn); % adjust rest of row
                    end
                end
                R(end,:) = []; % remove zero'ed out row
            end
        else
            R_k = R'\Gram(A(1:end-1),j);
            R_kk = sqrt(Gram(j,j)-R_k'*R_k);
            R = [R R_k; [zeros(1,size(R,2)) R_kk]]; % update R
        end
        GA1 = R\(R'\s);
        AA = 1/sqrt(sum(GA1.*s));
        w = AA*GA1;
    else
        GA1 = R\(R'\s);
        AA = 1/sqrt(sum(GA1.*s));
        w = AA*GA1;
    end
    u = X(:,A)*w; % equiangular direction (unit vector)

    if vars == nvars % if all variables active, go all the way to the lsq solution
        gamma = C/AA;
    else
        a = X'*u; % correlation between each variable and eqiangular vector
        temp = [(C - c(I))./(AA - a(I)); (C + c(I))./(AA + a(I))];
        gamma = min([temp(temp > 0); C/AA]);
    end

    % LASSO modification
    if lasso
        lassocond = 0;
        if isempty(Cardi)
            temp = -beta(A,k)./w;
        else
            temp = -betak(A)./w;
        end
        [gamma_tilde] = min([temp(temp > 0); gamma]);
        j = find(temp == gamma_tilde);
        if gamma_tilde < gamma,
            gamma = gamma_tilde;
            lassocond = 1;
        end
    end

    mu = mu + gamma*u;
    if isempty(Cardi)
        if size(beta,2) < k+1
            if bSparse
                beta = [beta sparse(p,size(beta,1))];
            else
                beta = [beta zeros(p,size(beta,1))];
            end
        end
        beta(A,k+1) = beta(A,k) + gamma*w;
    else
        tmpbetak = betak(A) + gamma*w;
        betak = zeros(p,1);
        betak(A) = tmpbetak;
        idx = find(Cardi==vars);
        if ~isempty(idx)
            beta(:,idx) = betak;
        end
    end

    % Early stopping at specified bound on L1 norm of beta
    if isempty(Cardi)
        if stop > 0
            t2 = sum(abs(beta(:,k+1)));
            if t2 >= stop
                t1 = sum(abs(beta(:,k)));
                s = (stop - t1)/(t2 - t1); % interpolation factor 0 < s < 1
                beta(:,k+1) = beta(:,k) + s*(beta(:,k+1) - beta(:,k));
                stopcond = 1;
            end
        end
    end

    % If LASSO condition satisfied, drop variable from active set
    if lassocond == 1
        if ~useGram
            R(:,j) = []; % remove column j
            tmpn = size(R,2);
            for tmpk = j:tmpn
                tmpp = tmpk:tmpk+1;
                [G,R(tmpp,tmpk)] = planerot(R(tmpp,tmpk)); % remove extra element in column
                if tmpk < tmpn
                    R(tmpp,tmpk+1:tmpn) = G*R(tmpp,tmpk+1:tmpn); % adjust rest of row
                end
            end
            R(end,:) = []; % remove zero'ed out row
        end
        rJ = j;
        I = [I A(j)];
        A(j) = [];
        vars = vars - 1;
        if trace
            disp(sprintf('%d\t\t\t\t%d\t\t\t%d', k, j, vars));
        end
    end

    % Early stopping at specified number of variables
    if stop < 0
        stopcond = vars >= -stop;
    end

%     TimeLoop(k) = cputime - tmpT;
%     tmpT = cputime;

%     if vars < 1000
%         if mod(vars,500) == 0
%             tmpT = cputime - tmpT;
%             disp(['LARS: ',num2str(vars),' features selected. Time: ',num2str(tmpT)]);
%             tmpT = cputime;
%         end
%     elseif vars < 2000
%         if mod(vars,200) == 0
%             tmpT = cputime - tmpT;
%             disp(['LARS: ',num2str(vars),' features selected. Time: ',num2str(tmpT)]);
%             tmpT = cputime;
%         end
%     elseif vars < 3000
%         if mod(vars,100) == 0
%             tmpT = cputime - tmpT;
%             disp(['LARS: ',num2str(vars),' features selected. Time: ',num2str(tmpT)]);
%             tmpT = cputime;
%         end
%     else
%         if mod(vars,50) == 0
%             tmpT = cputime - tmpT;
%             disp(['LARS: ',num2str(vars),' features selected. Time: ',num2str(tmpT)]);
%             tmpT = cputime;
%         end
%     end
end

if isempty(Cardi)
    % trim beta
    if size(beta,2) > k+1
        beta(:,k+2:end) = [];
    end
end

if k == maxk
    disp('LARS warning: Forced exit. Maximum number of iteration reached.');
end

%% To do
%
% There is a modification that turns least angle regression into stagewise
% (epsilon) regression. This has not been implemented.

author: ZJULearning

repo: MatlabFunc


text:
function []=olsss(Y,X,n,m,p,Bhat,stringdates1,decimaldates1,endo,pref)

% function []=olsss(Y,X,n,m,p,Bhat,stringdates1,decimaldates1,endo,datapath)
% calculates and displays the steady-state for the OLS VAR
% inputs:  - matrix 'Y': matrix of regressands for the VAR model (defined in 1.1.8)
%          - matrix 'X': matrix of regressors for the VAR model (defined in 1.1.8)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - matrix 'Bhat': OLS VAR coefficients, in non vectorised form (defined in 1.1.9)
%          - cell 'stringdates1': date strings for the sample period
%          - vector 'decimaldates1': dates converted into decimal values, for the sample period
%          - cell 'endo': list of endogenous variables of the model
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: none

% this function estimates the steady-state, using (a.7.6)

% first compute the steady-state values

% recover the coefficient matrices A1,...,Ap and C, as defined in (1.1.2)
% first, calculate B and take its transpose BT
BT=Bhat';
% estimate the summation term I-A1-...-Ap in
summation=eye(n);
for jj=1:p
    summation=summation-BT(:,(jj-1)*n+1:jj*n);
end
% recover C
C=BT(:,end-m+1:end);
% now calculate the product of the inverse of the summation with C
product=summation\C;
% keep only the exogenous regressor part of X
X_exo=X(:,end-m+1:end)';
% compute the steady-state values
ssvalues=product*X_exo;

% create steady-state figure
if pref.plot
    sstate=figure('Tag','BEARresults');
    set(sstate,'Color',[0.9 0.9 0.9]);
    set(sstate,'name','steady-state');
    ncolumns=ceil(n^0.5);
    nrows=ceil(n/ncolumns);
    for ii=1:n
        subplot(nrows,ncolumns,ii);
        hold on
        ss=plot(decimaldates1,ssvalues(ii,:),'Color',[0.4 0.4 1],'LineWidth',2);
        actual=plot(decimaldates1,Y(:,ii),'Color',[0 0 0],'LineWidth',2);
        plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
        hold off
        minband=min(min(ssvalues(ii,:)),min(Y(:,ii)));
        maxband=max(max(ssvalues(ii,:)),max(Y(:,ii)));
        space=maxband-minband;
        Ymin=minband-0.2*space;
        Ymax=maxband+0.2*space;
        set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'YLim',[Ymin,Ymax],'FontName','Times New Roman');
        title(endo{ii,1},'FontName','Times New Roman','FontSize',10,'interpreter','none');
        if ii==1
            plotlegend=legend([ss,actual],'steady-state','actual');
            set(plotlegend,'FontName','Times New Roman');
        end
    end
end

% finally, save in excel

% create the cell that will be saved on excel
sscell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},size(stringdates1,1)+3,1);
% loop over variables (horizontal dimension)
for ii=1:n
    % create cell of steady-state record for variable ii
    temp=['steady-state and actual: ' endo{ii,1}];
    ss_i=[temp {''} {''};{''} {''} {''};{''} {'sample'} {'median'};stringdates1 num2cell(Y(:,ii)) num2cell(ssvalues(ii,:))'];
    sscell=[sscell ss_i vertspace];
end
% trim
sscell=sscell(:,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),sscell,'steady state','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [S, mu, sigma2] = standardize(M, mu, sigma2)
% function S = standardize(M, mu, sigma2)
% Make each column of M be zero mean, std 1.
% Thus each row is scaled separately.
%
% If mu, sigma2 are omitted, they are computed from M

M = double(M);
if nargin < 2
  mu = mean(M,2);
  sigma2 = std(M,0,2);
  sigma2 = sigma2 + eps*(sigma2==0);
end

[nrows ncols] = size(M);
S = M - repmat(mu(:), [1 ncols]);
S = S ./ repmat(sigma2, [1 ncols]);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function this = parse_input_parameter(this, param)
%PARSE_INPUT_PARAMETER parses and set parameters in the struct PARAM
%
%   This is a member function of the class 'hs_optical_flow'.

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-11-30 $
%   $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

if mod(length(param), 2) ~=0
    error('Parse_input_parameter: Input parameters must be given in pairs (name and value)!');
end;

i = 1;

while (i <= length(param))

    if ischar(param{i+1})
        param{i+1} = str2num(param{i+1});
    end;

    switch lower(param{i})

        case 'lambda'
            this.lambda         = param{i+1};

        case 'pyramid_levels'
            this.pyramid_levels  = param{i+1};

        case 'pyramid_spacing'
            this.pyramid_spacing = param{i+1};

        case 'max_warping_iters'
            this.max_warping_iters       = param{i+1};

        case 'median_filter_size'
            this.median_filter_size = [param{i+1} param{i+1}];

        case 'image_filter'
            if strcmpi(param{i+1}, 'gaussian')
                this.image_filter = fspecial('gaussian', [15 15], 0.4);
            elseif strcmpi(param{i+1}, 'LoG')
                this.image_filter = - fspecial('gaussian', [15 15], 1);
                this.image_filter(8,8) = 1 + this.image_filter(8,8);
            elseif strcmpi(param{i+1}, 'DoG')
                this.image_filter =  fspecial('gaussian', [15 15], 0.4) - fspecial('gaussian', [15 15], 1);
            else
                this.image_filter = [];
            end;
    end;

    i = i+2;
end;

author: kristinbranson

repo: JAABA


text:
function [fal,dx0] = aliasing_frequency(x0,conf)
%ALIASING_FREQUENCY aliasing frequency for the given secondary sources
%
%   Usage: [fal,dx0] = aliasing_frequency([x0],conf)
%
%   Input parameters:
%       x0      - secondary sources / m
%       conf    - configuration struct (see SFS_config)
%
%   Output parameters:
%       fal     - aliasing frequency / Hz
%       dx0     - mean distance between secondary sources / m
%
%   ALIASING_FREQUENCY(x0,conf) returns the aliasing frequency for the given
%   secondary sources. First the mean distance dx0 between the secondary sources
%   is calculated, afterwards the aliasing frequency is calculated after Spors
%   (2009) as fal = c/(2*dx0). If no secondary sources x0 are provided, they are
%   first calculated by calling secondary_source_positions().
%   For a calculation that includes the dependency on the listener position have
%   a look at Start (1997).
%
%   See also: sound_field_mono_wfs, secondary_source_positions,
%       secondary_source_distance
%
%   References:
%       Spors and Ahrens (2009) - "Spatial sampling artifacts of wave field
%       synthesis for the reproduction of virtual point sources", 126th
%       Convention of the Audio Engineering Society, Paper 7744,
%       http://www.aes.org/e-lib/browse.cfm?elib=14940
%
%       Start (1997) - "Direct Sound Enhancement by Wave Field Synthesis",
%       PhD thesis, TU Delft,
%       http://resolver.tudelft.nl/uuid:c80d5b58-67d3-4d84-9e73-390cd30bde0d

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameters ====================================
nargmin = 1;
nargmax = 2;
narginchk(nargmin,nargmax);
if nargin<nargmax
    conf = x0;
    x0 = [];
end
isargstruct(conf);

%% ===== Configuration ==================================================
c = conf.c;

%% ===== Computation =====================================================
% If no explicit secondary source distribution is given, calculate one
if isempty(x0)
    x0 = secondary_source_positions(conf);
end
% Get average distance between secondary sources
dx0 = secondary_source_distance(x0);
% Calculate aliasing frequency
fal = c/(2*dx0);

author: sfstoolbox

repo: sfs-matlab


text:
function ind = fundamental_region(q,cs,ss)

if isempty(q), ind = []; return; end

c = {};

% eliminiate 3 fold symmetry axis of cubic symmetries
switch cs.LaueName

  case   {'m-3m','m-3'}

    c{end+1}.v = vector3d([1 1 1 1 -1 -1 -1 -1],[1 1 -1 -1 1 1 -1 -1],[1 -1 1 -1 1 -1 1 -1]);
    c{end}.h = sqrt(3)/3;

    if strcmp(cs.LaueName,'m-3m')
      c{end+1}.v = vector3d([1 -1 0 0 0 0],[0 0 1 -1 0 0],[0 0 0 0 1 -1]);
      c{end}.h = sqrt(2)-1;
    end
end

switch ss.LaueName
  case 'mmm'
   c{end+1}.v = vector3d([-1 0],[0 -1],[0 0]);
   c{end}.h = 0;
end

% find rotation not part of the fundamental region
rodrigues = Rodrigues(q); clear q;
ind = false(length(rodrigues),1);
for i = 1:length(c)
  for j = 1:length(c{i}.v)
    p = dot(rodrigues,1/norm(c{i}.v(j)) * c{i}.v(j));
    ind = ind | (p(:)>c{i}.h);
  end
end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [contact, fric_coef, geometry] = RightSole(robot)

    param = sys.GetExtraParams();

    r_foot_frame = robot.Joints(getJointIndices(robot, 'r_leg_akx'));
    contact = CoordinateFrame(...
        'Name','RightSole',...
        'Reference',r_foot_frame,...
        'Offset',[0, 0, param.hf],...
        'R',[0,0,0]... % z-axis is the normal axis, so no rotation required
        );

    fric_coef.mu = param.mu;
    fric_coef.gamma = param.gamma;

    geometry.la = param.wf/2;
    geometry.lb = param.wf/2;
    geometry.La = param.lt;
    geometry.Lb = param.lh;
end

author: ayonga

repo: frost-dev

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


text:
%%
function Metrics = ...
    performCalcu(datasetStruct,algStructArray)

evaluateSal = @(sMap,thresholds,gtMap) thresholdBased_HR_FR(sMap,thresholds,gtMap);

thresholds = [255:-1:0]/255;         % thresholding image from 0 - 255
GTfiles=dir([datasetStruct.GTdir '/*.png']); % gather the GT files in png
GTfiles = [GTfiles; dir([datasetStruct.GTdir '/*.jpg'])]; % gather the GT files in jpg
GTfiles = [GTfiles; dir([datasetStruct.GTdir '/*.bmp'])]; % gather the GT files in bmp

numOfFiles = size(GTfiles,1); % totalImgNum*1
numOfAlgos = length(algStructArray);  % num of algorithmos to be evaluated

[hitRate, falseAlarm] = deal(zeros(numOfFiles,length(thresholds),numOfAlgos));%

%% initial

[IoU, TPR, FPR, Pre, Recall] = deal(zeros(numOfFiles,length(thresholds),numOfAlgos));

[mean_Fmeasure] = deal(zeros(numOfFiles,3,numOfAlgos)); % totalImgNum*3*numOfAlgos [precision, recall, fmeasure]

[MAE] = deal(zeros(numOfFiles,1,numOfAlgos)); % totalImgNum*1*numOfAlgos

[S_measure] = deal(zeros(numOfFiles,1,numOfAlgos)); % totalImgNum*1*numOfAlgos

%Iterate over images
hwait  = waitbar(0,['Processing ', num2str(numOfFiles),' images on ',datasetStruct.datasetName,': 0%']);

for imIndx=1:numOfFiles
    % read the ground truth files
    [~,base_name,ext] = fileparts(GTfiles(imIndx).name);

    gtMap = im2double(imread([datasetStruct.GTdir base_name ext]));
    gtSize = size(gtMap);
    if (length(gtSize) == 3)
        gtMap = rgb2gray(gtMap);
        gtSize(3)= [];
    end

    gtMap = logical(gtMap>=0.1); % get binary mask
    totalNum = numOfFiles* ones(numOfAlgos,1); %3*1   ( [totalImgNum;totalImgNum;totalImgNum;totalImgNum] )
    for algIdx = 1:numOfAlgos % compute over each algorithms iteratively
        sMap = readSaliencyMap(algStructArray{algIdx},[base_name],gtSize);% read and resize saliency maps
        if sum(sum(sMap)) == 0
            totalNum(algIdx) = totalNum(algIdx) - 1;
        end

        [IoU(imIndx,:,algIdx), TPR(imIndx,:,algIdx), FPR(imIndx,:,algIdx), Pre(imIndx,:,algIdx), Recall(imIndx,:,algIdx), ...        % calculate IoU, TPR, FPR, Pre, Recall, hitRate, falseAlarm
         hitRate(imIndx,:,algIdx), falseAlarm(imIndx,:,algIdx)] ...
            = evaluateSal(sMap,thresholds,gtMap);

        [mean_Fmeasure(imIndx,:,algIdx)] = Fmeasure_calu(sMap,gtMap,gtSize); % calculate F-measure using 2 x MeanValueOfSalMap [precision,recall,fmeasure]

        MAE(imIndx,:,algIdx) = CalMAE(sMap, gtMap); % clculate MAE

        S_measure(imIndx,:,algIdx) = StructureMeasure(sMap, gtMap); % calculate S-measure

    end
    pct     = sprintf('%5.1f',imIndx/numOfFiles*100);
    waitbar(imIndx/numOfFiles,hwait,['Processing ', num2str(numOfFiles),...
        ' images on ',datasetStruct.datasetName,': ',pct,'%']);
end %End of image loop

close(hwait);
%Average across images -
mmHitRate = permute( sum(hitRate,1),[2 3 1] );
mmFalseAlarm = permute( sum(falseAlarm,1),[2 3 1]);

nanind = isnan(Pre);
Pre(nanind) = 0;
mmPre = permute( sum(Pre,1),[2 3 1]);
mmRecall = permute( sum(Recall,1),[2 3 1]);
mmTPR = permute( sum(TPR,1),[2 3 1]);
mmFPR = permute( sum(FPR,1),[2 3 1]);
mmIoU = permute( sum(IoU,1),[2 3 1]);

mmmean_Fmeasure = permute( sum(mean_Fmeasure,1),[2 3 1]);
mmMAE = zeros(1,numOfAlgos);
mmS_measure = zeros(1,numOfAlgos);

for j=1:numOfAlgos
    nanind = isnan(MAE(:,1,j));
    mmmMAE = MAE(:,1,j);
    mmmMAE(nanind) = [];

    nanind = isnan(S_measure(:,1,j));
    mmmS_measure = S_measure(:,1,j);
    mmmS_measure(nanind) = [];

    mmHitRate(:,j) = mmHitRate(:,j)./totalNum(j);
    mmFalseAlarm(:,j) = mmFalseAlarm(:,j)./totalNum(j);
    mmPre(:,j) = mmPre(:,j)./totalNum(j);
    mmRecall(:,j) = mmRecall(:,j)./totalNum(j);
    mmTPR(:,j) = mmTPR(:,j)./totalNum(j);
    mmFPR(:,j) = mmFPR(:,j)./totalNum(j);
    mmIoU(:,j) = mmIoU(:,j)./totalNum(j);

    mmmean_Fmeasure(:,j) = mmmean_Fmeasure(:,j)./totalNum(j);
    mmMAE(:,j) = mean(mmmMAE,1);
    mmS_measure(:,j) = mean(mmmS_measure,1);
end
Metrics.HitRate = mmHitRate;
Metrics.FalseAlarm = mmFalseAlarm;
Metrics.Pre = mmPre;
Metrics.Recall = mmRecall;
Metrics.TPR = mmTPR;
Metrics.FPR = mmFPR;
Metrics.IoU = mmIoU;

Metrics.mean_Fmeasure = mmmean_Fmeasure;
Metrics.MAE = mmMAE;
Metrics.S_measure = mmS_measure;
Metrics.AUC = nan(1,size(Metrics.FalseAlarm,2));
for algIdx=1:numOfAlgos
        Metrics.AUC(algIdx) = trapz(Metrics.FalseAlarm(:,algIdx),Metrics.HitRate(:,algIdx));
end

Metrics.mean_Fmeasure_image = mean_Fmeasure;
Metrics.MAE_image = MAE;
Metrics.S_measure_image = S_measure;
Metrics.TPR_image = TPR;
Metrics.FPR_image = FPR;
Metrics.Pre_image = Pre;
Metrics.Recall_image = Recall;

Metrics.Fmeasure_Curve                  = (1+0.3).*Metrics.Pre.*Metrics.Recall./...
                                            (0.3*Metrics.Pre+Metrics.Recall);
[Metrics.MaxFmeasure, ind]              = max(Metrics.Fmeasure_Curve, [], 1);
Metrics.IoU_at_maxF                     = Metrics.IoU(ind);
[Metrics.maxIoU, ~]                     = max(Metrics.IoU, [], 1);
Metrics.meanIoU                         = mean(Metrics.IoU, 1);
end

%% Read Saliency Maps

% Read and resize saliency map
function sMap = readSaliencyMap(algStruct,base_name,gtSize)
if strcmp(algStruct.dir(end-3:end),'DSR')
    algStruct.ext='jpg';
elseif strcmp(algStruct.dir(end-12:end-10),'DSR')
    algStruct.ext='jpg';
end
file_name = fullfile(algStruct.dir,[algStruct.prefix base_name algStruct.postfix '.' algStruct.ext]);
sMap = imresize(im2double(imread(file_name)),gtSize(1:2));
if (size(sMap,3)==3)
    sMap = rgb2gray(sMap);
end
sMap(sMap<0)=0;
maxnum = max(sMap(:));
if maxnum==0
    sMap = zeros(gtSize(1:2));
else
    sMap = sMap./maxnum;
end

end

%% Calculate the scores
function [IoU, TPR, FPR, Pre, Recall, hitRate, falseAlarm] ...
    = thresholdBased_HR_FR(sMap,thresholds,gtMap) %calculate the threshold based scores
    gtPxlNum = sum(gtMap(:));
    totalNum = size(gtMap,1)*size(gtMap,2);
    targetHist = histc(sMap(gtMap), fliplr(thresholds));
    nontargetHist = histc(sMap(~gtMap), fliplr(thresholds));

    targetHist = flipud(targetHist);
    nontargetHist = flipud(nontargetHist);

    targetHist = cumsum( targetHist );
    nontargetHist = cumsum( nontargetHist );

    Pre = targetHist ./ (targetHist + nontargetHist);
    Recall = targetHist / gtPxlNum;
    TPR = Recall;
    FPR = nontargetHist / (totalNum-gtPxlNum);
    IoU = targetHist ./ (gtPxlNum + nontargetHist);
    hitRate = targetHist / gtPxlNum;
    falseAlarm = 1 - (repmat(totalNum-gtPxlNum,size(nontargetHist)) - nontargetHist)/(totalNum-gtPxlNum);
end
%%

author: ArcherFMY


text:
function [tfr,t,f] = tfrcw(x,t,N,g,h,sigma,trace);
%TFRCW	Choi-Williams time-frequency distribution.
%	[TFR,T,F]=TFRCW(X,T,N,G,H,SIGMA,TRACE) computes the Choi-Williams
%	distribution of a discrete-time signal X, or the
%	cross Choi-Williams representation between two signals.
%
%	X     : signal if auto-CW, or [X1,X2] if cross-CW.
%	T     : time instant(s)          (default : 1:length(X)).
%	N     : number of frequency bins (default : length(X)).
%	G     : time smoothing window, G(0) being forced to 1.
%	                                 (default : Hamming(N/10)).
%	H     : frequency smoothing window, H(0) being forced to 1.
%	                                 (default : Hamming(N/4)).
%	SIGMA : kernel width             (default : 1)
%	TRACE : if nonzero, the progression of the algorithm is shown
%                                        (default : 0).
%	TFR   : time-frequency representation. When called without
%         output arguments, TFRCW runs TFRQVIEW.
%	F     : vector of normalized frequencies.
%
%	Example:
%	 sig=fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4);
%	 g=tftb_window(9,'Kaiser'); h=tftb_window(27,'Kaiser');
%	 t=1:128; tfrcw(sig,t,128,g,h,3.6,1);
%
%	See also all the time-frequency representations listed in
%	 the file CONTENTS (TFR*)

%	F. Auger, May-August 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least 1 parameter required');
end;
[xrow,xcol] = size(x);
if (xcol==0)|(xcol>2),
 error('X must have one or two columns');
end

if (nargin <= 2),
 N=xrow;
elseif (N<0),
 error('N must be greater than zero');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
end;

hlength=floor(N/4); hlength=hlength+1-rem(hlength,2);
glength=floor(N/10);glength=glength+1-rem(glength,2);

if (nargin == 1),
 t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); sigma = 1.0; trace = 0;
elseif (nargin == 2)|(nargin == 3),
 g = tftb_window(glength); h = tftb_window(hlength); sigma = 1.0; trace = 0;
elseif (nargin == 4),
 h = tftb_window(hlength); sigma = 1.0; trace = 0;
elseif (nargin == 5),
 sigma = 1.0; trace = 0;
elseif (nargin == 6),
 trace = 0;
end;

[trow,tcol] = size(t);
if (trow~=1),
 error('T must only have one row');
end;

[grow,gcol]=size(g); Lg=(grow-1)/2;
if (gcol~=1)|(rem(grow,2)==0),
 error('G must be a smoothing window with odd length');
end;

[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1);
if (hcol~=1)|(rem(hrow,2)==0),
  error('H must be a smoothing window with odd length');
end;

if (sigma<=0.0),
 error('SIGMA must be strictly positive');
end;

normfac = 16.0*pi/sigma; spreadfac = 16.0/sigma;
taumax = min([round(N/2),Lh]); tau = 1:taumax; points = -Lg:Lg;
CWKer = exp(-kron( points.' .^2, 1.0 ./ (spreadfac*tau.^2)));
CWKer = diag(g) * CWKer;

tfr= zeros (N,tcol) ;
if trace, disp('Choi-Williams distribution'); end;
for icol=1:tcol,
 ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]);
 if trace, disprog(icol,tcol,10); end;
 tfr(1,icol)= x(ti,1) .* conj(x(ti,xcol));

 for tau=1:taumax,
  points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]);
  g2 = CWKer(Lg+1+points,tau); g2=g2/sum(g2);
  R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)));
  tfr(  1+tau,icol)=h(Lh+tau+1)*R;
  R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)));
  tfr(N+1-tau,icol)=h(Lh-tau+1)*R;
 end;

 tau=round(N/2);
 if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh),
  points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]);
  g2 = CWKer(Lg+1+points,tau); g2=g2/sum(g2);
  tfr(tau+1,icol) = 0.5 * ...
   (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+...
    h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))));
 end;
end;

clear CWKer;

if trace, fprintf('\n'); end;
tfr= fft(tfr);
if (xcol==1), tfr=real(tfr); end ;

if (nargout==0),
 tfrqview(tfr,x,t,'tfrcw',g,h,sigma);
elseif (nargout==3),
 f=(0.5*(0:N-1)/N)';
end;

author: HeLiangHIT


text:
classdef Quadrature_Triangle < Quadrature

    methods
        function computeQuadrature(obj,order)
            computeQuadrature@Quadrature(obj,order);
            switch order
                case 'CONSTANT'
                    obj.ngaus = 1;
                    obj.weigp = 1/2;
                    obj.posgp = [1/3;1/3];

                case 'LINEAR'
                    obj.ngaus = 1;
                    obj.weigp = [1/2];
                    obj.posgp = [1/3;1/3];

                case 'QUADRATIC'
                    %obj.ngaus = 3;
                    %obj.weigp = [1/6;1/6;1/6];
                    %obj.posgp = [0,0.5;0.5,0;0.5,0.5]';
                    obj.ngaus = 3;
                    obj.weigp = [1/6;1/6;1/6];
                    obj.posgp = [2/3,1/6;1/6,1/6;1/6,2/3]';
                case 'CUBIC'
                    obj.ngaus = 4;
                    obj.weigp = [-27/96;25/96;25/96;25/96];
                    obj.posgp = [1/3,1/3;1/5,1/5;3/5,1/5;1/5,3/5]';

                case 'QUADRATIC2'
                    obj.ngaus = 2;
                    obj.weigp = [1/6;1/6;1/6];
                    obj.posgp = [1/6,1/6;2/3,1/6;1/6,2/3]';

                case 'QUADRATICMASS'
                    obj.ngaus = 3;
                    obj.weigp = [1/6;1/6;1/6];
                    obj.posgp = [2/3,1/6,1/6;1/6,2/3,1/6];
                otherwise
                    disp('Quadrature not implemented for triangle elements')
            end
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function shortcuts = GetShortcuts(category, labels)
%GETSHORTCUTS Returns details of existing toolbar shortcuts.
%
% SHORTCUTS = GETSHORTCUTS(CATEGORY) returns the labels, callback, icon
% path and whether or not the shortcut is editable, for each existing
% toolbar shortcut in the category; those in the 'Toolbar Shortcuts' by
% default.
%
% SHORTCUTS = GETSHORTCUTS(CATEGORY, LABELS) works as
% above, but only returns shortcuts with labels matching LABELS.
%
% Examples:
%
% shortcuts = GetShortcuts();
%
% See also: GetShortcutNames, AddShortcut

% $Author: rcotton $	$Date: 2010/08/23 14:35:09 $	$Revision: 1.1 $
% Copyright: Health and Safety Laboratory 2010

shortcutUtils = GetShortcutUtils();

SetDefaultValue(1, 'category', DefaultToolbarCategoryName());

filter = nargin >= 2 && ~isempty(labels);

shortcutArray = ...
   shortcutUtils.getShortcutsByCategory(category).toArray();

shortcuts = struct('label', {}, 'callback', {}, ...
   'icon', {}, 'category', {}, 'editable', {});

for i = 1:shortcutArray.length
   label = char(shortcutArray(i).getLabel());
   if ~filter || strcontains(label, labels)
      shortcuts(end + 1).label = label;           %#ok<AGROW>
      shortcuts(end).callback = ...
         char(shortcutArray(i).getCallback());
      shortcuts(end).icon = ...
         char(shortcutArray(i).getIconPath());
      shortcuts(end).category = category;
      shortcuts(end).editable = ...
         shortcutArray(i).isEditable;
   end
end

end

author: Sable

repo: mcbench-benchmarks


text:
function [data,units] = compute_angle_smallest_wing(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = -trx(fly).wing_anglel;
  idx = trx(fly).wing_arear_mm < trx(fly).wing_areal_mm;
  data{i}(idx) = trx(fly).wing_angler(idx);

end
units = parseunits('rad');


text:
function test_bug2773

% WALLTIME 00:20:00
% MEM 2gb
% DEPENDENCY ft_dipolefitting ft_movieplotER ft_prepare_sourcemodel ft_prepare_layout

orig = load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2773.mat'));
vol  = orig.cfg.vol;
elec = orig.cfg.elec;

cfg = [];
cfg.elec = elec;
cfg.channel =    {
  '1'
  '2'
  '3'
  '4'
  '5'
  '6'
  '7'
  '8'
  '9'
  '10'
  '11'
  '12'
  '13'
  '14'
  '15'
  '16'
  '17'
  '18'
  '19'
  '20'
  '21'
  '22'
  '23'
  '24'
  '25'
  '26'
  '27'
  '28'
  '29'
  '30'
  '31'
  '32'
  '33'
  '34'
  '35'
  '36'
  '37'
  '38'
  '39'
  '40'
  '41'
  '42'
  '43'
  '44'
  '45'
  '46'
  '47'
  '48'
  '49'
  '50'
  '51'
  '52'
  '53'
  '54'
  '55'
  '56'
  '57'
  '58'
  '59'
  '60'
  '61'
  };

layout = ft_prepare_layout(cfg); % it is not a very nice layout, but will do for the test script

figure; ft_plot_layout(layout);
figure; ft_plot_sens(elec);
figure; ft_plot_headmodel(vol);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function s = y_to_s_islamic ( y )

%*****************************************************************************80
%
%% Y_TO_S_ISLAMIC writes an Islamic year into a string.
%
%  Format:
%
%    YearNumber AH
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, string S, a representation of the year.
%
  s = sprintf ( 'AH %d', y );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function F = tag(F,text)
% TAG
%
% Sets the tag on a constraint.

F = flatten(F);
if nargin == 1
    F = F.clauses{1}.handle;
else
for i = 1:length(F.clauses)
    if ~isequal(F.clauses{i}.handle,'placeholder')
        F.clauses{i}.handle = text;
    end
end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP