From 61d1f40a1e4fd8b40f82aaafa909527b494515ba Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Thu, 1 Feb 2024 11:35:56 +1100 Subject: [PATCH 1/5] Handle change from `TP_data` -> `conn` TP_data (transition probability) implies the rows should sum to one (source to sink). Changing to `conn` (for connectivity) which indicates the proportion of larvae that survive and arrive/settle on a sink location. The remaining proportion indicates larvae "lost" out to sea, hence there is no expectation for rows to sum to 1.0 --- src/ExtInterface/ADRIA/Domain.jl | 4 ++-- src/ExtInterface/ReefMod/Domain.jl | 2 +- src/decision/dMCDA.jl | 2 +- src/decision/location_selection.jl | 4 ++-- src/ecosystem/connectivity.jl | 27 ++++++++++++--------------- src/io/result_io.jl | 4 ++-- test/connectivity.jl | 4 ++-- test/data_loading.jl | 6 +++--- 8 files changed, 25 insertions(+), 28 deletions(-) diff --git a/src/ExtInterface/ADRIA/Domain.jl b/src/ExtInterface/ADRIA/Domain.jl index 6aa885e96..2e4de5a99 100644 --- a/src/ExtInterface/ADRIA/Domain.jl +++ b/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::Σ # site 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 diff --git a/src/ExtInterface/ReefMod/Domain.jl b/src/ExtInterface/ReefMod/Domain.jl index 403aaf0ee..beb902213 100644 --- a/src/ExtInterface/ReefMod/Domain.jl +++ b/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 diff --git a/src/decision/dMCDA.jl b/src/decision/dMCDA.jl index 1105ed721..0ec0fbc9a 100644 --- a/src/decision/dMCDA.jl +++ b/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, diff --git a/src/decision/location_selection.jl b/src/decision/location_selection.jl index 83835a879..61c94e9c1 100644 --- a/src/decision/location_selection.jl +++ b/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. diff --git a/src/ecosystem/connectivity.jl b/src/ecosystem/connectivity.jl index b2bf3a92c..9e738e81a 100644 --- a/src/ecosystem/connectivity.jl +++ b/src/ecosystem/connectivity.jl @@ -26,7 +26,7 @@ NOTE: Transposes transitional probability matrix if `swap == true` # Returns NamedTuple: -- `TP_data` : Matrix, containing the transition probability for all sites +- `conn` : Matrix, containing the connectivity for all sites - `truncated` : ID of sites removed - `site_ids` : ID of sites kept """ @@ -41,16 +41,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 +129,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, diff --git a/src/io/result_io.jl b/src/io/result_io.jl index 6dc05b62e..f34164e6c 100644 --- a/src/io/result_io.jl +++ b/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 diff --git a/test/connectivity.jl b/test/connectivity.jl index b8248dbe7..d4864a6b1 100644 --- a/test/connectivity.jl +++ b/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 diff --git a/test/data_loading.jl b/test/data_loading.jl index d14c807d6..61d2400e3 100644 --- a/test/data_loading.jl +++ b/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 From fcdcf4f52db10773714e969bd25995ef9351455c Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Thu, 1 Feb 2024 11:50:39 +1100 Subject: [PATCH 2/5] Clean up variable names and usage for clarity --- src/ecosystem/connectivity.jl | 35 ++++++++++++++++++++--------------- src/scenario.jl | 17 +++++++++++------ 2 files changed, 31 insertions(+), 21 deletions(-) diff --git a/src/ecosystem/connectivity.jl b/src/ecosystem/connectivity.jl index 9e738e81a..f793540dc 100644 --- a/src/ecosystem/connectivity.jl +++ b/src/ecosystem/connectivity.jl @@ -159,16 +159,22 @@ 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: @@ -176,8 +182,8 @@ NamedTuple: - `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) @@ -206,17 +212,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 diff --git a/src/scenario.jl b/src/scenario.jl index f3e7e0e95..aab27bea6 100644 --- a/src/scenario.jl +++ b/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, From f795e83679a0e714cb5f6fe34b07ac8cd309e748 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Thu, 1 Feb 2024 11:52:52 +1100 Subject: [PATCH 3/5] Add notes to docstring Small simplification to column summing. --- src/ecosystem/corals/growth.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/ecosystem/corals/growth.jl b/src/ecosystem/corals/growth.jl index 7ae7871a2..222cd6b3c 100644 --- a/src/ecosystem/corals/growth.jl +++ b/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. From b03c234786a3e7940c333f0a59b86eb5bc09ba57 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Thu, 1 Feb 2024 13:52:51 +1100 Subject: [PATCH 4/5] Add note regarding connectivity --- src/ecosystem/connectivity.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/ecosystem/connectivity.jl b/src/ecosystem/connectivity.jl index f793540dc..ba8ab7eb0 100644 --- a/src/ecosystem/connectivity.jl +++ b/src/ecosystem/connectivity.jl @@ -26,9 +26,10 @@ NOTE: Transposes transitional probability matrix if `swap == true` # Returns NamedTuple: -- `conn` : Matrix, containing the connectivity 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, From 0ecc6b1a3f5e1ee94a90d5765e74adec8844d73d Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Thu, 1 Feb 2024 13:55:30 +1100 Subject: [PATCH 5/5] Remove reference to site where changes have been made --- src/ExtInterface/ADRIA/Domain.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ExtInterface/ADRIA/Domain.jl b/src/ExtInterface/ADRIA/Domain.jl index 2e4de5a99..3a574218f 100644 --- a/src/ExtInterface/ADRIA/Domain.jl +++ b/src/ExtInterface/ADRIA/Domain.jl @@ -14,7 +14,7 @@ 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 conn::Σ # 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