Add LeastCost and RandomWalk MovementModes

Project.toml

@@ -4,30 +4,41 @@ version = "0.3.0"
 
 [deps]
 ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
+ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
+[weakdeps]
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+
+[extensions]
+ConScapeLinearSolveExt = "LinearSolve"
+
 [compat]
 ArnoldiMethod = "0.0.4, 0.4"
+CommonSolve = "0.2"
+ConstructionBase = "1.5.8"
 DelimitedFiles = "1"
 Graphs = "1"
 LaTeXStrings = "1.1"
+LinearSolve = "2.38.0"
 Plots = "1.4"
 ProgressLogging = "0.1"
-Rasters = "0.12"
+Rasters = "0.14"
 SimpleWeightedGraphs = "1.1"
 julia = "1.10"
 
 [extras]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [targets]
-test = ["ArchGDAL"]
+test = ["ArchGDAL","Plots"]

examples/2_landmarks.jmd

@@ -1,6 +1,7 @@
 # Landmark function for testing purposes
 
 ```julia
+nothing
 using ConScape, Rasters, ArchGDAL, SparseArrays, LinearAlgebra, Plots
 ```
 
@@ -83,7 +84,7 @@ heatmap(qbetw)
 
 ```julia
 qbetw_coarse = ConScape.betweenness_qweighted(h_coarse);
-ConScape.heatmap(qbetw_coarse)
+heatmap(qbetw_coarse)
 ```
 
 ```julia
@@ -99,6 +100,80 @@ kbetw_coarse = @time ConScape.betweenness_kweighted(h_coarse, distance_transform
 cor(filter(!isnan, kbetw), filter(!isnan, kbetw_coarse))
 ```
 
+We can write this using a lazy problem definition:
+
+```julia
+# Put least cost, random walk, and rsp
+problem = ConScape.Problem(; 
+    graph_measures = (;
+        func=ConScape.ConnectedHabitat(),
+        qbetw=ConScape.BetweennessQweighted(),
+        kbetw=ConScape.BetweennessKweighted(),
+    ),
+    connectivity_measure=ConScape.ExpectedCost(
+        θ=1.0, 
+        distance_transformation = (exp=x -> exp(-x/75), oddsfor=ConScape.OddsFor())
+    ),
+    solver=ConScape.MatrixSolver(),
+    # solver=ConScape.VectorSolver(nothing),
+)
+```
+
+Then run it for all operations on both normal and coarse grids
+
+```julia
+rast = RasterStack((; affinities=mov_prob, qualities=hab_qual))
+rast_coarse = ConScape.coarse_graining(rast, 20)
+result = ConScape.solve(problem, rast)
+result_coarse = ConScape.solve(problem, rast_coarse)
+```
+
+We can plot these outputs:
+```julia
+heatmap(result.qbetw_oddsfor)
+heatmap(result_coarse.qbetw_oddsfor)
+heatmap(result.qbetw_exp)
+heatmap(result_coarse.qbetw_exp)
+heatmap(result.kbetw_exp)
+heatmap(result_coarse.kbetw_exp)
+heatmap(result.func_exp)
+heatmap(result.coarse_func_exp)
+```
+
+
+```julia
+windowed_problem = ConScape.WindowedProblem(problem; 
+    radius=20, overlap=5, threaded=true
+)
+@time result = ConScape.solve(windowed_problem, rast)
+plot(result)
+plot(result.func_exp)
+```
+
+
+```julia
+stored_problem = ConScape.BatchProblem(problem; 
+    path=".", radius=20, overlap=30, threaded=true
+)
+ConScape.solve(stored_problem, rast)
+result = mosaic(stored_problem; to=rast)
+plot(result)
+plot(result.func_exp)
+
+# Calculate job ids to run on a cluster
+jobs = ConScape.batch_ids(stored_problem, rast)
+```
+
+And we can check the corelation similarly to above, by getting
+layers from `res` and `res_coarse`
+```julia
+cor(filter(!isnan, res.qbetw.oddsfor), filter(!isnan, res_coarse.qbetw.oddsfor))
+cor(filter(!isnan, res.qbetw.oddsfor), filter(!isnan, res_coarse.qbetw.oddsfor))
+cor(filter(!isnan, res.kbetw.oddsfor), filter(!isnan, res_coarse.kbetw.oddsfor))
+cor(filter(!isnan, res.kbetw.exp), filter(!isnan, res_coarse.kbetw.exp))
+cor(filter(!isnan, res.func.exp), filter(!isnan, res_coarse.func.exp))
+```
+
 # Test landmark performance for amount of connected habitat
 
 ```julia

ext/ConScapeCUDSS.jl

@@ -0,0 +1,37 @@
+module ConScapeCUDSS
+
+using ConScape
+using CUDA
+using CUDA.CUSPARSE
+using CUDSS
+using SparseArrays
+using LinearAlgebra
+
+using ConScape: FundamentalMeasure, AbstractProblem, Grid, setup_sparse_problem
+
+struct CUDSSsolver <: Solver end
+
+function ConScape.solve(m::CUDSSsolver, cm::FundamentalMeasure, p::AbstractProblem, g::Grid) 
+    (; A, B, Pref, W) = setup_sparse_problem(g, cm)
+    Z = zeros(T, size(B))
+
+    A_gpu = CuSparseMatrixCSR(A |> tril)
+    Z_gpu = CuMatrix(Z)
+    B_gpu = CuMatrix(B)
+
+    solver = CudssSolver(A_gpu, "S", "L")
+
+    cudss("analysis", solver, Z_gpu, B_gpu)
+    cudss("factorization", solver, Z_gpu, B_gpu)
+    cudss("solve", solver, Z_gpu, B_gpu)
+
+    Z .= Z_gpu
+    # TODO: maybe graph measures can run on GPU as well?
+    grsp = GridRSP(g, cm.θ, Pref, W, Z)
+    results = map(p.graph_measures) do gm
+        compute(gm, p, grsp)
+    end 
+    return _merge_to_stack(results)
+end
+
+end

ext/ConScapeLinearSolveExt.jl

@@ -0,0 +1,53 @@
+module ConScapeLinearSolveExt
+
+using ConScape
+using LinearSolve
+import CommonSolve
+
+function CommonSolve.init(solver::LinearSolver, A::AbstractMatrix)
+    b = zeros(eltype(A), size(A, 2))
+    # Define and initialise the linear problem
+    linprob = LinearProblem(A, b)
+    linsolve = init(linprob, solver.args...; solver.keywords...)
+    # TODO what is needed here?
+    # Create a channel to store problem b vectors for threads
+    # see https://juliafolds2.github.io/OhMyThreads.jl/stable/literate/tls/tls/
+    if isthreaded(solver)
+        nbuffers = Threads.nthreads()
+        channel = Channel{Tuple{typeof(linsolve),Vector{Float64}}}(nbuffers)
+        for i in 1:nbuffers
+            # TODO fix this in LinearSolve.jl with batching
+            # We should not need to `deepcopy` the whole problem we 
+            # just need to replicate the specific workspace arrays 
+            # that will cause race conditions.
+            # But currently there is no parallel mode for LinearSolve.jl
+            # See https://github.com/SciML/LinearSolve.jl/issues/552
+            put!(channel, (deepcopy(linsolve), Vector{eltype(A)}(undef, size(A, 2))))
+        end
+        return channel
+    else
+        return linsolve
+    end
+end
+
+function LinearAlgebra.ldiv!(s::LinearSolver, B, init, B_copy)
+    # TODO: for now we define a Z matrix, but later modify ops 
+    # to run column by column without materialising Z
+    if isthreaded(s)
+        channel = init
+        # Get column memory from the channel
+        linsolve = take!(channel)
+        # Update solver with new b values
+        reinit!(linsolve; b=vec(B_copy), reuse_precs=true)
+        sol = LinearSolve.solve!(vec(B), linsolve, s.args...; s.keywords...)
+        vec(B) .= sol.u
+        put!(channel, linsolve)
+    else
+        linsolve = init
+        reinit!(linsolve; b=vec(B_copy), reuse_precs=true)
+        sol = LinearSolve.solve(linsolve, s.args...; s.keywords...)
+        vec(B) .= sol.u
+    end
+    return B
+end
+end

src/ConScape.jl

@@ -1,30 +1,59 @@
 module ConScape
 
-    using SparseArrays, LinearAlgebra
-    using Graphs, Plots, SimpleWeightedGraphs, ProgressLogging, ArnoldiMethod
-    using Rasters
-    using Rasters.DimensionalData
-
-    abstract type ConnectivityFunction <: Function end
-    abstract type DistanceFunction <: ConnectivityFunction end
-    abstract type ProximityFunction <: ConnectivityFunction end
-
-    struct least_cost_distance   <: DistanceFunction end
-    struct expected_cost         <: DistanceFunction end
-    struct free_energy_distance  <: DistanceFunction end
-
-    struct survival_probability  <: ProximityFunction end
-    struct power_mean_proximity  <: ProximityFunction end
-
-    # Randomized shortest path algorithms
-    include("randomizedshortestpath.jl")
-    # Grid struct and methods
-    include("grid.jl")
-    # GridRSP (randomized shortest path) struct and methods
-    include("gridrsp.jl")
-    # IO
-    include("io.jl")
-    # Utilities
-    include("utils.jl")
+using ArnoldiMethod
+using ConstructionBase
+using Graphs
+using LinearAlgebra
+using Rasters
+using SimpleWeightedGraphs
+using SparseArrays
+using Rasters.DimensionalData
+
+import CommonSolve
+import CommonSolve: solve, init
+
+export RandomisedShortestPath, LeastCost, RandomWalk
+
+export ExpectedCost, FreeEnergyDistance, SurvivalProbability, PowerMeanProximity, KullbackLeiblerDivergence
+
+export Betweenness, EdgeBetweenness, ConnectedHabitat, Criticality, EigMax, Sensitivity
+
+export QualityWeighted, QualityAndProximityWeighted, ProximityWeighted
+
+export VectorSolver, LinearSolver
+
+export MinusLog, MinusLogAlpha, Inv, OddsFor, OddsAgainst, ExpMinus, ExpMinusAlpha
+
+export solve, init, assess
+
+export WindowedProblem, BatchProblem
+
+"""
+    Solver
+
+Abstract supertype for ConScape solvers.
+"""
+abstract type AbstractProblem end
+abstract type Solver end
+
+# Randomized shortest path algorithms
+# Grid struct and methods
+include("transformations.jl")
+# Grid struct and methods
+include("grid.jl")
+# Utilities
+include("utils.jl")
+# Problems
+include("workspaces.jl")
+include("measures.jl")
+include("movement_modes.jl")
+include("problem.jl")
+include("initialisation.jl")
+include("return.jl")
+include("solvers.jl")
+include("compute_measures.jl")
+include("windows.jl")
+include("assessment.jl")
+include("simulations.jl")
 
 end