Add a Lowess smoother filter (from https://github.com/xKDR/Lowess.jl) (…

…#1684)

* Test/demo file for the whittaker smoothing function.

* Add lowess tests

* Include and export lowess smoother.

* Fix an Example in docs.

* Fix a case when plotting lines, symbols and a nested plot.

* File with lowess smother filter.

src/GMT.jl

@@ -178,7 +178,7 @@ export
 	findpeaks, makeDCWs, mksymbol, circfit,
 
 	gunique, sortslicesperm,
-	hampel, hampel!, whittaker,
+	hampel, hampel!, lowess, whittaker,
 
 	Ginnerjoin, Gouterjoin, Gleftjoin, Grightjoin, Gcrossjoin, Gsemijoin, Gantijoin, spatialjoin,
 	groupby, stats,
@@ -337,6 +337,7 @@ 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")

src/common_options.jl

@@ -4407,7 +4407,8 @@ function finish_PS_module(d::Dict, cmd::Vector{String}, opt_extra::String, K::Bo
 		CTRL.pocket_call[3] = nothing		# Reset
 		fi = 2		# First index, the one that the args... respect, start at 2
 	end
-	for k = fi:lastindex(cmd)
+	n_cmds = length(cmd)
+	for k = fi:n_cmds
 		is_psscale = (startswith(cmd[k], "psscale") || startswith(cmd[k], "colorbar"))
 		is_pscoast = (startswith(cmd[k], "pscoast") || startswith(cmd[k], "coast"))
 		is_basemap = (startswith(cmd[k], "psbasemap") || startswith(cmd[k], "basemap"))
@@ -4440,7 +4441,11 @@ function finish_PS_module(d::Dict, cmd::Vector{String}, opt_extra::String, K::Bo
 			continue
 		end
 		# Allow also plot data from a nested call to plot
-		P = !(k > fi && is_plot && (CTRL.pocket_call[1] !== nothing)) ? gmt(cmd[k], args...) : gmt(cmd[k], get_pocket_call())
+		cond = !(k > fi && is_plot && (CTRL.pocket_call[1] !== nothing))
+		# Next line means cases like 'plot(D, marker=:circ, ms="4p", lt=0.5, plot=(data=...))' where data in 'D'
+		# must be used twice before we pass to the next 'plot' (a nested call).
+		!cond && is_plot && (k == fi + 1) && (k < n_cmds) && contains(cmd[k], " -S") && (cond = true)
+		P = cond ? gmt(cmd[k], args...) : gmt(cmd[k], get_pocket_call())
 
 		# If we had a double frame to plot Geog on a Cartesian plot we must reset memory to original -J & -R so
 		# that appending other plots to same fig can continue to work and not fail because proj had become Geog.

src/lowess.jl

@@ -0,0 +1,254 @@
+# From https://github.com/xKDR/Lowess.jl with some (few) modifications.
+"""
+```julia
+lowess(x, y, span = 2 / 3, nsteps = 3, delta = 0.01 * (maximum(x) - minimum(x)))
+```
+
+Compute the smooth of a scatterplot of `y` against `x` using robust locally weighted regression. Input vectors `x` and `y` must contain either integers or floats. Parameters `span` and `delta` must be of type `T`, where `T <: AbstractFloat`. Returns a vector `ys`; `ys[i]` is the fitted value at `x[i]`. To get the smooth plot, `ys` must be plotted against `x`.
+
+# Arguments
+
+  - `x::Vector`: Abscissas of the points on the scatterplot. `x` must be ordered.
+  - `y::Vector`: Ordinates of the points in the scatterplot.
+  - `span::T`: The amount of smoothing.
+  - `nsteps::Integer`: Number of iterations in the robust fit.
+  - `delta::T`: A nonnegative parameter which may be used to save computations.
+
+# Example
+
+```julia
+x = sort(10 .* rand(100))
+y = sin.(x) .+ 0.5 * rand(100)
+ys = lowess(x, y, span=0.2)
+scatter(x, y)
+plot!(x, ys)
+```
+"""
+function lowess(D::GMTdataset; span=2/3, nsteps=3, delta=0.0)
+	(delta == 0.0) && (delta = 0.01 * (D.bbox[2] - D.bbox[1]))
+	lowess(view(D.data, :, 1), view(D.data, :, 2); span=span, nsteps=nsteps, delta=delta)
+end
+
+function lowess(x::AbstractVector{R}, y::AbstractVector{S}; span::T=2/3, nsteps::Integer=3,
+                delta::T = 0.01 * (maximum(x) - minimum(x))) where {R <: Real, S <: Real, T <: AbstractFloat}
+	lowess((R <: AbstractFloat) ? x : Vector{Float64}(x), (S <: AbstractFloat) ? y : Vector{Float64}(y); span=span, nsteps=nsteps, delta=delta)
+end
+
+function lowess(x::AbstractVector{T}, y::AbstractVector{T}; span::T=2/3, nsteps::Integer=3,
+                delta::T = 0.01 * (maximum(x) - minimum(x))) where {T <: AbstractFloat}
+	# defining needed variables
+
+	n::Int = length(x)
+	ys::Vector{T} = Vector{T}(undef, n)
+	rw::Vector{T} = Vector{T}(undef, n)
+	res::Vector{T} = Vector{T}(undef, n)
+
+	iter::Int = 0
+	ok::Vector{Int} = Vector{Int}(undef, 1)
+	# for safety, initialize ok to 0
+	ok[1] = 0
+
+	i::Int = 0
+	j::Int = 0
+	last::Int = 0
+	m1::Int = 0
+	m2::Int = 0
+	nleft::Int = 0
+	nright::Int = 0
+	ns::Int = 0
+	d1::T = 0.0
+	d2::T = 0.0
+	denom::T = 0.0
+	alpha::T = 0.0
+	cut::T = 0.0
+	cmad::T = 0.0
+	c9::T = 0.0
+	c1::T = 0.0
+	r::T = 0.0
+
+	if (n < 2)
+		ys[1] = y[1]
+		return ys
+	end
+
+	ns = max(min(floor(Int, span * n), n), 2)  # at least two, at most n points
+	for iter = 1:(nsteps + 1)  # robustness iterations
+		nleft = 0
+		nright = ns - 1
+		last = -1   # index of prev estimated point
+		i = 0   # index of current point
+
+		while true
+			while (nright < n - 1)
+				# move nleft, nright to right if radius decreases
+				d1 = x[i + 1] - x[nleft + 1]
+				d2 = x[nright + 2] - x[i + 1]
+				# if d1 <= d2 with x[nright + 2] == x[nright + 1], lowest fixes
+				(d1 <= d2) && break
+				# radius will not decrease by move right
+				nleft = nleft + 1
+				nright = nright + 1
+			end
+
+			lowest(x, y, n, x[i + 1], ys, i, nleft, nright, res, (iter > 1), rw, ok)
+
+			# fitted value at x[i + 1]
+			if (ok == 0)
+				ys[i + 1] = y[i + 1]
+			end
+
+			# all weights zero - copy over value (all rw==0)
+			if (last < i - 1)   # skipped points -- interpolate
+				denom = x[i + 1] - x[last + 1]  # non-zero - proof?
+				j = last + 1
+				for t = (last + 1):(i - 1) # t = j at all times
+					alpha = (x[j + 1] - x[last + 1]) / denom
+					ys[j + 1] = alpha * ys[i + 1] + (1.0 - alpha) * ys[last + 1]
+					j = j + 1
+				end
+			end
+
+			last = i    # last point actually estimated
+			cut = x[last + 1] + delta   # x coord of close points
+
+			# find close points
+			i = last + 1
+			for t = (last + 1):(n - 1) # t = i at all times
+				if (x[i + 1] > cut) # i one beyond last pt within cut
+					break
+				end
+
+				if (x[i + 1] == x[last + 1])
+					ys[i + 1] = ys[last + 1]
+					last = i
+				end
+				i = i + 1
+			end
+			i = max(last + 1, i - 1)
+
+			# back 1 point so interpolation within delta, but always go forward
+			# check do while loop condition
+			(last < n - 1) || break
+		end
+
+		# residuals
+		for i = 0:(n - 1)
+			res[i + 1] = y[i + 1] - ys[i + 1]
+		end
+
+		(iter > nsteps) && break  # compute robustness weights except last time
+
+		for i = 0:(n - 1)
+			rw[i + 1] = abs(res[i + 1])
+		end
+
+		sort!(rw)
+
+		m1 = floor(1 + n / 2)
+		m2 = n - m1 + 1
+		cmad = 3.0 * (rw[m1 + 1] + rw[m2 + 1])  # 6 median abs resid
+		c9 = 0.999 * cmad
+		c1 = 0.001 * cmad
+		for i = 0:(n - 1)
+			r = abs(res[i + 1])
+			if (r <= c1)    # near 0, avoid underflow
+				rw[i + 1] = 1.0
+			elseif (r > c9) # near 1, avoid underflow
+				rw[i + 1] = 0.0
+			else
+				rw[i + 1] = (1.0 - (r / cmad)^2)^2
+			end
+		end
+	end
+	return ys
+end
+
+# --------------------------------------------------------------------------------------------------
+function lowest(x::AbstractVector{T}, y::AbstractVector{T}, n::Integer, xs::T, ys::AbstractVector{T},
+	ys_pos::Integer, nleft::Integer, nright::Integer, w::AbstractVector{T}, userw::Bool, rw::AbstractVector{T},
+	ok::Vector{Int}) where {T <: AbstractFloat}
+	b::T = 0.0
+	c::T = 0.0
+	r::T = 0.0
+	nrt::Int = 0
+
+	# Julia indexing starts at 1, so add 1 to all indexes
+	range::T = x[n] - x[1]
+	h::T = max(xs - x[nleft + 1], x[nright + 1] - xs)
+	h9::T = 0.999 * h
+	h1::T = 0.001 * h
+
+	# compute weights (pick up all ties on right)
+	a::T = 0.0     # sum of weights
+	j::Int = nleft   # initialize j
+
+	for i = nleft:(n - 1)  # i = j at all times
+		w[j + 1] = 0.0
+		r = abs(x[j + 1] - xs)      # replaced fabs with abs
+		if (r <= h9)       # small enough for non-zero weight
+			if (r > h1)
+				w[j + 1] = (1.0 - (r / h)^3)^3
+			else
+				w[j + 1] = 1.0
+			end
+			if (userw)
+				w[j + 1] = rw[j + 1] * w[j + 1]
+			end
+			a += w[j + 1]
+		elseif (x[j + 1] > xs)      # get out at first zero wt on right
+			break
+		end
+		j = j + 1
+	end
+
+	nrt = j - 1     # rightmost pt (may be greater than nright because of ties)
+	if (a <= 0.0)
+		ok[1] = 0   # ok is a 1 length vector
+	else   # weighted least squares
+		ok[1] = 1
+
+		# make sum of w[j + 1] == 1
+		j = nleft
+		for i = nleft:nrt      # i = j at all times
+			w[j + 1] = w[j + 1] / a
+			j = j + 1
+		end
+
+		if (h > 0.0)    # use linear fit
+			# find weighted center of x values
+			j = nleft
+			a = 0.0
+			for i = nleft:nrt  # i = j at all times
+				a += w[j + 1] * x[j + 1]
+				j = j + 1
+			end
+
+			b = xs - a
+
+			j = nleft
+			c = 0.0
+			for i = nleft:nrt  # i = j at all times
+				c += w[j + 1] * (x[j + 1] - a) * (x[j + 1] - a)
+				j = j + 1
+			end
+
+			if (sqrt(c) > 0.001 * range)
+				# points are spread out enough to compute slope
+				b = b / c
+
+				j = nleft
+				for i = nleft:nrt  # i = j at all times
+					w[j + 1] = w[j + 1] * (1.0 + b * (x[j + 1] - a))
+					j = j + 1
+				end
+			end
+		end
+
+		j = nleft
+		ys[ys_pos + 1] = 0.0
+		for i = nleft:nrt  # i = j at all times
+			ys[ys_pos + 1] += w[j + 1] * y[j + 1]
+			j = j + 1
+		end
+	end
+end

src/whittaker.jl

@@ -38,7 +38,7 @@ t = 2001:0.003:2007;
 _v = 5*cospi.((t .- 2000)/2); v = _v + (5*rand(length(t)) .- 2.5);
 v[2002.6 .< t .< 2003.4] .= NaN;
 z = whittaker(t, v, 0.01, 3);
-plot(t, v, legend="Noisy", plot=(data=[t _v], lc=:green, lt=1, legend="Original"), show=1)
+plot(t, v, legend="Noisy", plot=(data=[t _v], lc=:green, lt=1, legend="Original"))
 plot!(t, z, lc=:red, lt=1, legend="Degree 3", show=1)
 ```
 """

test/runtests.jl

@@ -27,6 +27,7 @@ using FFTW
 	println("		Entering: test_proj4.jl")
 	include("test_proj4.jl")
 
+	include("test_lowess.jl")
 	include("test_whittaker.jl")
 	include("test_signalcorr.jl")
 	println("	MAGREF")
@@ -83,9 +84,7 @@ using FFTW
 	println("		Entering: test_solidos.jl")
 	include("test_solidos.jl")
 	include("test_statplots.jl")
-	println("		Entering: test_texture.jl")
 	include("test_texture.jl")
-	println("		Entering: test_pca.jl")
 	include("test_pca.jl")
 	println("		Entering: test_las.jl")
 	include("test_las.jl")

test/test_lowess.jl

@@ -0,0 +1,9 @@
+@testset "LOWESS" begin
+
+	println("	LOWESS")
+
+	x = sort(10 .* rand(100));
+	y = sin.(x) .+ 0.5 * rand(100);
+	ys = lowess(x, y, span=0.2);
+	ys = lowess(mat2ds([x y]), span=0.2);
+end

test/test_whittaker.jl

@@ -6,4 +6,6 @@
 	whittaker(D[:,2], 2e4);
 	whittaker(D[:,2], 2e4, 2; weights=D[:,3]);
 	whittaker(D, 2e4, 2);
+	t = 2001:0.003:2007; _v = 5*cospi.((t .- 2000)/2); v = _v + (5*rand(length(t)) .- 2.5); v[2002.6 .< t .< 2003.4] .= NaN;
+	whittaker(t, v, 0.01);
 end