Skip to content

Commit

Permalink
inversedistance fix
Browse files Browse the repository at this point in the history
  • Loading branch information
montyvesselinov committed Jul 14, 2020
1 parent 5dbf210 commit 7f6ecdb
Show file tree
Hide file tree
Showing 2 changed files with 23 additions and 7 deletions.
24 changes: 20 additions & 4 deletions examples/example.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,12 @@ import PyPlot
import Random

Random.seed!(0)
X = [3.0 7 1 5 9; 3 2 7 9 6]
X = [3 7 1 5 9; 3 2 7 9 6]
Z = exp.(randn(size(X, 2)))
fig, ax = PyPlot.subplots()
for i = 1:size(X, 2)
ax.plot([X[1, i]], [X[2, i]], "r.", ms=10)
ax.text([X[1, i] + 0.1], [X[2, i] + 0.1], "$(Z[i])"[1:3], fontsize=16)
ax.text([X[1, i] + 0.1], [X[2, i] + 0.1], "$(round(Z[i]; sigdigits=2))", fontsize=16)
end
ax.plot([5], [5], "r.", ms=10)
ax.text([5.1], [5.1], "?", fontsize=16)
Expand All @@ -26,8 +26,9 @@ 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]
end

fig, ax = PyPlot.subplots()
cax = ax.imshow(krigedfield', extent=[0, 10, 0, 10], origin="lower")
cax = ax.imshow(permutedims(krigedfield), extent=[0, 10, 0, 10], origin="lower")
for i = 1:size(X, 2)
ax.plot([X[1, i]], [X[2, i]], "r.", ms=10)
end
Expand All @@ -36,6 +37,21 @@ display(fig); println()
fig.savefig("kriging.pdf")
PyPlot.close(fig)

inversedistancefield = Array{Float64}(undef, length(xs), length(ys))
@time for (i, x) in enumerate(xs), (j, y) in enumerate(ys)
inversedistancefield[i, j] = Kriging.inversedistance(permutedims([x y]), X, Z, 1/2)[1]
end

fig, ax = PyPlot.subplots()
cax = ax.imshow(permutedims(inversedistancefield), extent=[0, 10, 0, 10], origin="lower")
for i = 1:size(X, 2)
ax.plot([X[1, i]], [X[2, i]], "r.", ms=10)
end
fig.colorbar(cax)
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)
Expand All @@ -51,7 +67,7 @@ for (i, x) in enumerate(xs), (j, y) in enumerate(ys)
k += 1
end
fig, ax = PyPlot.subplots()
cax = ax.imshow(condsimfield', extent=[0, 10, 0, 10], origin="lower")
cax = ax.imshow(permutedims(condsimfield), extent=[0, 10, 0, 10], origin="lower")
for i = 1:size(X, 2)
ax.plot([X[1, i]], [X[2, i]], "r.", ms=10)
end
Expand Down
6 changes: 3 additions & 3 deletions src/Kriging.jl
Original file line number Diff line number Diff line change
Expand Up @@ -119,15 +119,15 @@ end
function distance(a::AbstractArray, b::AbstractArray)
result = 0.0
for i = 1:length(a)
result += (a[i] - b[i])
result += abs.(a[i] .- b[i])
end
return result
end

function distsquared(a::AbstractArray, b::AbstractArray)
result = 0.0
for i = 1:length(a)
result += (a[i] - b[i]) ^ 2
result += abs.(a[i] .- b[i]) ^ 2
end
return result
end
Expand All @@ -137,7 +137,7 @@ function inversedistance(x0mat::AbstractMatrix, X::AbstractMatrix, Z::AbstractVe
weights = Array{Float64}(undef, size(X, 2))
for i = 1:size(x0mat, 2)
for j = 1:size(X, 2)
weights[j] = 1 / distance(x0mat[:, i], X[:, j]) ^ pow
weights[j] = 1. / distance(x0mat[:, i], X[:, j]) ^ pow
end
result[i] = LinearAlgebra.dot(weights, Z) / sum(weights)
end
Expand Down

0 comments on commit 7f6ecdb

Please sign in to comment.