Add comprehensive API documentation and restructure docs organization

src/FiniteDiff.jl

@@ -9,26 +9,126 @@ using LinearAlgebra, ArrayInterface
 
 import Base: resize!
 
+"""
+    _vec(x)
+
+Internal utility function to vectorize arrays while preserving scalars.
+
+# Arguments
+- `x`: Array or scalar
+
+# Returns
+- `vec(x)` for arrays, `x` unchanged for scalars
+"""
 _vec(x) = vec(x)
 _vec(x::Number) = x
 
+"""
+    _mat(x)
+
+Internal utility function to ensure matrix format.
+
+Converts vectors to column matrices while preserving existing matrices.
+Used internally to ensure consistent matrix dimensions for operations.
+
+# Arguments  
+- `x`: Matrix or vector
+
+# Returns
+- `x` unchanged if already a matrix
+- Reshaped column matrix if `x` is a vector
+"""
 _mat(x::AbstractMatrix) = x
 _mat(x::AbstractVector) = reshape(x, (axes(x, 1), Base.OneTo(1)))
 
 # Setindex overloads without piracy
 setindex(x...) = Base.setindex(x...)
 
+"""
+    setindex(x::AbstractArray, v, i...)
+
+Non-mutating setindex operation that returns a copy with modified elements.
+
+Creates a mutable copy of array `x`, sets the specified indices to value `v`,
+and returns the modified copy. This avoids type piracy while providing 
+setindex functionality for immutable arrays.
+
+# Arguments
+- `x::AbstractArray`: Array to modify (not mutated)
+- `v`: Value to set at the specified indices
+- `i...`: Indices where the value should be set
+
+# Returns
+- Modified copy of `x` with `x[i...] = v`
+
+# Examples
+```julia
+x = [1, 2, 3]
+y = setindex(x, 99, 2)  # y = [1, 99, 3], x unchanged
+```
+"""
 function setindex(x::AbstractArray, v, i...)
     _x = Base.copymutable(x)
     _x[i...] = v
     return _x
 end
 
+"""
+    setindex(x::AbstractVector, v, i::Int)
+
+Broadcasted setindex operation for vectors using boolean masking.
+
+Sets the i-th element of vector `x` to value `v` using broadcasting operations.
+This implementation uses boolean masks to avoid explicit copying and provide
+efficient vectorized operations.
+
+# Arguments
+- `x::AbstractVector`: Input vector  
+- `v`: Value to set at index `i`
+- `i::Int`: Index to modify
+
+# Returns
+- Vector with `x[i] = v`, computed via broadcasting
+
+# Examples
+```julia
+x = [1.0, 2.0, 3.0]
+y = setindex(x, 99.0, 2)  # [1.0, 99.0, 3.0]
+```
+"""
 function setindex(x::AbstractVector, v, i::Int)
     n = length(x)
     x .* (i .!== 1:n) .+ v .* (i .== 1:n)
 end
 
+"""
+    setindex(x::AbstractMatrix, v, i::Int, j::Int)
+
+Broadcasted setindex operation for matrices using boolean masking.
+
+Sets the (i,j)-th element of matrix `x` to value `v` using broadcasting operations.
+This implementation uses boolean masks to avoid explicit copying and provide
+efficient vectorized operations.
+
+# Arguments
+- `x::AbstractMatrix`: Input matrix
+- `v`: Value to set at position (i,j)
+- `i::Int`: Row index to modify
+- `j::Int`: Column index to modify
+
+# Returns
+- Matrix with `x[i,j] = v`, computed via broadcasting
+
+# Examples
+```julia
+x = [1.0 2.0; 3.0 4.0]
+y = setindex(x, 99.0, 1, 2)  # [1.0 99.0; 3.0 4.0]
+```
+
+# Notes
+The implementation uses transposed broadcasting `(j .!== i:m)'` which appears
+to be a typo - should likely be `(j .!== 1:m)'` for correct column masking.
+"""
 function setindex(x::AbstractMatrix, v, i::Int, j::Int)
     n, m = Base.size(x)
     x .* (i .!== 1:n) .* (j .!== i:m)' .+ v .* (i .== 1:n) .* (j .== i:m)'

src/epsilons.jl

@@ -2,28 +2,132 @@
 Very heavily inspired by Calculus.jl, but with an emphasis on performance and DiffEq API convenience.
 =#
 
-#=
-Compute the finite difference interval epsilon.
+"""
+    compute_epsilon(::Val{:forward}, x::T, relstep::Real, absstep::Real, dir::Real) where T<:Number
+
+Compute the finite difference step size (epsilon) for forward finite differences.
+
+The step size is computed as `max(relstep*abs(x), absstep)*dir`, which ensures 
+numerical stability by using a relative step scaled by the magnitude of `x` 
+when `x` is large, and an absolute step when `x` is small.
+
+# Arguments
+- `::Val{:forward}`: Finite difference type indicator for forward differences
+- `x::T`: Point at which to compute the step size
+- `relstep::Real`: Relative step size factor
+- `absstep::Real`: Absolute step size fallback
+- `dir::Real`: Direction multiplier (typically ±1)
+
+# Returns
+- Step size `ϵ` for forward finite difference: `f(x + ϵ)`
+
 Reference: Numerical Recipes, chapter 5.7.
-=#
-@inline function compute_epsilon(::Val{:forward}, x::T, relstep::Real, absstep::Real, dir::Real) where T<:Number
+"""
+@inline function compute_epsilon(
+        ::Val{:forward}, x::T, relstep::Real, absstep::Real, dir::Real) where {T <: Number}
     return max(relstep*abs(x), absstep)*dir
 end
 
-@inline function compute_epsilon(::Val{:central}, x::T, relstep::Real, absstep::Real, dir=nothing) where T<:Number
+"""
+    compute_epsilon(::Val{:central}, x::T, relstep::Real, absstep::Real, dir=nothing) where T<:Number
+
+Compute the finite difference step size (epsilon) for central finite differences.
+
+The step size is computed as `max(relstep*abs(x), absstep)`, which ensures 
+numerical stability by using a relative step scaled by the magnitude of `x` 
+when `x` is large, and an absolute step when `x` is small.
+
+# Arguments
+- `::Val{:central}`: Finite difference type indicator for central differences
+- `x::T`: Point at which to compute the step size
+- `relstep::Real`: Relative step size factor
+- `absstep::Real`: Absolute step size fallback
+- `dir`: Direction parameter (unused for central differences)
+
+# Returns
+- Step size `ϵ` for central finite difference: `(f(x + ϵ) - f(x - ϵ)) / (2ϵ)`
+"""
+@inline function compute_epsilon(::Val{:central}, x::T, relstep::Real,
+        absstep::Real, dir = nothing) where {T <: Number}
     return max(relstep*abs(x), absstep)
 end
 
-@inline function compute_epsilon(::Val{:hcentral}, x::T, relstep::Real, absstep::Real, dir=nothing) where T<:Number
+"""
+    compute_epsilon(::Val{:hcentral}, x::T, relstep::Real, absstep::Real, dir=nothing) where T<:Number
+
+Compute the finite difference step size (epsilon) for central finite differences in Hessian computations.
+
+The step size is computed as `max(relstep*abs(x), absstep)`, which ensures 
+numerical stability by using a relative step scaled by the magnitude of `x` 
+when `x` is large, and an absolute step when `x` is small.
+
+# Arguments
+- `::Val{:hcentral}`: Finite difference type indicator for Hessian central differences
+- `x::T`: Point at which to compute the step size
+- `relstep::Real`: Relative step size factor
+- `absstep::Real`: Absolute step size fallback
+- `dir`: Direction parameter (unused for central differences)
+
+# Returns
+- Step size `ϵ` for Hessian central finite differences
+"""
+@inline function compute_epsilon(::Val{:hcentral}, x::T, relstep::Real,
+        absstep::Real, dir = nothing) where {T <: Number}
     return max(relstep*abs(x), absstep)
 end
 
-@inline function compute_epsilon(::Val{:complex}, x::T, ::Union{Nothing,T}=nothing, ::Union{Nothing,T}=nothing, dir=nothing) where T<:Real
+"""
+    compute_epsilon(::Val{:complex}, x::T, ::Union{Nothing,T}=nothing, ::Union{Nothing,T}=nothing, dir=nothing) where T<:Real
+
+Compute the finite difference step size (epsilon) for complex step differentiation.
+
+For complex step differentiation, the step size is simply the machine epsilon `eps(T)`,
+which provides optimal accuracy since complex step differentiation doesn't suffer from
+subtractive cancellation errors.
+
+# Arguments
+- `::Val{:complex}`: Finite difference type indicator for complex step differentiation
+- `x::T`: Point at which to compute the step size (unused, type determines epsilon)
+- Additional arguments are unused for complex step differentiation
+
+# Returns
+- Machine epsilon `eps(T)` for complex step differentiation: `imag(f(x + iϵ)) / ϵ`
+
+# Notes
+Complex step differentiation computes derivatives as `imag(f(x + iϵ)) / ϵ` where `ϵ = eps(T)`.
+This method provides machine precision accuracy without subtractive cancellation.
+"""
+@inline function compute_epsilon(::Val{:complex}, x::T, ::Union{Nothing, T} = nothing,
+        ::Union{Nothing, T} = nothing, dir = nothing) where {T <: Real}
     return eps(T)
 end
 
-default_relstep(::Type{V}, T) where V = default_relstep(V(), T)
-@inline function default_relstep(::Val{fdtype}, ::Type{T}) where {fdtype,T<:Number}
+"""
+    default_relstep(fdtype, ::Type{T}) where T<:Number
+
+Compute the default relative step size for finite difference approximations.
+
+Returns optimal default step sizes based on the finite difference method and 
+numerical type, balancing truncation error and round-off error.
+
+# Arguments
+- `fdtype`: Finite difference type (`Val(:forward)`, `Val(:central)`, `Val(:hcentral)`, etc.)
+- `::Type{T}`: Numerical type for which to compute the step size
+
+# Returns
+- `sqrt(eps(real(T)))` for forward differences
+- `cbrt(eps(real(T)))` for central differences  
+- `eps(T)^(1/4)` for Hessian central differences
+- `one(real(T))` for other types
+
+# Notes
+These step sizes minimize the total error (truncation + round-off) for each method:
+- Forward differences have O(h) truncation error, optimal h ~ sqrt(eps)
+- Central differences have O(h²) truncation error, optimal h ~ eps^(1/3)
+- Hessian methods have O(h²) truncation error but involve more operations
+"""
+default_relstep(::Type{V}, T) where {V} = default_relstep(V(), T)
+@inline function default_relstep(::Val{fdtype}, ::Type{T}) where {fdtype, T <: Number}
     if fdtype==:forward
         return sqrt(eps(real(T)))
     elseif fdtype==:central
@@ -35,7 +139,20 @@ default_relstep(::Type{V}, T) where V = default_relstep(V(), T)
     end
 end
 
-function fdtype_error(::Type{T}=Float64) where T
+"""
+    fdtype_error(::Type{T}=Float64) where T
+
+Throw an informative error for unsupported finite difference type combinations.
+
+# Arguments
+- `::Type{T}`: Return type of the function being differentiated
+
+# Errors
+- For `Real` return types: suggests `Val{:forward}`, `Val{:central}`, `Val{:complex}`
+- For `Complex` return types: suggests `Val{:forward}`, `Val{:central}` (no complex step)
+- For other types: suggests the return type should be Real or Complex subtype
+"""
+function fdtype_error(::Type{T} = Float64) where {T}
     if T<:Real
         error("Unrecognized fdtype: valid values are Val{:forward}, Val{:central} and Val{:complex}.")
     elseif T<:Complex