add lazy compute operations

Project.toml

@@ -8,6 +8,7 @@ DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
@@ -20,6 +21,7 @@ ArnoldiMethod = "0.0.4, 0.4"
 DelimitedFiles = "1"
 Graphs = "1"
 LaTeXStrings = "1.1"
+LinearSolve = "2.38.0"
 Plots = "1.4"
 ProgressLogging = "0.1"
 Rasters = "0.12"

examples/2_landmarks.jmd

@@ -99,6 +99,41 @@ 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
+ops = (;
+    func=ConScape.ConnectedHabitat(; distance_transformation=x -> exp(-x/75)),
+    qbetw=ConScape.BetweennessQweighted(),
+    kbetw=ConScape.BetweennessKweighted(),
+)
+problem = ConScape.Problem(ops; θ=theta, mode=ConScape.MatrixSolve())
+problem = ConScape.Problem(ops; θ=theta, mode=ConScape.VectorSolve())
+```
+
+Then run it for all operations on both normal and coarse grids
+
+```julia
+res = ConScape.compute(problem, g)
+res_coarse = ConScape.compute(problem, g_coarse)
+```
+
+We can plot these outputs:
+```julia
+heatmap(res.qbetw)
+heatmap(res_coarse.qbetw)
+heatmap(res.kbetw)
+heatmap(res_coarse.kbetw)
+heatmap(res.func)
+heatmap(res_coarse.func)
+```
+
+And we can check the corelation similarly to above, by getting
+layers from `res` and `res_coarse`
+```julia
+cor(filter(!isnan, res.kbetw), filter(!isnan, res_coarse.kbetw))
+```
+
 # Test landmark performance for amount of connected habitat
 
 ```julia

src/ConScape.jl

@@ -3,6 +3,7 @@ module ConScape
     using SparseArrays, LinearAlgebra
     using Graphs, Plots, SimpleWeightedGraphs, ProgressLogging, ArnoldiMethod
     using Rasters
+    using LinearSolve
     using Rasters.DimensionalData
 
     abstract type ConnectivityFunction <: Function end
@@ -26,5 +27,6 @@ module ConScape
     include("io.jl")
     # Utilities
     include("utils.jl")
-
+    include("problem.jl")
+    include("tiles.jl")
 end

src/grid.jl

@@ -193,6 +193,12 @@ function _heatmap(canvas, g; kwargs...)
     end
 end
 
+function assess(g::Grid) 
+    targetidx, targetnodes = _targetidx_and_nodes(g)
+    # Calculate memory use and expected flops for 
+    # targetnodes, or something...
+end
+
 """
     is_strongly_connected(g::Grid)::Bool
 

src/gridrsp.jl

@@ -1,3 +1,4 @@
+
 struct GridRSP
     g::Grid
     θ::Float64
@@ -12,7 +13,6 @@ end
 Construct a GridRSP from a `g::Grid` based on the inverse temperature parameter `θ::Real`.
 """
 function GridRSP(g::Grid; θ=nothing)
-
     Pref = _Pref(g.affinities)
     W    = _W(Pref, θ, g.costmatrix)
 
@@ -48,9 +48,19 @@ function Base.show(io::IO, ::MIME"text/html", grsp::GridRSP)
     write(io, "<h4>$t</h4>")
     show(io, MIME"text/html"(), plot_outdegrees(grsp.g))
 end
-
 DimensionalData.dims(grsp::GridRSP) = dims(grsp.g)
 
+abstract type AbstractOperation end 
+abstract type RSPOperation <: AbstractOperation end
+
+function keywords(o::T) where T<:AbstractOperation
+    vals = map(f -> getfield(o, f), fieldnames(T))
+    return NamedTuple{fieldnames(T)}(vals) 
+end
+
+struct BetweennessQweighted <: RSPOperation end
+
+compute(r::BetweennessQweighted, grsp::GridRSP) = betweenness_qweighted(grsp)
 
 """
     betweenness_qweighted(grsp::GridRSP)::Matrix{Float64}
@@ -76,7 +86,8 @@ function betweenness_qweighted(grsp::GridRSP)
     return _maybe_raster(bet, grsp)
 end
 
-
+struct EdgeBetweennessQweighted <: RSPOperation end
+compute(r::EdgeBetweennessQweighted, grsp::GridRSP) = edge_betweenness_qweighted(grsp)
 
 """
     edge_betweenness_qweighted(grsp::GridRSP)::Matrix{Float64}
@@ -98,6 +109,14 @@ function edge_betweenness_qweighted(grsp::GridRSP)
     return betmatrix
 end
 
+@kwdef struct BetweennessKweighted{CV,DT,DV} <: RSPOperation 
+    connectivity_function::CV=expected_cost
+    distance_transformation::DT=nothing
+    diagvalue::DV=nothing
+end
+
+compute(o::BetweennessKweighted, grsp::GridRSP) = 
+    betweenness_kweighted(grsp; keywords(o)...)
 
 """
     betweenness_kweighted(grsp::GridRSP;
@@ -153,7 +172,13 @@ function betweenness_kweighted(grsp::GridRSP;
     return _maybe_raster(bet, grsp)
 end
 
+@kwdef struct EdgeBetweennessK{CV,DT,DV} <: RSPOperation 
+    distance_transformation::DT=inv(grsp.g.costfunction)
+    diagvalue::DV=nothing
+end
+
 
+compute(r::EdgeBetweennessK, grsp::GridRSP) = edge_betweenness_kweighted(grsp; keywords(r)...)
 
 """
     edge_betweenness_kweighted(grsp::GridRSP; [distance_transformation=inv(grsp.g.costfunction), diagvalue=nothing])::SparseMatrixCSC{Float64,Int}
@@ -189,6 +214,9 @@ function edge_betweenness_kweighted(grsp::GridRSP; distance_transformation=inv(g
 end
 
 
+struct ExpectedCost <: RSPOperation end
+
+compute(r::ExpectedCost, grsp::GridRSP) = expected_cost(grsp)
 
 """
     expected_cost(grsp::GridRSP)::Matrix{Float64}
@@ -200,23 +228,48 @@ function expected_cost(grsp::GridRSP)
     return RSP_expected_cost(grsp.W, grsp.g.costmatrix, grsp.Z, targetnodes)
 end
 
+
+struct FreeDnergyDistance <: RSPOperation end
+
+compute(r::FreeDnergyDistance, grsp::GridRSP) = free_energy_distance(grsp)
+
 function free_energy_distance(grsp::GridRSP)
     targetidx, targetnodes = _targetidx_and_nodes(grsp.g)
     return RSP_free_energy_distance(grsp.Z, grsp.θ, targetnodes)
+
 end
 
+struct SurvivalProbability <: RSPOperation end
+
+compute(r::SurvivalProbability, grsp::GridRSP) = survival_probability(grsp)
+
 function survival_probability(grsp::GridRSP)
     targetidx, targetnodes = _targetidx_and_nodes(grsp.g)
+
     return RSP_survival_probability(grsp.Z, grsp.θ, targetnodes)
 end
 
+struct PowerMeanProximity <: RSPOperation end
+
+compute(r::PowerMeanProximity, grsp::GridRSP) = power_mean_proximity(grsp)
+
 function power_mean_proximity(grsp::GridRSP)
+
     targetidx, targetnodes = _targetidx_and_nodes(grsp.g)
     return RSP_power_mean_proximity(grsp.Z, grsp.θ, targetnodes)
 end
 
+struct LeastCostDistance <: RSPOperation end
+
+
+compute(r::LeastCostDistance, grsp::GridRSP) = least_cost_distance(grsp)
+
 least_cost_distance(grsp::GridRSP) = least_cost_distance(grsp.g)
 
+struct MeanKullbackLeiblerDivergence <: RSPOperation end
+
+compute(r::MeanKullbackLeiblerDivergence, grsp::GridRSP) = mean_kl_divergence(grsp)
+
 """
     mean_kl_divergence(grsp::GridRSP)::Float64
 
@@ -229,6 +282,10 @@ function mean_kl_divergence(grsp::GridRSP)
     return qs'*(RSP_free_energy_distance(grsp.Z, grsp.θ, targetnodes) - expected_cost(grsp))*qt*grsp.θ
 end
 
+struct MeanLeastCostKullbackLeiblerDivergence <: RSPOperation end
+
+compute(r::MeanLeastCostKullbackLeiblerDivergence, grsp::GridRSP) = mean_kl_divergence(grsp)
+
 
 """
     mean_lc_kl_divergence(grsp::GridRSP)::Float64
@@ -281,12 +338,20 @@ function least_cost_kl_divergence(C::SparseMatrixCSC, Pref::SparseMatrixCSC, tar
         # Pointer swap
         tmp  = from
         from = to
+
         to   = tmp
     end
 
     return kl_div
 end
 
+@kwdef struct LeastCostKullbackLeiblerDivergence <: RSPOperation 
+    target::Tuple{Int,Int}
+end
+
+compute(r::LeastCostKullbackLeiblerDivergence, grsp::GridRSP) = 
+    least_cost_kl_divergence(grsp, r.target)
+
 """
     least_cost_kl_divergence(grsp::GridRSP, target::Tuple{Int,Int})
 
@@ -305,6 +370,18 @@ function least_cost_kl_divergence(grsp::GridRSP, target::Tuple{Int,Int})
     return reshape(div, grsp.g.nrows, grsp.g.ncols)
 end
 
+
+@kwdef struct ConnectedHabitat{CV,DT,DV} <: RSPOperation
+    # TODO not sure which kw to use here
+    connectivity_function::CV=expected_cost
+    distance_transformation::DT=nothing
+    diagvalue::DV=nothing
+    θ::Union{Nothing,Real}=nothing
+    approx::Bool=false
+end
+
+compute(r::ConnectedHabitat, grsp::GridRSP) = connected_habitat(grsp; keywords(r)...)
+
 """
     connected_habitat(grsp::Union{Grid,GridRSP};
         connectivity_function=expected_cost,
@@ -436,9 +513,19 @@ function connected_habitat(grsp::GridRSP,
 
     newh = GridRSP(newg, θ=grsp.θ)
 
+
     return connected_habitat(newh; diagvalue=diagvalue, distance_transformation=distance_transformation)
 end
 
+@kwdef struct EigMax{F,DT,DV,T} <: RSPOperation 
+    connectivity_function::F=expected_cost
+    Tdistance_transformation::DT=nothing
+    diagvalue::DV=nothing
+    tol::T=1e-14
+end
+
+compute(r::EigMax, grsp::GridRSP) = eigmax(grsp; keywords(r)...)
+
 """
     eigmax(grsp::GridRSP;
         connectivity_function=expected_cost,
@@ -564,6 +651,16 @@ function LinearAlgebra.eigmax(grsp::GridRSP;
     return vˡ, λ₀[1], vʳ
 end
 
+@kwdef struct Criticality{DT,DV,AV,QT,QS} <: RSPOperation 
+    distance_transformation::DT=inv(grsp.g.costfunction)
+    diagvalue::DV=nothing
+    avalue::AV=floatmin()
+    qˢvalue::QS=0.0
+    qᵗvalue::QT=0.0
+end
+
+compute(r::Criticality, grsp::GridRSP) = criticality(grsp; keywords(r)...)
+
 """
     criticality(grsp::GridRSP[;
                 distance_transformation=inv(grsp.g.costfunction),
@@ -605,3 +702,4 @@ function criticality(grsp::GridRSP;
 
     return _maybe_raster(landscape, grsp)
 end
+