Add moving_window function

Project.toml

@@ -13,6 +13,7 @@ ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
 SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 
 [compat]
 ArnoldiMethod = "0.0.4"

src/ConScape.jl

@@ -2,6 +2,7 @@ module ConScape
 
     using SparseArrays, LinearAlgebra
     using LightGraphs, Plots, SimpleWeightedGraphs, ProgressLogging, ArnoldiMethod
+	using ThreadsX
 
     abstract type ConnectivityFunction <: Function end
     abstract type DistanceFunction <: ConnectivityFunction end

src/grid.jl

@@ -36,7 +36,8 @@ end
               source_qualities::Matrix=qualities,
               target_qualities::AbstractMatrix=qualities,
               costs::Union{Transformation,SparseMatrixCSC{Float64,Int}}=MinusLog(),
-              prune=true)::Grid
+              prune=true,
+			  verbose=true)::Grid
 
 Construct a `Grid` from an `affinities` matrix of type `SparseMatrixCSC`. It is possible
 to also supply matrices of `source_qualities` and `target_qualities` as well as
@@ -51,7 +52,8 @@ function Grid(nrows::Integer,
               source_qualities::Matrix=qualities,
               target_qualities::AbstractMatrix=qualities,
               costs::Union{Transformation,SparseMatrixCSC{Float64,Int}}=MinusLog(),
-              prune=true)
+              prune=true,
+			  verbose=true)
 
     if affinities === nothing
         throw(ArgumentError("matrix of affinities must be supplied"))
@@ -106,7 +108,7 @@ function Grid(nrows::Integer,
     )
 
     if prune
-        return largest_subgraph(g)
+        return largest_subgraph(g, verbose=verbose)
     else
         return g
     end
@@ -199,20 +201,22 @@ false
 LightGraphs.is_strongly_connected(g::Grid) = is_strongly_connected(SimpleWeightedDiGraph(g.affinities))
 
 """
-    largest_subgraph(g::Grid)::Grid
+    largest_subgraph(g::Grid; verbose::Bool=true)::Grid
 
 Extract the largest fully connected subgraph of the `Grid`. The returned `Grid`
 will have the same size as the input `Grid` but only nodes associated with the
 largest subgraph of the affinities will be active.
 """
-function largest_subgraph(g::Grid)
+function largest_subgraph(g::Grid; verbose::Bool=true)
     # Convert cost matrix to graph
     graph = SimpleWeightedDiGraph(g.costmatrix, permute=false)
 
     # Find the subgraphs
     scc = strongly_connected_components(graph)
 
-    @info "cost graph contains $(length(scc)) strongly connected subgraphs"
+    if verbose
+		@info "cost graph contains $(length(scc)) strongly connected subgraphs"
+	end
 
     # Find the largest subgraph
     i = argmax(length.(scc))
@@ -221,10 +225,13 @@ function largest_subgraph(g::Grid)
     scci = sort(scc[i])
 
     ndiffnodes = size(g.costmatrix, 1) - length(scci)
-    if ndiffnodes > 0
-        @info "removing $ndiffnodes nodes from affinity and cost graphs"
-    end
-
+
+	if verbose
+		if ndiffnodes > 0
+			@info "removing $ndiffnodes nodes from affinity and cost graphs"
+		end
+	end
+
     # Extract the adjacency matrix of the largest subgraph
     affinities = g.affinities[scci, scci]
     # affinities = convert(SparseMatrixCSC{Float64,Int}, graph[scci])

src/utils.jl

@@ -190,3 +190,214 @@ function mapnz(f, A::SparseMatrixCSC)
     map!(f, B.nzval, A.nzval)
     return B
 end
+
+"""
+    make_landscape(pth::String; q::String, a::String, c::Union{String,Nothing})::Array
+
+Read asc files of raster data from a directory, into a 3D Array of quality, affinity, and optionally cost.
+The moving_window function takes the output from this function.
+"""
+
+function make_landscape(pth::String; q::String, a::String, c::Union{String,Nothing}=nothing)
+    mov_prob, meta_p = readasc(joinpath(pth, a))
+    hab_qual, meta_q = readasc(joinpath(pth, q))
+
+    if collect(values(meta_p))[1:end .!= 3] != collect(values(meta_q))[1:end .!= 3]
+        throw(ArgumentError("Maps of quality and permeability do not match"))
+    end
+
+    if c === nothing
+        lnd = reshape([hab_qual mov_prob], size(hab_qual)[1], size(hab_qual)[2], :)
+    else
+        mov_cost, meta_c = readasc(joinpath(pth, c))
+        if collect(values(meta_p))[1:end .!= 3] != collect(values(meta_c))[1:end .!= 3]
+            throw(ArgumentError("The cost map does not match with the other maps"))
+        end
+        lnd = reshape([hab_qual mov_prob mov_cost], size(hab_qual)[1], size(hab_qual)[2], :)
+    end    
+
+    return lnd
+end
+
+"""
+    moving_window_helper(win::Tuple{Int64, Int64}; dwin::Dict, win_cntr::Matrix, lnd::Array, cst_fun::Union{Transformation,Nothing}, fnc::String, args::Dict)::SparseMatrixCSC
+
+The serial version of the helper function. The function clips the window from the landscape and computes the desired function. The result is returned as a SparseMatrix, which are combined in the main moving_window function.  
+"""
+
+function moving_window_helper(win::Tuple{Int64, Int64}; dwin::Dict, win_cntr::Matrix, lnd::Array, cst_fun::Union{Transformation,Nothing}, fnc::String, args::Dict)
+    res = zeros(size(lnd)[1:2])
+
+    tgt_lms = [[((win[1]-1)*dwin["d"])+1, minimum([(win[1]*dwin["d"]), size(lnd)[1]])], 
+        [((win[2]-1)*dwin["d"])+1, minimum([(win[2]*dwin["d"]), size(lnd)[2]])]]
+    src_lms = [[maximum([tgt_lms[1][1]-dwin["buf"],1]), minimum([tgt_lms[1][2]+dwin["buf"], size(lnd)[1]])], 
+        [maximum([tgt_lms[2][1]-dwin["buf"],1]), minimum([tgt_lms[2][2]+dwin["buf"], size(lnd)[2]])]]
+
+    lnd_clp = lnd[src_lms[1][1]:src_lms[1][2],src_lms[2][1]:src_lms[2][2],:]
+    tgts = zeros(size(lnd_clp[:,:,1]))
+    tgts[(dwin["buf"]+1):minimum([dwin["buf"]+dwin["d"], size(tgts)[1]]),
+                        (dwin["buf"]+1):minimum([dwin["buf"]+dwin["d"], size(tgts)[2]])] = 
+                    map(t -> isnan(t) ? 0.0 : t, lnd[tgt_lms[1][1]:tgt_lms[1][2],tgt_lms[2][1]:tgt_lms[2][2],1]) .* 
+                                        win_cntr[1:(tgt_lms[1][2]-tgt_lms[1][1]+1),1:(tgt_lms[2][2]-tgt_lms[2][1]+1)]                                                            
+    tgts = dropzeros(sparse(tgts));
+
+    #skip empty affinities clips
+    if !((sum(map!(x -> isnan(x) ? 0 : x, lnd_clp[:,:,1], lnd_clp[:,:,1])) == 0) || 
+        (sum(map!(x -> isnan(x) ? 0 : x, lnd_clp[:,:,2], lnd_clp[:,:,2])) == 0) || 
+        length(nonzeros(tgts)) ==0)
+
+        adjacency_matrix = graph_matrix_from_raster(lnd_clp[:,:,2]);
+        # compute costs, if missing
+        if (size(lnd_clp)[3]==2)    
+            cost_matrix = mapnz(cst_fun, adjacency_matrix)
+        else 
+            cost_matrix = graph_matrix_from_raster(lnd_clp[:,:,3]);
+        end
+
+
+        h = GridRSP(Grid(size(lnd_clp[:,:,2])...,
+                    affinities=adjacency_matrix,
+                    source_qualities=lnd_clp[:,:,1],
+                    target_qualities=tgts,
+                    costs=cost_matrix, prune=true, verbose=false), θ=get(args, "θ", 1))
+
+        if fnc === "connected_habitat"
+            tmp = connected_habitat(h, 
+                connectivity_function=get(args, "connectivity_function", expected_cost), 
+                distance_transformation=get(args, "connectivity_function", ExpMinus()));
+        elseif fnc === "betweenness_kweighted"
+            tmp = betweenness_kweighted(h, 
+                connectivity_function=get(args, "connectivity_function", expected_cost), 
+                distance_transformation=get(args, "connectivity_function", ExpMinus()));
+        elseif fnc === "betweenness_qweighted"
+            tmp = betweenness_qweighted(h);
+        else 
+            throw(ArgumentError("Currently only the connected_habitat and both betweenness functions are supported"))
+        end
+
+        # return the result as a sparce matrix
+        res[src_lms[1][1]:src_lms[1][2],src_lms[2][1]:src_lms[2][2]] = map(t -> isnan(t) ? 0.0 : t, tmp)
+    end
+    res = sparse(res)
+    return(res)
+end
+
+"""
+    clip_landscape(win::Tuple{Int64, Int64}; dwin::Dict, win_cntr::Matrix, lnd::Array)::Dict
+
+Data preparation function for the parallel version of the moving_window_helper. It clips a window (win) from the landscape (lnd), which is returned as a Dictionary.
+The parameters for the window computation are provided in the window dictionary (dwin).
+The dictionary with the clipped landscape contains the boundaries or limits (i.e. upper and lower rows and columns) for the targets (tgt_lms) and sources (src_lms) together with the clipped landscape array (lnd_clp) and a sparce targets matrix (tgts). The targets matrix is the same size as the clipped landscape, but has only non-zero values on the non-zero values in the window center. 
+"""
+
+function clip_landscape(win::Tuple{Int64, Int64}; dwin::Dict, win_cntr::Matrix, lnd::Array)
+    tgt_lms = [[((win[1]-1)*dwin["d"])+1, minimum([(win[1]*dwin["d"]), size(lnd)[1]])], 
+        [((win[2]-1)*dwin["d"])+1, minimum([(win[2]*dwin["d"]), size(lnd)[2]])]]
+    src_lms = [[maximum([tgt_lms[1][1]-dwin["buf"],1]), minimum([tgt_lms[1][2]+dwin["buf"], size(lnd)[1]])], 
+        [maximum([tgt_lms[2][1]-dwin["buf"],1]), minimum([tgt_lms[2][2]+dwin["buf"], size(lnd)[2]])]]
+
+    lnd_clp = lnd[src_lms[1][1]:src_lms[1][2],src_lms[2][1]:src_lms[2][2],:]
+    tgts = zeros(size(lnd_clp[:,:,1]))
+    tgts[(dwin["buf"]+1):minimum([dwin["buf"]+dwin["d"], size(tgts)[1]]),
+                        (dwin["buf"]+1):minimum([dwin["buf"]+dwin["d"], size(tgts)[2]])] = 
+                    map(t -> isnan(t) ? 0.0 : t, lnd[tgt_lms[1][1]:tgt_lms[1][2],tgt_lms[2][1]:tgt_lms[2][2],1]) .* 
+                                        win_cntr[1:(tgt_lms[1][2]-tgt_lms[1][1]+1),1:(tgt_lms[2][2]-tgt_lms[2][1]+1)]                                                            
+    tgts = dropzeros(sparse(tgts));
+    res = Dict("tgt_lms" => tgt_lms, "src_lms" => src_lms, "clp" => lnd_clp, "tgts" => tgts)
+    return(res)
+end
+
+"""
+    moving_window_helper(lnd_clp::Dict; cst_fun::Union{Transformation,Nothing}, fnc::String, args::Dict, lnd_sz::Tuple{Int64, Int64})::SparseMatrixCSC
+
+Helper function for the parallel version of the moving_window, it takes the dictionary of the clipped landscape (lnd_clp) as the main input (this dictionary is the output from the clip_landscape function). 
+The size of the original landscape is provided as a Tupple (lnd_sz) to return the results to the main function. 
+"""
+function moving_window_helper(lnd_clp::Dict; cst_fun::Union{Transformation,Nothing}, fnc::String, args::Dict, lnd_sz::Tuple{Int64, Int64})
+    res = zeros(lnd_sz)
+
+    #skip empty affinities clips
+    if !((sum(map!(x -> isnan(x) ? 0 : x, lnd_clp["clp"][:,:,1], lnd_clp["clp"][:,:,1])) == 0) || 
+        (sum(map!(x -> isnan(x) ? 0 : x, lnd_clp["clp"][:,:,2], lnd_clp["clp"][:,:,2])) == 0) || 
+        length(nonzeros(lnd_clp["tgts"])) ==0)
+
+        adjacency_matrix = graph_matrix_from_raster(lnd_clp["clp"][:,:,2]);
+        # compute costs, if missing
+        if (size(lnd_clp["clp"])[3]==2)    
+            cost_matrix = mapnz(cst_fun, adjacency_matrix)
+        else 
+            cost_matrix = graph_matrix_from_raster(lnd_clp["clp"][:,:,3]);
+        end
+
+
+        h = GridRSP(Grid(size(lnd_clp["clp"][:,:,2])...,
+                    affinities=adjacency_matrix,
+                    source_qualities=lnd_clp["clp"][:,:,1],
+                    target_qualities=lnd_clp["tgts"],
+                    costs=cost_matrix, prune=true, verbose=false), θ=get(args, "θ", 1))
+
+        if fnc === "connected_habitat"
+            tmp = connected_habitat(h, 
+                connectivity_function=get(args, "connectivity_function", expected_cost), 
+                distance_transformation=get(args, "connectivity_function", ExpMinus()));
+        elseif fnc === "betweenness_kweighted"
+            tmp = betweenness_kweighted(h, 
+                connectivity_function=get(args, "connectivity_function", expected_cost), 
+                distance_transformation=get(args, "connectivity_function", ExpMinus()));
+        elseif fnc === "betweenness_qweighted"
+            tmp = betweenness_qweighted(h);
+        else 
+            throw(ArgumentError("Currently only the connected_habitat and both betweenness functions are supported"))
+        end
+
+        # return the result as a sparce matrix
+        res[lnd_clp["src_lms"][1][1]:lnd_clp["src_lms"][1][2],lnd_clp["src_lms"][2][1]:lnd_clp["src_lms"][2][2]] = 
+                                                                                        map(t -> isnan(t) ? 0.0 : t, tmp)
+    end
+    res = sparse(res)
+    return(res)
+end
+
+"""
+    moving_window(win_buf::Integer; win_cntr:::Matrix, lnd::Array, cst_fun::Union{Transformation,Nothing}, fnc::String, args::Dict, parallel::bool)::Array
+
+moving_window function to run connected_habitat or betweenness functions (fnc) as a moving window over the landscape (lnd, an output from make_landscape). 
+It relies on three helper functions:
+1) make_landscape: which creates an Array with a quality, affinity, and optionally a cost layer
+2) clip_landscape: which clips a window from the landscape and returns it as a dictionary
+3) moving_window_helper: which runs the actual computations for each window.
+
+The main moving_window function creates a window dictionary (dwin) with the window buffer (buf), the distance (d) between windows, a tupple (n) with the number of windows row- and columnwise, and the size (sz) of the window.
+Windows are identified as a tupple for the row and column direction. It dispatches the computation to the moving_window_helper based on whether computation should be parallelized or not. If parallel then the landscape will first be clipped to reduce memory useage and these clips will be processed in parallel. If serial (parallel=false) then the clipping will occur within the helper function. Thus, dispatch occurs based on whether the main argument is a dictionary with landscape clips or a tupple with the window identification. The helpers return a sparse array, which is then squashed along the window axis to return the final results of the computation.
+
+The moving window requires a center (win_cntr), which is a square matrix of zeros and ones, denoting the targets. The size of the window is given by the size of the win_cntr and the padding or buffer around the center: win_buf. The arguments for the fnc are to be provided in args. Finally, the function can be ran in parallel.
+"""
+
+function moving_window(win_buf::Integer=5; win_cntr::Matrix=[1], lnd::Array, cst_fun::Union{Transformation,Nothing}=nothing, fnc::String, args::Dict, parallel::Bool=false)
+    if length(size(win_cntr)) === 2
+        if size(win_cntr)[1] != size(win_cntr)[2]
+            throw(ArgumentError("Only square windows and window centers are currently supported"))
+        end
+    end
+
+    dwin = Dict("buf" => win_buf, "d" => size(win_cntr)[1], "n" => Int.(ceil.(size(lnd)[1:2] ./ size(win_cntr)[1])))
+    dwin["sz"] = (dwin["buf"]*2+dwin["d"])
+    wins = vec(collect(Iterators.product(1:dwin["n"][1], 1:dwin["n"][2])))
+
+    if !(fnc in ["connected_habitat", "betweenness_qweighted", "betweenness_kweighted"])
+        throw(ArgumentError("Currently only the connected_habitat and betweenness functions are supported"))
+    end
+
+    if parallel
+        lnd_clps = [clip_landscape(win; dwin=dwin, win_cntr=win_cntr, lnd=lnd) for win in wins]
+        res = ThreadsX.map(lnd_clp -> moving_window_helper(lnd_clp; cst_fun=cst_fun, fnc=fnc, args=args, lnd_sz=size(lnd)[1:2]), lnd_clps)
+    else
+        #lnd_clps = [clip_landscape(win; dwin=dwin, win_cntr=win_cntr, lnd=lnd) for win in wins]
+        #res = [moving_window_helper(lnd_clp; cst_fun=cst_fun, fnc=fnc, args=args) for lnd_clp in lnd_clps]
+        res = [moving_window_helper(win; dwin=dwin, win_cntr=win_cntr, lnd=lnd, cst_fun=cst_fun, fnc=fnc, args=args) for win in wins]
+    end
+
+    res = reshape(vcat(res...), size(lnd)[1], :, size(lnd)[2])
+    res = sum(res, dims=2)[:,1,:];
+    return(res)
+end