Update docstring references

As a follow-up (patch?) to #47, we remove the Markdown-style [^ABC] references from the docstrings and keep the '@cite' formatting of Documenter. However, we maintain the separate 'References' sections at the bottom of each docstring that actually provide the full citation for immediate readability in the REPL and in the source code.

Author: Luis-Varona

CHANGELOG.md

@@ -6,6 +6,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 
 ## [Unreleased]
 
+### Added
+
+- Added **References** sections to docstrings for immediate readability in the REPL and in the source code without needing to open the Documenter-generated website (#47, #48).
+
 ## [0.1.1] - 2025-08-05
 
 ### Changed

src/core.jl

@@ -8,7 +8,13 @@
     s_bandwidth(g::AbstractGraph, S) -> SBandMinimizationResult
     s_bandwidth(L::AbstractMatrix{<:Integer}, S) -> SBandMinimizationResult
 
-[TODO: Write here. Also, comment inline and cite [JP25](@cite).]
+[TODO: Write here. Also, comment inline and cite [JP25].]
+
+# References
+
+- [JP25](@cite): N. Johnston and S. Plosker. *Laplacian {−1,0,1}- and {−1,1}-diagonalizable
+    graphs*. *Linear Algebra and its Applications* **704**, 309–39 (2025).
+    https://doi.org/10.1016/j.laa.2024.10.016.
 """
 function s_bandwidth(g::AbstractGraph, S::Tuple{Vararg{Integer}})
     _assert_graph_has_defined_s_bandwidth(g)
@@ -75,7 +81,13 @@ end
     has_s_bandwidth_at_most_k(g::AbstractGraph, S, k) -> SBandRecognitionResult
     has_s_bandwidth_at_most_k(L::AbstractMatrix{<:Integer}, S, k) -> SBandRecognitionResult
 
-[TODO: Write here. Also, comment inline and cite [JP25](@cite).]
+[TODO: Write here. Also, comment inline and cite [JP25].]
+
+# References
+
+- [JP25](@cite): N. Johnston and S. Plosker. *Laplacian {−1,0,1}- and {−1,1}-diagonalizable
+    graphs*. *Linear Algebra and its Applications* **704**, 309–39 (2025).
+    https://doi.org/10.1016/j.laa.2024.10.016.
 """
 function has_s_bandwidth_at_most_k(g::AbstractGraph, S::Tuple{Vararg{Integer}}, k::Integer)
     _assert_graph_has_defined_s_bandwidth(g)

src/eigenvector_generation.jl

@@ -102,13 +102,18 @@ julia> hcat(SDiagonalizability._pot_kernel_01neg_eigvecs(3)...)
 
 # Notes
 The number of potential kernel eigenvectors (unique up to span) for an order ``n`` Laplacian
-matrix is given by ``(3ⁿ - 1) / 2``. See also the relevant OEIS sequence [Slo25](@cite).
+matrix is given by ``(3ⁿ - 1) / 2``. See also the relevant OEIS sequence [Slo25].
 
 Regrettably, the implementation here is rather clunky and unidiomatic, but it is worth
 noting that eigenvector generation is one of two major bottlenecks in the overall
 *S*-bandwidth minimization algorithm. Given how much potential there is for optimization in
 this piece of code, we thus prioritize performance over readability in this particular case,
 making every effort to include inline comments wherever clarification may be needed.
+
+# References
+
+- [Slo25](@cite): N. J. Sloane, *a(n) = (3^n - 1)/2*. Entry A003462 (2025). Accessed:
+    2025-05-22. https://oeis.org/A003462.
 """
 function _pot_kernel_01neg_eigvecs(n::Integer)
     # Cache to avoid redundant recomputations of the `leading` vector
@@ -197,13 +202,18 @@ julia> hcat(SDiagonalizability._pot_nonkernel_01neg_eigvecs(4)...)
 # Notes
 The number of potential non-kernel eigenvectors (unique up to span) for an order ``n``
 Laplacian matrix is, by non-trivial combinatorial arguments, equal to the number of humps in
-all Motzkin paths of length ``n``. See also the relevant OEIS sequence [Deu25](@cite).
+all Motzkin paths of length ``n``. See also the relevant OEIS sequence [Deu25].
 
 Regrettably, the implementation here is rather clunky and unidiomatic, but it is worth
 noting that eigenvector generation is one of two major bottlenecks in the overall
 *S*-bandwidth minimization algorithm. Given how much potential there is for optimization in
 this piece of code, we thus prioritize performance over readability in this particular case,
 making every effort to include inline comments wherever clarification may be needed.
+
+# References
+
+- [Deu25](@cite): E. Deutsch. *Number of humps in all Motzkin paths of length n*. Entry
+    A097861 (2025). Accessed: 2025-05-22. https://oeis.org/A097861.
 """
 function _pot_nonkernel_01neg_eigvecs(n::Integer)
     # Caches to avoid redundant recomputations of the `leading` and `entries` vectors

src/factories/orthogonality_factory.jl

@@ -14,9 +14,9 @@ Recall that an (ordered) collection of vectors ``v₁, v₂, ..., vₙ`` is said
 (i.e., if every pair of vectors at least ``k`` indices apart is orthogonal). This is
 equivalent to the vectors' Gram matrix having bandwidth at most `k`, where we define the
 bandwidth of a matrix ``A`` to be the minimum integer ``k ∈ \\{1, 2, …, n\\}`` such that
-``Aᵢⱼ = 0`` whenever ``|i - j| ≥ k`` [JP25; p. 313](@cite). (Note that many texts instead
-define matrix bandwidth using zero-based indexing—that is, with the condition
-``|i - j| > k`` [Maf14; p. 186](@cite).)
+``Aᵢⱼ = 0`` whenever ``|i - j| ≥ k`` [JP25; p. 313]. (Note that many texts instead define
+matrix bandwidth using zero-based indexing—that is, with the condition ``|i - j| > k``
+[Maf14; p. 186].)
 
 This type is used as a template for concretely defined properties corresponding to specific
 values of ``k``. In the context of the overarching *S*-bandwidth algorithm, we perform a
@@ -26,6 +26,15 @@ different depth-first search for each family of values of ``k`` on our "tree" of
 # Interface
 Concrete subtypes of `KOrthogonality` **must** implement the following fields:
 - `k::Int`: the ``k``-orthogonality parameter. Must be a positive integer.
+
+# References
+
+- [JP25](@cite): N. Johnston and S. Plosker. *Laplacian {−1,0,1}- and {−1,1}-diagonalizable
+    graphs*. *Linear Algebra and its Applications* **704**, 309–39 (2025).
+    https://doi.org/10.1016/j.laa.2024.10.016.
+- [Maf14](@cite): L. O. Mafteiu-Scai. *The Bandwidths of a Matrix. A Survey of Algorithms*.
+    *Annals of West University of Timisoara - Mathematics and Computer Science* **52**,
+    183–223 (2014). https://doi.org/10.2478/awutm-2014-0019.
 """
 abstract type KOrthogonality end
 
@@ -61,7 +70,7 @@ The property of quasi-orthogonality for a collection of vectors.
 Recall that an (ordered) collection of vectors ``v₁, v₂, ..., vₙ`` is said to be
 *quasi-orthogonal* if we have the inner product ``⟨vᵢ, vⱼ⟩ = 0`` whenever ``|i - j| ≥ 2``
 (i.e., if every pair of vectors at least ``2`` indices apart is orthogonal). This is
-equivalent to the vectors' Gram matrix being tridiagonal [JP25; p. 313](@cite).
+equivalent to the vectors' Gram matrix being tridiagonal [JP25; p. 313].
 
 # Fields
 - `k::Int`: the ``k``-orthogonality parameter; always necessarily ``2``.
@@ -71,6 +80,12 @@ equivalent to the vectors' Gram matrix being tridiagonal [JP25; p. 313](@cite).
 
 # Constructors
 - `QuasiOrthogonality()`: constructs a new `QuasiOrthogonality` object with `k = 2`.
+
+# References
+
+- [JP25](@cite): N. Johnston and S. Plosker. *Laplacian {−1,0,1}- and {−1,1}-diagonalizable
+    graphs*. *Linear Algebra and its Applications* **704**, 309–39 (2025).
+    https://doi.org/10.1016/j.laa.2024.10.016.
 """
 struct QuasiOrthogonality <: KOrthogonality
     k::Int

src/laplacian_s_spectra.jl

@@ -46,7 +46,7 @@ constructed as well.
 
 # Examples
 Confirm that the rotation matrix by ``π/2`` radians counterclockwise is not spectrum
-integral (rather, it has eigenvalues ``±i`` [Joy15; p. 1](@cite)):
+integral (rather, it has eigenvalues ``±i`` [Joy15; p. 1]):
 ```jldoctest
 julia> R = Int8.([0 -1; 1 0])
 2×2 Matrix{Int8}:
@@ -68,7 +68,7 @@ true
 ```
 
 Confirm that the adjacency matrix of the Petersen graph is spectrum integral, with correct
-eigenvalues and multiplicities of ``\\{3: 1, -2: 4, 1: 5\\}`` [Fox09; p. 2](@cite):
+eigenvalues and multiplicities of ``\\{3: 1, -2: 4, 1: 5\\}`` [Fox09; p. 2]:
 ```jldoctest
 julia> using Graphs
 
@@ -116,8 +116,19 @@ OrderedCollections.OrderedDict{Int64, Int64} with 3 entries:
 # Notes
 If an undirected graph with integer edge weights is ``\\{-1, 0, 1\\}``-diagonalizable (or,
 more restrictively, ``\\{-1, 1\\}``-diagonalizable), then its Laplacian matrix has integer
-eigenvalues [JP25; p. 312](@cite). Hence, validating Laplacian integrality serves as a
-useful screening step in this package's principal *S*-bandwidth minimization algorithm.
+eigenvalues [JP25; p. 312]. Hence, validating Laplacian integrality serves as a useful
+screening step in this package's principal *S*-bandwidth minimization algorithm.
+
+# References
+
+- [Fox09](@cite): J. Fox. *Lecture 19: The Petersen graph and Moore graphs*. Lecture notes,
+    MAT 307: Combinatorics (2009). Accessed: 2025-07-25.
+    https://math.mit.edu/~fox/MAT307.html.
+- [Joy15](@cite): D. Joyce. *Rotations and complex eigenvalues*. Lecture notes, Math 130:
+    Linear Algebra (2015). http://aleph0.clarku.edu/~ma130/complexeigen.pdf.
+- [JP25](@cite): N. Johnston and S. Plosker. *Laplacian {−1,0,1}- and {−1,1}-diagonalizable
+    graphs*. *Linear Algebra and its Applications* **704**, 309–39 (2025).
+    https://doi.org/10.1016/j.laa.2024.10.016.
 """
 function check_spectrum_integrality(A::AbstractMatrix{<:Integer})
     A_copy = Matrix{Int}(A) # Avoid shared mutability and cast to `Matrix{Int}`

src/types.jl

@@ -173,8 +173,14 @@ eigenvalues are indeed all integers. (Otherwise, the associated field is simply
 # Notes
 If an undirected graph with integer edge weights is ``\\{-1, 0, 1\\}``-diagonalizable (or,
 more restrictively, ``\\{-1, 1\\}``-diagonalizable), then its Laplacian matrix has integer
-eigenvalues [JP25; p. 312](@cite). Hence, validating Laplacian integrality serves as a
-useful screening step in this package's principal *S*-bandwidth minimization algorithm.
+eigenvalues [JP25; p. 312]. Hence, validating Laplacian integrality serves as a useful
+screening step in this package's principal *S*-bandwidth minimization algorithm.
+
+# References
+
+- [JP25](@cite): N. Johnston and S. Plosker. *Laplacian {−1,0,1}- and {−1,1}-diagonalizable
+    graphs*. *Linear Algebra and its Applications* **704**, 309–39 (2025).
+    https://doi.org/10.1016/j.laa.2024.10.016.
 """
 struct SpectrumIntegralResult{T<:Union{Nothing,OrderedDict{Int,Int}}}
     matrix::AbstractMatrix{<:Integer}

src/utils.jl

@@ -39,7 +39,7 @@ an SVD), while the pivots are used to extract a spanning set of independent colu
 
 The rank-revealing Businger–Golub QR algorithm is used for the pivoting strategy, appending
 the "most independent" column with respect to the current set of pivots at each step via
-Householder transformations [BG65; pp. 269--70](@cite).
+Householder transformations [BG65; pp. 269--70].
 
 # Arguments
 - `A::AbstractMatrix{T<:Integer}`: the matrix whose independent columns to extract.
@@ -75,9 +75,9 @@ julia> SDiagonalizability._extract_independent_cols(A)
 Since we already need a pivoted QR decomposition to identify independent columns of `A` (or,
 rather, to order the columns in such a way that the first `rank(A)` ones are guaranteed to
 be independent), it makes sense to use data from the resulting factorization object to
-compute the rank of `A` rather than compute a separate SVD. We thus count the nonzero scaling
-coefficients—that is, the diagonal entries of the `R` matrix in `A = QR`—to determine the
-rank, similarly to how we count the nonzero singular values in an SVD.
+compute the rank of `A` rather than compute a separate SVD. We thus count the nonzero
+scaling coefficients—that is, the diagonal entries of the `R` matrix in `A = QR`—to
+determine the rank, similarly to how we count the nonzero singular values in an SVD.
 
 It is worth noting that we manually specify a higher relative tolerance for this rank
 computation. Further discussion can be found in the [`_rank_rtol`](@ref) documentation, but
@@ -87,7 +87,18 @@ tall-and-skinny and short-and-fat matrices (precisely the type we expect to enco
 dealing with all ``\\{-1, 0, 1\\}``-eigenvectors of a Laplacian matrix, which is the
 intended use case of this helper function in this package). Our replacement tolerance, on
 the other hand, is a widely accepted standard in numerical analysis which uses the maximum
-dimension instead [PTVF07; p. 795](@cite).
+dimension instead [PTVF07; p. 795].
+
+# References
+
+- [BG65](@cite): P. Businger and G. H. Golub. *Linear Least Squares Solutions by Householder
+    Transformations*. *Numerische Mathematik* **7**, 269–76 (1965).
+    https://doi.org/10.1007/BF01436084.
+
+- [PTVF07](@cite): W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery.
+    *Numerical Recipes: The Art of Scientific Computing*. 3rd Edition (Cambridge University
+    Press, Cambridge, UK, 2007). ISBN: 978-0-521-88068-8.
+    https://dl.acm.org/doi/10.5555/1403886.
 """
 function _extract_independent_cols(A::AbstractMatrix{<:Integer})
     F = qr(A, ColumnNorm())
@@ -156,7 +167,7 @@ Matrix has nonzero row sums; cannot be an (undirected) Laplacian
 
 Both the in-degree and out-degree Laplacian matrices of this random tournament digraph have
 zero row sums but are not symmetric, so they fail the check. (These are the two standard
-ways of extending the concept of the Laplacian to directed graphs [VL20; p. 196](@cite).)
+ways of extending the concept of the Laplacian to directed graphs [VL20; p. 196].)
 ```jldoctest
 julia> using Graphs
 
@@ -204,6 +215,12 @@ At first blush, it may seem as though the choice of `DomainError` over something
 this is informed by the simple *ad hoc* use of this function to validate inputs for other
 functions requiring Laplacian matrices. Certainly, this function is never meant to be
 publicly exposed on its own.
+
+# References
+
+- [VL20](@cite): J. J. Veerman and R. Lyons. *A Primer on Laplacian Dynamics in Directed
+    Graphs*. *Nonlinear Phenomena in Complex Systems* **23**, 196–206 (2020).
+    https://doi.org/10.33581/1561-4085-2020-23-2-196-206.
 """
 function _assert_matrix_is_undirected_laplacian(L::AbstractMatrix{<:Integer})
     if !issymmetric(L)
@@ -266,11 +283,26 @@ type of `A` when `eltype(A)` is not an `AbstractFloat`.
 `LinearAlgebra.rank`'s default `rtol` of `min(m,n) * ϵ` for computing the rank of an
 ``m×n`` matrix may result in overestimating rank when ``|m - n| ≫ 0``, since condition
 number (which determines how numerically stable SVD and QRD are) grows with both dimensions
-[CD05; p. 603](@cite). Given that we often deal with short-and-fat matrices in this package
+[CD05; p. 603]. Given that we often deal with short-and-fat matrices in this package
 (particularly when processing all ``\\{-1, 0, 1\\}``-eigenvectors of a Laplacian matrix), we
 turn instead to the same relative tolerance used by NumPy's and MATLAB's rank
-functions—`max(m,n) * ϵ` [Num25, MAT25](@cite). (Indeed, this is a widely adopted standard
-across the field of numerical analysis [PTVF07; p. 795](@cite).)
+functions—`max(m,n) * ϵ` [Num25, MAT25]. (Indeed, this is a widely adopted standard across
+the field of numerical analysis [PTVF07; p. 795].)
+
+# References
+
+- [CD05](@cite): Z. Chen and J. Dongarra. *Condition Numbers of Gaussian Random Matrices*.
+    *SIAM Journal on Matrix Analysis and Applications* **27**, 603–20 (2005).
+    https://doi.org/10.1137/040616413.
+- [MAT25](@cite): MATLAB Developers, *rank*. MATLAB reference documentation – R2025a (2025).
+    Accessed: 2025-05-29. https://www.mathworks.com/help/matlab/ref/rank.html.
+- [Num25](@cite): NumPy Developers, *numpy.linalg.matrix_rank*. NumPy reference
+    documentation – v2.2 (2025). Accessed: 2025-05-22.
+    https://numpy.org/doc/stable/reference/generated/numpy.linalg.matrix_rank.html.
+- [PTVF07](@cite): W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery.
+    *Numerical Recipes: The Art of Scientific Computing*. 3rd Edition (Cambridge University
+    Press, Cambridge, UK, 2007). ISBN: 978-0-521-88068-8.
+    https://dl.acm.org/doi/10.5555/1403886.
 """
 function _rank_rtol(A::AbstractMatrix{T}) where {T}
     return maximum(size(A)) * eps(LinearAlgebra.eigtype(T))