module LocalHillResponses

using ForwardDiff

export local_sensitivity, respose_coefficient, local_hills

"""
    local_sensitivity(f, s)

Compute the local sensitivity (`df/ds`) of the function `f` for the singular input `s`.
"""
function local_sensitivity(f, s)
    return ForwardDiff.derivative(f, s)
end

"""
    respose_coefficient(f, s)

Compute the response coefficient (a.k.a. relative local sensitivity): 

```math
R^{x}_{s} = \frac{dx}{ds} \frac{s}{x}
```

where ``R^{x}_{s}`` is the response coefficient of the input ``s`` and the output ``x = f(s)``.

# Warning

The result may be fluctuated due to the singularity point.
"""
function respose_coefficient(f, s, y = f(s))
    return local_sensitivity(f, s) * s / y
end

const relative_local_sensitivity = respose_coefficient

"""
Check if function `f(s)` is monotonous between `smin` and `smax`. A positive `direction` checks increasing function; A negative `direction` checks decreasing function. If no `direction` is given, checks either increasing or decreasing.
"""
function is_monotonous(f, smin, smax, direction; Δs = (smax - smin) * 1e-3)
    for s in smin:Δs:smax
        if ForwardDiff.derivative(f, s) * direction < 0
            return false
        end
    end
    return true
end

function is_monotonous(f, smin, smax; Δs = (smax - smin) * 1e-3)
    is_monotonous(f, smin, smax, +1; Δs = Δs) || is_monotonous(f, smin, smax, -1; Δs = Δs)
end

"""
Activated fraction for reponse `y` bounded by minimal response `ymin` and maximal response `ymax`.
"""
activated_fraction(y, ymin, ymax) = (y - ymin) / (ymax - ymin)

"""
Local Hill coefficient.
"""
function local_hills(f, s, ymin, ymax)
    responseCoef = respose_coefficient(f, s)
    actFrac = activated_fraction(y, ymin, ymax)
    nh = responseCoef / (one(actFrac) - actFrac)
    return nh
end



end # module
