v1.3.0

Project.toml

@@ -1,12 +1,13 @@
 name = "Kriging"
 uuid = "6133562b-a8e8-5ea3-954b-37642b766dad"
-version = "1.2.1"
+version = "1.3.0"
 
 [deps]
 DocumentFunction = "e1f3b4f0-2dc4-57d3-83f7-d4b7faf3b05b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 DocumentFunction = "1"

examples/example.jl

@@ -25,6 +25,8 @@ PyPlot.close(fig)
 covfun(h) = Kriging.expcov(h, 0.1, 3)
 xs = collect(range(0; stop=10, length=100))
 ys = collect(range(0, stop=10, length=100))
+x0mat = Kriging.getgridpoints(xs, ys)
+
 krigedfield = Array{Float64}(undef, length(xs), length(ys))
 @time for (i, x) in enumerate(xs), (j, y) in enumerate(ys)
 	krigedfield[i, j] = Kriging.krige(permutedims([x y]), X, Z, covfun)[1]
@@ -41,20 +43,8 @@ display(fig); println()
 fig.savefig("kriging.pdf")
 PyPlot.close(fig)
 
-x0mat = Array{Float64}(undef, 2, length(xs) * length(ys))
-global k = 1
-for (i, x) in enumerate(xs), (j, y) in enumerate(ys)
-	x0mat[1, k] = x
-	x0mat[2, k] = y
-	global k += 1
-end
-@time z0 = Kriging.interpolate_neighborhood(x0mat, X, Z, covfun, 20; neighborsearch=50)
-kriging_neighborhood_field = Array{Float64}(undef, length(xs), length(ys))
-global k = 1
-for (i, x) in enumerate(xs), (j, y) in enumerate(ys)
-	kriging_neighborhood_field[i, j] = z0[k]
-	global k += 1
-end
+@time z0 = Kriging.interpolate_neighborhood(x0mat, X, Z; cov=covfun, neighborsearch=50)
+kriging_neighborhood_field = Kriging.putgridpoints(xs, ys, z0)
 fig, ax = PyPlot.subplots()
 cax = ax.imshow(permutedims(kriging_neighborhood_field), extent=[0, 10, 0, 10], origin="lower")
 for i = 1:size(X, 2)
@@ -82,20 +72,8 @@ display(fig); println()
 fig.savefig("inversedistance.pdf")
 PyPlot.close(fig)
 
-x0mat = Array{Float64}(undef, 2, length(xs) * length(ys))
-global k = 1
-for (i, x) in enumerate(xs), (j, y) in enumerate(ys)
-	x0mat[1, k] = x
-	x0mat[2, k] = y
-	global k += 1
-end
-@time z0 = Kriging.interpolate_neighborhood(x0mat, X, Z, covfun, 20; neighborsearch=50, interpolate=Kriging.inversedistance, pow=2)
-inversedistance_neighborhood_field = Array{Float64}(undef, length(xs), length(ys))
-global k = 1
-for (i, x) in enumerate(xs), (j, y) in enumerate(ys)
-	inversedistance_neighborhood_field[i, j] = z0[k]
-	global k += 1
-end
+@time z0 = Kriging.interpolate_neighborhood(x0mat, X, Z; neighborsearch=50, interpolate=Kriging.inversedistance, pow=2)
+inversedistance_neighborhood_field = Kriging.putgridpoints(xs, ys, z0)
 fig, ax = PyPlot.subplots()
 cax = ax.imshow(permutedims(inversedistance_neighborhood_field), extent=[0, 10, 0, 10], origin="lower")
 for i = 1:size(X, 2)
@@ -107,20 +85,8 @@ display(fig); println()
 fig.savefig("inversedistance_neighborhood.pdf")
 PyPlot.close(fig)
 
-x0mat = Array{Float64}(undef, 2, length(xs) * length(ys))
-global k = 1
-for (i, x) in enumerate(xs), (j, y) in enumerate(ys)
-	x0mat[1, k] = x
-	x0mat[2, k] = y
-	global k += 1
-end
 @time z0 = Kriging.condsim(x0mat, X, Z, covfun, 20; neighborsearch=50)
-condsimfield = Array{Float64}(undef, length(xs), length(ys))
-global k = 1
-for (i, x) in enumerate(xs), (j, y) in enumerate(ys)
-	condsimfield[i, j] = z0[k]
-	global k += 1
-end
+condsimfield = Kriging.putgridpoints(xs, ys, z0)
 fig, ax = PyPlot.subplots()
 cax = ax.imshow(permutedims(condsimfield), extent=[0, 10, 0, 10], origin="lower")
 for i = 1:size(X, 2)

src/Kriging.jl

@@ -5,6 +5,11 @@ import DocumentFunction
 import LinearAlgebra
 import Random
 
+"Test Kriging"
+function test()
+	include(joinpath(Base.pkgdir(Kriging), "test", "runtests.jl"))
+end
+
 """
 Gaussian spatial covariance function
 
@@ -171,7 +176,7 @@ Returns:
 """
 function simplekrige(mu, x0mat::AbstractMatrix, X::AbstractMatrix, Z::AbstractVector, cov::Function)
 	@assert size(X, 2) == length(Z) "Number of data points ($(size(X, 2))) and observations ($(length(Z))) do not match!"
-	@assert size(x0mat, 1) == size(X, 1) "Dimensions of data points ($(size(x0mat, 1))) and observations ($(size(X, 1))) do not match!"
+	@assert size(x0mat, 2) == size(X, 1) "Dimensions of data points ($(size(x0mat, 2))) and observations ($(size(X, 1))) do not match!"
 	result = fill(mu, size(x0mat, 2))
 	resultvariance = fill(cov(0), size(x0mat, 2))
 	covmat = getcovmat(X, cov)
@@ -278,7 +283,7 @@ Returns:
 function condsim(x0mat::AbstractMatrix, X::AbstractMatrix, Z::AbstractVector, cov::Function, numneighbors, numobsneighbors=min(1000, size(X, 2)); neighborsearch=min(1000, size(x0mat, 2)))
 	@assert size(X, 2) == length(Z) "Number of data points ($(size(X, 2))) and observations ($(length(Z))) do not match!"
 	@assert size(x0mat, 1) == size(X, 1) "Dimensions of data points ($(size(x0mat, 1))) and observations ($(size(X, 1))) do not match!"
-	nnindices, nnindices_obs = kdtree_indices(x0mat, X; neighborsearch=neighborsearch, numobsneighbors=numobsneighbors)
+	nnindices, nnindices_obs = kdtree_indices(x0mat, X; numpredneighbors=neighborsearch, numobsneighbors=numobsneighbors)
 	z0 = Vector{Float64}(undef, size(x0mat, 2))
 	filledin = fill(false, size(x0mat, 2))
 	perm = Random.randperm(size(x0mat, 2))
@@ -309,10 +314,10 @@ function condsim(x0mat::AbstractMatrix, X::AbstractMatrix, Z::AbstractVector, co
 	return z0
 end
 
-function interpolate_neighborhood(x0mat::AbstractMatrix, X::AbstractMatrix, Z::AbstractVector, cov::Function, numneighbors, numobsneighbors=min(1000, size(X, 2)); neighborsearch=min(1000, size(x0mat, 2)), interpolate=krigevariance, return_variance::Bool=false, kw...)
+function interpolate_neighborhood(x0mat::AbstractMatrix, X::AbstractMatrix, Z::AbstractVector; cov::Function=cov(h)=Kriging.expcov(h, 0.1, 3), numobsneighbors=min(1000, size(X, 2)), neighborsearch=min(1000, size(x0mat, 2)), interpolate=krigevariance, return_variance::Bool=false, cutoff::Number=0, kw...)
 	@assert size(X, 2) == length(Z) "Number of data points ($(size(X, 2))) and observations ($(length(Z))) do not match!"
 	@assert size(x0mat, 1) == size(X, 1) "Dimensions of data points ($(size(x0mat, 1))) and observations ($(size(X, 1))) do not match!"
-	_, nnindices_obs = kdtree_indices(x0mat, X; neighborsearch=neighborsearch, numobsneighbors=numobsneighbors)
+	_, nnindices_obs = kdtree_indices(x0mat, X; numpredneighbors=neighborsearch, numobsneighbors=numobsneighbors, cutoff_obs=cutoff)
 	z0 = Vector{Float64}(undef, size(x0mat, 2))
 	if return_variance
 		var0 = Vector{Float64}(undef, size(x0mat, 2))
@@ -332,17 +337,31 @@ function interpolate_neighborhood(x0mat::AbstractMatrix, X::AbstractMatrix, Z::A
 	end
 end
 
-function kdtree_indices(x0mat::AbstractMatrix{T}, X::AbstractMatrix=Matrix{T}(undef, 0, 0); neighborsearch=min(1000, size(x0mat, 2)), numobsneighbors=min(1000, size(X, 2))) where T <: Real
-	kdtree = NearestNeighbors.KDTree(x0mat)
-	nnindices, _ = NearestNeighbors.knn(kdtree, x0mat, neighborsearch, true)
-	if numobsneighbors > 0 && size(X, 2) > 0
-		@assert size(x0mat, 1) == size(X, 1) "Dimensions of data points ($(size(x0mat, 1))) and observations ($(size(X, 1))) do not match!"
-		obs_kdtree = NearestNeighbors.KDTree(X)
-		nnindices_obs, _ = NearestNeighbors.knn(obs_kdtree, x0mat, numobsneighbors, true)
-		return nnindices, nnindices_obs
+function kdtree_indices(x_pred::AbstractMatrix{T}, x_obs::AbstractMatrix=Matrix{T}(undef, 0, 0); numpredneighbors=min(1000, size(x_pred, 2)), numobsneighbors=min(1000, size(x_obs, 2)), cutoff_pred::Number=0, cutoff_obs::Number=0) where T <: Real
+	if numpredneighbors > 0
+		kdtree = NearestNeighbors.KDTree(x_pred)
+		nnindices_pred, nndistances_pred = NearestNeighbors.knn(kdtree, x_pred, numpredneighbors, true)
+		if cutoff_pred > 0
+			for i = 1:length(nnindices_pred)
+				nnindices_pred[i] = nnindices_pred[i][nndistances_pred[i] .<= cutoff_pred]
+			end
+		end
+	else
+		nnindices_pred = []
+	end
+	if numobsneighbors > 0 && size(x_obs, 2) > 0
+		@assert size(x_pred, 1) == size(x_obs, 1) "Dimensions of data points ($(size(x_pred, 1))) and observations ($(size(x_obs, 1))) do not match!"
+		obs_kdtree = NearestNeighbors.KDTree(x_obs)
+		nnindices_obs, nndistances_obs = NearestNeighbors.knn(obs_kdtree, x_pred, numobsneighbors, true)
+		if cutoff_obs > 0
+			for i = 1:length(nnindices_obs)
+				nnindices_obs[i] = nnindices_obs[i][nndistances_obs[i] .<= cutoff_obs]
+			end
+		end
 	else
-		return nnindices, []
+		nndistances_obs = []
 	end
+	return nnindices_pred, nnindices_obs
 end
 
 """
@@ -396,7 +415,7 @@ end
 
 function estimationerror(w::AbstractVector, x0::AbstractVector, X, cov::Function)
 	covmat = getcovmat(X, cov)
-	covvec = Array{Float64}(undef, size(X, 2))
+	covvec = Vector{Float64}(undef, size(X, 2))
 	getcovvec!(covvec, x0, X, cov)
 	cov0 = cov(0.)
 	return estimationerror(w, x0, X, covmat, covvec, cov0)
@@ -423,8 +442,8 @@ Returns:
 - estimation kriging error
 """ estimationerror
 
-function getgridpoints(xs::AbstractVector, ys::AbstractVector)
-	gridxy = Array{Float64}(undef, 2, length(xs) * length(ys))
+function getgridpoints(xs::Union{AbstractVector,AbstractRange}, ys::Union{AbstractVector,AbstractRange})
+	gridxy = Matrix{Float64}(undef, 2, length(xs) * length(ys))
 	local i = 1
 	for x in xs
 		for y in ys
@@ -435,8 +454,8 @@ function getgridpoints(xs::AbstractVector, ys::AbstractVector)
 	end
 	return gridxy
 end
-function getgridpoints(xs::AbstractVector, ys::AbstractVector, zs::AbstractVector)
-	gridxyz = Array{Float64}(undef, 3, length(xs) * length(ys) * length(zs))
+function getgridpoints(xs::Union{AbstractVector,AbstractRange}, ys::Union{AbstractVector,AbstractRange}, zs::Union{AbstractVector,AbstractRange})
+	gridxyz = Matrix{Float64}(undef, 3, length(xs) * length(ys) * length(zs))
 	local i = 1
 	for x in xs
 		for y in ys
@@ -464,6 +483,19 @@ Returns:
 - grid points
 """ getgridpoints
 
+function putgridpoints(xs::Union{AbstractVector,AbstractRange}, ys::Union{AbstractVector,AbstractRange}, zs::AbstractVector)
+	@assert length(xs) * length(ys) == length(zs) "Number of grid points ($(length(xs) * length(ys))) and values ($(length(zs))) do not match!"
+	gridxy = Matrix{Float64}(undef, length(xs), length(ys))
+	local k = 1
+	for (i, _) in enumerate(xs)
+		for (j, _) in enumerate(ys)
+			gridxy[i, j] = zs[k]
+			k += 1
+		end
+	end
+	return gridxy
+end
+
 function grid2layers(obs::AbstractVector, xs::AbstractVector, ys::AbstractVector, zs::AbstractVector)
 	layers = Array{Array{Float64, 2}}(undef, length(zs))
 	for k = 1:length(zs)

test/runtests.jl

@@ -15,8 +15,8 @@ mycov1 = h->Kriging.gaussiancov(h, 2, 300)
 mycov2 = h->Kriging.expcov(h, 2, 300)
 mycov3 = h->Kriging.sphericalcov(h, 2, 300)
 
-x0 = [.5 .5; .49 .49; .01 .01; .99 1.; 0. 1.; 1. 0.]
-xs = [0. 0.; 0. 1.; 1. 0.; 1. 1.; 0.500001 0.5]'
+x0 = permutedims([.5 .5; .49 .49; .01 .01; .99 1.; 0. 1.; 1. 0.])
+xs = permutedims([0. 0.; 0. 1.; 1. 0.; 1. 1.; 0.500001 0.5])
 zs = [-20., .6, .4, 1., 20.]
 
 krige_results_1 = Kriging.krige(x0, xs, zs, mycov1)
@@ -48,9 +48,9 @@ estimation_error_3 = Kriging.estimationerror(ones(size(xs, 2)), zs, xs, mycov3)
 	end
 
 	@Test.testset "Krige" begin
-		@Test.test isapprox(krige_results_1, [19.8891, 19.8891], atol=0.1)
-		@Test.test isapprox(krige_results_2, [19.4586, 19.4586], atol=0.1)
-		@Test.test isapprox(krige_results_3, [19.4586, 19.4586], atol=0.1)
+		@Test.test isapprox(krige_results_1, [19.99993512367364, 19.778132898057812, -18.621832575427835, 1.4811239442322404, 0.6000381594407371, 0.39998541082022715], atol=0.1)
+		@Test.test isapprox(krige_results_2, [19.999945901520327, 19.172644554834168, -19.219551440274383, 1.2150507275634252, 0.5999999999977101, 0.3999999999961652], atol=0.1)
+		@Test.test isapprox(krige_results_3, [19.999945901706013, 19.17264693896976, -19.219551074375026, 1.2151032013981584, 0.600000000000781, 0.39999999999867536], atol=0.1)
 	 end
 
 	# Testing Kriging.estimationerror()