From f5ae274ee40723ca190f6a4f7424a3b3a7ebdeee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Tue, 7 Jul 2026 09:55:38 -0300 Subject: [PATCH 1/2] Choose different factorization in the presence of missing data --- .github/workflows/Benchmarks.yml | 20 ++++++++++++++++++++ benchmark/Project.toml | 4 ++++ benchmark/benchmarks.jl | 25 +++++++++++++++++++++++++ src/krig.jl | 32 +++++++++++++++++--------------- 4 files changed, 66 insertions(+), 15 deletions(-) create mode 100644 .github/workflows/Benchmarks.yml create mode 100644 benchmark/Project.toml create mode 100644 benchmark/benchmarks.jl diff --git a/.github/workflows/Benchmarks.yml b/.github/workflows/Benchmarks.yml new file mode 100644 index 0000000..45dd538 --- /dev/null +++ b/.github/workflows/Benchmarks.yml @@ -0,0 +1,20 @@ +name: Benchmarks +on: + pull_request_target: + branches: [ main ] +permissions: + pull-requests: write + +jobs: + bench: + runs-on: ubuntu-latest + strategy: + matrix: + version: + - 'lts' + - '1' + steps: + - uses: MilesCranmer/AirspeedVelocity.jl@action-v1 + with: + julia-version: ${{ matrix.version }} + bench-on: ${{ github.event.pull_request.head.sha }} \ No newline at end of file diff --git a/benchmark/Project.toml b/benchmark/Project.toml new file mode 100644 index 0000000..62ff9a1 --- /dev/null +++ b/benchmark/Project.toml @@ -0,0 +1,4 @@ +[deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +GeoTables = "e502b557-6362-48c1-8219-d30d308dcdb0" +Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl new file mode 100644 index 0000000..11925dd --- /dev/null +++ b/benchmark/benchmarks.jl @@ -0,0 +1,25 @@ +using BenchmarkTools +using GeoStatsModels +using GeoTables +using Meshes + +# auxiliary variables +table1 = (; z=[1.0, 2.0, 3.0]) +table2 = (; z=[missing, 2.0, 3.0]) +coords = [(25.0, 25.0), (50.0, 75.0), (75.0, 50.0)] +data1 = georef(table1, coords) +data2 = georef(table2, coords) +igrid = CartesianGrid(100, 100) +model = Kriging(GaussianVariogram(range=35.0)) + +# initialize benchmark suite +const SUITE = BenchmarkGroup() + +# -------- +# KRIGING +# -------- + +SUITE["kriging"] = BenchmarkGroup() + +SUITE["kriging"]["full"] = @benchmarkable GeoStatsModels.fitpredict($model, $data1, $igrid, neighbors=true) +SUITE["kriging"]["miss"] = @benchmarkable GeoStatsModels.fitpredict($model, $data2, $igrid, neighbors=true) diff --git a/src/krig.jl b/src/krig.jl index 142e516..8b50d95 100644 --- a/src/krig.jl +++ b/src/krig.jl @@ -67,7 +67,7 @@ function fit(model::KrigingModel, data) nfun, miss = setlhs!(model, LHS, data) # factorize LHS - FHS = lhsfactorize(model, LHS) + FHS = lhsfactorize(model, LHS, miss) # record Kriging state state = KrigingState(data, LHS, RHS, FHS, nfun, miss) @@ -139,24 +139,26 @@ function setlhs!(model::KrigingModel, LHS, data) end # choose appropriate factorization of LHS -lhsfactorize(model::GeoStatsModel, LHS) = _lhsfactorize(model.fun, LHS) - -# enforce Bunch-Kaufman factorization for symmetric functions -_lhsfactorize(::Variogram, LHS) = bunchkaufman!(Symmetric(LHS), check=false) -_lhsfactorize(::Covariance, LHS) = bunchkaufman!(Symmetric(LHS), check=false) +lhsfactorize(model::GeoStatsModel, LHS, miss) = _lhsfactorize(model.fun, LHS, miss) # enforce SVD factorization for rank-deficient matrices # use QRIteration to avoid LAPACK bug: https://github.com/Reference-LAPACK/lapack/issues/672 -_lhsfactorize(::Transiogram, LHS) = svd!(LHS, alg=QRIteration()) - -# choose appropriate factorization for other functions (e.g., CompositeFunction) -function _lhsfactorize(fun::GeoStatsFunction, LHS) - if issymmetric(fun) - # enforce Bunch-Kaufman factorization - bunchkaufman!(Symmetric(LHS), check=false) +_lhsfactorize(::Transiogram, LHS, miss) = svd!(LHS, alg=QRIteration()) + +# fallback method for all geostatistical functions (e.g., CompositeFunction) +function _lhsfactorize(fun::GeoStatsFunction, LHS, miss) + if isempty(miss) + if issymmetric(fun) + # enforce Bunch-Kaufman factorization + bunchkaufman!(Symmetric(LHS), check=false) + else + # fallback to LU factorization + lu!(LHS, check=false) + end else - # fallback to LU factorization - lu!(LHS, check=false) + # fallback to SVD factorization for rank-deficient matrices + # use QRIteration to avoid LAPACK bug: https://github.com/Reference-LAPACK/lapack/issues/672 + svd!(LHS, alg=QRIteration()) end end From 6a5c16461f2811369a0ff1ee84939a672a611a36 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Tue, 7 Jul 2026 09:59:33 -0300 Subject: [PATCH 2/2] Undo changes in krig.jl --- src/krig.jl | 32 +++++++++++++++----------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/src/krig.jl b/src/krig.jl index 8b50d95..142e516 100644 --- a/src/krig.jl +++ b/src/krig.jl @@ -67,7 +67,7 @@ function fit(model::KrigingModel, data) nfun, miss = setlhs!(model, LHS, data) # factorize LHS - FHS = lhsfactorize(model, LHS, miss) + FHS = lhsfactorize(model, LHS) # record Kriging state state = KrigingState(data, LHS, RHS, FHS, nfun, miss) @@ -139,26 +139,24 @@ function setlhs!(model::KrigingModel, LHS, data) end # choose appropriate factorization of LHS -lhsfactorize(model::GeoStatsModel, LHS, miss) = _lhsfactorize(model.fun, LHS, miss) +lhsfactorize(model::GeoStatsModel, LHS) = _lhsfactorize(model.fun, LHS) + +# enforce Bunch-Kaufman factorization for symmetric functions +_lhsfactorize(::Variogram, LHS) = bunchkaufman!(Symmetric(LHS), check=false) +_lhsfactorize(::Covariance, LHS) = bunchkaufman!(Symmetric(LHS), check=false) # enforce SVD factorization for rank-deficient matrices # use QRIteration to avoid LAPACK bug: https://github.com/Reference-LAPACK/lapack/issues/672 -_lhsfactorize(::Transiogram, LHS, miss) = svd!(LHS, alg=QRIteration()) - -# fallback method for all geostatistical functions (e.g., CompositeFunction) -function _lhsfactorize(fun::GeoStatsFunction, LHS, miss) - if isempty(miss) - if issymmetric(fun) - # enforce Bunch-Kaufman factorization - bunchkaufman!(Symmetric(LHS), check=false) - else - # fallback to LU factorization - lu!(LHS, check=false) - end +_lhsfactorize(::Transiogram, LHS) = svd!(LHS, alg=QRIteration()) + +# choose appropriate factorization for other functions (e.g., CompositeFunction) +function _lhsfactorize(fun::GeoStatsFunction, LHS) + if issymmetric(fun) + # enforce Bunch-Kaufman factorization + bunchkaufman!(Symmetric(LHS), check=false) else - # fallback to SVD factorization for rank-deficient matrices - # use QRIteration to avoid LAPACK bug: https://github.com/Reference-LAPACK/lapack/issues/672 - svd!(LHS, alg=QRIteration()) + # fallback to LU factorization + lu!(LHS, check=false) end end