MODULE Network_Geometry
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Random_Numbers
      USE Sorting_Ranking
      USE Space_Filling_Curves
      USE Network_Data_Structures
      USE Lattice_Geometry
      USE Graph_Algorithms
      USE Support_Trees
      IMPLICIT NONE
      PUBLIC::InitializeNetwork,DestroyNetwork,CreateNetwork,AssignSuppliesDemands
      PRIVATE
      
      CHARACTER,DIMENSION(n_dim),PUBLIC,SAVE::&
      source_arcs_cost='L',sink_arcs_cost='L',&
      sources_status='R',sinks_status='R'
      CHARACTER,PUBLIC,SAVE::plate_arcs_cost='L'
      REAL(KIND=r_wp),PUBLIC,DIMENSION(n_dim),SAVE::&
      sources_flow_capacity=0.0_r_wp,sinks_flow_capacity=0.0_r_wp
      REAL(KIND=r_wp),PUBLIC,SAVE::injection_flow_capacity=1.0_r_wp,outflow_flow_capacity=-1.0_r_wp
      CHARACTER,PUBLIC,SAVE::injection_flow_distribution='F',outflow_flow_distribution='F'
      REAL(KIND=r_wp),PUBLIC,DIMENSION(2)::injection_flow_parameters=(/1.0_r_wp,0.0_r_wp/),outflow_flow_parameters=(/1.0_r_wp,0.0_&
     &r_wp/)
      CHARACTER,PUBLIC,SAVE::SFC='H'
      
      
      INTEGER(KIND=i_byte),DIMENSION(n_dim),SAVE::source_arcs_cost_,sink_arcs_cost_,&
      sources_status_,sinks_status_
      INTEGER(KIND=i_byte),SAVE::plate_arcs_cost_
      INTEGER,SAVE::n_random_matchings=0
      
      NAMELIST/ProblemOptions/sources_status,sinks_status,source_arcs_cost,sink_arcs_cost,&
      plate_arcs_cost,sources_flow_capacity,sinks_flow_capacity,&
      injection_flow_capacity,outflow_flow_capacity,&
      injection_flow_distribution,outflow_flow_distribution,&
      injection_flow_parameters,outflow_flow_parameters,SFC,n_random_matchings
      
      CONTAINS
      
      
      SUBROUTINE InitializeNetwork()
      
      REWIND(UNIT=program_options_unit)
      READ(UNIT=program_options_unit,NML=ProblemOptions,IOSTAT=error_status)
      IF(error_status/=0)&
      CALL CriticalError(message="NAMELIST ProblemOptions was not read successfully"//&
      " from file "//TRIM(options_file),caller="InitializeNetwork")
      WRITE(UNIT=message_log_unit,NML=ProblemOptions)
      
      END SUBROUTINE InitializeNetwork
      
      
      
      SUBROUTINE DestroyNetwork()
      IMPLICIT NONE
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_mask),KIND=i_wp),&
      mold=.TRUE._l_wp,caller="DestroyNetwork")
      DEALLOCATE(nodes_mask)
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_mask),KIND=i_wp),&
      mold=.TRUE._l_wp,caller="DestroyNetwork")
      DEALLOCATE(arcs_mask)
      
      CALL RecordAllocation(n_elements=-INT(SIZE(special_nodes_status),KIND=i_wp),&
      mold=1_i_byte,caller="DestroyNetwork")
      DEALLOCATE(special_nodes_status)
      CALL RecordAllocation(n_elements=-INT(SIZE(special_arcs_status),KIND=i_wp),&
      mold=1_i_byte,caller="DestroyNetwork")
      DEALLOCATE(special_arcs_status)
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_potentials),KIND=i_wp),&
      mold=1.0_r_wp,caller="DestroyNetwork")
      DEALLOCATE(nodes_potentials)
      CALL RecordAllocation(n_elements=-INT(SIZE(supplies_demands),KIND=i_wp),&
      mold=1.0_r_wp,caller="DestroyNetwork")
      DEALLOCATE(supplies_demands)
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_flows),KIND=i_wp),&
      mold=1.0_r_wp,caller="DestroyNetwork")
      DEALLOCATE(arcs_flows)
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_voltages),KIND=i_wp),&
      mold=1.0_r_wp,caller="DestroyNetwork")
      DEALLOCATE(arcs_voltages)
      
      END SUBROUTINE DestroyNetwork
      
      
      
      SUBROUTINE CreateNetwork(reordering_timer)
      INTEGER,INTENT(IN),OPTIONAL::reordering_timer
      
      INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::nodes_reordering,nodes_renumbering,&
      arcs_reordering
      INTEGER(KIND=i_wp)::node,arc,special_node,special_arc,source,sink
      INTEGER::dim
      
      
      DO dim=1,n_dim
      
      SELECT CASE(source_arcs_cost(dim))
      CASE('D','d')
      source_arcs_cost_(dim)=dummy_arc
      CASE('L','l')
      source_arcs_cost_(dim)=low_cost_arc
      CASE('H','h')
      source_arcs_cost_(dim)=high_cost_arc
      CASE DEFAULT
      source_arcs_cost_(dim)=regular_cost_arc
      ENDSELECT
      
      SELECT CASE(sink_arcs_cost(dim))
      CASE('D','d')
      sink_arcs_cost_(dim)=dummy_arc
      CASE('L','l')
      sink_arcs_cost_(dim)=low_cost_arc
      CASE('H','h')
      sink_arcs_cost_(dim)=high_cost_arc
      CASE DEFAULT
      sink_arcs_cost_(dim)=regular_cost_arc
      ENDSELECT
      
      SELECT CASE(sources_status(dim))
      CASE('D','d')
      sources_status_(dim)=dummy_node
      CASE('F','f')
      sources_status_(dim)=fixed_potential_node
      CASE DEFAULT
      sources_status_(dim)=regular_node
      ENDSELECT
      
      SELECT CASE(sinks_status(dim))
      CASE('D','d')
      sinks_status_(dim)=dummy_node
      CASE('F','f')
      sinks_status_(dim)=fixed_potential_node
      CASE DEFAULT
      sinks_status_(dim)=regular_node
      ENDSELECT
      END DO
      
      SELECT CASE(plate_arcs_cost)
      CASE('D','d')
      plate_arcs_cost_=dummy_arc
      CASE('L','l')
      plate_arcs_cost_=low_cost_arc
      CASE('H','h')
      plate_arcs_cost_=high_cost_arc
      CASE DEFAULT
      plate_arcs_cost_=regular_cost_arc
      ENDSELECT
      
      WRITE(message_log_unit,*)"Creating a connected lattice network with:"
      WRITE(message_log_unit,*)"n_arcs=",n_arcs," n_nodes=",n_nodes,&
      " n_special_nodes=",n_special_nodes," n_special_arcs=",n_special_arcs
      
      
      ALLOCATE(special_nodes_status(-n_special_nodes:0),STAT=error_status)
      CALL RecordAllocation(n_elements=(n_special_nodes+1),mold=1_i_byte,&
      caller="CreateNetwork",alloc_status=error_status)
      ALLOCATE(special_arcs_status(-n_special_arcs:0),STAT=error_status)
      CALL RecordAllocation(n_elements=(n_special_arcs+1),mold=1_i_byte,&
      caller="CreateNetwork",alloc_status=error_status)
      
      
      DO dim=1,n_dim
      IF(boundaries(bc_source_sink,dim))THEN
      source=sources(dim)
      special_nodes_status(source)=source_arcs_cost_(dim)
      sink=sinks(dim)
      special_nodes_status(sink)=sink_arcs_cost_(dim)
      END IF
      END DO
      
      n_regular_cost_arcs=0;
      n_high_cost_arcs=0;
      n_low_cost_arcs=0
      special_arcs_status(0)=dummy_arc
      
      DO special_arc=-n_special_arcs,-1
      
      SELECT CASE(arcs_status(special_arc))
      CASE(plate_arc)
      special_arcs_status(special_arc)=plate_arcs_cost_
      CASE(source_arc)
      source=heads_tails(1,special_arc)
      special_arcs_status(special_arc)=special_nodes_status(source)
      CASE(sink_arc)
      sink=heads_tails(2,special_arc)
      special_arcs_status(special_arc)=special_nodes_status(sink)
      CASE DEFAULT
      special_arcs_status(special_arc)=regular_cost_arc
      END SELECT
      
      SELECT CASE(special_arcs_status(special_arc))
      CASE(low_cost_arc)
      n_low_cost_arcs=n_low_cost_arcs+1
      CASE(high_cost_arc)
      n_high_cost_arcs=n_high_cost_arcs+1
      CASE(regular_cost_arc)
      n_regular_cost_arcs=n_regular_cost_arcs+1
      END SELECT
      
      END DO
      
      special_nodes_status(0)=dummy_arc
      special_nodes_status(-n_special_nodes:(-1))=regular_node
      DO dim=1,n_dim
      IF(boundaries(bc_source_sink,dim))THEN
      source=sources(dim)
      special_nodes_status(source)=sources_status_(dim)
      sink=sinks(dim)
      special_nodes_status(sink)=sinks_status_(dim)
      END IF
      END DO
      
      n_fixed_potential_nodes=0;
      n_regular_nodes=0
      DO special_node=-n_special_nodes,-1
      SELECT CASE(special_nodes_status(special_node))
      CASE(fixed_potential_node)
      n_fixed_potential_nodes=n_fixed_potential_nodes+1
      CASE(regular_node)
      n_regular_nodes=n_regular_nodes+1
      END SELECT
      END DO
      
      
      IF(PRESENT(reordering_timer))CALL StartTimer(reordering_timer)
      
      ALLOCATE(nodes_reordering(-n_special_nodes:n_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=(n_special_nodes+n_nodes+1),mold=1_i_wp,&
      caller="CreateNetwork",alloc_status=error_status)
      ALLOCATE(nodes_renumbering(-n_special_nodes:n_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=(n_special_nodes+n_nodes+1),mold=1_i_wp,&
      caller="CreateNetwork",alloc_status=error_status)
      nodes_reordering(0)=0
      
      SELECT CASE(SFC)
      CASE('h','H','m','M')
      CALL SFCOrder(points_coords=nodes_coords(:,1:n_nodes),&
      SFC_order=nodes_reordering(1:n_nodes),SFC=SFC,&
      coords_min=(/(1.0,dim=1,n_dim)/),coords_max=REAL(lengths))
      
      CASE('C','c')
      CALL ST_MatchingNodesMapping(node_offset=-1,arc_offset=-1,&
      heads_tails=heads_tails(:,1:),nodes_mapping=nodes_renumbering(1:n_nodes),&
      mapping_offset=1,heavy_matching=.TRUE.,&
      n_random_matchings=n_random_matchings)
      DO node=1,n_nodes
      nodes_reordering(nodes_renumbering(node))=node
      END DO
      CASE('R','r')
      CALL DisorderPermutation(disorder=1.0,permutation=nodes_reordering(1:n_nodes))
      
      CASE DEFAULT
      DO node=1,n_nodes
      nodes_reordering(node)=node
      END DO
      END SELECT
      
      IF(n_special_nodes>0)THEN
      CALL QuickRank(array=special_nodes_status(-n_special_nodes:(-1)),&
      permutation=nodes_reordering(-n_special_nodes:(-1)),pivot_selection='R')
      nodes_reordering(-n_special_nodes:(-1))=&
      nodes_reordering(-n_special_nodes:(-1))-(n_special_nodes+1)
      special_nodes_status(-n_special_nodes:(-1))=&
      special_nodes_status(nodes_reordering(-n_special_nodes:(-1)))
      END IF
      
      nodes_coords=nodes_coords(:,nodes_reordering)
      nodes_status=nodes_status(nodes_reordering)
      
      DO node=-n_special_nodes,0
      nodes_renumbering(nodes_reordering(node))=node
      END DO
      IF((SFC/='C').OR.(SFC/='c'))THEN
      DO node=1,n_nodes
      nodes_renumbering(nodes_reordering(node))=node
      END DO
      END IF
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_reordering),KIND=i_wp),&
      mold=1_i_wp,caller="CreateNetwork")
      DEALLOCATE(nodes_reordering)
      
      DO dim=1,n_dim
      IF(boundaries(bc_source_sink,dim))THEN
      source=sources(dim)
      sources(dim)=nodes_renumbering(source)
      sink=sinks(dim)
      sinks(dim)=nodes_renumbering(sink)
      END IF
      END DO
      
      DO arc=-n_special_arcs,n_arcs
      heads_tails(:,arc)=nodes_renumbering(heads_tails(:,arc))
      END DO
      
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_renumbering),KIND=i_wp),&
      mold=1_i_wp,caller="CreateNetwork")
      DEALLOCATE(nodes_renumbering)
      
      
      ALLOCATE(arcs_reordering(-n_special_arcs:n_arcs),STAT=error_status)
      CALL RecordAllocation(n_elements=(n_special_arcs+n_arcs+1),mold=1_i_wp,&
      caller="CreateNetwork",alloc_status=error_status)
      arcs_reordering(0)=0
      
      CALL QuickRank(array=heads_tails(1,1:n_arcs),permutation=arcs_reordering(1:n_arcs),&
      pivot_selection='U')
      CALL ShellInsertionRank(array=heads_tails(:,1:n_arcs),&
      permutation=arcs_reordering(1:n_arcs),partially_ranked=.TRUE.,method="Insertion")
      
      IF(n_special_arcs>0)THEN
      CALL QuickRank(array=special_arcs_status(-n_special_arcs:(-1)),&
      permutation=arcs_reordering(-n_special_arcs:(-1)),pivot_selection='R')
      arcs_reordering(-n_special_arcs:(-1))=&
      arcs_reordering(-n_special_arcs:(-1))-(n_special_arcs+1)
      special_arcs_status(-n_special_arcs:(-1))=&
      special_arcs_status(arcs_reordering(-n_special_arcs:(-1)))
      END IF
      
      heads_tails=heads_tails(:,arcs_reordering)
      arcs_status=arcs_status(arcs_reordering)
      
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_reordering),KIND=i_wp),&
      mold=1_i_wp,caller="CreateNetwork")
      DEALLOCATE(arcs_reordering)
      
      
      IF(PRESENT(reordering_timer))CALL StopTimer(reordering_timer)
      
      ALLOCATE(nodes_mask(-n_special_nodes:n_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1,mold=.TRUE._l_wp,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      nodes_mask=.TRUE._l_wp;
      nodes_mask(0)=.FALSE._l_wp
      ALLOCATE(arcs_mask(-n_special_arcs:n_arcs),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=.TRUE._l_wp,&
      caller="CreateLatticeNetwork",alloc_status=error_status)
      arcs_mask=.TRUE._l_wp;
      arcs_mask(0)=.FALSE._l_wp
      
      ALLOCATE(nodes_potentials(-n_special_nodes:n_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1,mold=1.0_r_wp,&
      caller="CreateNetwork",alloc_status=error_status)
      ALLOCATE(supplies_demands(-n_special_nodes:n_nodes),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1,mold=1.0_r_wp,&
      caller="CreateNetwork",alloc_status=error_status)
      ALLOCATE(arcs_flows(-n_special_arcs:n_arcs),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=1.0_r_wp,&
      caller="CreateNetwork",alloc_status=error_status)
      ALLOCATE(arcs_voltages(-n_special_arcs:n_arcs),STAT=error_status)
      CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=1.0_r_wp,&
      caller="CreateNetwork",alloc_status=error_status)
      
      
      
      END SUBROUTINE CreateNetwork
      
      
      
      SUBROUTINE AssignSuppliesDemands()
      IMPLICIT NONE
      
      INTEGER(KIND=i_wp)::node,arc,special_node,special_arc,source,sink
      INTEGER::dim
      REAL(KIND=r_wp)::injection_flow_scaling,injection_flow_mean,injection_flow_std,&
      injection_flow_lb,injection_flow_ub,injection_flow,&
      outflow_flow_scaling,outflow_flow_mean,outflow_flow_std,&
      outflow_flow_lb,outflow_flow_ub,outflow_flow,&
      total_injection_flow,total_outflow_flow,total_flow,mean_flow
      
      supplies_demands=0.0_r_wp
      
      DO dim=1,n_dim
      IF(boundaries(bc_source_sink,dim))THEN
      source=sources(dim)
      supplies_demands(source)=sources_flow_capacity(dim)
      sink=sinks(dim)
      supplies_demands(sink)=sinks_flow_capacity(dim)
      END IF
      END DO
      
      
      SuppliesDemands_injection:IF(n_injection_nodes>0)THEN
      
      injection_flow=injection_flow_capacity/REAL(n_injection_nodes,KIND=r_wp)
      
      SELECT CASE(injection_flow_distribution)
      CASE('U','u')
      mean_flow=0.5_r_dp*(injection_flow_parameters(1)+injection_flow_parameters(2))
      IF(ABS(mean_flow)>100.0_r_wp*EPSILON(1.0_r_wp))THEN
      injection_flow_scaling=injection_flow/mean_flow
      ELSE
      injection_flow_scaling=1.0_r_wp
      END IF
      injection_flow_lb=injection_flow_scaling*injection_flow_parameters(1)
      
      injection_flow_ub=injection_flow_scaling*injection_flow_parameters(2)
      
      CASE('N','n','G','g')
      mean_flow=injection_flow_parameters(1)
      IF(ABS(mean_flow)>10.0_r_wp*EPSILON(1.0_r_wp))THEN
      injection_flow_scaling=injection_flow/mean_flow
      ELSE
      injection_flow_scaling=1.0_r_wp
      END IF
      injection_flow_mean=injection_flow_scaling*injection_flow_parameters(1)
      injection_flow_std=injection_flow_scaling*injection_flow_parameters(2)
      ENDSELECT
      
      total_injection_flow=0.0_r_wp
      DO node=1,n_nodes
      IF(nodes_status(node)==injection_node)THEN
      SELECT CASE(injection_flow_distribution)
      CASE('U','u')
      CALL RandomUniform(supplies_demands(node),&
      range=(/injection_flow_lb,injection_flow_ub/))
      CASE('N','n','G','g')
      CALL RandomNormal(supplies_demands(node),&
      mean_std=(/injection_flow_mean,injection_flow_std/))
      CASE DEFAULT
      supplies_demands(node)=injection_flow
      ENDSELECT
      total_injection_flow=total_injection_flow+supplies_demands(node)
      END IF
      END DO
      
      CorrectFlow_injection:DO node=1,n_nodes
      IF(nodes_status(node)==injection_node)THEN
      supplies_demands(node)=supplies_demands(node)&
      +injection_flow_capacity-total_injection_flow
      EXIT CorrectFlow_injection
      END IF
      END DO CorrectFlow_injection
      
      END IF SuppliesDemands_injection
      
      SuppliesDemands_outflow:IF(n_outflow_nodes>0)THEN
      
      outflow_flow=outflow_flow_capacity/REAL(n_outflow_nodes,KIND=r_wp)
      
      SELECT CASE(outflow_flow_distribution)
      CASE('U','u')
      mean_flow=0.5_r_dp*(outflow_flow_parameters(1)+outflow_flow_parameters(2))
      IF(ABS(mean_flow)>100.0_r_wp*EPSILON(1.0_r_wp))THEN
      outflow_flow_scaling=outflow_flow/mean_flow
      ELSE
      outflow_flow_scaling=1.0_r_wp
      END IF
      outflow_flow_lb=outflow_flow_scaling*outflow_flow_parameters(1)
      
      outflow_flow_ub=outflow_flow_scaling*outflow_flow_parameters(2)
      
      CASE('N','n','G','g')
      mean_flow=outflow_flow_parameters(1)
      IF(ABS(mean_flow)>10.0_r_wp*EPSILON(1.0_r_wp))THEN
      outflow_flow_scaling=outflow_flow/mean_flow
      ELSE
      outflow_flow_scaling=1.0_r_wp
      END IF
      outflow_flow_mean=outflow_flow_scaling*outflow_flow_parameters(1)
      outflow_flow_std=outflow_flow_scaling*outflow_flow_parameters(2)
      ENDSELECT
      
      total_outflow_flow=0.0_r_wp
      DO node=1,n_nodes
      IF(nodes_status(node)==outflow_node)THEN
      SELECT CASE(outflow_flow_distribution)
      CASE('U','u')
      CALL RandomUniform(supplies_demands(node),&
      range=(/outflow_flow_lb,outflow_flow_ub/))
      CASE('N','n','G','g')
      CALL RandomNormal(supplies_demands(node),&
      mean_std=(/outflow_flow_mean,outflow_flow_std/))
      CASE DEFAULT
      supplies_demands(node)=outflow_flow
      ENDSELECT
      total_outflow_flow=total_outflow_flow+supplies_demands(node)
      END IF
      END DO
      
      CorrectFlow_outflow:DO node=1,n_nodes
      IF(nodes_status(node)==outflow_node)THEN
      supplies_demands(node)=supplies_demands(node)&
      +outflow_flow_capacity-total_outflow_flow
      EXIT CorrectFlow_outflow
      END IF
      END DO CorrectFlow_outflow
      
      END IF SuppliesDemands_outflow
      
      total_flow=SUM(supplies_demands(-(n_special_nodes-n_fixed_potential_nodes):n_nodes))
      WRITE(message_log_unit,*)"The total supply/demand is: ",total_flow
      
      IF(ABS(total_flow)>MAX(100.0_r_wp*EPSILON(1.0_r_wp),1E-10_r_wp))THEN
      
      CALL CriticalError(message=&
      "The sum of supplies and demands was not zero within loose precision:",&
      caller="CreateNetwork")
      END IF
      
      END SUBROUTINE AssignSuppliesDemands
      
      
      
      END MODULE Network_Geometry