MODULE Network_Spanning_Trees
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Sorting_Ranking
      USE Graph_Algorithms
      USE Network_Data_Structures
      IMPLICIT NONE
      PUBLIC::InitializeSpanningTree,CreateSpanningTree,UpdateSpanningTree,DestroySpanningTree
      PRIVATE
      
      PUBLIC::any_tree,random_tree,min_cost_tree,max_cost_tree
      INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE,TARGET,PUBLIC::tree_nodes_parents,&
      tree_nodes_cardinalities,tree_nodes_ordering
      REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE,TARGET,PUBLIC::tree_path_labels
      INTEGER(KIND=i_byte),DIMENSION(:),ALLOCATABLE,TARGET,PUBLIC::&
      tree_nodes_orientations
      LOGICAL(KIND=l_wp),DIMENSION(:),ALLOCATABLE,TARGET,PUBLIC::tree_arcs_mask
      
      CONTAINS
      
      
      SUBROUTINE InitializeSpanningTree()
      
      INTEGER::alloc_status
      
      ALLOCATE(tree_nodes_parents(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_nodes+n_special_nodes+1,mold=1_i_wp,&
      caller="InitializeSpanningTree",alloc_status=alloc_status)
      ALLOCATE(tree_nodes_orientations(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_nodes+n_special_nodes+1,mold=1_i_byte,&
      caller="InitializeSpanningTree",alloc_status=alloc_status)
      ALLOCATE(tree_nodes_cardinalities(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_nodes+n_special_nodes+1,mold=1_i_wp,&
      caller="InitializeSpanningTree",alloc_status=alloc_status)
      ALLOCATE(tree_nodes_ordering(0:n_nodes+n_special_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_nodes+n_special_nodes+1,mold=1_i_wp,&
      caller="InitializeSpanningTree",alloc_status=alloc_status)
      ALLOCATE(tree_path_labels(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_nodes+n_special_nodes+1,mold=1.0_r_wp,&
      caller="InitializeSpanningTree",alloc_status=alloc_status)
      ALLOCATE(tree_arcs_mask(-n_special_arcs:n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=.TRUE._l_wp,&
      caller="InitializeSpanningTree",alloc_status=alloc_status)
      
      tree_nodes_orientations=no_parent
      
      END SUBROUTINE InitializeSpanningTree
      
      
      
      SUBROUTINE CreateSpanningTree(tree_type,arcs_weights,&
      total_weight,weights_distribution,tree_timer,thread_timer)
      INTEGER(KIND=i_byte),INTENT(IN)::tree_type
      REAL(KIND=r_wp),DIMENSION(:),INTENT(IN),OPTIONAL::arcs_weights
      REAL(KIND=r_wp),INTENT(OUT),OPTIONAL::total_weight
      CHARACTER(LEN=*),INTENT(IN),OPTIONAL::weights_distribution
      INTEGER,INTENT(IN),OPTIONAL::tree_timer,thread_timer
      
      INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::arcs_ordering
      INTEGER(KIND=i_wp)::arc,node,tree_arc
      INTEGER::alloc_status
      LOGICAL::order_arcs
      CHARACTER::distribution
      
      distribution='R'
      IF(PRESENT(weights_distribution))distribution=weights_distribution(1:1)
      
      order_arcs=(PRESENT(arcs_weights).AND.(tree_type/=any_tree))&
      .OR.(tree_type==random_tree)
      
      OrderArcs:IF(order_arcs)THEN
      
      ALLOCATE(arcs_ordering(-n_special_arcs:n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=1_i_wp,&
      caller="CreateSpanningTree",alloc_status=alloc_status)
      
      SELECT CASE(tree_type)
      CASE(min_cost_tree,max_cost_tree)
      SinglePrecision:IF(r_wp==r_sp)
      SELECT CASE(distribution)
      CASE('U','u','N','n','G','g')
      CALL HashRank(array=arcs_weights,permutation=arcs_ordering,&
      distribution=distribution,partially_ranked=.FALSE.)
      CASE DEFAULT
      CALL RadixRank(array=arcs_weights,permutation=arcs_ordering)
      END SELECT
      ELSE
      CALL QuickRank(array=arcs_weights,permutation=arcs_ordering,&
      pivot_selection='R',partially_ranked=.FALSE.)
      END IF SinglePrecision
      CASE DEFAULT
      CALL DisorderPermutation(disorder=1.0,permutation=arcs_ordering)
      END SELECT
      arcs_ordering=arcs_ordering-n_special_arcs-1
      
      IF(PRESENT(tree_timer))THEN
      CALL ResetTimer(tree_timer)
      CALL StartTimer(tree_timer)
      END IF
      IF(tree_type==max_cost_tree)THEN
      CALL BuildSpanningTree(arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      arcs_list=arcs_ordering(n_arcs:(-n_special_arcs):(-1)),heads_tails=heads_tails,&
      orientations=tree_nodes_orientations,parents=tree_nodes_parents)
      ELSE
      CALL BuildSpanningTree(arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      arcs_list=arcs_ordering,heads_tails=heads_tails,&
      orientations=tree_nodes_orientations,parents=tree_nodes_parents)
      END IF
      IF(PRESENT(tree_timer))THEN
      CALL StopTimer(tree_timer)
      END IF
      
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_ordering),KIND=i_wp),&
      mold=1_i_wp,caller="CreateSpanningTree")
      DEALLOCATE(arcs_ordering)
      
      ELSE
      
      CALL BuildSpanningTree(arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      arcs_list=(/(arc,arc=-n_special_arcs,n_arcs)/),heads_tails=heads_tails,&
      orientations=tree_nodes_orientations,parents=tree_nodes_parents)
      
      END IF OrderArcs
      
      IF(PRESENT(thread_timer))THEN
      CALL ResetTimer(thread_timer)
      CALL StartTimer(thread_timer)
      END IF
      CALL BuildTreeStructures(arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      heads_tails=heads_tails,orientations=tree_nodes_orientations,&
      parents=tree_nodes_parents,cardinalities=tree_nodes_cardinalities,&
      level_ordering=tree_nodes_ordering,tree_mask=tree_arcs_mask)
      IF(PRESENT(thread_timer))THEN
      CALL StopTimer(thread_timer)
      END IF
      
      IF(PRESENT(total_weight).AND.PRESENT(arcs_weights))THEN
      
      total_weight=0.0_r_wp
      DO node=-n_special_nodes,n_nodes
      IF(tree_nodes_orientations(node)/=no_parent)THEN
      tree_arc=tree_nodes_parents(node)
      total_weight=total_weight+arcs_weights(tree_arc+n_special_arcs+1)
      END IF
      END DO
      
      
      END IF
      
      END SUBROUTINE CreateSpanningTree
      
      
      
      SUBROUTINE UpdateSpanningTree(tree_type,arcs_weights,&
      total_weight,n_examined_arcs,n_pivots,tree_timer,thread_timer)
      INTEGER(KIND=i_byte),INTENT(IN)::tree_type
      REAL(KIND=r_wp),DIMENSION(:),INTENT(IN),OPTIONAL::arcs_weights
      REAL(KIND=r_wp),INTENT(OUT),OPTIONAL::total_weight
      INTEGER(KIND=i_wp),INTENT(OUT),OPTIONAL::n_examined_arcs,n_pivots
      INTEGER,INTENT(IN),OPTIONAL::tree_timer,thread_timer
      
      INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::candidate_arcs,candidate_order
      INTEGER(KIND=i_wp)::arc,arc_index,tree_arc,node,head,tail,n_candidate_arcs
      INTEGER::alloc_status
      LOGICAL::tree_arc_test,path_labels_test
      REAL::percent_pivoted
      
      ALLOCATE(candidate_arcs(0:n_special_arcs+n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=1_i_wp,&
      caller="UpdateSpanningTree",alloc_status=alloc_status)
      
      CALL ResetTimer(tree_timer)
      CALL StartTimer(tree_timer)
      
      CALL CalculatePathLabels(arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      heads_tails=heads_tails,level_ordering=tree_nodes_ordering,&
      orientations=tree_nodes_orientations,parents=tree_nodes_parents,&
      path_labels=tree_path_labels,arcs_weights=arcs_weights,&
      tree_type=tree_type)
      
      n_candidate_arcs=0
      DO arc=1,n_arcs+n_special_arcs+1
      arc_index=arc-n_special_arcs-1
      head=heads_tails(1,arc_index)
      tail=heads_tails(2,arc_index)
      
      tree_arc_test=.NOT.(&
      ((tree_nodes_orientations(head)==tail_is_parent).AND.(tree_nodes_parents(head)==arc_index))&
      .OR.&
      ((tree_nodes_orientations(tail)==head_is_parent).AND.(tree_nodes_parents(tail)==arc_index)))
      IF(tree_type==min_cost_tree)THEN
      path_labels_test=(arcs_weights(arc)<=MAX(tree_path_labels(head),tree_path_labels(tail)))
      ELSE
      path_labels_test=(arcs_weights(arc)>=MIN(tree_path_labels(head),tree_path_labels(tail)))
      END IF
      IF(tree_arc_test.AND.path_labels_test.AND.(arc_index/=0))THEN
      n_candidate_arcs=n_candidate_arcs+1
      candidate_arcs(n_candidate_arcs)=arc
      END IF
      END DO
      CALL ReBuildSpanningTree(arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      arcs_list=candidate_arcs(1:n_candidate_arcs),&
      heads_tails=heads_tails,arcs_weights=arcs_weights,&
      orientations=tree_nodes_orientations,parents=tree_nodes_parents,&
      cardinalities=tree_nodes_cardinalities,tree_type=tree_type,&
      percent_pivoted=percent_pivoted)
      IF(PRESENT(n_examined_arcs))n_examined_arcs=n_candidate_arcs
      IF(PRESENT(n_pivots))n_pivots=INT(percent_pivoted*REAL(n_candidate_arcs),i_wp)
      
      CALL StopTimer(tree_timer)
      
      CALL ResetTimer(thread_timer)
      CALL StartTimer(thread_timer)
      CALL BuildTreeStructures(arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      heads_tails=heads_tails,orientations=tree_nodes_orientations,parents=tree_nodes_parents,&
      cardinalities=tree_nodes_cardinalities,level_ordering=tree_nodes_ordering,&
      rebuild_cardinalities=.FALSE.)
      CALL StopTimer(thread_timer)
      
      CALL RecordAllocation(n_elements=-INT(SIZE(candidate_arcs),KIND=i_wp),&
      mold=1_i_wp,caller="UpdateSpanningTree")
      DEALLOCATE(candidate_arcs)
      
      IF(PRESENT(total_weight))THEN
      
      total_weight=0.0_r_wp
      DO node=-n_special_nodes,n_nodes
      IF(tree_nodes_orientations(node)/=no_parent)THEN
      tree_arc=tree_nodes_parents(node)
      total_weight=total_weight+arcs_weights(tree_arc+n_special_arcs+1)
      END IF
      END DO
      
      
      END IF
      
      END SUBROUTINE UpdateSpanningTree
      
      
      
      SUBROUTINE DestroySpanningTree()
      
      INTEGER::alloc_status
      
      CALL RecordAllocation(n_elements=-INT(SIZE(tree_nodes_parents),KIND=i_wp),&
      mold=1_i_wp,caller="DestroySpanningTree")
      DEALLOCATE(tree_nodes_parents)
      CALL RecordAllocation(n_elements=-INT(SIZE(tree_nodes_orientations),KIND=i_wp),&
      mold=1_i_byte,caller="DestroySpanningTree")
      DEALLOCATE(tree_nodes_orientations)
      CALL RecordAllocation(n_elements=-INT(SIZE(tree_nodes_cardinalities),KIND=i_wp),&
      mold=1_i_wp,caller="DestroySpanningTree")
      DEALLOCATE(tree_nodes_cardinalities)
      CALL RecordAllocation(n_elements=-INT(SIZE(tree_nodes_ordering),KIND=i_wp),&
      mold=1_i_wp,caller="DestroySpanningTree")
      DEALLOCATE(tree_nodes_ordering)
      CALL RecordAllocation(n_elements=-INT(SIZE(tree_path_labels),KIND=i_wp),&
      mold=1.0_r_wp,caller="DestroySpanningTree")
      DEALLOCATE(tree_path_labels)
      CALL RecordAllocation(n_elements=-INT(SIZE(tree_arcs_mask),KIND=i_wp),&
      mold=.TRUE._l_wp,caller="DestroySpanningTree")
      DEALLOCATE(tree_arcs_mask)
      
      END SUBROUTINE DestroySpanningTree
      
      
      
      END MODULE Network_Spanning_Trees