Skip to content

Commit

Permalink
fix important bug reading external anisotropy data
Browse files Browse the repository at this point in the history
  • Loading branch information
@rmcaixeta
rmcaixeta committed Dec 2, 2024
1 parent 6213ec4 commit c2faff9
Show file tree
Hide file tree
Showing 8 changed files with 102 additions and 32 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "LocalAnisotropies"
uuid = "cb698b3a-2f1d-48b4-a1e9-ddc0b8f648e2"
authors = ["Rafael Caixeta and contributors"]
version = "0.6.2"
version = "0.7.0"

[deps]
Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
Expand Down
6 changes: 4 additions & 2 deletions ext/LocalAnisotropiesMakieExt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,13 +17,15 @@ Makie.plottype(::LocalAnisotropy, ::LocalAnisotropies.SpatialData) =
function pre_ellipsoid(lpars, geom, f)
# pre process
coords = LocalAnisotropies.coords_(geom)
mag = lpars.magnitude
mx = maximum
D3 = size(coords, 1) >= 3
x = coords[1, :]
y = coords[2, :]
z = D3 ? coords[3, :] : map(i -> 0, x)
s =
D3 ? [f * Makie.Vec3f(i[1], i[2], i[3]) for i in eachcol(lpars.magnitude)] :
[f * Makie.Vec3f(i[1], i[2], i[2]) for i in eachcol(lpars.magnitude)]
D3 ? [f * Makie.Vec3f(i[1]/mx(i), i[2]/mx(i), i[3]/mx(i)) for i in eachcol(mag)] :
[f * Makie.Vec3f(i[1]/mx(i), i[2]/mx(i), i[2]/mx(i)) for i in eachcol(mag)]
q = map(q -> Makie.Quaternion(q.q1, q.q2, q.q3, q.q0), lpars.rotation)

x, y, z, s, q
Expand Down
7 changes: 6 additions & 1 deletion src/estimators/kriging.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,12 @@ krig_estimator(model::KC_SKModel, localpar = nothing) = SimpleKriging(model.γ,
neighs_localaniso(model::MWModels, m) = nothing
neighs_localaniso(model::KCModels, m) = view(model.hdlocalaniso, m)

function neighs_localaniso(localaniso::LocalAnisotropy, dom, neigh; method::Symbol=:KernelConvolution)
function neighs_localaniso(
localaniso::LocalAnisotropy,
dom,
neigh;
method::Symbol = :KernelConvolution,
)
method != :KernelConvolution ? nothing : qmat(nnpars(localaniso, dom, neigh))
end

Expand Down
16 changes: 8 additions & 8 deletions src/estimators/sgs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,21 +146,21 @@ end

function local_probmodel(probmodel, localaniso, hdlocalaniso)
isnothing(hdlocalaniso) ? MW_SKModel(localaniso, probmodel.γ, probmodel.μ) :
KC_SKModel(localaniso, probmodel.γ, probmodel.μ, hdlocalaniso)
KC_SKModel(localaniso, probmodel.γ, probmodel.μ, hdlocalaniso)
end

function Meshes._pboxes(::Type{𝔼{N}}, points) where {N}
p = first(points)
ℒ = lentype(p)
cmin = fill(typemax(ℒ), N)
cmax = fill(typemin(ℒ), N)

for p in points
c = getfield(coords(p), :coords)
for i in 1:N
cmin[i] = min(c[i], cmin[i])
cmax[i] = max(c[i], cmax[i])
end
c = getfield(coords(p), :coords)
for i = 1:N
cmin[i] = min(c[i], cmin[i])
cmax[i] = max(c[i], cmax[i])
end
end
Box(withcrs(p, Tuple(cmin)), withcrs(p, Tuple(cmax)))
end
end
10 changes: 2 additions & 8 deletions src/estimators/variograms.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,8 @@ end

qmat(L::LocalAnisotropy, i::Int) = qmat(localpair(L, i)...)
qmat(L::LocalGeoData, i::Int) = qmat(localpair(L, i)...)
qmat(L::Union{LocalAnisotropy,Nothing}) = isnothing(L) ? nothing : [qmat(L, i) for i = 1:nvals(L)]
qmat(L::Union{LocalAnisotropy,Nothing}) =
isnothing(L) ? nothing : [qmat(L, i) for i = 1:nvals(L)]

function mw_estimator(μ, γ, localpar)
# get local Mahalanobis matrix
Expand Down Expand Up @@ -69,10 +70,3 @@ function kccov(γ::Variogram, xi, xj, Qi::AbstractMatrix, Qj::AbstractMatrix)
Cij = (sill(γl) - γl(xi, xj))
(det(Qi)^0.25) * (det(Qj)^0.25) * (det(Qij)^-0.5) * Cij
end

function localball(radii, rot, convention)
P, Λ = rotmat(radii, rot, convention)
Q = P' * Λ * P
Qs = ustrip.(Q ./ collect(radii)' .^ 2)
MetricBall(1.0, Mahalanobis(Symmetric(Qs)))
end
18 changes: 17 additions & 1 deletion src/parametrization/conventions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,22 @@ function rotmat(

# rotation matrix
P = angle_to_dcm(angles..., rule.order)[SOneTo(N), SOneTo(N)]
!rule.extrinsic && (P = P')
rule.extrinsic && (P = P')
P, Λ
end

function anisodistance(
semiaxes::AbstractVector,
angles::AbstractVector,
convention = :TaitBryanExtr;
rev = false,
)
P, Λ = rotmat(semiaxes, angles, convention; rev)
Q = P' * Λ * P
Mahalanobis(Symmetric(Q))
end

function localball(radii, rot, convention)
metric_ = anisodistance(radii, rot, convention)
MetricBall(1.0, metric_)
end
13 changes: 8 additions & 5 deletions src/parametrization/conversions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,18 +51,21 @@ function localanisotropies(

if transf
rule = convention isa Symbol ? rules[convention] : convention
rule.main == :y && (ranges = ranges[reverse(1:dim, 1, 2)])
if rule.main == :y
ranges = ranges[reverse(1:dim, 1, 2)]
rule = RotationRule(rule.order, rule.motion, rule.radian, :x, rule.extrinsic)
end
end

for i = 1:len
xranges =
[Tables.getcolumn(tab, x)[i] for x in ranges] ./ Tables.getcolumn(tab, ranges[1])[i]
vranges = [Tables.getcolumn(tab, x)[i] for x in ranges]
xranges = vranges ./ maximum(vranges)
xrot = [Tables.getcolumn(tab, x)[i] for x in rotation]

quat = if isquat
Quaternion(xrot)
else
P, Λ = rotmat(xranges, xrot, rule; rev = false)
P, Λ = rotmat(xranges, xrot, rule)
size(P, 1) == 2 && (P = DCM([P[1, 1] P[1, 2] 0; P[2, 1] P[2, 2] 0; 0 0 1]))
dcm_to_quat(P)
end
Expand Down Expand Up @@ -196,7 +199,7 @@ function rotmat2angles(dcm::AbstractMatrix, convention::RotConvention)
P = N == 2 ? DCM([P[1, 1] P[1, 2] 0; P[2, 1] P[2, 2] 0; 0 0 1]) : dcm

rule = convention isa Symbol ? rules[convention] : convention
!rule.extrinsic && (P = P')
rule.extrinsic && (P = P')
preangs = dcm_to_angle(DCM(P), rule.order)
angles = [preangs.a1, preangs.a2, preangs.a3]

Expand Down
62 changes: 56 additions & 6 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
@@ -1,18 +1,66 @@
using LocalAnisotropies
using Distances
using GeoStats
using LocalAnisotropies
using Test
import LocalAnisotropies: rotmat
import LocalAnisotropies: rotmat, anisodistance


@testset "LocalAnisotropies.jl" begin

d₁ = anisodistance([1.0, 1.0], [0.0])
d₂ = anisodistance([1.0, 2.0], [0.0])
@test evaluate(d₁, [1.0, 0.0], [0.0, 0.0]) == evaluate(d₁, [0.0, 1.0], [0.0, 0.0])
@test evaluate(d₂, [1.0, 0.0], [0.0, 0.0]) != evaluate(d₂, [0.0, 1.0], [0.0, 0.0])

d₃ = anisodistance([1.0, 0.5, 0.5], [π / 4, 0.0, 0.0])
@test evaluate(d₃, [1.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈ √2
@test evaluate(d₃, [-1.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈ √8

# intrinsic conventions
gslib = anisodistance([50.0, 25.0, 5.0], [30.0, -30.0, 30.0], :GSLIB)
tait = anisodistance([25.0, 50.0, 5.0], [-π / 6, -π / 6, π / 6], :TaitBryanIntr)
euler = anisodistance(
[50.0, 25.0, 5.0],
[-deg2rad(78), -deg2rad(41), -deg2rad(50)],
:EulerIntr,
)
lpf = anisodistance([50.0, 25.0, 5.0], [78.0, 41.0, 50.0], :Leapfrog)
dm = anisodistance([50.0, 25.0, 5.0], [78.0, 41.0, 50.0], :Datamine)

@test evaluate(gslib, [1.0, 0.0, 0.0], [0.0, 0.0, 0.0]) ≈ 0.1325707358356285
@test evaluate(gslib, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈ 0.039051248379533283
@test evaluate(gslib, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈ 0.15132745950421558

@test evaluate(gslib, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈
evaluate(tait, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0])
@test evaluate(euler, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈
evaluate(dm, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0])
@test evaluate(euler, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈
evaluate(lpf, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0])
@test evaluate(euler, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) -
evaluate(gslib, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) < 10^-3

# extrinsic conventions
xtait = anisodistance([50.0, 25.0, 5.0], [π, 0, π / 2], :TaitBryanExtr)
xeuler = anisodistance([50.0, 25.0, 5.0], [-π / 2, -π / 2, -π / 2], :EulerExtr)

@test evaluate(xtait, [1.0, 0.0, 0.0], [0.0, 0.0, 0.0]) ≈ 0.20
@test evaluate(xtait, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈ 0.04
@test evaluate(xtait, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈ 0.02

@test evaluate(xtait, [1.0, 0.0, 0.0], [0.0, 0.0, 0.0]) ≈
evaluate(xeuler, [1.0, 0.0, 0.0], [0.0, 0.0, 0.0])
@test evaluate(xtait, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈
evaluate(xeuler, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0])
@test evaluate(xtait, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈
evaluate(xeuler, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0])

# convert data to LocalAnisotropy
dummy = georef(
(az = 1:10, r1 = 1:10, r2 = 1:10),
PointSet([(i / 2, (i + 1) / 2) for i = 1:2:19]),
)
pars = localanisotropies(dummy, [:az], [:r1, :r2], :GSLIB)
@test round(pars.rotation[1][4], digits = 4) ≈ 0.0087

tabpars_quat = to_table(pars)
tabpars_angs = to_table(pars, :GSLIB)
Expand All @@ -30,9 +78,11 @@ import LocalAnisotropies: rotmat
# convert between different rotation conventions
angs1 = convertangles([30, 30, 30], :GSLIB, :Datamine)
angs2 = convertangles(angs1, :Datamine, :GSLIB)
rmat = rotmat([1, 1, 1], [30, 30, 30], :GSLIB)
@test all(rmat .≈ rotmat([1, 1, 1], angs1, :Datamine))
@test all(rmat .≈ rotmat([1, 1, 1], angs2, :GSLIB))
rmat, _ = rotmat([1, 1, 1], [30, 30, 30], :GSLIB)
dmat, _ = rotmat([1, 1, 1], angs1, :Datamine)
gmat, _ = rotmat([1, 1, 1], angs2, :GSLIB)
@test all(rmat .≈ dmat)
@test all(rmat .≈ gmat)
angs_ = convertangles.(pars.rotation, :GSLIB)
@test all([x[1] for x in angs_] .≈ 1:10)

Expand Down

2 comments on commit c2faff9

@rmcaixeta
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/120557

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.7.0 -m "<description of version>" c2faff991451c943a5c330d910ac769c8c8a0ef3
git push origin v0.7.0

Please sign in to comment.