SciML · ChrisRackauckas · Jun 4, 2026 · Jun 3, 2026 · Jun 4, 2026
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -3,4 +3,5 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Evolutionary = "86b6b26d-c046-49b6-aa0b-5f0f74682bd6"
 
 [compat]
-Documenter = "0.26"
+Documenter = "0.26, 1.17"
+Evolutionary = "0.11.2"
diff --git a/docs/src/crossover.md b/docs/src/crossover.md
@@ -68,19 +68,3 @@ SSX
 ```@docs
 Evolutionary.crosstree
 ```
-
-## References
-
-[^1]: H. Mühlenbein, D. Schlierkamp-Voosen, "Predictive Models for the Breeder Genetic Algorithm: I. Continuous Parameter Optimization". Evolutionary Computation, 1 (1), pp. 25-49, 1993.
-
-[^2]: K. V. Price and R. M. Storn and J. A. Lampinen, "Differential evolution: A practical approach to global optimization", Springer, 2005.
-
-[^3]: Z. Michalewicz, T. Logan,  S. Swaminathan. "Evolutionary operators for continuous convex parameter spaces." Proceedings of the 3rd Annual conference on Evolutionary Programming, 1994.
-
-[^4]: K. Deep, M. Thakur, "A new crossover operator for real coded genetic algorithms", Applied Mathematics and Computation 188, 2007, 895–912
-
-[^5]: K. Deep, K. P. Singh, M. L. Kansal, and C. Mohan, "A real coded  genetic algorithm for solving integer and mixed integer optimization problems.", Appl. Math. Comput. 212, 505-518, 2009
-
-[^6]: K. Deb, R. B. Agrawal, "Simulated Binary Crossover for Continuous Search Space", Complex Syst., 9., 1995
-
-[^7]: M. A. Wolters, “A Genetic Algorithm for Selection of Fixed-Size Subsets with Application to Design Problems”, J. Stat. Soft., vol. 68, no. 1, pp. 1–18, Nov. 2015.
diff --git a/docs/src/gp.md b/docs/src/gp.md
@@ -1,5 +1,9 @@
 # Genetic Programming
 
+```@meta
+CurrentModule = Evolutionary
+```
+
 ```@docs
 TreeGP
 ```
@@ -14,6 +18,24 @@ The [Genetic Programming](https://en.wikipedia.org/wiki/Genetic_programming) is
 Evolutionary.Expression
 Evolutionary.randterm
 Evolutionary.simplify!
+Base.rand(::AbstractRNG, ::TreeGP)
+```
+
+## Protected Functions
+
+The following protected operators guard against domain errors (e.g. division by
+zero, logarithm of zero) so that arbitrary expression trees can be evaluated safely.
+
+```@docs
+Evolutionary.pdiv
+Evolutionary.aq
+Evolutionary.pexp
+Evolutionary.plog
+Evolutionary.psqrt
+Evolutionary.psin
+Evolutionary.pcos
+Evolutionary.ppow
+Evolutionary.cond
 ```
 
 ## References

diff --git a/docs/src/mutation.md b/docs/src/mutation.md
@@ -23,13 +23,20 @@ cauchy(::AbstractVector, ::IsotropicStrategy)
 List of ES strategy mutation operations:
 
 ```@docs
+Evolutionary.nop(::AbstractStrategy)
 gaussian(::IsotropicStrategy)
 gaussian(::AnisotropicStrategy)
 ```
 
 
 ## Genetic Algorithm
 
+The default (no-op) mutation operator simply returns the individual unchanged.
+
+```@docs
+Evolutionary.genop(recombinant)
+```
+
 ### Binary Mutations
 
 ```@docs
@@ -75,23 +82,3 @@ point
 hoist
 shrink
 ```
-
-## References
-
-[^1]: Mühlenbein, H. and Schlierkamp-Voosen, D., "Predictive Models for the Breeder Genetic Algorithm: I. Continuous Parameter Optimization", Evolutionary Computation, 1 (1), 25-49, 1993.
-
-[^2]: Yao, Xin, and Yong Liu, "Fast evolution strategies", In International Conference on Evolutionary Programming, 149-161, Springer, 1997.
-
-[^3]: K. Deep, M. Thakur, "A new crossover operator for real coded genetic algorithms", Applied Mathematics and Computation 188, 895-912, 2007.
-
-[^4]: K. Deep, K. P. Singh, M. L. Kansal, and C. Mohan, "A real coded  genetic algorithm for solving integer and mixed integer optimization problems", Appl. Math. Comput. 212, 505-518, 2009
-
-[^5]: K. E. Kinnear, Jr., "Evolving a sort: Lessons in genetic programming", In Proceedings of the 1993 International Conference on Neural Networks, vol 2, 881-888, IEEE Press, 1993.
-
-[^6]: B. McKay, M. J. Willis, and G. W. Barton., "Using a tree structured genetic algorithm to perform symbolic regression", GALESIA, vol 414, 487-492, 1995.
-
-[^7]: K. E. Kinnear, Jr., "Fitness landscapes and difficulty in genetic programming", In Proceedings of the 1994 IEEE World Conference on Computational Intelligence, vol 1, 142-147, IEEE Press, 1994.
-
-[^8]: P. J. Angeline, "An investigation into the sensitivity of genetic programming to the frequency of leaf selection during subtree crossover", Genetic Programming 1996: Proceedings of the First Annual Conference, 21–29, 1996.
-
-[^9]: K. Deb, R. B. Agrawal, "Simulated Binary Crossover for Continuous Search Space", Complex Syst., 9., 1995
diff --git a/docs/src/selection.md b/docs/src/selection.md
@@ -32,12 +32,4 @@ random
 permutation
 randomoffset
 best
-```
-
-## References
-
-[^1]: Baker J.E., Adaptive selection methods for genetic algorithms, In Proceedings of International Conference on Genetic Algorithms and Their Applications, pp. 100-111, 1985.
-
-[^2]: Schwefel H.P., Evolution and Optimum Seeking, Wiley, New York, 1995.
-
-[^3]: Baker, J. E., Reducing Bias and Inefficiency in the Selection Algorithm. In [ICGA2], pp. 14-21, 1987.
+```
diff --git a/examples/Project.toml b/examples/Project.toml
@@ -10,3 +10,16 @@ PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
+
+[compat]
+Distributions = "0.25.125"
+Evolutionary = "0.11.2"
+Flux = "0.16.10"
+IJulia = "1.34.4"
+MLDatasets = "0.7.21"
+NLSolversBase = "8.0.0"
+Plots = "1.41.6"
+PyCall = "1.96.4"
+Revise = "3.14.5"
+StableRNGs = "1.0.4"
+SymPy = "2.3.3"
diff --git a/src/mutations.jl b/src/mutations.jl
@@ -55,6 +55,9 @@ Performs isotropic mutation of the recombinant `x` given the strategy `s`  by ad
 
 where ``\\delta`` is a Cauchy random variable with the scale parameter ``t = 1`` [^2].
 
+# References
+
+[^2]: Yao, Xin, and Yong Liu, "Fast evolution strategies", In International Conference on Evolutionary Programming, 149-161, Springer, 1997.
 """
 function cauchy(
         recombinant::AbstractVector, s::IsotropicStrategy;
@@ -155,6 +158,9 @@ The mutated individual is given by
 
 - ``x_i^\\prime = x_i + z_i``
 
+# References
+
+[^1]: Mühlenbein, H. and Schlierkamp-Voosen, D., "Predictive Models for the Breeder Genetic Algorithm: I. Continuous Parameter Optimization", Evolutionary Computation, 1 (1), 25-49, 1993.
 """
 function uniform(r::Real = 1.0)
     function mutation(
@@ -228,6 +234,10 @@ scheme within `lower` and `upper` bound, and an index of
 mutation `p`[^3].
 
 *Note:* The implementation is a degenerate case of Mixed Integer Power Mutation ([`MIPM`](@ref))
+
+# References
+
+[^3]: K. Deep, M. Thakur, "A new crossover operator for real coded genetic algorithms", Applied Mathematics and Computation 188, 895-912, 2007.
 """
 function PM(lower::Vector, upper::Vector, p::Float64 = 5.0) # index of distribution p
     return mipmmutation(lower, upper, p)
@@ -237,6 +247,10 @@ end
     MIPM(lower, upper, p_real = 10, p_int = 4)
 
 Returns an in-place real valued mutation function that performs the Mixed Integer Power Mutation (MI-PM) scheme within `lower` and `upper` bound, and an index of mutation `p_real` for real value and `p_int` for integer values[^4].
+
+# References
+
+[^4]: K. Deep, K. P. Singh, M. L. Kansal, and C. Mohan, "A real coded  genetic algorithm for solving integer and mixed integer optimization problems", Appl. Math. Comput. 212, 505-518, 2009
 """
 function MIPM(lowerBounds::Vector, upperBounds::Vector, p_real::Float64 = 10.0, p_int::Float64 = 4.0) # index of distribution p
     return mipmmutation(lowerBounds, upperBounds, p_real, p_int)
@@ -284,6 +298,10 @@ end
 
 Returns an in-place real valued mutation function that performs the Polynomial Mutation (PLM) scheme
 within `lower` and `upper` bounds, and a mutation distribution index `η`[^9].
+
+# References
+
+[^9]: K. Deb, R. B. Agrawal, "Simulated Binary Crossover for Continuous Search Space", Complex Syst., 9., 1995
 """
 function PLM(Δ::Union{Real, Vector} = 1.0; η = 2, pm::Real = NaN) # index of distribution p
     function mutation(recombinant::T; rng::AbstractRNG = default_rng()) where {T <: AbstractVector}
@@ -457,6 +475,10 @@ Returns an in-place expression mutation function that performs mutation of an ar
 Parameters:
 - `growth`: Growth restriction on the offspring in comparison to the parent.
             The offspring cannot be more than `growth`% deeper than its parent.  (default: `0.0`)
+
+# References
+
+[^5]: K. E. Kinnear, Jr., "Evolving a sort: Lessons in genetic programming", In Proceedings of the 1993 International Conference on Neural Networks, vol 2, 881-888, IEEE Press, 1993.
 """
 function subtree(method::TreeGP; growth::Real = 0.0)
     function mutation(recombinant::Expr; rng::AbstractRNG = default_rng())
@@ -481,6 +503,10 @@ end
 Returns an in-place expression mutation function that replaces an arbitrary node in the tree by the randomly selected one.
 Node replacement mutation is similar to bit string mutation in that it randomly changes a point in the individual.
 To ensure the tree remains legal, the replacement node has the same number of arguments as the node it is replacing [^6].
+
+# References
+
+[^6]: B. McKay, M. J. Willis, and G. W. Barton., "Using a tree structured genetic algorithm to perform symbolic regression", GALESIA, vol 414, 487-492, 1995.
 """
 function point(method::TreeGP)
     function mutation(recombinant::Expr; rng::AbstractRNG = default_rng())
@@ -506,6 +532,10 @@ end
 Returns an in-place expression mutation function that creates a new offspring individual which is copy
 of a randomly chosen subtree of the parent. Thus, the offspring will be smaller than the parent
 and will have a different root node [^7].
+
+# References
+
+[^7]: K. E. Kinnear, Jr., "Fitness landscapes and difficulty in genetic programming", In Proceedings of the 1994 IEEE World Conference on Computational Intelligence, vol 1, 142-147, IEEE Press, 1994.
 """
 function hoist(method::TreeGP)
     function mutation(recombinant::Expr; rng::AbstractRNG = default_rng())
@@ -528,6 +558,10 @@ end
 Returns an in-place expression mutation function that replaces a randomly chosen subtree with a randomly
 created terminal. This is a special case of subtree mutation where the replacement tree is a terminal.
 As with hoist mutation, it is motivated by the desire to reduce program size [^8].
+
+# References
+
+[^8]: P. J. Angeline, "An investigation into the sensitivity of genetic programming to the frequency of leaf selection during subtree crossover", Genetic Programming 1996: Proceedings of the First Annual Conference, 21–29, 1996.
 """
 function shrink(method::TreeGP)
     function mutation(recombinant::Expr; rng::AbstractRNG = default_rng())

diff --git a/src/recombinations.jl b/src/recombinations.jl
@@ -135,6 +135,10 @@ end
 Returns a uniform (binomial) crossover function, see [Recombination Interface](@ref), function with the probability `Cr` [^2].
 
 The crossover probability value must be in unit interval, ``Cr \\in [0,1]``.
+
+# References
+
+[^2]: K. V. Price and R. M. Storn and J. A. Lampinen, "Differential evolution: A practical approach to global optimization", Springer, 2005.
 """
 function BINX(Cr::Real = 0.5)
     function binxvr(v1::T, v2::T; rng::AbstractRNG = default_rng()) where {T <: AbstractVector}
@@ -184,8 +188,12 @@ end
     BSX(k::Int)
 
 Binary Subset Crossover[^7]. Produces an offspring by first pooling the unique
-items of the two parents, and then creating each offspring by sampling without 
+items of the two parents, and then creating each offspring by sampling without
 replacement at most `k` elements from the pool of items.
+
+# References
+
+[^7]: M. A. Wolters, “A Genetic Algorithm for Selection of Fixed-Size Subsets with Application to Design Problems”, J. Stat. Soft., vol. 68, no. 1, pp. 1–18, Nov. 2015.
 """
 function BSX(k::Int)
     function BSX(v1::T, v2::T; rng::AbstractRNG = default_rng()) where {T <: AbstractVector{Bool}}
@@ -298,6 +306,10 @@ Heuristic crossover (HX) recombination operation[^3] generates offspring `u` and
 - ``v = y + r (y - x)``
 
 where ``r`` is chosen uniform randomly in the interval ``[0;1)``.
+
+# References
+
+[^3]: Z. Michalewicz, T. Logan,  S. Swaminathan. "Evolutionary operators for continuous convex parameter spaces." Proceedings of the 3rd Annual conference on Evolutionary Programming, 1994.
 """
 function HX(v1::T, v2::T; rng::AbstractRNG = default_rng()) where {T <: AbstractVector}
     c1 = v1 .+ rand(rng) * (v1 .- v2)
@@ -309,6 +321,10 @@ end
     LX(μ::Real = 0.0, b::Real = 0.2)
 
 Returns a Laplace crossover (LX) recombination operation[^4], see [Recombination Interface](@ref).
+
+# References
+
+[^4]: K. Deep, M. Thakur, "A new crossover operator for real coded genetic algorithms", Applied Mathematics and Computation 188, 2007, 895–912
 """
 function LX(μ::Real = 0.0, b::Real = 0.2) # location μ, scale b > 0
     function lxxvr(v1::T, v2::T; rng::AbstractRNG = default_rng()) where {T <: AbstractVector}
@@ -326,6 +342,10 @@ end
     MILX(μ::Real = 0.0, b_real::Real = 0.15, b_int::Real = 0.35)
 
 Returns a mixed integer Laplace crossover (MI-LX) recombination operation[^5], see [Recombination Interface](@ref).
+
+# References
+
+[^5]: K. Deep, K. P. Singh, M. L. Kansal, and C. Mohan, "A real coded  genetic algorithm for solving integer and mixed integer optimization problems.", Appl. Math. Comput. 212, 505-518, 2009
 """
 function MILX(μ::Real = 0.0, b_real::Real = 0.15, b_int::Real = 0.35) # location μ, scale b > 0
     function milxxvr(v1::T, v2::T; rng::AbstractRNG = default_rng()) where {T <: AbstractVector}
@@ -347,6 +367,10 @@ end
 
 Returns a Simulated Binary Crossover (SBX) recombination operation, see [Recombination Interface](@ref),
 with the mutation probability `pm` of the recombinant component, and is the crossover distribution index `η`[^6].
+
+# References
+
+[^6]: K. Deb, R. B. Agrawal, "Simulated Binary Crossover for Continuous Search Space", Complex Syst., 9., 1995
 """
 function SBX(pm::Real = 0.5, η::Integer = 2)
     function sbxv(v1::T, v2::T; rng::AbstractRNG = default_rng()) where {T <: AbstractVector}

diff --git a/src/selections.jl b/src/selections.jl
@@ -14,6 +14,9 @@ Consider ``M`` the number of individuals in the population, ``P`` the position o
 
 Linear ranking allows values of selective pressure in [1.0, 2.0].
 
+# References
+
+[^1]: Baker J.E., Adaptive selection methods for genetic algorithms, In Proceedings of International Conference on Genetic Algorithms and Their Applications, pp. 100-111, 1985.
 """
 function ranklinear(sp::Real)
     @assert 1.0 <= sp <= 2.0 "Selective pressure has to be in range [1.0, 2.0]."
@@ -38,6 +41,10 @@ end
 Returns a (μ, λ)-uniform ranking selection function, see [Selection Interface](@ref) with the best individuals parameter `μ`.
 
 In uniform ranking, the best ``\\mu`` individuals are assigned a selection probability of ``1/\\mu`` while the rest them are discarded [^2].
+
+# References
+
+[^2]: Schwefel H.P., Evolution and Optimum Seeking, Wiley, New York, 1995.
 """
 function uniformranking(μ::Int)
     function uniformrank(
@@ -93,6 +100,9 @@ Consider ``N`` the number of individuals to be selected, then the distance betwe
 
 *Note*: Best used in maximization context.
 
+# References
+
+[^3]: Baker, J. E., Reducing Bias and Inefficiency in the Selection Algorithm. In [ICGA2], pp. 14-21, 1987.
 """
 function sus(fitness::Vector{<:Real}, N::Int; rng::AbstractRNG = default_rng())
     F = sum(abs, fitness)