Skip to content

Add smooth_clamp function #2586

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Add smooth_clamp function #2586

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions core/src/solve.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -436,7 +436,7 @@ function linear_resistance_flow(

Δh = h_a - h_b
q_unlimited = Δh / resistance[node_id.idx]
q = clamp(q_unlimited, -max_flow_rate[node_id.idx], max_flow_rate[node_id.idx])
q = smooth_clamp(q_unlimited, -max_flow_rate[node_id.idx], max_flow_rate[node_id.idx])
return q * low_storage_factor_resistance_node(p, q_unlimited, inflow_id, outflow_id)
end

Expand Down Expand Up @@ -698,7 +698,7 @@ function formulate_pump_or_outlet_flow!(
end
end
end
q = clamp(q, lower_bound, upper_bound)
q = smooth_clamp(q, lower_bound, upper_bound)

# Special case for outlet: check level difference
if reduce_Δlevel
Expand Down
55 changes: 55 additions & 0 deletions core/src/util.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -314,6 +314,60 @@ function reduction_factor(x::T, threshold::Real)::T where {T<:Real}
end
end

"""
Function that goes smoothly from lo to hi in the interval [lo - δ, hi + δ],
where δ = α * (hi - lo) / 2 with smoothing parameter α ∈ [0, 1]. This function is identical
to the non-smooth clamp function outside the intervals [lo - δ, lo + δ] and [hi - δ, hi + δ].
"""
function smooth_clamp(x, lo, hi; α = 0.2)
# Zero length interval
Copy link
Contributor

@Huite Huite Oct 1, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we want to check for lo > hi, or 0 <= α <= 1?

...I just tested, the clamp function doesn't care about lo and hi, so it's probably okay to leave unchecked, see also its source:

function clamp(x::X, lo::L, hi::H) where {X,L,H}
    T = promote_type(X, L, H)
    return (x > hi) ? convert(T, hi) : (x < lo) ? convert(T, lo) : convert(T, x)
end

The behavior isn't the same though if you fill in values "wrongly":

julia> smooth_clamp(0.01, 1.0, 0.0,α=0.01)
1.0
julia> clamp(0.01, 1.0, 0.0)
0.0

Might be best in some sense to have it mimick the lo vs hi behavior.

if lo == hi
return lo
end

lo_inf = isinf(lo)
hi_inf = isinf(hi)

# No clamping
if lo_inf && hi_inf
return x
end

δ = α * (hi - lo) / 2
Copy link
Contributor

@Huite Huite Oct 1, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If either hi or lo are inf, we just get the hard clamp. Ideally though, I think you'd want the clamping one on side not to be influenced on what happens on the other side, i.e.:

smooth_clamp(x_near_lo, lo, inf) == smooth_clamp(x_near_lo, lo, hi)

Thinking about this, why α instead of δ? δ works regardless of finite lo or hi. Arguably, δ is more intuitive as well, since your docstring explicitly starts off with δ too:

Function that goes smoothly from lo to hi in the interval [lo - δ, hi + δ]

And then defines it in terms of α. So in use, I have to do the α arithmetic first.

I think in terms of general use, we'd mostly have an absolute range in mind? If a relative zone is useful, it can be either provided in the call; alternatively you could support both arguments in an abstol, reltol kind of way.

Having α makes it easier to choose a default value, but... we probably don't want a default value here! The width of the clamping can make a big difference. I.e. in this line:

q = smooth_clamp(q_unlimited, -max_flow_rate[node_id.idx], max_flow_rate[node_id.idx])

The smoothing width should probably be explicit.


# Happens for one-sided clamping
if isinf(δ)
return x
end

lower_lower_bound = lo - δ

if x < lower_lower_bound
return lo
end

lower_upper_bound = lo + δ
denom = 4δ

if x < lower_upper_bound
return lo + ((x - lower_lower_bound)^2) / denom
end

upper_lower_bound = hi - δ

if x < upper_lower_bound
return x
end

upper_upper_bound = hi + δ

return if x < upper_upper_bound
hi - ((x - upper_upper_bound)^2) / denom
else
hi
end
end

function get_low_storage_factor(p::Parameters, id::NodeID)
(; current_low_storage_factor) = p.state_time_dependent_cache
if id.type == NodeType.Basin
Expand Down Expand Up @@ -684,6 +738,7 @@ end

# Overloads for SparseConnectivityTracer
reduction_factor(x::GradientTracer, ::Real) = x
smooth_clamp(x::GradientTracer, ::Real, ::Real) = x
low_storage_factor_resistance_node(::Parameters, q::GradientTracer, ::NodeID, ::NodeID) = q
relaxed_root(x::GradientTracer, threshold::Real) = x
get_level_from_storage(basin::Basin, state_idx::Int, storage::GradientTracer) = storage
Expand Down
8 changes: 8 additions & 0 deletions core/test/utils_test.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -336,6 +336,14 @@ end
@test reduction_factor(-Inf, 2.0) === 0.0
end

@testitem "smooth_clamp" begin
using Ribasim: smooth_clamp
@test smooth_clamp(2.5, 2.0, 3.0) == 2.5
@test smooth_clamp(3.14, -Inf, Inf) == 3.14
@test smooth_clamp(1.0, 1.41, 1.41) == 1.41
@test smooth_clamp(1.95, 1.0, 2.0) ≈ 1.94375
end

@testitem "Node types" begin
using Ribasim:
node_types,
Expand Down
Loading