add_edge Module Subroutine

  module subroutine add_edge( &
       this, edge, index, weight, feature, directed, id, update_adjacency &
  )
    !! Add an edge to the graph.
    implicit none

    ! Arguments
    class(graph_type), intent(inout) :: this
    !! Parent. Instance of the graph structure.
    type(edge_type), intent(in), optional :: edge
    !! Edge to be added.
    integer, dimension(2), intent(in), optional :: index
    !! Vertex indices of the edge.
    real(real32), intent(in), optional :: weight
    !! Weight of the edge.
    real(real32), dimension(:), intent(in), optional :: feature
    !! Feature vector of the edge.
    logical, intent(in), optional :: directed
    !! Boolean whether the edge is directed. Default is False.
    integer, intent(in), optional :: id
    !! Identifier of the edge.
    logical, intent(in), optional :: update_adjacency
    !! Boolean whether to update the adjacency matrix. Default is True.

    ! Local variables
    type(edge_type) :: edge_
    !! Initialised edge.
    real(real32) :: weight_
    !! Weight of the edge.
    logical :: directed_
    real(real32), dimension(:,:), allocatable :: edge_features
    !! Feature vectors of the edges.
    logical :: update_adjacency_
    !! Boolean whether to update the adjacency matrix.


    directed_ = .false.
    
    ! Validate input: either edge object OR individual parameters, not both
    if(present(edge) .and. (present(index) .or. present(weight) .or. &
                             present(directed) .or. present(feature)))then
       call stop_program('Specify either edge object or parameters, not both')
    elseif(.not. present(edge) .and. .not. present(index) .and. &
             .not. present(weight) .and. .not. present(directed) .and. &
             .not. present(feature))then
       call stop_program('Must specify either edge object or parameters')
    end if

    ! Construct edge from provided inputs
    if(present(edge))then
       edge_ = edge
    else
       if(.not. present(index))then
          call stop_program('Index must be specified when creating edge from parameters')
       end if
       
       weight_ = 0._real32
       if(present(weight)) weight_ = weight
       
       if(present(directed))then
          directed_ = directed
       else
          directed_ = any(index .lt. 0)
       end if
       
       if(directed_ .and. (directed_ .neqv. this%directed))then
          call stop_program('Edge direction must match graph direction')
       end if
       
       if(present(feature))then
          edge_ = edge_type_init(index, weight_, feature, directed_)
       else
          edge_ = edge_type_init(index, weight_, directed=directed_)
       end if
    end if

    ! Handle edge feature allocation and validation
    if(.not. allocated(edge_%feature))then
       allocate(edge_%feature(this%num_edge_features), source=0.0_real32)
    elseif(this%num_edge_features .eq. 0)then
       this%num_edge_features = size(edge_%feature)
    elseif(size(edge_%feature) .ne. this%num_edge_features)then
       call stop_program('Edge feature size does not match graph')
    end if

    if(present(id)) edge_%id = id

    if(.not. allocated(this%edge)) allocate(this%edge(0))
    this%num_edges = this%num_edges + 1
    
    ! Add edge to appropriate storage
    if(this%is_sparse)then
       allocate(edge_features(this%num_edge_features, this%num_edges))
       edge_features(:, 1:this%num_edges-1) = this%edge_features
       edge_features(:, this%num_edges) = edge_%feature
       if(allocated(this%edge_features)) deallocate(this%edge_features)
       call move_alloc(edge_features, this%edge_features)
       deallocate(edge_%feature)
       this%edge = [this%edge, edge_]
    else
       this%edge = [this%edge, edge_]
    end if

    ! Update adjacency matrix and degrees if requested
    update_adjacency_ = .true.
    if(present(update_adjacency)) update_adjacency_ = update_adjacency
    if(update_adjacency_)then
       call this%generate_adjacency()
       if(edge_%index(1) .gt. 0)then
          this%vertex(edge_%index(1))%degree = this%vertex(edge_%index(1))%degree + 1
       end if
       if(.not. directed_)then
          this%vertex(abs(edge_%index(2)))%degree = &
               this%vertex(abs(edge_%index(2)))%degree + 1
       end if
    end if

  end subroutine add_edge