Non-linear rate functions for ContinuousObservable

[EwoutH]

Building on Discussion #2529 and PR #2851 which adds ContinuousObservable to mesa_signals, we need to address how to handle non-linear rate functions with threshold detection.

The Problem

The current implementation assumes linear rates, which allows exact threshold crossing calculation: if value_new = value_old + rate * elapsed_time, we can solve algebraically for when value_old + rate * t = threshold to get t = (threshold - value_old) / rate.

However, many real-world phenomena follow non-linear dynamics. Energy doesn't always deplete at a constant rate - it might decay exponentially (proportional to current energy level), or follow more complex patterns like logistic growth for populations or oscillatory patterns for circadian rhythms. With arbitrary rate functions rate_func(value, elapsed, agent), there's generally no closed-form solution for threshold crossings.

This creates a fundamental tension: we want to support realistic non-linear dynamics while maintaining the efficiency and precision of lazy evaluation with exact threshold detection.

Common Non-linear Patterns

After linear rates, the most common patterns in agent-based modeling are:

Exponential decay/growth appears in radioactive decay, population growth, chemical reactions, drug metabolism, and temperature equilibration. The differential equation dv/dt = k*v has the solution v(t) = v₀ * e^(k*t), and threshold crossings can be solved exactly:

t_cross = ln(threshold / v₀) / k

This pattern is probably the single most important non-linear case to support.

Logistic growth models populations with carrying capacity, disease spread, market saturation, and learning curves. Given dv/dt = r*v*(1 - v/K), the solution is v(t) = K / (1 + A*e^(-r*t)) where A = (K - v₀)/v₀. Threshold crossings require solving:

t_cross = -ln((K/threshold - 1) / A) / r

Still analytical, but more complex than exponential.

Power laws appear in allometric scaling, diffusion, and network growth. For dv/dt = k*v^n, the solution depends on the exponent n. Most values have analytical solutions, though they require careful case handling:

# For n ≠ 1:
t_cross = (threshold^(1-n) - v₀^(1-n)) / (k*(1-n))

Harmonic/oscillatory patterns like dv/dt = A*sin(ωt + φ) + offset model circadian rhythms, seasonal effects, and economic cycles. These generally require numerical root-finding to locate threshold crossings.

Together, linear, exponential, and logistic patterns likely cover around 90% of realistic ABM use cases where exact threshold detection matters.

Possible Approaches

1. Hand-coded Analytical Solutions

We could implement analytical solutions for the most common patterns as distinct classes. A RateFunction base class would define methods for both calculating rates and finding threshold crossings:

class RateFunction:
    def __call__(self, value, elapsed, agent):
        """Calculate the rate at the given value."""
        raise NotImplementedError
    
    def integrate(self, value_old, elapsed, agent):
        """Calculate new value after elapsed time."""
        raise NotImplementedError
    
    def find_threshold_crossing(self, value_old, threshold, agent):
        """Return time when threshold is crossed, or None."""
        raise NotImplementedError

class ExponentialDecay(RateFunction):
    def __init__(self, k):
        self.k = k
    
    def __call__(self, value, elapsed, agent):
        return -self.k * value
    
    def integrate(self, value_old, elapsed, agent):
        return value_old * np.exp(-self.k * elapsed)
    
    def find_threshold_crossing(self, value_old, threshold, agent):
        if value_old <= 0 or threshold <= 0 or self.k <= 0:
            return None
        if threshold > value_old:
            return None  # Decaying, can't reach higher threshold
        return np.log(threshold / value_old) / (-self.k)

Users would then write:

energy = ContinuousObservable(
    initial_value=100.0,
    rate=ExponentialDecay(k=0.05)
)

Advantages: No additional dependencies, exact threshold detection for common cases, clear semantics, ~200-300 lines of documented mathematical code.

Disadvantages: Only covers predefined patterns, users need to learn the rate function classes, doesn't handle truly arbitrary functions.

2. Numerical Integration with Adaptive Sampling

When checking for threshold crossings, simulate the continuous evolution by sampling at intervals:

def _check_thresholds_numerical(self, old_value, new_value, elapsed_time, agent):
    n_samples = min(max(int(elapsed_time * 10), 10), 1000)
    dt = elapsed_time / n_samples
    
    current_value = old_value
    for i in range(n_samples):
        rate = self.rate_func(current_value, i * dt, agent)
        next_value = current_value + rate * dt
        
        for threshold in self._thresholds:
            if self._crosses_threshold(current_value, next_value, threshold):
                # Linear interpolation for approximate crossing time
                t_cross = i * dt + self._interpolate_crossing_time(...)
                crossed.append((threshold, direction, t_cross))
        
        current_value = next_value
    
    return crossed

Advantages: Works with any rate function, configurable precision, can use more sophisticated integration methods (RK4, adaptive step size).

Disadvantages: Approximation error (might miss rapid crossings between samples), performance overhead with many agents, complexity in choosing sampling resolution, not truly "exact" timing.

3. Hybrid Approach with Tiered Complexity

Combine analytical and numerical methods based on what the user provides:

# Tier 1: Simple constant rate (current implementation)
energy = ContinuousObservable(initial_value=100.0, rate=-0.5)

# Tier 2: Common patterns with analytical solutions
energy = ContinuousObservable(
    initial_value=100.0,
    rate=ExponentialDecay(k=0.05)
)

population = ContinuousObservable(
    initial_value=10.0,
    rate=LogisticGrowth(r=0.3, K=100.0)
)

# Tier 3: Arbitrary functions with numerical integration
energy = ContinuousObservable(
    initial_value=100.0,
    rate_func=lambda v, e, a: -a.base_rate * v * (1 + 0.1 * np.sin(a.model.time)),
    integration_samples=100
)

The implementation would auto-detect which tier to use:

class ContinuousObservable(Observable):
    def __init__(self, initial_value=0.0, rate=None, rate_func=None, integration_samples=50):
        if isinstance(rate, (int, float)):
            self._method = 'linear'
        elif isinstance(rate, RateFunction):
            self._method = 'analytical'
        elif rate_func is not None:
            self._method = 'numerical'

Advantages: Best of both worlds, optimizes common cases, provides escape hatch for complex scenarios, clear performance/accuracy trade-offs.

Disadvantages: More complex implementation, two code paths to maintain, users need to understand when to use which tier.

4. SciPy Integration for Numerical Cases

For arbitrary rate functions, leverage SciPy's excellent ODE solvers:

from scipy.integrate import solve_ivp
from scipy.optimize import brentq

def integrate_numerical(self, value_old, elapsed, agent):
    result = solve_ivp(
        lambda t, v: self.rate_func(v[0], t, agent),
        (0, elapsed),
        [value_old],
        method='RK45',
        dense_output=True
    )
    return result.y[0][-1]

def find_crossing_numerical(self, value_old, threshold, agent, max_time):
    result = solve_ivp(..., dense_output=True)
    try:
        return brentq(lambda t: result.sol(t)[0] - threshold, 0, max_time)
    except ValueError:
        return None

Advantages: SciPy is already a transitive dependency for some Mesa features, excellent numerical solvers, adaptive step size, handles stiff equations.

Disadvantages: Heavier dependency if made required, still numerical approximation, need to solve forward in time then search for crossing, overhead for simple cases.

5. Lazy with Safety Warnings

Accept the limitations of lazy evaluation but add safeguards:

def __get__(self, instance, owner):
    elapsed_time = self._get_elapsed_time(instance)
    
    if elapsed_time > self.max_time_between_checks:
        warnings.warn(
            f"ContinuousObservable '{self.name}' hasn't been accessed "
            f"for {elapsed_time} time units. Thresholds may have been missed."
        )
    
    # Use numerical integration if elapsed time is large
    if elapsed_time > NUMERICAL_INTEGRATION_THRESHOLD:
        return self._calculate_with_numerical_integration(...)
    else:
        return self._calculate_with_linear_approximation(...)

Advantages: Simple API, automatically adapts to usage patterns, warns about potential issues.

Disadvantages: Still not perfect, warnings might be annoying, feels like a patch rather than a solution.

Existing Libraries

I surveyed existing Python libraries to see if any could help:

SymPy can solve ODEs symbolically and find roots, but it's a heavy dependency (~100MB installed) with significant performance overhead. The symbolic manipulation is powerful but overkill for Mesa's runtime needs. It might be useful as an optional advanced feature, but not as a required dependency.

SciPy provides excellent numerical ODE solvers (solve_ivp) and root finding (brentq), but no symbolic solutions. Mesa already has SciPy as a transitive dependency for some features. This would be useful for the numerical fallback case but doesn't help with analytical solutions.

NumPy has polynomial root finding (np.roots) which could help for polynomial rate functions, but doesn't solve the general ODE problem. Since Mesa already depends on NumPy, this could be useful for specific cases.

There doesn't appear to be a perfect existing library that provides lightweight analytical solutions for common ODE patterns. Most libraries are either too heavy (SymPy) or purely numerical (SciPy).

Questions for Discussion

• How important is exact threshold detection vs approximate? Are there use cases where approximate timing creates problems?
• Should we prioritize simplicity (fewer options) or flexibility (support more patterns)?
• What's the right balance between "batteries included" (built-in rate functions) and "bring your own" (user-provided functions)?
• Is the performance overhead of numerical integration acceptable for the "long tail" of complex rate functions?
• How do we make the distinction between analytical and numerical methods clear to users without overwhelming them?
• Should this be in Mesa core or start as an experimental feature with room for API changes?

[tpike3]

This is awesome @EwoutH! Although more cumbersome, I like option 3. Critically, to me this is phase transition of complex systems is a key hard problem. By allowing option 3, we can crowd source the widest net of measures that find indicators for phase transitions. Whether this be specific measures or trained models. (I am a particular fan of Martin Scheffer's work in this area.)

To me the goal is find universal signals of phase transition of complex systems, and Mesa is serves as a knowledge bank of different measurements, so having options to increase awareness and an extensible format to let users add new ones to me is the goal.

Per your questions for discussion, I would add the base architecture of option 3, as experimental to tweak to the optimal format of the architecture, but then I would argue to let this component of Mesa be at "edge of chaos", by having a low threshold for what measures get in.

My thoughts, and excited for the discussion.

[quaquel]

Some quick responses from my side to the questions

• How important is exact threshold detection vs approximate? Are there use cases where approximate timing creates problems?

This directly interacts with how time is represented. In a classic ABM with fixed time steps, it's highly unlikely that the crossing of a treshold coincides exactly with a timestep. Most often, a threshold crossing will occur between two timesteps.

• Should we prioritize simplicity (fewer options) or flexibility (support more patterns)?

I would favor flexibility.

• What's the right balance between "batteries included" (built-in rate functions) and "bring your own" (user-provided functions)?

I would favor "bring your own" because of the wide variety of use cases. I doubt it is possible to cover them all. The main thing is to ensure that we have an easy to use and highly flexible structure for the signaling of threshold crossings.

• Is the performance overhead of numerical integration acceptable for the "long tail" of complex rate functions?

Given that I favor bring your own, this largely becomes a non-issue for MESA.

• How do we make the distinction between analytical and numerical methods clear to users without overwhelming them?

Most modelers will have a basic understanding of maths and should be familiar with the distinction so I would not worry about this too much.

• Should this be in Mesa core or start as an experimental feature with room for API changes?

Experimental, clearly.