How to use the wntr.metrics function in wntr

def population_impacted_metrics(wn):
    # Compute population per node
    pop = wntr.metrics.population(wn)
    total_population = pop.sum()
    print("Total population: " + str(total_population))
    wntr.graphics.plot_network(wn, node_attribute=pop, node_range = [0,400], node_size=40,
                          title='Population, Total = ' + str(total_population))

    # Find population and nodes impacted by pressure less than 40 m
    sim = wntr.sim.EpanetSimulator(wn)
    results = sim.run_sim()
    pressure = results.node['pressure'].loc[:,wn.junction_name_list]
    pop_impacted = wntr.metrics.population_impacted(pop, pressure, np.less, 40)
    plt.figure()
    pop_impacted.sum(axis=1).plot(title='Total population with pressure < 40 m')
    nodes_impacted = wntr.metrics.query(pressure, np.less, 40)
    wntr.graphics.plot_network(wn, node_attribute=nodes_impacted.any(axis=0), node_size=40,
                          title='Nodes impacted')

sim = wntr.sim.EpanetSimulator(wn)
results = sim.run_sim()
junctions = [name for name, node in wn.junctions()]
pop_impacted = wntr.metrics.population_impacted(pop, results.node['pressure',:,junctions], np.less, 40)
plt.figure()
pop_impacted.sum(axis=1).plot(title='Total population with pressure < 40 m') 
nodes_impacted = wntr.metrics.query(results.node['pressure',:,junctions], np.less, 40)
wntr.network.draw_graph(wn, node_attribute=nodes_impacted.any(axis=0), node_size=40,
                      title='Nodes impacted')

# Copute network cost
network_cost = wntr.metrics.cost(wn)
print "Network cost: $" + str(round(network_cost,2))

# COmpute green house gas emissions
network_ghg = wntr.metrics.ghg_emissions(wn)
print "Network GHG emissions: " + str(round(network_ghg,2))

# Isolate node results at consumer nodes (nzd = non-zero demand)
nzd_junctions = wn.query_node_attribute('base_demand', np.greater, 0, 
                                        node_type=wntr.network.Junction).keys()
node_results = results.node.loc[:, :, nzd_junctions]

# Compute population per node
pop = wntr.metrics.population(wn)

# Compute FDV, FDD, and FDQ, change average_times and average_nodes to average results
quality_upper_bound = 0.0035 # kg/m3 (3.5 mg/L)                   
demand_factor = 0.9
average_times = False
average_nodes = False
fdv = wntr.metrics.fdv(node_results, average_times, average_nodes)
fdd = wntr.metrics.fdd(node_results, demand_factor, average_times, average_nodes)
#fdq = wntr.metrics.fdq(node_results, quality_upper_bound, average_times, average_nodes)

# Plot results
if average_times == False and average_nodes == False:
    fdv.plot(ylim=(-0.05, 1.05), legend=False, title='FDV for each node and time')
    
    fdd.plot(ylim=(-0.05, 1.05),legend=False,  title='FDD for each node and time')
    # Fraction of nodes with reduced demand
    fdd_fraction_impacted = 1-fdd.sum(axis=1)/fdd.shape[1]
    plt.figure()
    fdd_fraction_impacted.plot(ylim=(-0.05, 1.05), title='Fraction of nodes not receiving adaquate demand')        
    # Population impacted by reduced demand
    fdd_pop_impacted = wntr.metrics.population_impacted(pop, fdd, np.less, 1)
    plt.figure()
    fdd_pop_impacted.plot(legend=False, title='Population not receiving adaquate demand')        
    # Timeseries of fraction of population not meeting demand

import wntr
import matplotlib.pyplot as plt

# Create a water network model
inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)

# Define WQ scenarios
WQscenario = wntr.scenario.Waterquality('CHEM', '121', 'SETPOINT', 1000, 2*3600, 15*3600)

# Simulate hydraulics and water quality for each scenario
sim = wntr.sim.EpanetSimulator(wn)
results_CHEM = sim.run_sim(WQscenario)

MC = wntr.metrics.mass_contaminant_consumed(results_CHEM.node)
VC = wntr.metrics.volume_contaminant_consumed(results_CHEM.node, 0.001)
EC = wntr.metrics.extent_contaminant(results_CHEM.node, results_CHEM.link, wn, 0.001)

wntr.network.draw_graph(wn, node_attribute=MC.sum(axis=0), node_range = [0,400], node_size=40,
                      title='Total mass consumed')
    
plt.figure()           
EC.sum(axis=1).plot(title='Extent of contamination')

str(time_of_failure) + ', End Time: ' + \
                str(time_of_failure+duration_of_failure))
    results[i] = sim.run_sim()
    
    # Reload the water network model
    f=open('wn.pickle','rb')
    wn = pickle.load(f)
    f.close()

# Plot water service availability and tank water level for each realization
for i in results.keys():
    
    # Water service availability at each junction and time
    expected_demand = wntr.metrics.expected_demand(wn)
    demand = results[i].node['demand'].loc[:,wn.junction_name_list]
    wsa_nt = wntr.metrics.water_service_availability(expected_demand, demand)
    
    # Average water service availability at each time
    wsa_t = wntr.metrics.water_service_availability(expected_demand.sum(axis=1), 
                                                  demand.sum(axis=1))
                               
    # Tank water level
    tank_level = results[i].node['pressure'].loc[:,wn.tank_name_list]
    
    # Plot results
    plt.figure()
    
    plt.subplot(2,1,1)
    wsa_nt.plot(ax=plt.gca(), legend=False)
    wsa_t.plot(ax=plt.gca(), label='Average', color='k', linewidth=3.0, legend=False)
    plt.ylim( (-0.05, 1.05) )
    plt.ylabel('Water service availability')

sim = wntr.sim.WNTRSimulator(wn, mode='PDD')
    print('Pipe Breaks: ' + str(pipes_to_fail) + ', Start Time: ' + \
                str(time_of_failure) + ', End Time: ' + \
                str(time_of_failure+duration_of_failure))
    results[i] = sim.run_sim()
    
    # Reload the water network model
    f=open('wn.pickle','rb')
    wn = pickle.load(f)
    f.close()

# Plot water service availability and tank water level for each realization
for i in results.keys():
    
    # Water service availability at each junction and time
    expected_demand = wntr.metrics.expected_demand(wn)
    demand = results[i].node['demand'].loc[:,wn.junction_name_list]
    wsa_nt = wntr.metrics.water_service_availability(expected_demand, demand)
    
    # Average water service availability at each time
    wsa_t = wntr.metrics.water_service_availability(expected_demand.sum(axis=1), 
                                                  demand.sum(axis=1))
                               
    # Tank water level
    tank_level = results[i].node['pressure'].loc[:,wn.tank_name_list]
    
    # Plot results
    plt.figure()
    
    plt.subplot(2,1,1)
    wsa_nt.plot(ax=plt.gca(), legend=False)
    wsa_t.plot(ax=plt.gca(), label='Average', color='k', linewidth=3.0, legend=False)

for t in np.arange(0, 24*3600+1,3600):
        attr = results.link['flowrate'].loc[t,:]
        G_flowrate_t.weight_graph(link_attribute=attr)
        entropy = wntr.metrics.entropy(G_flowrate_t)
        shat.append(entropy[1])
    plt.figure()
    plt.plot(shat)
    plt.ylabel('System Entropy')
    plt.xlabel('Time, hr')
    print("Entropy")
    print("  Mean: " + str(np.mean(shat)))
    print("  Max: " + str(np.nanmax(shat)))
    print("  Min: " + str(np.nanmin(shat)))

    # Compute water service availability, for each junction
    expected_demand = wntr.metrics.expected_demand(wn)
    demand = results.node['demand'].loc[:,wn.junction_name_list]
    wsa = wntr.metrics.water_service_availability(expected_demand.sum(axis=0), 
                                                  demand.sum(axis=0))
    wntr.graphics.plot_network(wn, node_attribute=wsa, node_size=40, node_range=[0,1], 
                               title='Water service availability, averaged over times')

# Simulate hydraulics and water quality
inp_file = 'networks/Net3.inp'
wn = wntr.network.WaterNetworkModel(inp_file)
wn.options.duration = 48*3600
for name, node in wn.junctions():
    node.nominal_pressure = 60
sim = wntr.sim.ScipySimulator(wn, pressure_driven=True)
results = sim.run_sim()
       
# Isolate node results at consumer nodes (nzd = non-zero demand)
nzd_junctions = wn.query_node_attribute('base_demand', np.greater, 0, 
                                        node_type=wntr.network.Junction).keys()
node_results = results.node.loc[:, :, nzd_junctions]

# Compute population per node
pop = wntr.metrics.population(wn)

# Compute FDV, FDD, and FDQ, change average_times and average_nodes to average results
quality_upper_bound = 0.0035 # kg/m3 (3.5 mg/L)                   
demand_factor = 0.9
average_times = False
average_nodes = False
fdv = wntr.metrics.fdv(node_results, average_times, average_nodes)
fdd = wntr.metrics.fdd(node_results, demand_factor, average_times, average_nodes)
#fdq = wntr.metrics.fdq(node_results, quality_upper_bound, average_times, average_nodes)

# Plot results
if average_times == False and average_nodes == False:
    fdv.plot(ylim=(-0.05, 1.05), legend=False, title='FDV for each node and time')
    
    fdd.plot(ylim=(-0.05, 1.05),legend=False,  title='FDD for each node and time')
    # Fraction of nodes with reduced demand

def cost_ghg_metrics(wn):
    # Compute network cost
    network_cost = wntr.metrics.annual_network_cost(wn)
    print("Network cost: $" + str(round(network_cost,2)))

    # Compute green house gas emissions
    network_ghg = wntr.metrics.annual_ghg_emissions(wn)
    print("Network GHG emissions: " + str(round(network_ghg,2)))

junctions = [name for name, node in wn.junctions()]

# Define pressure lower bound
P_lower = 21.09 # m (30 psi)

# Query pressure
pressure = results.node.loc['pressure', :, junctions]
mask = wntr.metrics.query(pressure, np.greater, P_lower)
pressure_regulation = mask.all(axis=0).sum() # True over all time
print "Fraction of nodes > 30 psi: " + str(pressure_regulation)
print "Average node pressure: " +str(pressure.mean().mean()) + " m"
wntr.network.draw_graph(wn, node_attribute=pressure.min(axis=0), node_size=40, 
                      title= 'Min pressure')
              
# Compute todini index
todini = wntr.metrics.todini(results.node,results.link,wn, P_lower)
plt.figure()
plt.plot(todini)
plt.ylabel('Todini Index')
plt.xlabel('Time, hr')
print "Todini Index"
print "  Mean: " + str(np.mean(todini))
print "  Max: " + str(np.max(todini))
print "  Min: " + str(np.min(todini))

# Create a weighted graph for flowrate at time 36 hours
t = 36*3600
attr = results.link.loc['flowrate', t, :]   
G_flowrate_36hrs = wn.get_graph_deep_copy()
G_flowrate_36hrs.weight_graph(link_attribute=attr)     
 
# Compute betweenness-centrality time 36 hours