MODULE FF_Cost_Parameters
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Random_Numbers
      USE Network_Data_Structures
      IMPLICIT NONE
      PUBLIC::InitializeCostParameters,CreateCostParameters,DestroyCostParameters,&
      AssignCostParameters,CalculateFeasibleFlow
      PRIVATE
      
      CHARACTER,PUBLIC,SAVE::critical_current_distribution='F',width_distribution='F',&
      resistance_distribution='U'
      REAL(KIND=r_wp),PUBLIC,DIMENSION(2),SAVE::&
      critical_current_parameters=(/2.0_r_wp,0.0_r_wp/),&
      width_parameters=(/0.1_r_wp,0.0_r_wp/),resistance_parameters=(/0.5_r_wp,1.0_r_wp/)
      REAL(KIND=r_wp)::resistance_low,resistance_high,critical_current_low,critical_current_high
      
      NAMELIST/FF_CostParametersOptions/critical_current_distribution,width_distribution,&
      resistance_distribution,critical_current_parameters,&
      width_parameters,resistance_parameters,resistance_low,&
      resistance_high,critical_current_low,critical_current_high
      
      CONTAINS
      
      SUBROUTINE InitializeCostParameters()
      IMPLICIT NONE
      
      REWIND(UNIT=program_options_unit)
      READ(UNIT=program_options_unit,NML=FF_CostParametersOptions,IOSTAT=error_status)
      IF(error_status/=0)&
      CALL CriticalError(message="NAMELIST FF_CostParametersOptions was not read successfully"//&
      " from file "//TRIM(options_file),caller="InitializeCostsParameters")
      WRITE(UNIT=message_log_unit,NML=FF_CostParametersOptions)
      
      END SUBROUTINE InitializeCostParameters
      
      
      
      SUBROUTINE CreateCostParameters()
      n_cost_parameters=5
      
      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)
      
      END SUBROUTINE CreateCostParameters
      
      
      
      SUBROUTINE AssignCostParameters()
      IMPLICIT NONE
      
      INTEGER(KIND=i_wp)::node,arc,special_node,special_arc
      REAL(KIND=r_wp)::critical_current,width,resistance
      
      
      DO special_arc=0,n_special_arcs
      
      SELECT CASE(width_distribution)
      CASE('U','u')
      CALL RandomUniform(width,&
      range=(/width_parameters(1),width_parameters(2)/))
      CASE('N','n','G','g')
      CALL RandomNormal(width,&
      mean_std=(/width_parameters(1),width_parameters(2)/))
      CASE DEFAULT
      width=width_parameters(1)
      END SELECT
      SELECT CASE(special_arcs_status(-special_arc))
      CASE(regular_cost_arc)
      
      SELECT CASE(critical_current_distribution)
      CASE('U','u')
      CALL RandomUniform(critical_current,&
      range=(/critical_current_parameters(1),critical_current_parameters(2)/))
      CASE('N','n','G','g')
      CALL RandomNormal(critical_current,&
      mean_std=(/critical_current_parameters(1),critical_current_parameters(2)/))
      CASE DEFAULT
      critical_current=critical_current_parameters(1)
      END SELECT
      
      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
      arcs_cost_parameters(1,-special_arc)=critical_current
      arcs_cost_parameters(2,-special_arc)=width
      arcs_cost_parameters(3,-special_arc)=resistance
      arcs_cost_parameters(4,-special_arc)=SQRT(arcs_cost_parameters(1,-special_arc)**2&
      +arcs_cost_parameters(2,-special_arc)**2)
      CASE(low_cost_arc)
      arcs_cost_parameters(1,-special_arc)=critical_current_high
      arcs_cost_parameters(2,-special_arc)=width
      arcs_cost_parameters(3,-special_arc)=resistance_low
      arcs_cost_parameters(4,-special_arc)=SQRT(arcs_cost_parameters(1,-special_arc)**2&
      +arcs_cost_parameters(2,-special_arc)**2)
      CASE(high_cost_arc)
      arcs_cost_parameters(1,-special_arc)=critical_current_low
      arcs_cost_parameters(2,-special_arc)=width
      arcs_cost_parameters(3,-special_arc)=resistance_high
      arcs_cost_parameters(4,-special_arc)=SQRT(arcs_cost_parameters(1,-special_arc)**2&
      +arcs_cost_parameters(2,-special_arc)**2)
      CASE DEFAULT
      arcs_cost_parameters(1,-special_arc)=1.0_r_wp
      arcs_cost_parameters(2,-special_arc)=0.01_r_wp
      arcs_cost_parameters(3,-special_arc)=1.0_r_wp
      arcs_cost_parameters(4,-special_arc)=SQRT(arcs_cost_parameters(1,-special_arc)**2&
      +arcs_cost_parameters(2,-special_arc)**2)
      END SELECT
      
      END DO
      
      
      DO arc=1,n_arcs
      
      SELECT CASE(critical_current_distribution)
      CASE('U','u')
      CALL RandomUniform(critical_current,&
      range=(/critical_current_parameters(1),critical_current_parameters(2)/))
      CASE('N','n','G','g')
      CALL RandomNormal(critical_current,&
      mean_std=(/critical_current_parameters(1),critical_current_parameters(2)/))
      CASE DEFAULT
      critical_current=critical_current_parameters(1)
      END SELECT
      
      SELECT CASE(width_distribution)
      CASE('U','u')
      CALL RandomUniform(width,&
      range=(/width_parameters(1),width_parameters(2)/))
      CASE('N','n','G','g')
      CALL RandomNormal(width,&
      mean_std=(/width_parameters(1),width_parameters(2)/))
      CASE DEFAULT
      width=width_parameters(1)
      END SELECT
      
      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
      arcs_cost_parameters(1,arc)=critical_current
      arcs_cost_parameters(2,arc)=width
      arcs_cost_parameters(3,arc)=resistance
      arcs_cost_parameters(4,arc)=SQRT(arcs_cost_parameters(1,arc)**2&
      +arcs_cost_parameters(2,arc)**2)
      END DO
      END SUBROUTINE AssignCostParameters
      
      
      
      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=max_cost_tree,arcs_weights=ABS(arcs_cost_parameters(1,:)),&
      weights_distribution=critical_current_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_nodes_orientations,parents=tree_nodes_parents,&
      level_ordering=tree_nodes_ordering,&
      arcs_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 FF_Cost_Parameters