Skip to content

Commit 12f2bb7

Browse files
rafaqzrafaqz
committed
move sensitivity code to its own file
1 parent 3f6559a commit 12f2bb7

File tree

3 files changed

+161
-161
lines changed

3 files changed

+161
-161
lines changed

‎src/ConScape.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -62,6 +62,7 @@ include("solve.jl")
6262
include("return.jl")
6363
include("solvers.jl")
6464
include("compute.jl")
65+
include("sensitivity.jl")
6566
include("windows.jl")
6667
include("assessment.jl")
6768

‎src/compute.jl‎

Lines changed: 0 additions & 161 deletions
Original file line numberDiff line numberDiff line change
@@ -436,167 +436,6 @@ _weight(::QualityAndProximityWeighted, ti::TargetInit) = ti.M
436436
_weight(::QualityWeighted, ti::TargetInit) = ti.Q
437437
_weight(w::CustomWeighted, ti::TargetInit) = w.weight
438438

439-
######################################################################################
440-
# Sensitivity
441-
function compute(m::SensitivityAnalysis{<:SourceQuality}, ti::TargetInit{<:Union{RSP,RandomWalk}})
442-
(; qˢ, qᵗ, K, workspace) = ti
443-
if sensitivitytype(m) isa Elasticity
444-
target_sensitivity = workspace .*= qˢ .* K .* qᵗ[target(ti).node]
445-
else # sensitivitytype(m) isa Sensitivity
446-
target_sensitivity = workspace .= K .* qᵗ[target(ti).node]
447-
end
448-
return target_sensitivity
449-
end
450-
function compute(m::SensitivityAnalysis{<:TargetQuality}, ti::TargetInit{<:Union{RSP,RandomWalk}})
451-
(; qˢ, qᵗ, K, workspace) = ti
452-
if sensitivitytype(m) isa Elasticity
453-
target_sensitivity = workspace .*= qˢ .* K .* qᵗ[target(ti).node]
454-
else # sensitivitytype(m) isa Sensitivity
455-
target_sensitivity = workspace .= K .* qˢ
456-
end
457-
return target_sensitivity
458-
end
459-
compute(m::SensitivityAnalysis{<:Permeability}, ti::TargetInit{<:Union{RSP,RandomWalk}}) =
460-
compute(proximity_measure(ti), m, ti)
461-
function compute(::ExpectedCost, ::SensitivityAnalysis{<:Permeability}, ti::TargetInit)
462-
(; K, qˢ, qᵗ, Y, CW, IW_adj_factorization, Z, Zⁱ, Zrows) = ti
463-
node = target(ti).node
464-
Kd = ti.workspace .= _diff_KD(distance_transformation(ti)).(K)
465-
Md = ti.workspace .= Kd .* qˢ .* qᵗ
466-
m = sum(Md)
467-
MdZⁱ = Md .* Zⁱ
468-
MdᵀZ = ti.workspace .= MdZⁱ
469-
ldiv!(ti, IW_adj_factorization, MdᵀZ) # MdᵀZ = MdZⁱ' / IW
470-
k̂diagZⁱ = m * Zⁱ[target(ti).node]
471-
472-
X3 = ti.workspace .= k̂diagZⁱ .* Zrows
473-
X6 = ti.workspace .= MdᵀZ .- X3
474-
C̄ᵣ = ti.workspace .= Y .* Zⁱ
475-
476-
k̂diagC̄Zⁱ = k̂diagZⁱ * C̄ᵣ[node]
477-
RHS = vec(MdZⁱ .* C̄ᵣ) .- vec((MdᵀZ' * CW) .+ (X3' * CW))
478-
X5 = ldiv!(ti, IW_adj_factorization, RHS) .- k̂diagC̄Zⁱ .* Zrows
479-
480-
return (; Z, Y, X5, X6)
481-
end
482-
function compute(pmp::PowerMeanProximity, ::SensitivityAnalysis{<:Permeability}, ti::TargetInit)
483-
(; Z, θ) = ti
484-
weight = ti.workspace .= Z ./= Z[target(ti).node] .^ θ
485-
bet_node_k = compute(Betweenness(CustomWeighted(weight)), ti)
486-
bet_edge_k = compute(EdgeBetweenness(CustomWeighted(weight)), ti)
487-
return (; bet_node_k, bet_edge_k)
488-
end
489-
490-
491-
# kB and kΣ are set up as zeroed-out W matrices
492-
allocate_output(l::GridGraphLevel, m::SensitivityAnalysis{<:Permeability}, p::ConScapeProblem, args...) =
493-
allocate_output(l, ReturnDenseSpatialSum(), p, args...)
494-
allocate_output(l::Union{ConnectedGraphLevel,TargetLevel}, m::SensitivityAnalysis{<:Permeability}, p::ConScapeProblem, args...) =
495-
allocate_output(l, proximity_measure(p), m, p, args...)
496-
function allocate_output(l::Level, ::PowerMeanProximity, m::SensitivityAnalysis{<:Permeability}, p::ConScapeProblem, gg::GridGraph, cg::ConnectedGraph, precalculation)
497-
(; W) = precalculation
498-
bet_edge_k = allocate_output(l, EdgeBetweenness(CustomWeighted(nothing)), p, gg, cg, precalculation)
499-
bet_node_k = allocate_output(l, Betweenness(CustomWeighted(nothing)), p, gg, cg, precalculation)
500-
result = mapnz(_ -> 0.0, W)
501-
resultrows = zeros(size(W, 1))
502-
return l => (; bet_node_k, bet_edge_k, result, resultrows)
503-
end
504-
# kB and kΣ are set up as zeroed-out W matrices
505-
function allocate_output(l::Level, ::ExpectedCost, m::SensitivityAnalysis{<:Permeability}, ::ConScapeProblem, ::GridGraph, ::ConnectedGraph, precalculation)
506-
(; W) = precalculation
507-
kB = mapnz(_ -> 0.0, W)
508-
kΣ = mapnz(_ -> 0.0, W)
509-
result = mapnz(_ -> 0.0, W)
510-
resultrows = zeros(size(W, 1))
511-
return l => (; kB, kΣ, result, resultrows)
512-
end
513-
514-
# And we store into them for each target
515-
update_output!(output::NamedTuple, l::ConnectedGraphLevel, m::SensitivityAnalysis{<:Permeability}, ti::TargetInit, v) =
516-
update_output!(output, l, proximity_measure(ti), m, ti, v)
517-
function update_output!(output, l::ConnectedGraphLevel, ::ExpectedCost, m::SensitivityAnalysis{<:Permeability}, ti::TargetInit, v)
518-
(; W, C) = ti
519-
(; kB, kΣ) = output
520-
(; Z, Y, X5, X6) = v
521-
522-
foreachnz(W) do i, j, n
523-
kB.nzval[n] += W.nzval[n] * Z[j] * X6[i]
524-
kΣ.nzval[n] += Z[j] * X5[i] - Y[j] * X6[i] - C.nzval[n] * kB.nzval[n] / W.nzval[n]
525-
end
526-
return output
527-
end
528-
function update_output!(output, l::ConnectedGraphLevel, ::PowerMeanProximity, m::SensitivityAnalysis{<:Permeability}, ti::TargetInit, v)
529-
(; bet_edge_k, bet_node_k) = v
530-
# Just call update_output! on the component parts
531-
update_output!(output.bet_edge_k, EdgeBetweenness(CustomWeighted(nothing)), ti, bet_edge_k)
532-
update_output!(output.bet_node_k, Betweenness(CustomWeighted(nothing)), ti, bet_node_k)
533-
return nothing
534-
end
535-
536-
# And after all targets run, finalize them
537-
finalize_output!(output::Pair, m::SensitivityAnalysis{<:Permeability}, sgi::ConnectedGraphInit) =
538-
finalize_output!(output[2], proximity_measure(sgi), m, sgi)
539-
function finalize_output!(output::NamedTuple, ::ExpectedCost, m::SensitivityAnalysis{<:Permeability}, sgi::ConnectedGraphInit)
540-
(; kB, kΣ, result) = output
541-
(; A_rowsums, Aⁱ) = precalculation(sgi)
542-
543-
kΣ_node = sum(kΣ, dims=2)
544-
foreachnz(Aⁱ) do i, j, n
545-
S_cost = kB.nzval[n] + theta(sgi) * kΣ.nzval[n]
546-
S_likelihood = (kΣ_node[j] / A_rowsums[j]) * 1 - kΣ.nzval[n] * Aⁱ.nzval[n]
547-
S_e_likelihood = _maybe_scale(S_likelihood, sensitivitytype(m), wrt(m), sgi, n)
548-
S_e_cost_scaled = _maybe_scale(S_cost, sensitivitytype(m), wrt(m), sgi, n)
549-
result.nzval[n] = _combine_sensitivity(wrt(m), S_e_likelihood, S_e_cost_scaled, sgi, n)
550-
end
551-
552-
return output
553-
end
554-
function finalize_output!(output::NamedTuple, ::PowerMeanProximity, m::SensitivityAnalysis{<:Permeability}, sgi::ConnectedGraphInit)
555-
(; bet_edge_k, bet_node_k, result, resultrows) = output
556-
(; A_rowsums, Aⁱ) = precalculation(sgi)
557-
A = steplikelihood(sgi)
558-
foreachnz(Aⁱ) do i, j, n
559-
S_cost = bet_edge_k[2].nzval[n]
560-
S_likelihood = bet_edge_k[2].nzval[n] * Aⁱ.nzval[n] - bet_node_k[2][j] / A_rowsums[j] * A.nzval[n] * theta(sgi)
561-
S_e_likelihood = _maybe_scale(S_likelihood, sensitivitytype(m), wrt(m), sgi, n)
562-
S_e_cost_scaled = _maybe_scale(S_cost, sensitivitytype(m), wrt(m), sgi, n)
563-
result.nzval[n] = _combine_sensitivity(wrt(m), S_e_likelihood, S_e_cost_scaled, sgi, n)
564-
end
565-
# TODO dont allocate
566-
return resultrows .= vec(sum(result; dims=1))
567-
end
568-
569-
transfer_output!(dest::AbstractMatrix, source::NamedTuple, ::SensitivityAnalysis, cgi::ConnectedGraphInit) =
570-
dest[sourceids(cgi)] .= source.resultrows
571-
572-
_combine_sensitivity(::StepLikelihood, S_e_likelihood, S_e_cost, ti, n) = S_e_likelihood
573-
_combine_sensitivity(::StepCost, S_e_likelihood, S_e_cost, ti, n) = S_e_cost
574-
function _combine_sensitivity(::StepCostToLikelihood, S_e_likelihood, S_e_cost, ti, n)
575-
f = _diff_CA(costfunction(ti))
576-
C = stepcost(ti)
577-
S_e_likelihood + S_e_cost * f(C.nzval[n])
578-
end
579-
function _combine_sensitivity(::StepLikelihoodToCost, S_e_likelihood, S_e_cost, ti, n)
580-
f = _diff_AC(costfunction(ti))
581-
L = steplikelihood(ti)
582-
S_e_cost + S_e_likelihood * f(L.nzval[n])
583-
end
584-
585-
_maybe_scale(a, ::Elasticity, ::Union{StepLikelihood,StepLikelihoodToCost}, ti, n) =
586-
a * steplikelihood(ti).nzval[n]
587-
_maybe_scale(a, ::Elasticity, ::Union{StepCost,StepCostToLikelihood}, ti, n) =
588-
a * stepcost(ti).nzval[n]
589-
_maybe_scale(a, ::Sensitivity, ::Permeability, ti, n) = a
590-
591-
_diff_CA(::MinusLog) = x -> -inv(x)
592-
_diff_AC(::MinusLog) = x -> -(x)
593-
_diff_CA(::Inv) = x -> -inv(x^2)
594-
_diff_AC(::Inv) = x -> -inv(x^2)
595-
596-
_diff_KD(x::ExpMinusAlpha) = k -> -k * x.alpha
597-
_diff_KD(::ExpMinus) = k -> -k
598-
_diff_KD(::Inv) = k -> -k ^ 2
599-
600439
# TODO: handle self connectivity for single isolated nodes
601440
# fill_isolated_node(::FunctionalHabitat, init::Initalisation, target::CartesianIndex) =
602441
# diagvalue(init) * sourcequality_spatial(init)[target] * targetquality_spatial(init)[target]

‎src/sensitivity.jl‎

Lines changed: 160 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,160 @@
1+
######################################################################################
2+
# Sensitivity
3+
function compute(m::SensitivityAnalysis{<:SourceQuality}, ti::TargetInit{<:Union{RSP,RandomWalk}})
4+
(; qˢ, qᵗ, K, workspace) = ti
5+
if sensitivitytype(m) isa Elasticity
6+
target_sensitivity = workspace .*= qˢ .* K .* qᵗ[target(ti).node]
7+
else # sensitivitytype(m) isa Sensitivity
8+
target_sensitivity = workspace .= K .* qᵗ[target(ti).node]
9+
end
10+
return target_sensitivity
11+
end
12+
function compute(m::SensitivityAnalysis{<:TargetQuality}, ti::TargetInit{<:Union{RSP,RandomWalk}})
13+
(; qˢ, qᵗ, K, workspace) = ti
14+
if sensitivitytype(m) isa Elasticity
15+
target_sensitivity = workspace .*= qˢ .* K .* qᵗ[target(ti).node]
16+
else # sensitivitytype(m) isa Sensitivity
17+
target_sensitivity = workspace .= K .* qˢ
18+
end
19+
return target_sensitivity
20+
end
21+
compute(m::SensitivityAnalysis{<:Permeability}, ti::TargetInit{<:Union{RSP,RandomWalk}}) =
22+
compute(proximity_measure(ti), m, ti)
23+
function compute(::ExpectedCost, ::SensitivityAnalysis{<:Permeability}, ti::TargetInit)
24+
(; K, qˢ, qᵗ, Y, CW, IW_adj_factorization, Z, Zⁱ, Zrows) = ti
25+
node = target(ti).node
26+
Kd = ti.workspace .= _diff_KD(distance_transformation(ti)).(K)
27+
Md = ti.workspace .= Kd .* qˢ .* qᵗ
28+
m = sum(Md)
29+
MdZⁱ = Md .* Zⁱ
30+
MdᵀZ = ti.workspace .= MdZⁱ
31+
ldiv!(ti, IW_adj_factorization, MdᵀZ) # MdᵀZ = MdZⁱ' / IW
32+
k̂diagZⁱ = m * Zⁱ[target(ti).node]
33+
34+
X3 = ti.workspace .= k̂diagZⁱ .* Zrows
35+
X6 = ti.workspace .= MdᵀZ .- X3
36+
C̄ᵣ = ti.workspace .= Y .* Zⁱ
37+
38+
k̂diagC̄Zⁱ = k̂diagZⁱ * C̄ᵣ[node]
39+
RHS = vec(MdZⁱ .* C̄ᵣ) .- vec((MdᵀZ' * CW) .+ (X3' * CW))
40+
X5 = ldiv!(ti, IW_adj_factorization, RHS) .- k̂diagC̄Zⁱ .* Zrows
41+
42+
return (; Z, Y, X5, X6)
43+
end
44+
function compute(pmp::PowerMeanProximity, ::SensitivityAnalysis{<:Permeability}, ti::TargetInit)
45+
(; Z, θ) = ti
46+
weight = ti.workspace .= Z ./= Z[target(ti).node] .^ θ
47+
bet_node_k = compute(Betweenness(CustomWeighted(weight)), ti)
48+
bet_edge_k = compute(EdgeBetweenness(CustomWeighted(weight)), ti)
49+
return (; bet_node_k, bet_edge_k)
50+
end
51+
52+
53+
# kB and kΣ are set up as zeroed-out W matrices
54+
allocate_output(l::GridGraphLevel, m::SensitivityAnalysis{<:Permeability}, p::ConScapeProblem, args...) =
55+
allocate_output(l, ReturnDenseSpatialSum(), p, args...)
56+
allocate_output(l::Union{ConnectedGraphLevel,TargetLevel}, m::SensitivityAnalysis{<:Permeability}, p::ConScapeProblem, args...) =
57+
allocate_output(l, proximity_measure(p), m, p, args...)
58+
function allocate_output(l::Level, ::PowerMeanProximity, m::SensitivityAnalysis{<:Permeability}, p::ConScapeProblem, gg::GridGraph, cg::ConnectedGraph, precalculation)
59+
(; W) = precalculation
60+
bet_edge_k = allocate_output(l, EdgeBetweenness(CustomWeighted(nothing)), p, gg, cg, precalculation)
61+
bet_node_k = allocate_output(l, Betweenness(CustomWeighted(nothing)), p, gg, cg, precalculation)
62+
result = mapnz(_ -> 0.0, W)
63+
resultrows = zeros(size(W, 1))
64+
return l => (; bet_node_k, bet_edge_k, result, resultrows)
65+
end
66+
# kB and kΣ are set up as zeroed-out W matrices
67+
function allocate_output(l::Level, ::ExpectedCost, m::SensitivityAnalysis{<:Permeability}, ::ConScapeProblem, ::GridGraph, ::ConnectedGraph, precalculation)
68+
(; W) = precalculation
69+
kB = mapnz(_ -> 0.0, W)
70+
kΣ = mapnz(_ -> 0.0, W)
71+
result = mapnz(_ -> 0.0, W)
72+
resultrows = zeros(size(W, 1))
73+
return l => (; kB, kΣ, result, resultrows)
74+
end
75+
76+
# And we store into them for each target
77+
update_output!(output::NamedTuple, l::ConnectedGraphLevel, m::SensitivityAnalysis{<:Permeability}, ti::TargetInit, v) =
78+
update_output!(output, l, proximity_measure(ti), m, ti, v)
79+
function update_output!(output, l::ConnectedGraphLevel, ::ExpectedCost, m::SensitivityAnalysis{<:Permeability}, ti::TargetInit, v)
80+
(; W, C) = ti
81+
(; kB, kΣ) = output
82+
(; Z, Y, X5, X6) = v
83+
84+
foreachnz(W) do i, j, n
85+
kB.nzval[n] += W.nzval[n] * Z[j] * X6[i]
86+
kΣ.nzval[n] += Z[j] * X5[i] - Y[j] * X6[i] - C.nzval[n] * kB.nzval[n] / W.nzval[n]
87+
end
88+
return output
89+
end
90+
function update_output!(output, l::ConnectedGraphLevel, ::PowerMeanProximity, m::SensitivityAnalysis{<:Permeability}, ti::TargetInit, v)
91+
(; bet_edge_k, bet_node_k) = v
92+
# Just call update_output! on the component parts
93+
update_output!(output.bet_edge_k, EdgeBetweenness(CustomWeighted(nothing)), ti, bet_edge_k)
94+
update_output!(output.bet_node_k, Betweenness(CustomWeighted(nothing)), ti, bet_node_k)
95+
return nothing
96+
end
97+
98+
# And after all targets run, finalize them
99+
finalize_output!(output::Pair, m::SensitivityAnalysis{<:Permeability}, sgi::ConnectedGraphInit) =
100+
finalize_output!(output[2], proximity_measure(sgi), m, sgi)
101+
function finalize_output!(output::NamedTuple, ::ExpectedCost, m::SensitivityAnalysis{<:Permeability}, sgi::ConnectedGraphInit)
102+
(; kB, kΣ, result) = output
103+
(; A_rowsums, Aⁱ) = precalculation(sgi)
104+
105+
kΣ_node = sum(kΣ, dims=2)
106+
foreachnz(Aⁱ) do i, j, n
107+
S_cost = kB.nzval[n] + theta(sgi) * kΣ.nzval[n]
108+
S_likelihood = (kΣ_node[j] / A_rowsums[j]) * 1 - kΣ.nzval[n] * Aⁱ.nzval[n]
109+
S_e_likelihood = _maybe_scale(S_likelihood, sensitivitytype(m), wrt(m), sgi, n)
110+
S_e_cost_scaled = _maybe_scale(S_cost, sensitivitytype(m), wrt(m), sgi, n)
111+
result.nzval[n] = _combine_sensitivity(wrt(m), S_e_likelihood, S_e_cost_scaled, sgi, n)
112+
end
113+
114+
return output
115+
end
116+
function finalize_output!(output::NamedTuple, ::PowerMeanProximity, m::SensitivityAnalysis{<:Permeability}, sgi::ConnectedGraphInit)
117+
(; bet_edge_k, bet_node_k, result, resultrows) = output
118+
(; A_rowsums, Aⁱ) = precalculation(sgi)
119+
A = steplikelihood(sgi)
120+
foreachnz(Aⁱ) do i, j, n
121+
S_cost = bet_edge_k[2].nzval[n]
122+
S_likelihood = bet_edge_k[2].nzval[n] * Aⁱ.nzval[n] - bet_node_k[2][j] / A_rowsums[j] * A.nzval[n] * theta(sgi)
123+
S_e_likelihood = _maybe_scale(S_likelihood, sensitivitytype(m), wrt(m), sgi, n)
124+
S_e_cost_scaled = _maybe_scale(S_cost, sensitivitytype(m), wrt(m), sgi, n)
125+
result.nzval[n] = _combine_sensitivity(wrt(m), S_e_likelihood, S_e_cost_scaled, sgi, n)
126+
end
127+
# TODO dont allocate
128+
return resultrows .= vec(sum(result; dims=1))
129+
end
130+
131+
transfer_output!(dest::AbstractMatrix, source::NamedTuple, ::SensitivityAnalysis, cgi::ConnectedGraphInit) =
132+
dest[sourceids(cgi)] .= source.resultrows
133+
134+
_combine_sensitivity(::StepLikelihood, S_e_likelihood, S_e_cost, ti, n) = S_e_likelihood
135+
_combine_sensitivity(::StepCost, S_e_likelihood, S_e_cost, ti, n) = S_e_cost
136+
function _combine_sensitivity(::StepCostToLikelihood, S_e_likelihood, S_e_cost, ti, n)
137+
f = _diff_CA(costfunction(ti))
138+
C = stepcost(ti)
139+
S_e_likelihood + S_e_cost * f(C.nzval[n])
140+
end
141+
function _combine_sensitivity(::StepLikelihoodToCost, S_e_likelihood, S_e_cost, ti, n)
142+
f = _diff_AC(costfunction(ti))
143+
L = steplikelihood(ti)
144+
S_e_cost + S_e_likelihood * f(L.nzval[n])
145+
end
146+
147+
_maybe_scale(a, ::Elasticity, ::Union{StepLikelihood,StepLikelihoodToCost}, ti, n) =
148+
a * steplikelihood(ti).nzval[n]
149+
_maybe_scale(a, ::Elasticity, ::Union{StepCost,StepCostToLikelihood}, ti, n) =
150+
a * stepcost(ti).nzval[n]
151+
_maybe_scale(a, ::Sensitivity, ::Permeability, ti, n) = a
152+
153+
_diff_CA(::MinusLog) = x -> -inv(x)
154+
_diff_AC(::MinusLog) = x -> -(x)
155+
_diff_CA(::Inv) = x -> -inv(x^2)
156+
_diff_AC(::Inv) = x -> -inv(x^2)
157+
158+
_diff_KD(x::ExpMinusAlpha) = k -> -k * x.alpha
159+
_diff_KD(::ExpMinus) = k -> -k
160+
_diff_KD(::Inv) = k -> -k ^ 2

0 commit comments

Comments
 (0)