So I started trying to reverse engineer this comparison done in R that performs a very close implementation of INLA from scratch. Though, I’m having difficulty reproducing the results of the “from-scratch” implementation due to differences in how Julia can perform computations dealing with infinite values versus R.

Here’s the link to my example.

What I’m trying to reproduce is essentially this example done in R

Where I am trying to just bring in the res object into Julia and then use it for comparison against the julia implementation.

I could definitely use a second set of eyes with this but if we can reproduce the example done in R. I think we can have a good standing point for what to implement for an INLA.jl package. We just need to get the Julia version scratch going and then we can build from there.

Very interesting tutorial on INLA @bparbhu , thanks for sharing. I didn’t have time to read it carefully, but INLA is something that we are pursuing in our GeoStats.jl stack for spatial processes. We have implemented the Laplacian operators for arbitrary meshes already and the underlying Gaussian solvers for high-dimensional simulation. Any help is appreciated

Thanks again for the feedback I appreciate it! I didn’t know that JuliaEarth was working on an INLA solution and its very cool that you guys were able to replicate the Laplacian operators, that’s awesome. If that’s the case. Would you want to move this under JuliaEarth or do you have a branch in GeoStats.jl that we can work this out. Also, curious just about identifying differences in implementation through this kind of exercise as well. So if you have an implementation this can also serve as a test to see how close it is to R-INLA. So happy to help, just want to finish this up for my own sake.

using GeoStats
using GeoStatsViz
import GLMakie as Mke
mesh = simplexify(Sphere((0,0,0), 1))
# ask for 3 realizations of the process over the sphere
problem = SimulationProblem(mesh, :Z=>Float64, 3)
# latent process of INLA
solver = SPDEGS()
# perform simulation
ensemble = solve(problem, solver)
# view first realization
viewer(ensemble[1])