! Copyright (c) 2022-2026 Jason Christopherson
! SPDX-License-Identifier: MIT
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The Software is provided "as is", without warranty of any kind, express or
! implied, including but not limited to the warranties of merchantability,
! fitness for a particular purpose and noninfringement.
module dynamics_graph
    !! Provides a lightweight, undirected multigraph container along with the
    !! traversal utilities required to analyze the topology of closed-loop
    !! mechanisms.  Vertices are represented by integer indices, and edges are
    !! stored in a simple allocatable array.
    use iso_fortran_env
    use dynamics_error_handling
    implicit none
    private
    public :: graph_edge
    public :: graph_path
    public :: graph_loop
    public :: spanning_tree
    public :: graph

    type graph_edge
        !! Defines an undirected edge connecting two vertices.
        integer(int32), public :: vertex_1 = 0
            !! The index of the first vertex.
        integer(int32), public :: vertex_2 = 0
            !! The index of the second vertex.
    end type

    type graph_path
        !! Describes a path through a graph.
        integer(int32), allocatable, public, dimension(:) :: vertices
            !! An N+1 element array containing the vertices visited along the
            !! path, in order.
        integer(int32), allocatable, public, dimension(:) :: edges
            !! An N-element array containing the edges traversed along the path,
            !! in order.
        logical, allocatable, public, dimension(:) :: forward
            !! An N-element array that is true if the corresponding edge was
            !! traversed from vertex_1 to vertex_2, and false if the edge was
            !! traversed from vertex_2 to vertex_1.
    end type

    type graph_loop
        !! Describes an independent loop within the graph.  The loop is defined
        !! by an edge that is not part of the spanning tree; the loop itself is
        !! closed by the tree paths leading to each of the edge's vertices.
        integer(int32), public :: cut_edge = 0
            !! The index of the edge closing the loop.
        integer(int32), public :: vertex_1 = 0
            !! The first vertex of the cut edge.
        integer(int32), public :: vertex_2 = 0
            !! The second vertex of the cut edge.
    end type

    type spanning_tree
        !! Describes a spanning tree of a graph as generated by a breadth-first
        !! traversal.
        integer(int32), public :: root = 0
            !! The vertex from which the traversal was started.
        integer(int32), allocatable, public, dimension(:) :: parent_vertex
            !! An array, one entry per vertex, containing the index of the
            !! parent vertex.  The root vertex and any unreachable vertices are
            !! assigned a value of zero.
        integer(int32), allocatable, public, dimension(:) :: parent_edge
            !! An array, one entry per vertex, containing the index of the edge
            !! connecting the vertex to its parent.  The root vertex and any
            !! unreachable vertices are assigned a value of zero.
        logical, allocatable, public, dimension(:) :: parent_edge_forward
            !! An array, one entry per vertex, that is true if the edge
            !! connecting the vertex to its parent is traversed from vertex_1 to
            !! vertex_2 when moving from the parent to the vertex.
        integer(int32), allocatable, public, dimension(:) :: depth
            !! An array, one entry per vertex, containing the number of edges
            !! between the vertex and the root.  Unreachable vertices are
            !! assigned a value of -1.
        logical, allocatable, public, dimension(:) :: edge_in_tree
            !! An array, one entry per edge, that is true if the edge is part of
            !! the spanning tree.
        integer(int32), allocatable, public, dimension(:) :: visit_order
            !! An array containing the reachable vertices in the order in which
            !! they were visited.
    contains
        procedure, public :: get_path => st_get_path
        procedure, public :: get_cut_edges => st_get_cut_edges
        procedure, public :: is_connected => st_is_connected
    end type

    type graph
        !! Defines an undirected multigraph.  Multiple edges between the same
        !! pair of vertices, as well as self-loops, are permitted.
        integer(int32), private :: m_vertices = 0
        integer(int32), private :: m_count = 0
        type(graph_edge), allocatable, private, dimension(:) :: m_edges
    contains
        procedure, public :: initialize => gr_initialize
        procedure, public :: get_vertex_count => gr_vertex_count
        procedure, public :: get_edge_count => gr_edge_count
        procedure, public :: get_edge => gr_get_edge
        procedure, public :: add_edge => gr_add_edge
        procedure, public :: get_adjacent_edges => gr_adjacent_edges
        procedure, public :: build_spanning_tree => gr_spanning_tree
        procedure, public :: find_independent_loops => gr_find_loops
        procedure, public :: get_independent_loop_count => gr_loop_count
    end type

contains
! ******************************************************************************
! GRAPH MEMBERS
! ------------------------------------------------------------------------------
    subroutine gr_initialize(this, nvertices, ncapacity)
        !! Initializes the graph with the requested number of vertices and no
        !! edges.
        class(graph), intent(inout) :: this
            !! The graph object.
        integer(int32), intent(in) :: nvertices
            !! The number of vertices in the graph.  This value must be at least
            !! one.
        integer(int32), intent(in), optional :: ncapacity
            !! An optional estimate of the number of edges the graph will
            !! contain.  This value is used only to size the initial storage.

        ! Local Variables
        integer(int32) :: n, flag

        ! Input Checking
        if (nvertices < 1) error stop DYN_INVALID_INPUT_ERROR

        ! Process
        n = 10
        if (present(ncapacity)) n = max(ncapacity, 1)
        this%m_vertices = nvertices
        this%m_count = 0
        if (allocated(this%m_edges)) deallocate(this%m_edges)
        allocate(this%m_edges(n), stat = flag)
        if (flag /= 0) error stop DYN_MEMORY_ERROR
    end subroutine

! ------------------------------------------------------------------------------
    pure function gr_vertex_count(this) result(rst)
        !! Gets the number of vertices in the graph.
        class(graph), intent(in) :: this
            !! The graph object.
        integer(int32) :: rst
            !! The vertex count.

        rst = this%m_vertices
    end function

! ------------------------------------------------------------------------------
    pure function gr_edge_count(this) result(rst)
        !! Gets the number of edges in the graph.
        class(graph), intent(in) :: this
            !! The graph object.
        integer(int32) :: rst
            !! The edge count.

        rst = this%m_count
    end function

! ------------------------------------------------------------------------------
    function gr_get_edge(this, i) result(rst)
        !! Gets the requested edge.
        class(graph), intent(in) :: this
            !! The graph object.
        integer(int32), intent(in) :: i
            !! The index of the edge to retrieve (1 = first edge).
        type(graph_edge) :: rst
            !! The requested edge.

        if (i < 1 .or. i > this%m_count) error stop DYN_INDEX_OUT_OF_RANGE
        rst = this%m_edges(i)
    end function

! ------------------------------------------------------------------------------
    subroutine gr_add_edge(this, v1, v2)
        !! Adds an edge to the graph.
        class(graph), intent(inout) :: this
            !! The graph object.
        integer(int32), intent(in) :: v1
            !! The index of the first vertex.
        integer(int32), intent(in) :: v2
            !! The index of the second vertex.

        ! Local Variables
        integer(int32) :: flag
        type(graph_edge), allocatable, dimension(:) :: temp

        ! Input Checking
        if (this%m_vertices < 1) error stop DYN_INVALID_INPUT_ERROR
        if (v1 < 1 .or. v1 > this%m_vertices) error stop DYN_INDEX_OUT_OF_RANGE
        if (v2 < 1 .or. v2 > this%m_vertices) error stop DYN_INDEX_OUT_OF_RANGE

        ! Ensure sufficient capacity exists
        if (this%m_count >= size(this%m_edges)) then
            allocate(temp(2 * size(this%m_edges)), stat = flag)
            if (flag /= 0) error stop DYN_MEMORY_ERROR
            temp(1:this%m_count) = this%m_edges(1:this%m_count)
            call move_alloc(temp, this%m_edges)
        end if

        ! Store the edge
        this%m_count = this%m_count + 1
        this%m_edges(this%m_count)%vertex_1 = v1
        this%m_edges(this%m_count)%vertex_2 = v2
    end subroutine

! ------------------------------------------------------------------------------
    function gr_adjacent_edges(this, v) result(rst)
        !! Gets the indices of all edges connected to the requested vertex.
        class(graph), intent(in) :: this
            !! The graph object.
        integer(int32), intent(in) :: v
            !! The index of the vertex of interest.
        integer(int32), allocatable, dimension(:) :: rst
            !! An array containing the indices of the connected edges.

        ! Local Variables
        integer(int32) :: i, n
        integer(int32), allocatable, dimension(:) :: buffer

        ! Input Checking
        if (v < 1 .or. v > this%m_vertices) error stop DYN_INDEX_OUT_OF_RANGE

        ! Process
        allocate(buffer(this%m_count))
        n = 0
        do i = 1, this%m_count
            if (this%m_edges(i)%vertex_1 == v .or. &
                this%m_edges(i)%vertex_2 == v) &
            then
                n = n + 1
                buffer(n) = i
            end if
        end do
        rst = buffer(1:n)
    end function

! ------------------------------------------------------------------------------
    function gr_spanning_tree(this, root) result(rst)
        !! Constructs a spanning tree of the graph by means of a breadth-first
        !! traversal.
        class(graph), intent(in) :: this
            !! The graph object.
        integer(int32), intent(in), optional :: root
            !! The index of the vertex from which to start the traversal.  If
            !! not supplied, the first vertex is used.
        type(spanning_tree) :: rst
            !! The resulting spanning_tree object.

        ! Local Variables
        integer(int32) :: i, k, nv, ne, iv, ihead, itail, other
        integer(int32), allocatable, dimension(:) :: queue

        ! Initialization
        nv = this%m_vertices
        ne = this%m_count
        if (nv < 1) error stop DYN_INVALID_INPUT_ERROR
        rst%root = 1
        if (present(root)) rst%root = root
        if (rst%root < 1 .or. rst%root > nv) error stop DYN_INDEX_OUT_OF_RANGE

        allocate(rst%parent_vertex(nv), source = 0)
        allocate(rst%parent_edge(nv), source = 0)
        allocate(rst%parent_edge_forward(nv), source = .true.)
        allocate(rst%depth(nv), source = -1)
        allocate(rst%edge_in_tree(ne), source = .false.)
        allocate(queue(nv))

        ! Traverse the graph
        ihead = 1
        itail = 1
        queue(1) = rst%root
        rst%depth(rst%root) = 0
        do while (ihead <= itail)
            iv = queue(ihead)
            ihead = ihead + 1
            do i = 1, ne
                if (this%m_edges(i)%vertex_1 == iv) then
                    other = this%m_edges(i)%vertex_2
                else if (this%m_edges(i)%vertex_2 == iv) then
                    other = this%m_edges(i)%vertex_1
                else
                    cycle
                end if
                if (rst%depth(other) >= 0) cycle
                rst%depth(other) = rst%depth(iv) + 1
                rst%parent_vertex(other) = iv
                rst%parent_edge(other) = i
                rst%parent_edge_forward(other) = &
                    this%m_edges(i)%vertex_1 == iv
                rst%edge_in_tree(i) = .true.
                itail = itail + 1
                queue(itail) = other
            end do
        end do

        ! Store the visitation order
        allocate(rst%visit_order(itail))
        do k = 1, itail
            rst%visit_order(k) = queue(k)
        end do
    end function

! ------------------------------------------------------------------------------
    function gr_find_loops(this, tree) result(rst)
        !! Determines the set of independent loops within the graph.  The number
        !! of loops is the cyclomatic number of the graph.
        class(graph), intent(in) :: this
            !! The graph object.
        type(spanning_tree), intent(in) :: tree
            !! A spanning tree of this graph.
        type(graph_loop), allocatable, dimension(:) :: rst
            !! An array of the independent loops.

        ! Local Variables
        integer(int32) :: i, n
        integer(int32), allocatable, dimension(:) :: cuts

        ! Process
        cuts = tree%get_cut_edges()
        n = size(cuts)
        allocate(rst(n))
        do i = 1, n
            rst(i)%cut_edge = cuts(i)
            rst(i)%vertex_1 = this%m_edges(cuts(i))%vertex_1
            rst(i)%vertex_2 = this%m_edges(cuts(i))%vertex_2
        end do
    end function

! ------------------------------------------------------------------------------
    pure function gr_loop_count(this) result(rst)
        !! Gets the number of independent loops in the graph assuming the graph
        !! is connected.  The value is the cyclomatic number \(E - V + 1\).
        class(graph), intent(in) :: this
            !! The graph object.
        integer(int32) :: rst
            !! The number of independent loops.

        rst = this%m_count - this%m_vertices + 1
    end function

! ******************************************************************************
! SPANNING_TREE MEMBERS
! ------------------------------------------------------------------------------
    function st_get_path(this, v) result(rst)
        !! Gets the path from the root of the tree to the requested vertex.
        class(spanning_tree), intent(in) :: this
            !! The spanning_tree object.
        integer(int32), intent(in) :: v
            !! The index of the destination vertex.
        type(graph_path) :: rst
            !! The resulting path.  The path contains no edges if the requested
            !! vertex is the root vertex.

        ! Local Variables
        integer(int32) :: i, n, iv, ie

        ! Input Checking
        if (.not.allocated(this%depth)) error stop DYN_INVALID_INPUT_ERROR
        if (v < 1 .or. v > size(this%depth)) error stop DYN_INDEX_OUT_OF_RANGE
        if (this%depth(v) < 0) error stop DYN_INVALID_INPUT_ERROR

        ! Walk back to the root and store the results in forward order
        n = this%depth(v)
        allocate(rst%vertices(n + 1), rst%edges(n), rst%forward(n))
        iv = v
        rst%vertices(n + 1) = v
        do i = n, 1, -1
            ie = this%parent_edge(iv)
            rst%edges(i) = ie
            rst%vertices(i) = this%parent_vertex(iv)
            rst%forward(i) = this%parent_edge_forward(iv)
            iv = this%parent_vertex(iv)
        end do
    end function

! ------------------------------------------------------------------------------
    function st_get_cut_edges(this) result(rst)
        !! Gets the indices of all edges that are not part of the spanning tree.
        class(spanning_tree), intent(in) :: this
            !! The spanning_tree object.
        integer(int32), allocatable, dimension(:) :: rst
            !! An array containing the indices of the cut edges.

        ! Local Variables
        integer(int32) :: i, n
        integer(int32), allocatable, dimension(:) :: buffer

        if (.not.allocated(this%edge_in_tree)) error stop DYN_INVALID_INPUT_ERROR
        allocate(buffer(size(this%edge_in_tree)))
        n = 0
        do i = 1, size(this%edge_in_tree)
            if (.not.this%edge_in_tree(i)) then
                n = n + 1
                buffer(n) = i
            end if
        end do
        rst = buffer(1:n)
    end function

! ------------------------------------------------------------------------------
    pure function st_is_connected(this) result(rst)
        !! Determines if every vertex was reached by the traversal used to
        !! construct the tree.
        class(spanning_tree), intent(in) :: this
            !! The spanning_tree object.
        logical :: rst
            !! True if the graph is connected; else, false.

        rst = .false.
        if (.not.allocated(this%depth)) return
        rst = all(this%depth >= 0)
    end function

! ------------------------------------------------------------------------------
end module