Move the isoutlier function/methods to a separate file. (#1710)

Move more files to "extras"

Author: joa-quim

src/GMT.jl

@@ -264,7 +264,6 @@ include("grdview.jl")
 include("grdvolume.jl")
 include("greenspline.jl")
 include("gridit.jl")
-include("hampel_outliers.jl")
 include("img_funs.jl")
 include("imgtiles.jl")
 include("imshow.jl")
@@ -323,10 +322,14 @@ include("xyz2grd.jl")
 include("utils.jl")
 include("utils_project.jl")
 include("choropleth_utils.jl")
-include("webmapserver.jl")
+include("extras/hampel_outliers.jl")
+include("extras/isoutlier.jl")
+include("extras/lowess.jl")
 include("extras/seismicity.jl")
 include("extras/okada.jl")
 include("extras/weather.jl")
+include("extras/webmapserver.jl")
+include("extras/whittaker.jl")
 include("seis/psmeca.jl")
 include("seis/gmtisf.jl")
 include("geodesy/psvelo.jl")
@@ -343,8 +346,6 @@ include("potential/gravfft.jl")
 include("potential/grdseamount.jl")
 include("spotter/grdrotater.jl")
 include("windbarbs/windbarbs.jl")
-include("lowess.jl")
-include("whittaker.jl")
 include("zscale.jl")
 include("drawing.jl")
 include("get_enums.jl")

src/common_options.jl

@@ -3789,7 +3789,7 @@ function _read_data(d::Dict, cmd::String, arg, opt_R::String="", is3D::Bool=fals
 			(t_col == 0) ? (s[1] *= cT; s[2] *= cT) : (t_col == 1) ? (s[3] *= cT; s[4] *= cT) : (t_col == 2 && length(s) == 6) ? (s[5] *= cT; s[6] *= cT) : error("-f option not cmpatible with -R")
 			opt_R = join(s, "/")				# Reconstruct the -R with the T's in the time axis
 			# Remove the --TIME_UNIT from the command due to a GMT bug. But do it only when using the default -B. VERY RISKY. If fails => GB'
-			#(contains(cmd, DEF_FIG_AXES_BAK) && ((TU = scan_opt(cmd, "--TIME_UNIT=")) != "")) && (cmd = replace(cmd, " --TIME_UNIT="*TU => ""))
+			(contains(cmd, DEF_FIG_AXES_BAK) && ((TU = scan_opt(cmd, "--TIME_UNIT=")) == "y")) && (cmd = replace(cmd, " --TIME_UNIT="*TU => ""))
 			#((TU = scan_opt(cmd, "--TIME_UNIT=")) != "") && (cmd = replace(cmd, " --TIME_UNIT="*TU => ""))
 			#((TU = scan_opt(cmd, "--TIME_EPOCH=")) != "") && (cmd = replace(cmd, " --TIME_EPOCH="*TU => ""))
 		end

src/hampel_outliers.jl src/extras/hampel_outliers.jlsrc/hampel_outliers.jl renamed to src/extras/hampel_outliers.jl

@@ -0,0 +1,106 @@
+"""
+	bv = isoutlier(x; critic=3.0, method::Symbol=:median, threshold=nothing, width=0.0) -> Union{BitVector, BitMatrix}
+
+Return a logical array whose elements are true when an outlier is detected in the corresponding element of `x`.
+
+If `x` is a matrix, and `width` is not used, then ``isoutlier`` operates on each column of `x` separately.
+By default, an outlier is a value that is more than three median absolute deviations (MAD) from the median.
+But that can be changed via the `critic` option.
+
+- `critic`: A number that sets the threshold for detecting outliers. The default is 3.0, which means that
+   a value is considered an outlier if it is more than 3 MAD from the median (when `method` is `:median`),
+   or more than 3 standard deviations (when `method` is `:mean`).
+- `method`: The method used to calculate the threshold for outliers. It can be one of the following:
+  - `:median`: Uses the median absolute deviation (MAD) method. Outliers are defined as elements more than
+   `critic` MAD from the median.
+  - `:mean`: Uses the mean and standard deviation method. Outliers are defined as elements more than `critic`
+   standard deviations from the mean. This method is faster but less robust than "median".
+  - `:quartiles`: Uses the interquartile range (IQR) method. Outliers are defined as elements more than 1.5
+   interquartile ranges above the upper quartile (75 percent) or below the lower quartile (25 percent).
+   This method is useful when the data in `x` is not normally distributed.
+- `threshold`: Is an alternative to the `method` option. It specifies the percentile thresholds, given as a
+   two-element array whose elements are in the interval [0, 100]. The first element indicates the lower percentile
+   threshold, and the second element indicates the upper percentile threshold. It can also be a single number
+   between (0, 100), which is interpreted as the percentage of end member points that are considered outliers.
+   For example, `threshold=1` means that the lower and upper thresholds are the 1th and 99th percentiles.
+- `width`: If this option is used (only when `x` is a Matrix or a GMTdataset) we detect local outliers using
+   a moving window method with window length ``width``. This length is given in the same units as the input
+   data stored in first column of `x`.
+
+### Example:
+```julia
+x = [57, 59, 60, 100, 59, 58, 57, 58, 300, 61, 62, 60, 62, 58, 57];
+findall(isoutlier(x))
+2-element Vector{Int64}:
+ 4
+ 9
+```
+
+```julia
+x = -50.0:50;	y = x / 50 .+ 3 .+ 0.25 * rand(length(x));
+y[[30,50,60]] = [4,-3,6];	# Add 3 outliers
+findall(isoutlier([x y], width=5))
+```
+"""
+function isoutlier(x::AbstractVector{<:Real}; critic=3.0, method::Symbol=:median, threshold=nothing)
+	method in [:median, :mean, :quartiles] || throw(ArgumentError("Unknown method: $method. Use :median, :mean or :quartiles"))
+	if (threshold !== nothing)
+		!(isa(threshold, VecOrMat{<:Real}) || isa(threshold, Real)) &&
+			throw(ArgumentError("Threshold: $threshold must be a scalar or a 2 elements array."))
+		if (isa(threshold, VecOrMat{<:Real}))
+			@assert !(length(threshold) == 2 && threshold[1] >= 0.0 && threshold[2] <= 100) "Threshold must be a 2 elements array in the [0 100] interval"
+			_thresh = [Float64(threshold[1])/100, Float64(threshold[2])/100]	# Make it [0 1]
+		else
+			@assert threshold > 0.0 && threshold < 100 "When threshold is a scalar, it must be in the ]0 100[ interval."
+			_thresh = [Float64(threshold)/100, Float64(100 - threshold)/100]	# Make it [0 1]
+		end
+		method = :threshold
+	end
+
+	bv = BitVector(undef, length(x))
+	if (method == :median)
+		_mad, _med = mad(x);	_mad *= critic
+		@inbounds for k in eachindex(x)
+			bv[k] = abs(x[k] - _med)  > _mad
+		end
+	elseif (method == :mean)
+		_mean, _std = nanmean(x), nanstd(x) * critic;
+		@inbounds for k in eachindex(x)
+			bv[k] = abs(x[k] - _mean)  > _std
+		end
+	elseif (method == :threshold)
+		q_l, q_u = quantile(x, _thresh)
+		@inbounds for k in eachindex(x)
+			bv[k] = (x[k] < q_l) || (x[k] > q_u)
+		end
+	else			# method == :quartiles
+		q25, q75 = quantile(x, [0.25, 0.75])
+		crit = (q75 - q25) * 1.5
+		q_minus, q_plus  = q25 - crit, q75 + crit
+		@inbounds for k in eachindex(x)
+			bv[k] = (x[k] < q_minus) || (x[k] > q_plus)
+		end
+	end
+	return bv
+end
+
+function isoutlier(mat::AbstractMatrix{<:Real}; critic=3.0, method::Symbol=:median, threshold=nothing, width=0.0)
+	width > 0 && return isoutlier(mat2ds(mat), critic=critic, method=method, threshold=threshold, width=width)
+	bm = BitMatrix(undef, size(mat))
+	for k = 1:size(mat,2)
+		bm[:,k] .= isoutlier(view(mat, :, k), critic=critic, method=method, threshold=threshold)
+	end
+	return bm
+end
+
+function isoutlier(D::GMTdataset{<:Real,2}; critic=3.0, method::Symbol=:median, threshold=nothing, width=0.0)
+	(threshold === nothing && width <= 0) && error("Must provide the window length of via the `width` option or use the `threshold` option.")
+	if (width > 0)
+		method in [:median, :mean] || throw(ArgumentError("Unknown method: $method. Here use only one of :median or :mean"))
+		(method == :mean) && (method = boxcar)		# :boxcar is the name in GMT for 'mean'
+		Dres = filter1d(D, filter=(type=method, width=width, highpass=true), E=true)
+		isoutlier(view(Dres.data, :, 2), critic=critic, method=method, threshold=threshold)
+	else
+		isoutlier(view(D.data, :, 2), critic=critic, method=method, threshold=threshold)
+	end
+end

src/extras/isoutlier.jl

@@ -43,6 +43,7 @@ function trend1d(cmd0::String="", arg1=nothing; kwargs...)
 
 	cmd, = parse_common_opts(d, "", [:V_params :b :d :e :f :h :i :w :yx])
 	cmd  = parse_these_opts(cmd, d, [[:C :condition_number], [:I :conf_level :confidence_level], [:F :out :output], [:W :weights]])
+	opt_F = scan_opt(cmd, "-F")
 	((val = find_in_dict(d, [:N :model], false)[1]) === nothing) && error("The option 'model' must be specified")
 	if (isa(val, Tuple) && isa(val[1], NamedTuple))
 		# Complicated case here -- A mixed model ---. So input must be a Tuple of NamedTuples, and the +l+o+r options
@@ -68,7 +69,27 @@ function trend1d(cmd0::String="", arg1=nothing; kwargs...)
 	(tryparse(Int, opt_N[4:end]) !== nothing) && (opt_N = " -Np" * opt_N[4:end])	# -N2 is old syntax and will error GMT
 	cmd *= opt_N
 
-	common_grd(d, cmd0, cmd, "trend1d ", arg1)		# Finish build cmd and run it
+	R = common_grd(d, cmd0, cmd, "trend1d ", arg1)		# Finish build cmd and run it
+	!isa(R, GDtype) && return R							# Should be the output of Vd=2
+
+	if (opt_F != "")		# Extract column names from opt_F
+		if (opt_F[1] == 'p' || opt_F[1] == 'P' || opt_F[1] == 'c')
+			colnames = ["a$i" for i=0:size(D,2)-1]
+		else
+			s = collect(opt_F)
+			colnames = Vector{String}(undef, length(s))
+			for k = 1:numel(s)
+				if     (s[k] == 'x')  colnames[k] = "x"
+				elseif (s[k] == 'y')  colnames[k] = "y"
+				elseif (s[k] == 'm')  colnames[k] = "model"
+				elseif (s[k] == 'r')  colnames[k] = "residues"
+				elseif (s[k] == 'w')  colnames[k] = "weights"
+				end
+			end
+		end
+		isa(R, GMTdataset) ? (R.colnames = colnames) : (R[1].colnames = colnames)
+	end
+	return R
 end
 
 # ---------------------------------------------------------------------------------------------------

src/lowess.jl src/extras/lowess.jlsrc/lowess.jl renamed to src/extras/lowess.jl

@@ -1439,111 +1439,3 @@ end
 #r = Ref(tuple(5.0, 2.0, 1.0, 6.0))
 #p = Base.unsafe_convert(Ptr{Float64}, r)
 #u = unsafe_wrap(Array, p, 4)
-
-# ------------------------------------------------------------------------------------------------------
-"""
-	bv = isoutlier(x; critic=3.0, method::Symbol=:median, threshold=nothing, width=0.0) -> BitVector or BitMatrix
-
-Return a logical array whose elements are true when an outlier is detected in the corresponding element of `x`.
-
-If `x` is a matrix, and `width` is not used, then ``isoutlier`` operates on each column of `x` separately.
-By default, an outlier is a value that is more than three median absolute deviations (MAD) from the median.
-But that can be changed via the `critic` option.
-
-- `critic`: A number that sets the threshold for detecting outliers. The default is 3.0, which means that
-   a value is considered an outlier if it is more than 3 MAD from the median (when `method` is `:median`),
-   or more than 3 standard deviations (when `method` is `:mean`).
-- `method`: The method used to calculate the threshold for outliers. It can be one of the following:
-  - `:median`: Uses the median absolute deviation (MAD) method. Outliers are defined as elements more than
-   `critic` MAD from the median.
-  - `:mean`: Uses the mean and standard deviation method. Outliers are defined as elements more than `critic`
-   standard deviations from the mean. This method is faster but less robust than "median".
-  - `:quartiles`: Uses the interquartile range (IQR) method. Outliers are defined as elements more than 1.5
-   interquartile ranges above the upper quartile (75 percent) or below the lower quartile (25 percent).
-   This method is useful when the data in `x` is not normally distributed.
-- `threshold`: Is an alternative to the `method` option. It specifies the percentile thresholds, given as a
-   two-element array whose elements are in the interval [0, 100]. The first element indicates the lower percentile
-   threshold, and the second element indicates the upper percentile threshold. It can also be a single number
-   between (0, 100), which is interpreted as the percentage of end member points that are considered outliers.
-   For example, `threshold=1` means that the lower and upper thresholds are the 1th and 99th percentiles.
-- `width`: If this option is used (only when `x` is a Matrix or a GMTdataset) we detect local outliers using
-   a moving window method with window length ``width``. This length is given in the same units as the input
-   data stored in first column of `x`.
-
-### Example:
-```julia
-x = [57, 59, 60, 100, 59, 58, 57, 58, 300, 61, 62, 60, 62, 58, 57];
-findall(isoutlier(x))
-2-element Vector{Int64}:
- 4
- 9
-```
-
-```julia
-x = -50.0:50;	y = x / 50 .+ 3 .+ 0.25 * rand(length(x));
-y[[30,50,60]] = [4,-3,6];	# Add 3 outliers
-findall(isoutlier([x y], width=5))
-```
-"""
-function isoutlier(x::AbstractVector{<:Real}; critic=3.0, method::Symbol=:median, threshold=nothing)
-	method in [:median, :mean, :quartiles] || throw(ArgumentError("Unknown method: $method. Use :median, :mean or :quartiles"))
-	if (threshold !== nothing)
-		!(isa(threshold, VecOrMat{<:Real}) || isa(threshold, Real)) &&
-			throw(ArgumentError("Threshold: $threshold must be a scalar or a 2 elements array."))
-		if (isa(threshold, VecOrMat{<:Real}))
-			@assert !(length(threshold) == 2 && threshold[1] >= 0.0 && threshold[2] <= 100) "Threshold must be a 2 elements array in the [0 100] interval"
-			_thresh = [Float64(threshold[1])/100, Float64(threshold[2])/100]	# Make it [0 1]
-		else
-			@assert threshold > 0.0 && threshold < 100 "When threshold is a scalar, it must be in the ]0 100[ interval."
-			_thresh = [Float64(threshold)/100, Float64(100 - threshold)/100]	# Make it [0 1]
-		end
-		method = :threshold
-	end
-
-	bv = BitVector(undef, length(x))
-	if (method == :median)
-		_mad, _med = mad(x);	_mad *= critic
-		@inbounds for k in eachindex(x)
-			bv[k] = abs(x[k] - _med)  > _mad
-		end
-	elseif (method == :mean)
-		_mean, _std = nanmean(x), nanstd(x) * critic;
-		@inbounds for k in eachindex(x)
-			bv[k] = abs(x[k] - _mean)  > _std
-		end
-	elseif (method == :threshold)
-		q_l, q_u = quantile(x, _thresh)
-		@inbounds for k in eachindex(x)
-			bv[k] = (x[k] < q_l) || (x[k] > q_u)
-		end
-	else			# method == :quartiles
-		q25, q75 = quantile(x, [0.25, 0.75])
-		crit = (q75 - q25) * 1.5
-		q_minus, q_plus  = q25 - crit, q75 + crit
-		@inbounds for k in eachindex(x)
-			bv[k] = (x[k] < q_minus) || (x[k] > q_plus)
-		end
-	end
-	return bv
-end
-
-function isoutlier(mat::AbstractMatrix{<:Real}; critic=3.0, method::Symbol=:median, threshold=nothing, width=0.0)
-	width > 0 && return isoutlier(mat2ds(mat), critic=critic, method=method, threshold=threshold, width=width)
-	bm = BitMatrix(undef, size(mat))
-	for k = 1:size(mat,2)
-		bm[:,k] .= isoutlier(view(mat, :, k), critic=critic, method=method, threshold=threshold)
-	end
-	return bm
-end
-
-function isoutlier(D::GMTdataset{<:Real,2}; critic=3.0, method::Symbol=:median, threshold=nothing, width=0.0)
-	(threshold === nothing && width <= 0) && error("Must provide the window length of via the `width` option or use the `threshold` option.")
-	if (width > 0)
-		method in [:median, :mean] || throw(ArgumentError("Unknown method: $method. Here use only one of :median or :mean"))
-		(method == :mean) && (method = boxcar)		# :boxcar is the name in GMT for 'mean'
-		Dres = filter1d(D, filter=(type=method, width=width, highpass=true), E=true)
-		isoutlier(view(Dres.data, :, 2), critic=critic, method=method, threshold=threshold)
-	else
-		isoutlier(view(D.data, :, 2), critic=critic, method=method, threshold=threshold)
-	end
-end