MODULE Power_Cost_Parameters
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Random_Numbers
      USE Network_Data_Structures
      IMPLICIT NONE
      PUBLIC::InitializeCostParameters,CreateCostParameters,DestroyCostParameters,&
      CalculateFeasibleFlow
      PRIVATE
      
      REAL(KIND=r_wp),PUBLIC,SAVE::alpha=2.0,inv_alpha_minus_1
      REAL(KIND=r_wp),PUBLIC,SAVE::resistance_low=0.1E-3_r_wp,capacity_low=1.0_r_wp,&
      resistance_high=0.1E3_r_wp,capacity_high=1.0_r_wp
      CHARACTER,PUBLIC,SAVE::resistance_distribution='F',capacity_distribution='F'
      REAL(KIND=r_wp),PUBLIC,DIMENSION(2),SAVE::resistance_parameters=(/1.0_r_wp,0.0_r_wp/),capacity_parameters=(/1.0_r_wp,0.0_r_w&
     &p/)
      CHARACTER,PUBLIC,SAVE::costs_distribution='R'
      
      NAMELIST/CostParametersOptions/alpha,&
      resistance_low,capacity_low,resistance_high,capacity_high,&
      resistance_distribution,capacity_distribution,&
      resistance_parameters,capacity_parameters
      
      CONTAINS
      
      
      SUBROUTINE InitializeCostParameters()
      IMPLICIT NONE
      
      REWIND(UNIT=program_options_unit)
      READ(UNIT=program_options_unit,NML=CostParametersOptions,IOSTAT=error_status)
      IF(error_status/=0)&
      CALL CriticalError(message="NAMELIST CostParametersOptions was not read successfully"//&
      " from file "//TRIM(options_file),caller="InitializeCostsParameters")
      WRITE(UNIT=message_log_unit,NML=CostParametersOptions)
      
      inv_alpha_minus_1=1.0_r_wp/(alpha-1.0_r_wp)
      
      END SUBROUTINE InitializeCostParameters
      
      
      
      SUBROUTINE CreateCostParameters()
      
      INTEGER(KIND=i_wp)::node,arc,special_node,special_arc
      REAL(KIND=r_wp)::resistance,capacity
      
      n_cost_parameters=1
      
      ALLOCATE(arcs_cost_parameters(n_cost_parameters,-n_special_arcs:n_arcs),&
      STAT=error_status)
      CALL RecordAllocation(n_elements=n_cost_parameters*(n_special_arcs+n_arcs+1),&
      mold=1.0_r_wp,caller="CreateCostParameters",alloc_status=error_status)
      
      DO special_arc=0,n_special_arcs
      SELECT CASE(special_arcs_status(-special_arc))
      CASE(regular_cost_arc)
      
      SELECT CASE(resistance_distribution)
      CASE('U','u')
      CALL RandomUniform(resistance,&
      range=(/resistance_parameters(1),resistance_parameters(2)/))
      CASE('N','n','G','g')
      CALL RandomNormal(resistance,&
      mean_std=(/resistance_parameters(1),resistance_parameters(2)/))
      CASE DEFAULT
      resistance=resistance_parameters(1)
      END SELECT
      
      SELECT CASE(capacity_distribution)
      CASE('U','u')
      CALL RandomUniform(capacity,&
      range=(/capacity_parameters(1),capacity_parameters(2)/))
      CASE('N','n','G','g')
      CALL RandomNormal(capacity,&
      mean_std=(/capacity_parameters(1),capacity_parameters(2)/))
      CASE DEFAULT
      capacity=capacity_parameters(1)
      END SELECT
      arcs_cost_parameters(1,-special_arc)=ABS(resistance/capacity**alpha)
      CASE(low_cost_arc)
      arcs_cost_parameters(1,-special_arc)=resistance_low/capacity_low**alpha
      CASE(high_cost_arc)
      arcs_cost_parameters(1,-special_arc)=resistance_high/capacity_high**alpha
      CASE DEFAULT
      arcs_cost_parameters(1,-special_arc)=0.0_r_wp
      END SELECT
      END DO
      
      DO arc=1,n_arcs
      
      SELECT CASE(resistance_distribution)
      CASE('U','u')
      CALL RandomUniform(resistance,&
      range=(/resistance_parameters(1),resistance_parameters(2)/))
      CASE('N','n','G','g')
      CALL RandomNormal(resistance,&
      mean_std=(/resistance_parameters(1),resistance_parameters(2)/))
      CASE DEFAULT
      resistance=resistance_parameters(1)
      END SELECT
      
      SELECT CASE(capacity_distribution)
      CASE('U','u')
      CALL RandomUniform(capacity,&
      range=(/capacity_parameters(1),capacity_parameters(2)/))
      CASE('N','n','G','g')
      CALL RandomNormal(capacity,&
      mean_std=(/capacity_parameters(1),capacity_parameters(2)/))
      CASE DEFAULT
      capacity=capacity_parameters(1)
      END SELECT
      arcs_cost_parameters(1,arc)=ABS(resistance/capacity**alpha)
      END DO
      
      IF(capacity_distribution=='F')costs_distribution=resistance_distribution
      
      END SUBROUTINE CreateCostParameters
      
      
      
      SUBROUTINE DestroyCostParameters()
      
      CALL RecordAllocation(n_elements=-INT(SIZE(arcs_cost_parameters),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(arcs_cost_parameters)
      
      END SUBROUTINE DestroyCostParameters
      
      
      
      SUBROUTINE CalculateFeasibleFlow()
      USE Graph_Algorithms
      USE Network_Spanning_Trees
      IMPLICIT NONE
      
      REAL(KIND=r_wp)::flow_infeasibility
      
      CALL CreateSpanningTree(tree_type=min_cost_tree,arcs_weights=ABS(arcs_cost_parameters(1,:)),&
      weights_distribution=costs_distribution)
      arcs_flows=0.0_r_wp
      CALL PropagateArcsFlows(arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      heads_tails=heads_tails,supplies_demands=-supplies_demands,&
      orientations=tree_arcs_orientations,parents=tree_nodes_parents,&
      level_ordering=tree_nodes_ordering,&
      tree_flows=arcs_flows,flow_imbalance=flow_infeasibility)
      WRITE(message_log_unit,*)"Calculated initial flow with flow excess of: ",flow_infeasibility
      
      END SUBROUTINE CalculateFeasibleFlow
      
      
      
      END MODULE Power_Cost_Parameters