removing Matlab Toolbox dependencies (e.g., quantile); changing defau…

…lt option of batch processing to no parallelization

CODE/COMMON/calc_quantile.m

@@ -0,0 +1,61 @@
+function [quant] = calc_quantile(M, p)
+% This function calculates the p quantile for a matrix (each column) or 
+% vector and avoids the Statistics and Machine Learning Toolbox.
+%
+% Uses the same algorithm as Matlabm, described here:
+% https://de.mathworks.com/help/matlab/ref/quantile.html (accessed May 5, 2023)
+% For an n-element vector A, quantile computes quantiles by using a sorting-based algorithm:
+% The sorted elements in A are taken as the (0.5/n), (1.5/n), ..., ([n – 0.5]/n) quantiles. For example:
+% For a data vector of five elements such as {6, 3, 2, 10, 1}, the sorted elements {1, 2, 3, 6, 10} respectively correspond to the 0.1, 0.3, 0.5, 0.7, and 0.9 quantiles.
+% For a data vector of six elements such as {6, 3, 2, 10, 8, 1}, the sorted elements {1, 2, 3, 6, 8, 10} respectively correspond to the (0.5/6), (1.5/6), (2.5/6), (3.5/6), (4.5/6), and (5.5/6) quantiles.
+% quantile uses Linear Interpolation to compute quantiles for probabilities between (0.5/n) and ([n – 0.5]/n).
+% For the quantiles corresponding to the probabilities outside that range, quantile assigns the minimum or maximum values of the elements in A.
+% quantile treats NaNs as missing values and removes them.
+% 
+% INPUT:
+%   M       matrix or vector with data
+%   p       [0; 1], defines the quantile
+% OUTPUT:
+%	quant   vector or double, p quantile of input matrix or vector
+%
+% Revision:
+%   ...
+%
+% This function belongs to raPPPid, Copyright (c) 2023, M.F. Glaner
+% *************************************************************************
+
+
+% make sure that input vector is a column vector
+[n, m] = size(M);
+if (n == 1 || m == 1) && m > n
+    M = M';     % transpose row to column vector
+end
+
+% initialize ouput variable
+[~, m] = size(M);
+quant = NaN(1,m);
+
+for i = 1:m
+    vec = M(:,i);           % get current column of data
+    vec = sortrows(vec);    % sort values
+    vec(isnan(vec)) = [];  	% ignore NaNs
+    n = numel(vec);
+
+    % sorted elements in A are taken as the (0.5/n), (1.5/n), ..., ([n – 0.5]/n) quantiles
+    vec_quantiles = (0.5 : 1 : (n-0.5)) / n; 
+
+    % for quantiles outside that range, assign the minimum or maximum value
+    if p < 0.5/n
+        q = min(vec);
+    elseif p > (n-0.5)/n
+        q = max(vec);
+    else
+        % inside range -> use linear interpolation
+        q = lininterp1(vec_quantiles, vec, p);
+    end
+
+     % save quantile of current column
+     if ~isempty(q)
+         quant(i) = q;          
+     end
+end

CODE/FIXING/CheckSatellitesFixable.m

@@ -56,7 +56,7 @@
 excl_prn = settings.AMBFIX.exclude_sats_fixing;
 if ~isempty(excl_prn)
     [~,ind] = ismember(Epoch.sats, excl_prn);
-    Epoch.fixable(boolean(ind),:) = false;
+    Epoch.fixable(logical(ind),:) = false;
 end
 
 
@@ -188,13 +188,13 @@
 
 
 
-function boolean = frequency_convert(boolean, settings)
+function bool = frequency_convert(bool, settings)
 % this function checks which ambiguities can not be fixed depending on the
 % processed PPP model and number of frequencies
 if settings.INPUT.proc_freqs == 1           % e.g. 2-frequency IF LC
     % only one frequency is used in the fixing process, therefore exclude
     % if any of the observations is unfixable
-    boolean = any(boolean,2);  
+    bool = any(bool,2);  
 end
 
 % ||| extend for other PPP-AR models

CODE/OBSERVATIONS/cycleSlip_Doppler.m

@@ -154,7 +154,7 @@
 % L_diff_ = abs(L_now - L_pred_);
 % new_cs_found_ = L_diff_ > thresh;
 % % take the right prediction (+ or -), this is different for each receiver
-% if nansum(L_diff_) < nansum(L_diff)
+% if nansum(L_diff_) < nansum(L_diff)			% nansum should not be used
 %     new_cs_found = new_cs_found_;
 %     L_diff = L_diff_;
 % end

CODE/VISUAL/GUI/GUI_PPP.m

@@ -50,8 +50,7 @@ function GUI_PPP_OpeningFcn(hObject, eventdata, handles, varargin)  %#ok<*INUSL>
 axis off
 
 % Set Copyright and Version in the lower right
-set(handles.text_version, 'String', ['Version 1.5 ', char(169), ' TUW 2023']);
-set(handles.checkbox_realtime, 'Visible', 'Off');
+set(handles.text_version, 'String', ['Version 1.6 ', char(169), ' TUW 2023']);
 
 % load default filter settings for selected filter
 handles = LoadDefaultFilterSettings(handles);

CODE/VISUAL/SinglePlotting.m

@@ -178,7 +178,6 @@
 % -+-+-+-+- Figure: Wet Troposphere Plot -+-+-+-+-
 if settings.PLOT.wet_tropo && settings.TROPO.estimate_ZWD
     TropoPlot(hours, label_x_h, storeData, reset_h, obs.startdate, obs.station_long)
-%     vis_plotTroposphere(hours, label_x_h, storeData, reset_h);    % OLD
 end
 if STOP_CALC; return; end
 
@@ -333,7 +332,7 @@
 
 %     -+-+-+-+- Figures: Multipath Detection  -+-+-+-+-
 if settings.PLOT.mp
-    if (isfield(settings.OTHER, 'mp_detection') || settings.OTHER.mp_detection)
+    if (isfield(settings.OTHER, 'mp_detection') && settings.OTHER.mp_detection)
         C1_diff = zero2nan(storeData.mp_C1_diff_n);
         PlotObsDiff(epochs, C1_diff, label_x_epc, rgb, 'C1 difference', settings, satellites.obs, settings.OTHER.mp_thresh, settings.OTHER.mp_degree, '', false);
     else

CODE/VISUAL/StationResultPlot.m

@@ -126,20 +126,20 @@
         if PlotStruct.fixed
             [TTFF, TTCF] = prepTTFF(TTFF, TTCF, d.FIXED,  PlotStruct.thresh_2D, d.dT, d.E, d.N, i);
         end
-        % calculate station accuracy (0.95 quantile)
+        % calculate station accuracy
         d2D = sqrt(d.N.^2 + d.E.^2);
         d3D = sqrt(d.N.^2 + d.E.^2 + d.H.^2);
-        ACC_2D(i)  = quantile(d2D(:), .95);
-        ACC_3D(i)  = quantile(d3D(:), .95);
         % calculate mean station convergence
         CONV_2D_mean(i) = mean(conv_2D, 'omitnan');
         CONV_3D_mean(i) = mean(conv_3D, 'omitnan');
         % calculate median station convergence
         CONV_2D_medi(i) = median(conv_2D, 'omitnan');
         CONV_3D_medi(i) = median(conv_3D, 'omitnan');        
-        % calculate quantiles of ZTD
-        ZTD_68(i) = quantile(abs(d.ZTD(:)), .68);
-        ZTD_95(i) = quantile(abs(d.ZTD(:)), .95);
+        % calculate quantiles
+        ACC_2D(i) = calc_quantile(d2D(:), .95);
+        ACC_3D(i) = calc_quantile(d3D(:), .95);
+        ZTD_68(i) = calc_quantile(abs(d.ZTD(:)), .68);
+        ZTD_95(i) = calc_quantile(abs(d.ZTD(:)), .95);
 
     end         % end of loop over files of current label
 

CODE/VISUAL/ThreeCoordinatesPlot.m

@@ -21,8 +21,13 @@ function ThreeCoordinatesPlot(interval, seconds, dN, dE, dH, resets, str_xAxis,
 
 
 seconds = seconds(1:length(dN));    % if last epochs of processing delivered no results
-no_eps = numel(seconds);       % number of epochs;
-RMS_dN = rms(dN(dN~=0), 'omitnan'); RMS_dE = rms(dE(dE~=0), 'omitnan'); RMS_dH = rms(dH(dH~=0), 'omitnan');
+
+% RMS_dN = rms(dN(dN~=0), 'omitnan'); RMS_dE = rms(dE(dE~=0), 'omitnan'); RMS_dH = rms(dH(dH~=0), 'omitnan');
+% calculate rms with own implementation to avoid ToolBox Dependency
+RMS_dN = calculate_rms(dN);
+RMS_dE = calculate_rms(dE);
+RMS_dH = calculate_rms(dH);
+
 isnan_dN = isnan(dN); isnan_dE = isnan(dE); isnan_dH = isnan(dH);
 dN(isnan_dN) = 0;     dE(isnan_dE) = 0;     dH(isnan_dH) = 0;
 
@@ -123,3 +128,11 @@ function ThreeCoordinatesPlot(interval, seconds, dN, dE, dH, resets, str_xAxis,
     end
 end
 end
+
+
+
+function rms_calc = calculate_rms(vec)
+% calculate rms
+vec = vec(vec~=0 & ~isnan(vec));        % ignore zeros and NaNs
+rms_calc = sqrt(mean(vec.*vec));
+end

CODE/VISUAL/calcPerformanceIndicators.m

@@ -92,7 +92,8 @@
 ZTD_50 = median(ZTD, 'omitnan') * 100;        % [cm]
 
 % 68% quantile of the ZTD difference
-ZTD_68 = quantile(ZTD, .68) * 100;    % [cm]
+ZTD_68 = calc_quantile(ZTD, .68) * 100;    % [cm]
+
 
 
 %% Output to command window

CODE/VISUAL/prepMultiPlot.m

@@ -77,19 +77,20 @@
 % find quantiles of position differences
 if n ~= 1
     % calculate 68% quantile
-    dN_68   = quantile(dN_all , 0.68);       % quantile() is slow
-    dE_68   = quantile(dE_all , 0.68);
-    dH_68   = quantile(dH_all , 0.68);
-    d2D_68  = quantile(d2D_all, 0.68);
-    d3D_68  = quantile(d3D_all, 0.68);
-    dZTD_68 = quantile(dZTD_all, 0.68);
+    dN_68   = calc_quantile(dN_all , 0.68);       
+    dE_68   = calc_quantile(dE_all , 0.68);
+    dH_68   = calc_quantile(dH_all , 0.68);
+    d2D_68  = calc_quantile(d2D_all, 0.68);
+    d3D_68  = calc_quantile(d3D_all, 0.68);
+    dZTD_68 = calc_quantile(dZTD_all, 0.68);
     % calculate 95% quantile
-    dN_95   = quantile(dN_all , 0.95);
-    dE_95   = quantile(dE_all , 0.95);
-    dH_95   = quantile(dH_all , 0.95);
-    d2D_95  = quantile(d2D_all, 0.95);
-    d3D_95  = quantile(d3D_all, 0.95);
-    dZTD_95 = quantile(dZTD_all, 0.95);
+    dN_95   = calc_quantile(dN_all , 0.95);
+    dE_95   = calc_quantile(dE_all , 0.95);
+    dH_95   = calc_quantile(dH_all , 0.95);
+    d2D_95  = calc_quantile(d2D_all, 0.95);
+    d3D_95  = calc_quantile(d3D_all, 0.95);
+    dZTD_95 = calc_quantile(dZTD_all, 0.95);
+
 else    % only one convergence period
     dN_68  = NaN(1,m);  dE_68 = NaN(1,m); dH_68  = NaN(1,m);
     d2D_68 = NaN(1,m); d3D_68 = NaN(1,m); dZTD_68 = NaN(1,m);

CODE/VISUAL/vis_BoxPlot.m

@@ -37,9 +37,13 @@
 
 % Plot
 figure('Name', txt, 'NumberTitle','off')
-boxplot(cell2mat(BOX), grp, 'Notch','on','Whisker',0, 'Symbol', 'o', 'OutlierSize',5, 'jitter', 1)
-% whiskers = 0, values outside the 75% and 25% quantile are outliers
-% jitter = 1, outliers are maximally distributed 
+try  % requires Statistics and Machine Learning Toolbox
+    boxplot(cell2mat(BOX), grp, 'Notch','on','Whisker',0, 'Symbol', 'o', 'OutlierSize',5, 'jitter', 1)
+    % whiskers = 0, values outside the 75% and 25% quantile are outliers
+    % jitter = 1, outliers are maximally distributed 
+catch
+    fprintf('BoxPlot will work in a future version of raPPPid.\n')
+end
 
 
 % Style

CODE/VISUAL/vis_plotSatConstellation.m

@@ -43,15 +43,15 @@ function vis_plotSatConstellation(hours, epochs, strXAxis, satellites, cutoff, i
 sv_GAL_cutoff = sv_GAL .* cutoff(epochs, idx_gal);
 sv_BDS_cutoff = sv_BDS .* cutoff(epochs, idx_bds);
 % number of observed satellites and number of satellites under cut-off
-noAllGPS = nansum(observ_sats(:, idx_gps),2);
-noAllGLO = nansum(observ_sats(:, idx_glo),2);
-noAllGAL = nansum(observ_sats(:, idx_gal),2);
-noAllBDS = nansum(observ_sats(:, idx_bds),2);
+noAllGPS = sum(observ_sats(:, idx_gps),2,'omitnan');
+noAllGLO = sum(observ_sats(:, idx_glo),2,'omitnan');
+noAllGAL = sum(observ_sats(:, idx_gal),2,'omitnan');
+noAllBDS = sum(observ_sats(:, idx_bds),2,'omitnan');
 noAllGNSS = noAllGPS+noAllGLO+noAllGAL+noAllBDS;
-notUsedGPS = nansum( cutoff(:, idx_gps),2 );
-notUsedGLO = nansum( cutoff(:, idx_glo),2 );
-notUsedGAL = nansum( cutoff(:, idx_gal),2 );   
-notUsedBDS = nansum( cutoff(:, idx_bds),2 );   
+notUsedGPS = sum( cutoff(:, idx_gps),2,'omitnan');
+notUsedGLO = sum( cutoff(:, idx_glo),2,'omitnan');
+notUsedGAL = sum( cutoff(:, idx_gal),2,'omitnan');   
+notUsedBDS = sum( cutoff(:, idx_bds),2,'omitnan');   
 UsedGPS = noAllGPS-notUsedGPS;
 UsedGLO = noAllGLO-notUsedGLO;
 UsedGAL = noAllGAL-notUsedGAL;
@@ -107,8 +107,8 @@ function vis_plotSatConstellation(hours, epochs, strXAxis, satellites, cutoff, i
     legend_txt{end+1} = 'GPS visible';
     legend_txt{end+1} = 'GPS used';
     if ~isGLO && ~isGAL && ~isBDS       % plot number of GPS L5 satellites
-        noGPS_L5 = nansum(observ_sats(:, DEF.PRNS_GPS_L5),2);
-        notUsed_L5 = nansum( cutoff(:, DEF.PRNS_GPS_L5),2 );
+        noGPS_L5 = sum(observ_sats(:, DEF.PRNS_GPS_L5),2,'omitnan');
+        notUsed_L5 = sum( cutoff(:, DEF.PRNS_GPS_L5),2,'omitnan');
         UsedGPS_L5 = noGPS_L5-notUsed_L5;
         plot(hours, noGPS_L5,    '.', 'Color', [1 .65 0.8])
         plot(hours, UsedGPS_L5,  '.', 'Color', [1 .65 0])