Merge pull request #667 from open-AIMS/tp-conn

Clean up use of TP matrices to align with conceptualisation of connectivity

src/ExtInterface/ADRIA/Domain.jl

@@ -14,12 +14,12 @@ mutable struct ADRIADomain{
     RCP::String  # RCP scenario represented
     env_layer_md::EnvLayer  # Layers used
     scenario_invoke_time::String  # time latest set of scenarios were run
-    const TP_data::Σ  # site connectivity data
+    const conn::Σ  # connectivity data
     const in_conn::Vector{Float64}  # sites ranked by incoming connectivity strength (i.e., number of incoming connections)
     const out_conn::Vector{Float64}  # sites ranked by outgoing connectivity strength (i.e., number of outgoing connections)
     const strong_pred::Vector{Int64}  # strongest predecessor
     site_data::D  # table of site data (depth, carrying capacity, etc)
-    const site_id_col::String  # column to use as site ids, also used by the connectivity dataset (indicates order of `TP_data`)
+    const site_id_col::String  # column to use as site ids, also used by the connectivity dataset (indicates order of `conn`)
     const unique_site_id_col::String  # column of unique site ids
     init_coral_cover::M  # initial coral cover dataset
     const coral_growth::CoralGrowth  # coral

src/ExtInterface/ReefMod/Domain.jl

@@ -11,7 +11,7 @@ mutable struct ReefModDomain <: Domain
     RCP::String
     env_layer_md
     scenario_invoke_time::String  # time latest set of scenarios were run
-    const TP_data
+    const conn
     const in_conn
     const out_conn
     const strong_pred

src/decision/dMCDA.jl

@@ -106,7 +106,7 @@ function DMCDA_vars(
         domain.sim_constants.priority_zones,
         site_d.zone_type,
         domain.strong_pred,
-        domain.TP_data .* site_k_area(domain),
+        domain.conn .* site_k_area(domain),
         waves,
         dhws,
         site_d.depth_med,

src/decision/location_selection.jl

@@ -32,9 +32,9 @@ function _location_selection(
     # Determine connectivity strength
     # Account for cases where no coral cover
     in_conn, out_conn, strong_pred = connectivity_strength(
-        domain.TP_data .* site_k_area(domain),
+        domain.conn .* site_k_area(domain),
         vec(sum_cover),
-        similar(domain.TP_data)  # tmp cache matrix
+        similar(domain.conn)  # tmp cache matrix
     )
 
     # Perform location selection for seeding and shading.

src/ecosystem/connectivity.jl

@@ -26,9 +26,10 @@ NOTE: Transposes transitional probability matrix if `swap == true`
 
 # Returns
 NamedTuple:
-- `TP_data` : Matrix, containing the transition probability for all sites
-- `truncated` : ID of sites removed
-- `site_ids` : ID of sites kept
+- `conn` : Matrix, containing the connectivity for all locations.
+           Accounts for larvae which do not settle, so rows are not required to sum to 1
+- `truncated` : ID of locations removed
+- `site_ids` : ID of locations kept
 """
 function site_connectivity(
     file_loc::String,
@@ -41,16 +42,13 @@ function site_connectivity(
         error("Could not find location: $(file_loc)")
     end
 
-    local extracted_TP::Matrix{Float64}
-    if isfile(file_loc)
-        con_files::Vector{String} = String[file_loc]
-    elseif isdir(file_loc)
-        # Get connectivity years available in data store
-        years::Vector{String} = getindex(first(walkdir(file_loc)), 2)
-        year_conn_fns = NamedTuple{Tuple(Symbol.(years))}([
-            [joinpath.(first(fl), last(fl)) for fl in walkdir(joinpath(file_loc, yr))][1]
-            for yr in years
-        ])
+    local extracted_conn::Matrix{Float64}
+    if isfile(file_path)
+        conn_files::Vector{String} = String[file_path]
+        first_file = conn_files[1]
+    elseif isdir(file_path)
+        conn_fns = readdir(file_path)
+        conn_fns = String[fn for fn in conn_fns if endswith(fn, ".csv")]
 
         con_files = vcat([x for x in values(year_conn_fns)]...)
 
@@ -132,12 +130,12 @@ function site_connectivity(
         extracted_TP[extracted_TP .< con_cutoff] .= 0.0
     end
 
-    TP_base = NamedDimsArray(
-        extracted_TP; Source=unique_site_ids, Sink=unique_site_ids
+    conn = NamedDimsArray(
+        extracted_conn; Source=unique_site_ids, Sink=unique_site_ids
     )
-    @assert all(0.0 .<= TP_base .<= 1.0) "Connectivity data not scaled between 0 - 1"
+    @assert all(0.0 .<= conn .<= 1.0) "Connectivity data not scaled between 0 - 1"
 
-    return (TP_base=TP_base, truncated=invalid_ids, site_ids=unique_site_ids)
+    return (conn=conn, truncated=invalid_ids, site_ids=loc_ids)
 end
 function site_connectivity(
     file_loc::String,
@@ -162,25 +160,31 @@ function site_connectivity(
 end
 
 """
-    connectivity_strength(TP_base::AbstractArray)::NamedTuple
-    connectivity_strength(area_weighted_TP::AbstractMatrix{Float64}, cover::Vector{Float64}, TP_cache::AbstractMatrix{Float64})::NamedTuple
+    connectivity_strength(area_weighted_conn::AbstractMatrix{Float64}, cover::Vector{Float64}, conn_cache::AbstractMatrix{Float64})::NamedTuple
 
 Create in/out degree centralities for all nodes, and vector of their strongest predecessors.
 
 # Arguments
-- `TP_base` : Base transfer probability matrix to create Directed Graph from.
-- `area_weighted_TP` : Transfer probability matrix weighted by location `k` area
+- `area_weighted_conn` : Transfer probability matrix weighted by location `k` area
 - `cover` : Total relative coral cover at location
-- `TP_cache` : Cache matrix of same size as TP_base to hold intermediate values
+- `conn_cache` : Cache matrix of same size as TP_base to hold intermediate values
+
+
+    connectivity_strength(conn::AbstractArray)::NamedTuple
+
+Create in/out degree centralities for all nodes, and vector of their strongest predecessors.
+
+# Arguments
+- `conn` : Base connectivity matrix to create Directed Graph from.
 
 # Returns
 NamedTuple:
 - `in_conn` : sites ranked by incoming connectivity
 - `out_conn` : sites ranked by outgoing connectivity
 - `strongest_predecessor` : strongest predecessor for each site
 """
-function connectivity_strength(TP_base::AbstractMatrix{Float64})::NamedTuple
-    g = SimpleDiGraph(TP_base)
+function connectivity_strength(conn::AbstractMatrix{Float64})::NamedTuple
+    g = SimpleDiGraph(conn)
 
     # Measure centrality based on number of incoming connections
     C1 = indegree_centrality(g)
@@ -209,17 +213,16 @@ function connectivity_strength(TP_base::AbstractMatrix{Float64})::NamedTuple
     return (in_conn=C1, out_conn=C2, strongest_predecessor=strong_pred)
 end
 function connectivity_strength(
-    area_weighted_TP::AbstractMatrix{Float64},
+    area_weighted_conn::AbstractMatrix{Float64},
     cover::Vector{<:Union{Float32,Float64}},
-    TP_cache::AbstractMatrix{Float64},
+    conn_cache::AbstractMatrix{Float64},
 )::NamedTuple
-
     # Accounts for cases where there is no coral cover
-    TP_cache .= (area_weighted_TP .* cover)
-    max_TP = maximum(TP_cache)
-    if max_TP > 0.0
-        TP_cache .= TP_cache ./ max_TP
+    conn_cache .= (area_weighted_conn .* cover)
+    max_conn = maximum(conn_cache)
+    if max_conn > 0.0
+        conn_cache .= conn_cache ./ max_conn
     end
 
-    return connectivity_strength(TP_cache)
+    return connectivity_strength(conn_cache)
 end

src/ecosystem/corals/growth.jl

@@ -441,6 +441,7 @@ locations.
 - `k_area` : Absolute coral habitable area
 - `tp` : Transition Probability matrix indicating the proportion of larvae that reaches a
     sink location (columns) from a given source location (rows)
+    Columns should sum to 1
 - `settlers` : recruited corals for each sink location
 - `fec_params_per_m²` : Fecundity parameters for each species/size class combination
 - `h²` : narrow-sense heritability
@@ -472,7 +473,7 @@ function settler_DHW_tolerance!(
 
         # Calculate contribution to cover to determine weights for each species/group
         w::NamedDimsArray = @views settlers[:, sink_loc]' .* tp[source_locs, sink_loc]
-        w_per_group = w ./ sum.(eachcol(w))'
+        w_per_group = w ./ sum(w, dims=1)
         replace!(w_per_group, NaN=>0.0)
 
         # Determine new distribution mean for each species at all locations
@@ -684,7 +685,7 @@ function settler_cover(
     valid_sinks::BitVector = vec(sum(conn, dims=1) .> 0.0)
 
     # Send larvae out into the world (reuse potential_settlers to reduce allocations)
-    # Note, conn rows need not sum to 1.0 as this missing probability accounts for larvae 
+    # Note, conn rows need not sum to 1.0 as this missing probability accounts for larvae
     # which do not settle. Pers comm with C. Ani (2023-01-29 13:24 AEST).
     # [Larval pool for each location in larvae/m²] * [survival rate]
     # this is known as in-water mortality.

src/io/result_io.jl

@@ -99,7 +99,7 @@ Retrieve connectivity matrices from Domain for storage.
 Produce connectivity for a particular RCP saved to a Zarr data store.
 
 # Arguments
-- `conn_data` : connectivity data (e.g. `domain.TP_data`)
+- `conn_data` : connectivity data (e.g. `domain.conn`)
 - `file_loc` : path for Zarr data store
 - `rcp`: RCP associated with connectivity data.
 
@@ -478,7 +478,7 @@ function setup_result_store!(domain::Domain, scen_spec::DataFrame)::Tuple
         push!(
             connectivity,
             store_conn(
-                domain.TP_data,
+                domain.conn,
                 joinpath(z_store.folder, "connectivity"),
                 rcp,
                 COMPRESSOR

src/scenario.jl

@@ -409,7 +409,12 @@ function run_model(domain::Domain, param_set::NamedDimsArray)::NamedTuple
     fec_params_per_m²::Vector{Float64} = corals.fecundity  # number of larvae produced per m²
 
     # Caches
-    TP_data = domain.TP_data
+    conn = domain.conn
+
+    # Determine contribution of each source to a sink location
+    # i.e., columns should sum to 1!
+    TP_data = conn ./ sum(conn, dims=1)
+
     # sf = cache.sf  # unused as it is currently deactivated
     fec_all = cache.fec_all
     fec_scope = cache.fec_scope
@@ -560,8 +565,8 @@ function run_model(domain::Domain, param_set::NamedDimsArray)::NamedTuple
     growth::ODEProblem = ODEProblem{true}(growthODE, ode_u, tspan, p)
     alg_hint = Symbol[:nonstiff]
 
-    area_weighted_TP = TP_data .* site_k_area(domain)
-    TP_cache = similar(area_weighted_TP)
+    area_weighted_conn = conn .* site_k_area(domain)
+    conn_cache = similar(area_weighted_conn)
 
     # basal_area_per_settler is the area in m^2 of a size class one coral
     basal_area_per_settler = colony_mean_area(corals.mean_colony_diameter_m[corals.class_id.==1])
@@ -590,7 +595,7 @@ function run_model(domain::Domain, param_set::NamedDimsArray)::NamedTuple
         # Recruitment/settlement occurs after the full moon in October/November
         recruitment[:, valid_locs] .= settler_cover(
             fec_scope,
-            TP_data,
+            conn,
             leftover_space_m²,
             sim_params.max_settler_density,
             sim_params.max_larval_density,
@@ -602,7 +607,7 @@ function run_model(domain::Domain, param_set::NamedDimsArray)::NamedTuple
             c_mean_t_1,
             c_mean_t,
             site_k_area(domain),
-            TP_data,
+            TP_data,  # ! IMPORTANT: Pass in transition probability matrix, not connectivity!
             recruitment,
             fec_params_per_m²,
             param_set("heritability")
@@ -649,7 +654,7 @@ function run_model(domain::Domain, param_set::NamedDimsArray)::NamedTuple
 
             # Determine connectivity strength
             # Account for cases where there is no coral cover
-            in_conn, out_conn, strong_pred = connectivity_strength(area_weighted_TP, vec(loc_coral_cover), TP_cache)
+            in_conn, out_conn, strong_pred = connectivity_strength(area_weighted_conn, vec(loc_coral_cover), conn_cache)
             (seed_locs, fog_locs, rankings) = guided_site_selection(
                 mcda_vars,
                 MCDA_approach,

test/connectivity.jl

@@ -14,6 +14,6 @@ import ADRIA.GeoDataFrames as GDF
 
     conn_details = ADRIA.site_connectivity(conn_files, unique_site_ids)
 
-    TP_data = conn_details.TP_base
-    @test all(names(TP_data, 2) .== site_data.reef_siteid) || "Sites do not match expected order!"
+    conn = conn_details.conn
+    @test all(names(conn, 2) .== site_data.reef_siteid) || "Sites do not match expected order!"
 end

test/data_loading.jl

@@ -29,9 +29,9 @@ end
 
     conn_details = ADRIA.site_connectivity(conn_files, unique_site_ids)
 
-    TP_data = conn_details.TP_base
-    @test all(axiskeys(TP_data, 1) .== axiskeys(TP_data, 2)) || "Site order does not match between rows/columns."
-    @test all(axiskeys(TP_data, 2) .== site_data.reef_siteid) || "Sites do not match expected order."
+    conn = conn_details.conn
+    @test all(axiskeys(conn, 1) .== axiskeys(conn, 2)) || "Site order does not match between rows/columns."
+    @test all(axiskeys(conn, 2) .== site_data.reef_siteid) || "Sites do not match expected order."
     @test all(unique_site_ids .== conn_details.site_ids) || "Included site ids do not match length/order in geospatial file."
 end