From 13e931b9a2a8eac2ff9d04fa1f4e61b4f5d989ec Mon Sep 17 00:00:00 2001
From: Kevin Zhang <68650998+kzhang81@users.noreply.github.com>
Date: Fri, 13 Jan 2023 12:06:51 -0500
Subject: [PATCH 1/4] Add files via upload

---
 changelog.md | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/changelog.md b/changelog.md
index 378b2d0..70581b5 100644
--- a/changelog.md
+++ b/changelog.md
@@ -1,5 +1,16 @@
 # FTOT Change Log
 
+## v2022_4
+The FTOT 2022.4 public release includes updates related to candidate generation using network density reduction (NDR), Tableau and reporting enhancements, detailed reporting on shortest paths between origin and destination facilities, and other minor bug fixes. The following changes have been made:
+- Implements a new method for candidate generation scenarios through network density reduction. The new method incorporates NetworkX shortest path algorithms into (1) selection of candidate processor locations and (2) calculation of the optimal routing solution. Running candidate generation using NDR significantly decreases runtime by pre-solving the network for shortest paths; see documentation for details on methodology and known issues.
+- Adds metrics for processor capacity and utilization to the text report and CSV report.
+- Updates the Tableau workbook with new information on facility sizes, processor efficiencies and utilizations, and other minor aesthetic improvements.
+- Expands reporting on shortest path routes between facilities in the "all_routes" CSV report (formerly "optimal_routes") and a new Tableau routes dashboard when NDR is enabled. Reporting allows users to analyze routes across different scenarios, modes, facility pairs, and optimal vs non-optimal routes. The scenario comparison dashboard has also been updated to reflect these enhancements.
+- Corrects a minor bug on converting candidate processor build costs to default units.
+- Replaces "kgal" with "thousand_gallon" as default liquid units in the template and example scenario XML files.
+
+See documentation files for additional details.
+
 ## v2022_3
 The FTOT 2022.3 public release includes updates related to the FTOT multimodal GIS network, Tableau and reporting enhancements, user experience (including video tutorials and supplemental tools), and resilience analysis tools. The following changes have been made:
 - Refreshes the default continental United States multimodal network accompanying FTOT for 2022. Key highlights include a more detailed road network, updated pipeline tariff data, and rail and waterway networks have also been updated. For capacity-constrained scenarios, continue to use the older (2019) version of the network.

From 9d13906ef4fd19c81fee567a1bab12505addf9ad Mon Sep 17 00:00:00 2001
From: Kevin Zhang <68650998+kzhang81@users.noreply.github.com>
Date: Fri, 13 Jan 2023 12:07:20 -0500
Subject: [PATCH 2/4] Add files via upload

---
 program/ftot.py                           |   7 +-
 program/ftot_facilities.py                |  36 +-
 program/ftot_networkx.py                  | 508 +++++++++++++++++++---
 program/ftot_postprocess.py               | 126 +++++-
 program/ftot_processor.py                 |   6 +-
 program/ftot_pulp.py                      | 210 +++++----
 program/ftot_pulp_candidate_generation.py | 289 ++++++++----
 program/ftot_pulp_sourcing.py             |  12 +-
 program/ftot_report.py                    | 114 +++--
 program/ftot_scenario.py                  |   6 +-
 program/ftot_supporting.py                |  20 +-
 11 files changed, 1024 insertions(+), 310 deletions(-)

diff --git a/program/ftot.py b/program/ftot.py
index 54143ad..d9ffedb 100644
--- a/program/ftot.py
+++ b/program/ftot.py
@@ -19,6 +19,7 @@
 ureg = UnitRegistry()
 Q_ = ureg.Quantity
 ureg.define('usd = [currency]')  # add the US Dollar, "USD" to the unit registry
+ureg.define('thousand_gallon = kgal')
 # solves issue in pint 0.9
 if pint.__version__ == 0.9:
     ureg.define('short_hundredweight = short_hunderdweight')
@@ -26,9 +27,9 @@
     ureg.define('us_ton = US_ton')
 
 
-FTOT_VERSION = "2022.3"
-SCHEMA_VERSION = "6.0.3"
-VERSION_DATE = "10/5/2022"
+FTOT_VERSION = "2022.4"
+SCHEMA_VERSION = "6.0.4"
+VERSION_DATE = "1/13/2023"
 
 # ===================================================================================================
 
diff --git a/program/ftot_facilities.py b/program/ftot_facilities.py
index 4141151..b2b074c 100644
--- a/program/ftot_facilities.py
+++ b/program/ftot_facilities.py
@@ -286,7 +286,10 @@ def db_report_commodity_potentials(the_scenario, logger):
 
     # -----------------------------------
     with sqlite3.connect(the_scenario.main_db) as db_con:
-        sql = """ select c.commodity_name, fti.facility_type, fc.io, case when sum(fc.scaled_quantity) is null then 'Unconstrained' else sum(fc.scaled_quantity) end as scaled_quantity, fc.units
+        sql = """ select c.commodity_name, fti.facility_type, fc.io,
+                  (case when sum(case when fc.scaled_quantity is null then 1 else 0 end) = 0 then sum(fc.scaled_quantity)
+                  else 'Unconstrained' end) as scaled_quantity,
+                  fc.units
                   from facility_commodities fc
                   join commodities c on fc.commodity_id = c.commodity_id
                   join facilities f on f.facility_id = fc.facility_id
@@ -312,7 +315,10 @@ def db_report_commodity_potentials(the_scenario, logger):
         # Add the ignored processing capacity.
         # Note this doesn't happen until the bx step.
         # -------------------------------------------
-        sql = """ select c.commodity_name, fti.facility_type, fc.io, case when sum(fc.scaled_quantity) is null then 'Unconstrained' else sum(fc.scaled_quantity) end as scaled_quantity, fc.units, f.ignore_facility
+        sql = """ select c.commodity_name, fti.facility_type, fc.io,
+                  (case when sum(case when fc.scaled_quantity is null then 1 else 0 end) = 0 then sum(fc.scaled_quantity)
+                  else 'Unconstrained' end) as scaled_quantity,
+                  fc.units, f.ignore_facility
                   from facility_commodities fc
                   join commodities c on fc.commodity_id = c.commodity_id
                   join facilities f on f.facility_id = fc.facility_id
@@ -339,7 +345,10 @@ def db_report_commodity_potentials(the_scenario, logger):
 
             # Report out net quantities with ignored facilities removed from the query
             # -------------------------------------------------------------------------
-            sql = """ select c.commodity_name, fti.facility_type, fc.io, case when sum(fc.scaled_quantity) is null then 'Unconstrained' else sum(fc.scaled_quantity) end as scaled_quantity, fc.units
+            sql = """ select c.commodity_name, fti.facility_type, fc.io,
+                      (case when sum(case when fc.scaled_quantity is null then 1 else 0 end) = 0 then sum(fc.scaled_quantity)
+                      else 'Unconstrained' end) as scaled_quantity,
+                      fc.units
                       from facility_commodities fc
                       join commodities c on fc.commodity_id = c.commodity_id
                       join facilities f on f.facility_id = fc.facility_id
@@ -553,7 +562,7 @@ def load_facility_commodities_input_data(the_scenario, commodity_input_file, log
                 logger.debug('the CSV file has a blank in the first column. Skipping this line: {}'.format(
                     list(row.values())))
                 continue
-            # {'units': 'kgal', 'facility_name': 'd:01053', 'phase_of_matter': 'liquid', 'value': '9181.521484',
+            # {'units': 'thousand_gallon', 'facility_name': 'd:01053', 'phase_of_matter': 'liquid', 'value': '9181.521484',
             # 'commodity': 'diesel', 'io': 'o', 'share_max_transport_distance'; 'Y'}
             io                  = row["io"]
             facility_name       = str(row["facility_name"])
@@ -664,11 +673,22 @@ def load_facility_commodities_input_data(the_scenario, commodity_input_file, log
             if commodity_phase.lower() == 'solid':
                 commodity_unit = the_scenario.default_units_solid_phase
 
-            if commodity_name == 'cost_formula':
-                pass
-            else:
+            if commodity_name == 'cost_formula': # handle cost formula units
+                commodity_unit = commodity_unit.split("/")[-1]  # get the denominator from string version
+                if str(ureg(commodity_unit).dimensionality) == '[length] ** 3' : # if denominator unit looks like a volume
+                    commodity_phase = 'liquid' # then phase is liquid
+                    commodity_unit = ureg.usd/the_scenario.default_units_liquid_phase # and cost unit phase should use default liquid
+                elif str(ureg(commodity_unit).dimensionality) == '[mass]': # if denominator unit looks like a mass
+                    commodity_phase = 'solid' # then phase is solid
+                    commodity_unit = ureg.usd/the_scenario.default_units_solid_phase # and cost unit phase should use default solid
+            
+            # now that we process cost_formula, all properties can use this (with try statement in case)
+            try:
                 commodity_quantity = commodity_quantity_and_units.to(commodity_unit).magnitude
-
+            except Exception as e:
+                logger.error("FAIL: {} ".format(e))
+                raise Exception("FAIL: {}".format(e))
+            
             if max_processor_input != 'Null':
                 max_processor_input = max_input_quantity_and_units.to(commodity_unit).magnitude
             if min_processor_input != 'Null':
diff --git a/program/ftot_networkx.py b/program/ftot_networkx.py
index 42b678d..08912ae 100644
--- a/program/ftot_networkx.py
+++ b/program/ftot_networkx.py
@@ -15,6 +15,8 @@
 import os
 import multiprocessing
 import math
+from heapq import heappush, heappop
+from itertools import count
 from ftot_pulp import commodity_mode_setup
 from ftot import Q_
 
@@ -189,6 +191,9 @@ def presolve_network(the_scenario, G, logger):
     # Make subgraphs for combination of permitted modes
     commodity_subgraph_dict = make_mode_subgraphs(the_scenario, G, logger)
 
+    # if MTD, then make rmp-specific subgraphs with only nodes reachable in MTD
+    commodity_subgraph_dict = make_max_transport_distance_subgraphs(the_scenario, logger, commodity_subgraph_dict)
+
     # Create a dictionary of edge_ids from the database which is used later to uniquely identify edges
     edge_id_dict = find_edge_ids(the_scenario, logger)
 
@@ -208,12 +213,24 @@ def presolve_network(the_scenario, G, logger):
     # by commodity and destination
     for commodity_id in od_pairs:
         phase_of_matter = od_pairs[commodity_id]['phase_of_matter']
-        allowed_modes = commodity_subgraph_dict[commodity_id]['modes']
-        for a_target in od_pairs[commodity_id]['targets'].keys():
-            stuff_to_pass.append([commodity_subgraph_dict[commodity_id]['subgraph'],
-                                  od_pairs[commodity_id]['targets'][a_target],
-                                  a_target, all_route_edges, no_path_pairs,
-                                  edge_id_dict, phase_of_matter, allowed_modes])
+        allowed_modes = commodity_subgraph_dict[commodity_id]['modes'] 
+
+        if 'targets' in od_pairs[commodity_id].keys():
+            for a_target in od_pairs[commodity_id]['targets'].keys():
+                stuff_to_pass.append([commodity_subgraph_dict[commodity_id]['subgraph'],
+                                    od_pairs[commodity_id]['targets'][a_target],
+                                    a_target, all_route_edges, no_path_pairs,
+                                    edge_id_dict, phase_of_matter, allowed_modes, 'target'])
+        else:    
+            for a_source in od_pairs[commodity_id]['sources'].keys():
+                if 'facility_subgraphs' in commodity_subgraph_dict[commodity_id].keys():
+                    subgraph = commodity_subgraph_dict[commodity_id]['facility_subgraphs'][a_source]
+                else:
+                    subgraph = commodity_subgraph_dict[commodity_id]['subgraph']
+                stuff_to_pass.append([subgraph,
+                                    od_pairs[commodity_id]['sources'][a_source],
+                                    a_source, all_route_edges, no_path_pairs,
+                                    edge_id_dict, phase_of_matter, allowed_modes, 'source'])
 
     # Allow multiprocessing, with no more than 75% of cores to be used, rounding down if necessary
     logger.info("start: the multiprocessing route solve.")
@@ -271,7 +288,7 @@ def presolve_network(the_scenario, G, logger):
 # -----------------------------------------------------------------------------
 
 
-# Ensure shortest paths not calculated in presence of candidate generation, capacity, max transport distance
+# Ensure shortest paths not calculated in presence of capacity
 def update_ndr_parameter(the_scenario, logger):
     logger.debug("start: check NDR conditions")
     new_ndr_parameter = the_scenario.ndrOn
@@ -281,21 +298,6 @@ def update_ndr_parameter(the_scenario, logger):
         new_ndr_parameter = False
         logger.debug("NDR de-activated due to capacity enforcement")
 
-    # if candidate generation is active, then skip NDR
-    if the_scenario.processors_candidate_slate_data != 'None':
-        new_ndr_parameter = False
-        logger.debug("NDR de-activated due to candidate generation")
-
-    # if max_transport_distance field is used in commodities table, then skip NDR
-    with sqlite3.connect(the_scenario.main_db) as main_db_con:
-        # get count of commodities with a specified max_transport_distance
-        sql = "select count(commodity_id) from commodities where max_transport_distance is not null;"
-        db_cur = main_db_con.execute(sql)
-        count = db_cur.fetchone()[0]
-        if count > 0:
-            new_ndr_parameter = False
-            logger.debug("NDR de-activated due to use of max transport distance")
-
     the_scenario.ndrOn = new_ndr_parameter
     logger.debug("finish: check NDR conditions")
 
@@ -307,39 +309,72 @@ def update_ndr_parameter(the_scenario, logger):
 # network that are a part of the shortest path connecting an origin to a destination
 # for each commodity
 def multi_shortest_paths(stuff_to_pass):
+
     global all_route_edges, no_path_pairs
-    G, sources, target, all_route_edges, no_path_pairs, edge_id_dict, phase_of_matter, allowed_modes = stuff_to_pass
-    t = target
-
-    shortest_paths_to_t = nx.shortest_path(G, target=t, weight='{}_weight'.format(phase_of_matter))
-    for a_source in sources:
-        s = a_source
-        # This accounts for when a_source may not be connected to t,
-        # as is the case when certain modes may not be permitted
-        if a_source not in shortest_paths_to_t:
+    if stuff_to_pass[8] == 'target':
+        G, sources, target, all_route_edges, no_path_pairs, edge_id_dict, phase_of_matter, allowed_modes, st_dummy = stuff_to_pass
+        t = target
+        shortest_paths_to_t = nx.shortest_path(G, target=t, weight='{}_weight'.format(phase_of_matter))
+        for a_source in sources:
+            s = a_source
+            # This accounts for when a_source may not be connected to t,
+            # as is the case when certain modes may not be permitted
+            if a_source not in shortest_paths_to_t:
+                for i in sources[a_source]:
+                    rt_id = i
+                    no_path_pairs.append((s, t, rt_id))
+                continue
             for i in sources[a_source]:
                 rt_id = i
-                no_path_pairs.append((s, t, rt_id))
-            continue
-        for i in sources[a_source]:
-            rt_id = i
-            for index, from_node in enumerate(shortest_paths_to_t[s]):
-                if index < (len(shortest_paths_to_t[s]) - 1):
-                    to_node = shortest_paths_to_t[s][index + 1]
-                    # find the correct edge_id on the shortest path
-                    min_route_cost = 999999999
-                    min_edge_id = None
-                    for j in edge_id_dict[to_node][from_node]:
-                        edge_id, mode_source, route_cost = j
-                        if mode_source in allowed_modes:
-                            if route_cost < min_route_cost:
-                                min_route_cost = route_cost
-                                min_edge_id = edge_id
-                    if min_edge_id is None:
-                        error = """something went wrong finding the edge_id from node {} to node {}
-                                for scenario_rt_id {} in shortest path algorithm""".format(from_node, to_node, rt_id)
-                        raise Exception(error)
-                    all_route_edges.append((from_node, to_node, min_edge_id, rt_id, index + 1))
+                for index, from_node in enumerate(shortest_paths_to_t[s]):
+                    if index < (len(shortest_paths_to_t[s]) - 1):
+                        to_node = shortest_paths_to_t[s][index + 1]
+                        # find the correct edge_id on the shortest path
+                        min_route_cost = 999999999
+                        min_edge_id = None
+                        for j in edge_id_dict[to_node][from_node]:
+                            edge_id, mode_source, route_cost = j
+                            if mode_source in allowed_modes:
+                                if route_cost < min_route_cost:
+                                    min_route_cost = route_cost
+                                    min_edge_id = edge_id
+                        if min_edge_id is None:
+                            error = """something went wrong finding the edge_id from node {} to node {}
+                                    for scenario_rt_id {} in shortest path algorithm""".format(from_node, to_node, rt_id)
+                            raise Exception(error)
+                        all_route_edges.append((from_node, to_node, min_edge_id, rt_id, index + 1))
+    else:
+        G, targets, source, all_route_edges, no_path_pairs, edge_id_dict, phase_of_matter, allowed_modes, st_dummy = stuff_to_pass
+        s = source
+        shortest_paths_from_s = nx.shortest_path(G, source=s, weight='{}_weight'.format(phase_of_matter))
+        for a_target in targets:
+            t = a_target
+            # This accounts for when a_target may not be connected to s,
+            # as is the case when certain modes may not be permitted
+            if a_target not in shortest_paths_from_s:
+                for i in targets[a_target]:
+                    rt_id = i
+                    no_path_pairs.append((s, t, rt_id))
+                continue
+            for i in targets[a_target]:
+                rt_id = i
+                for index, from_node in enumerate(shortest_paths_from_s[t]):
+                    if index < (len(shortest_paths_from_s[t]) - 1):
+                        to_node = shortest_paths_from_s[t][index + 1]
+                        # find the correct edge_id on the shortest path
+                        min_route_cost = 999999999
+                        min_edge_id = None
+                        for j in edge_id_dict[to_node][from_node]:
+                            edge_id, mode_source, route_cost = j
+                            if mode_source in allowed_modes:
+                                if route_cost < min_route_cost:
+                                    min_route_cost = route_cost
+                                    min_edge_id = edge_id
+                        if min_edge_id is None:
+                            error = """something went wrong finding the edge_id from node {} to node {}
+                                    for scenario_rt_id {} in shortest path algorithm""".format(from_node, to_node, rt_id)
+                            raise Exception(error)
+                        all_route_edges.append((from_node, to_node, min_edge_id, rt_id, index + 1))
 
 
 # -----------------------------------------------------------------------------
@@ -587,6 +622,55 @@ def make_od_pairs(the_scenario, logger):
         '''
         db_cur.execute(sql)
 
+        if not os.path.exists(the_scenario.processor_candidates_commodity_data) and the_scenario.processors_candidate_slate_data != 'None':
+            # Create temp table for RMP-endcap
+            sql = '''
+            create table tmp_rmp_to_endcap as
+            select en.node_id as endcap_node,
+            en.commodity_id,
+            origin.location_id as rmp_location,
+            origin.facility_id as rmp_facility,
+            origin.node_id as rmp_node,
+            origin.phase_of_matter
+            from
+            endcap_nodes en
+            join tmp_connected_facilities_with_commodities as origin on
+            origin.node_id = en.source_node_id
+            '''
+            db_cur.execute(sql)
+
+            # Create temp table for endcap-proc or endcap-dest
+            # Join endcap_nodes with candidate_process_commodities for process_id
+            # then endcap_nodes with facility_commodities for facilities that take commodity as input
+
+            sql = '''
+            create table tmp_endcap_to_other as
+            select
+            en.node_id as endcap_node,
+            fc.facility_id as facility_id,
+            fc.location_id as location_id,
+            cpc.io,
+            cpc.commodity_id,
+            tmp.node_id as dest_node,
+            tmp.phase_of_matter
+            from
+            endcap_nodes en
+            join candidate_process_commodities cpc on en.process_id = cpc.process_id
+            join facility_commodities fc on
+            case
+            when cpc.io = 'o'
+            then
+            cpc.commodity_id = fc.commodity_id
+            and fc.io= 'i'
+            end
+            join tmp_connected_facilities_with_commodities tmp on
+            fc.facility_id = tmp.facility_id
+            and tmp.io = 'i'
+            and tmp.location_1 like '%_in'
+            ;
+            '''
+            db_cur.execute(sql)
+
         sql = '''
         insert into od_pairs (from_location_id, to_location_id, from_facility_id, to_facility_id, commodity_id, 
         phase_of_matter, from_node_id, to_node_id, from_location_1, to_location_1)
@@ -619,12 +703,45 @@ def make_od_pairs(the_scenario, logger):
         '''
         db_cur.execute(sql)
 
+        if not os.path.exists(the_scenario.processor_candidates_commodity_data) and the_scenario.processors_candidate_slate_data != 'None':
+            sql = '''
+            insert into od_pairs (from_location_id, from_facility_id, commodity_id, phase_of_matter, from_node_id, to_node_id)
+            select
+            tmp.rmp_location,
+            tmp.rmp_facility,
+            tmp.commodity_id,
+            tmp.phase_of_matter,
+            tmp.rmp_node,
+            tmp.endcap_node
+            from tmp_rmp_to_endcap as tmp
+            '''
+            db_cur.execute(sql)
+
+            sql = '''
+            insert into od_pairs (to_location_id, to_facility_id, commodity_id, phase_of_matter, from_node_id, to_node_id)
+            select
+            tmp.location_id,
+            tmp.facility_id,
+            tmp.commodity_id,
+            tmp.phase_of_matter,
+            tmp.endcap_node,
+            tmp.dest_node
+            from tmp_endcap_to_other as tmp
+            '''
+            db_cur.execute(sql)
+
         logger.debug("drop the tmp_connected_facilities_with_commodities table")
         db_cur.execute("drop table if exists tmp_connected_facilities_with_commodities;")
 
         logger.debug("drop the tmp_od_pairs table")
         db_cur.execute("drop table if exists tmp_od_pairs;")
 
+        logger.debug("drop the tmp_rmp_to_endcap table")
+        db_cur.execute("drop table if exists tmp_rmp_to_endcap")
+
+        logger.debug("drop the tmp_endcap_to_other table")
+        db_cur.execute("drop table if exists tmp_endcap_to_other")
+
         logger.info("end: create o-d pairs table")
 
         # Fetch all od-pairs, ordered by target
@@ -634,14 +751,38 @@ def make_od_pairs(the_scenario, logger):
         '''
         sql_list = db_cur.execute(sql).fetchall()
 
-        # Loop through the od_pairs
-        od_pairs = {}
-        for row in sql_list:
-            target = row[0]
-            source = row[1]
-            scenario_rt_id = row[2]
-            commodity_id = row[3]
-            phase_of_matter = row[4]
+        # If commodity associated with a max transport distance, always make shortest paths from source
+        sql = '''
+        select odp.commodity_id, count(distinct odp.from_node_id), count(distinct odp.to_node_id), cm.max_transport_distance
+        from od_pairs odp
+        join commodities cm on odp.commodity_id = cm.commodity_id
+        group by odp.commodity_id
+        '''
+        od_count = db_cur.execute(sql)
+
+        commodity_st = {}
+        for row in od_count:
+            commodity_id = row[0]
+            count_source = row[1]
+            count_target = row[2]
+            MTD = row[3]
+            if MTD is not None:
+                commodity_st[commodity_id] = 'source'
+            else:
+                if count_source < count_target:
+                    commodity_st[commodity_id] = 'source'
+                else:
+                    commodity_st[commodity_id] = 'target'
+
+    # Loop through the od_pairs
+    od_pairs = {}
+    for row in sql_list:
+        target = row[0]
+        source = row[1]
+        scenario_rt_id = row[2]
+        commodity_id = row[3]
+        phase_of_matter = row[4]
+        if commodity_st[commodity_id] == 'target':
             if commodity_id not in od_pairs:
                 od_pairs[commodity_id] = {}
                 od_pairs[commodity_id]['phase_of_matter'] = phase_of_matter
@@ -651,13 +792,256 @@ def make_od_pairs(the_scenario, logger):
             if source not in od_pairs[commodity_id]['targets'][target].keys():
                 od_pairs[commodity_id]['targets'][target][source] = []
             od_pairs[commodity_id]['targets'][target][source].append(scenario_rt_id)
-
+        else:
+            if commodity_id not in od_pairs:
+                od_pairs[commodity_id] = {}
+                od_pairs[commodity_id]['phase_of_matter'] = phase_of_matter
+                od_pairs[commodity_id]['sources'] = {}
+            if source not in od_pairs[commodity_id]['sources'].keys():
+                od_pairs[commodity_id]['sources'][source] = {}
+            if target not in od_pairs[commodity_id]['sources'][source].keys():
+                od_pairs[commodity_id]['sources'][source][target] = []
+            od_pairs[commodity_id]['sources'][source][target].append(scenario_rt_id)
+    
     return od_pairs
 
 
 # -----------------------------------------------------------------------------
 
 
+def make_max_transport_distance_subgraphs(the_scenario, logger, commodity_subgraph_dict):
+    
+    # Get facility, commodity, and MTD info from the db
+    logger.info("start: pull facility/commodity MTD from SQL")
+    with sqlite3.connect(the_scenario.main_db) as db_cur:
+        sql = """select cm.commodity_id, networkx_nodes.node_id, c.max_transport_distance, f.facility_type_id, fc.io
+               from commodity_mode cm
+               join facility_commodities fc on cm.commodity_id = fc.commodity_id 
+               join facilities f on fc.facility_id = f.facility_id  --likely need facility's node id from od_pairs instead
+               join networkx_nodes on networkx_nodes.location_id = f.location_id
+               join commodities c on cm.commodity_id = c.commodity_id
+               where cm.allowed_yn like 'Y' and fc.io = 'o' and networkx_nodes.location_1 like '%_out'
+			   group by networkx_nodes.node_id;"""
+        commodity_mode_data = db_cur.execute(sql).fetchall()
+    logger.info("end: pull commodity mode from SQL")
+
+    # Find which facilities have max transport distance defined
+    mtd_count = 0
+    for row in commodity_mode_data:
+        commodity_id = int(row[0])
+        facility_node_id = int(row[1])
+        commodity_mtd = row[2]
+        if commodity_mtd is not None:
+            mtd_count = mtd_count + 1
+            commodity_subgraph_dict[commodity_id]['MTD'] = commodity_mtd
+            if 'facilities' not in commodity_subgraph_dict[commodity_id]:
+                commodity_subgraph_dict[commodity_id]['facilities'] = []
+            commodity_subgraph_dict[commodity_id]['facilities'].append(facility_node_id)
+
+    if mtd_count == 0:
+        return commodity_subgraph_dict
+    
+    # Store nodes that can be reached from an RMP with MTD
+    ends = {}
+    for commodity_id, commodity_dict in commodity_subgraph_dict.items():
+        # For facility, MTD in commodities_with_mtd[commodity_id]
+        if 'MTD' in commodity_dict:
+            MTD = commodity_dict['MTD']
+            for facility_node_id in commodity_dict['facilities']:
+                # If 'facility_subgraphs' dictionary for commodity_subgraph_dict[commodity_id] doesn't exist, add it
+                if 'facility_subgraphs' not in commodity_subgraph_dict[commodity_id]:
+                    commodity_subgraph_dict[commodity_id]['facility_subgraphs'] = {}
+                
+                # Use shortest path to find the nodes that are within MTD
+                logger.info("start: dijkstra for facility node ID " + str(facility_node_id))
+                G = commodity_subgraph_dict[commodity_id]['subgraph']
+                # distances: key = node ids within cutoff, value = length of paths
+                # endcaps: key = node ids within cutoff, value = list of nodes labeled as endcaps
+                # modeled after NetworkX single_source_dijkstra:
+                # distances, paths = nx.single_source_dijkstra(G, facility_node_id, cutoff = MTD, weight = 'MILES')
+                fn_miles = lambda u, v, d: min(attr.get('MILES', 1) for attr in d.values())
+                distances, endcaps = dijkstra(G, facility_node_id, fn_miles, cutoff=MTD)
+
+                logger.info("start: distances/paths for facility node ID " + str(facility_node_id))
+                
+                # Creates a subgraph of the nodes and edges that are reachable from the facility
+                commodity_subgraph_dict[commodity_id]['facility_subgraphs'][facility_node_id] = G.subgraph(distances.keys()).copy()
+
+                # If in G1 step for candidate generation, find endcaps
+                if not os.path.exists(the_scenario.processor_candidates_commodity_data) and the_scenario.processors_candidate_slate_data != 'None':
+                    ends[facility_node_id] = {}
+                    ends[facility_node_id]['ends'] = endcaps
+                    ends[facility_node_id]['commodity_id'] = commodity_id     
+
+                logger.info("end: distances/paths for facility node ID " + str(facility_node_id))
+
+    # If in G1 step for candidate generation, add endcaps to endcap_nodes table
+    if not os.path.exists(the_scenario.processor_candidates_commodity_data) and the_scenario.processors_candidate_slate_data != 'None':                 
+        with sqlite3.connect(the_scenario.main_db) as db_cur:
+            # Create temp table to hold endcap nodes
+            sql = """
+                drop table if exists tmp_endcap_nodes;"""
+            db_cur.execute(sql)
+
+            sql = """
+                create table tmp_endcap_nodes(
+                node_id integer primary key,
+                location_id default null,
+                mode_source default null,
+                source_node_id integer,
+                commodity_id integer,
+                destination_yn default null)
+            ;"""
+            db_cur.execute(sql)
+
+            for facility_node_id in ends:
+                logger.info("Updating endcap_nodes for facility node ID " + str(facility_node_id))
+                commodity_id = ends[facility_node_id]['commodity_id']
+                end_list = ends[facility_node_id]['ends']
+                for i in range(len(end_list)):
+                    sql = """insert or replace into tmp_endcap_nodes (node_id, location_id, mode_source, source_node_id, commodity_id, destination_yn)
+                             values({}, NULL, NULL, {}, {}, NULL);""".format(end_list[i], facility_node_id, commodity_id)
+                    db_cur.execute(sql)
+                    db_cur.commit()
+            
+            sql = """alter table tmp_endcap_nodes
+                     add column source_facility_id;"""
+            db_cur.execute(sql)
+
+            sql = """update tmp_endcap_nodes
+                     set source_facility_id = (select temp.source_facility_id
+                                               from (select en.source_node_id, f.facility_id as source_facility_id
+									                 from tmp_endcap_nodes en
+									                 join networkx_nodes n on en.source_node_id = n.node_id
+                                                     join facilities f on f.location_id = n.location_id) temp
+									           where tmp_endcap_nodes.source_node_id = temp.source_node_id)"""
+            db_cur.execute(sql)
+            
+            sql = """drop table if exists endcap_nodes;"""
+            db_cur.execute(sql)
+
+            sql = """
+                create table endcap_nodes(
+                node_id integer primary key,
+                location_id default null,
+                mode_source default null,
+                source_node_id integer,
+                source_facility_id integer,
+                commodity_id integer,
+                process_id integer,
+                shape_x integer,
+                shape_y integer,
+                destination_yn default null)
+            ;"""
+            db_cur.execute(sql)
+
+            # Join temp endcap nodes with process id and edge shapes and insert into endcap nodes table
+            sql = """insert into endcap_nodes (node_id, location_id, mode_source, source_node_id, source_facility_id, commodity_id, process_id, shape_x, shape_y, destination_yn)
+                        SELECT
+                        tmp.node_id,
+                        tmp.location_id,
+                        tmp.mode_source,
+                        tmp.source_node_id,
+                        tmp.source_facility_id,
+                        tmp.commodity_id,
+                        cpc.process_id,
+                        nwx.shape_x,
+                        nwx.shape_y,
+                        tmp.destination_yn
+                        from tmp_endcap_nodes tmp
+                        join candidate_process_commodities cpc on tmp.commodity_id = cpc.commodity_id
+                        join networkx_nodes nwx on tmp.node_id = nwx.node_id"""
+            db_cur.execute(sql)
+
+            sql = "drop table if exists tmp_endcap_nodes;"
+            db_cur.execute(sql)
+
+    return commodity_subgraph_dict
+
+
+# -----------------------------------------------------------------------------
+
+
+def dijkstra(G, source, get_weight, pred=None, paths=None, cutoff=None,
+             target=None):
+    """Implementation of Dijkstra's algorithm
+    Parameters
+    ----------
+    G : NetworkX graph
+    source : node label
+        Starting node for path.
+    get_weight : function
+        Function for getting edge weight.
+    pred : list, optional (default=None)
+        List of predecessors of a node.
+    paths : dict, optional (default=None)
+        Path from the source to a target node.
+    cutoff : integer or float, optional (default=None)
+        Depth to stop the search. Only paths of length <= cutoff are returned.
+    target : node label, optional (default=None)
+        Ending node for path.
+    Returns
+    -------
+    distance, endcaps : dictionaries
+        Returns a tuple of two dictionaries keyed by node.
+        The first dictionary stores distance from the source.
+        The second stores all endcap nodes for that node.
+    """
+    G_succ = G.succ if G.is_directed() else G.adj
+
+    push = heappush
+    pop = heappop
+    dist = {}  # dictionary of final distances
+    seen = {source: 0}
+    c = count()
+    fringe = []  # use heapq with (distance,label) tuples
+    endcaps = []  # MARK MOD
+    push(fringe, (0, next(c), source))
+
+    while fringe:
+        (d, _, v) = pop(fringe)
+        if v in dist:
+            continue  # already searched this node.
+        dist[v] = d
+        if v == target:
+            break
+        added = 0  # MARK MOD
+        too_far = 0  # MARK MOD
+        for u, e in G_succ[v].items():
+            cost = get_weight(v, u, e)
+            if cost is None:
+                continue
+            vu_dist = dist[v] + get_weight(v, u, e)
+            if cutoff is not None:
+                if vu_dist > cutoff:
+                    too_far += 1  # MARK MOD
+                    # endcaps.append(v)  # MARK MOD
+                    continue
+            if u in dist:
+                if vu_dist < dist[u]:
+                    raise ValueError('Contradictory paths found:',
+                                     'negative weights?')
+            elif u not in seen or vu_dist < seen[u]:
+                added += 1  # MARK MOD
+                seen[u] = vu_dist
+                push(fringe, (vu_dist, next(c), u))
+                if paths is not None:
+                    paths[u] = paths[v] + [u]
+                if pred is not None:
+                    pred[u] = [v]
+            elif vu_dist == seen[u]:
+                if pred is not None:
+                    pred[u].append(v)
+        
+        if added == 0 and too_far >= 1:  # MARK MOD
+            endcaps.append(v)  # MARK MOD
+
+    return (dist, endcaps)  # MARK MOD
+
+
+# -----------------------------------------------------------------------------
+
+
 def export_fcs_from_main_gdb(the_scenario, logger):
     # export fcs from the main.GDB to individual shapefiles
     logger.info("start: export_fcs_from_main_gdb")
@@ -803,7 +1187,7 @@ def clean_networkx_graph(the_scenario, G, logger):
                         converted_pipeline_tariff_cost = "usd/{}".format(the_scenario.default_units_liquid_phase)
                         route_cost_scaling = Q_(pipeline_tariff_cost).to(converted_pipeline_tariff_cost).magnitude
                     else:
-                        # Below was hardcoded conversion of US cents per barrel to dollars per kgal. Updated code uses pint and we now provide the base_rates in dollars like other costs
+                        # Below was hardcoded conversion of US cents per barrel to dollars per thousand_gallon. Updated code uses pint and we now provide the base_rates in dollars like other costs
                         # this if statement maintains compatibility with older network which is still in cents.
                         route_cost_scaling = (((float(G.edges[u, v, keys]['base_rate']) / 100) / 42.0) * 1000.0)
 
diff --git a/program/ftot_postprocess.py b/program/ftot_postprocess.py
index 500042e..7c18337 100644
--- a/program/ftot_postprocess.py
+++ b/program/ftot_postprocess.py
@@ -1030,6 +1030,45 @@ def make_optimal_scenario_results_db(the_scenario, logger):
                             order by facility_name
                             ;"""
 
+        # Processor capacity report (input)
+        sql_proc_input_capacity = """insert into optimal_scenario_results
+                                     select
+                                     "facility_summary",
+                                     c.commodity_name,
+                                     f.facility_name,
+                                     "processor_input_capacity",
+                                     "total",
+                                     IFNULL(fc.scaled_quantity, "Unconstrained"),
+                                     fc.units,
+                                     ''
+                                     from facility_commodities fc
+                                     join facilities f on f.facility_id = fc.facility_id
+                                     join commodities c on c.commodity_id = fc.commodity_id
+                                     join facility_type_id fti on fti.facility_type_id  = f.facility_type_id
+                                     where fti.facility_type = 'processor' and fc.io = 'i'
+                                     order by facility_name
+                                     ;"""
+
+        # Processor capacity report (output)
+        sql_proc_output_capacity = """insert into optimal_scenario_results
+                                      select
+                                      "facility_summary",
+                                      c.commodity_name,
+                                      f.facility_name,
+                                      "processor_output_capacity",
+                                      "total",
+                                      IFNULL(fc.scaled_quantity, "Unconstrained"),
+                                      fc.units,
+                                      ''
+                                      from facility_commodities fc
+                                      join facilities f on f.facility_id = fc.facility_id
+                                      join commodities c on c.commodity_id = fc.commodity_id
+                                      join facility_type_id fti on fti.facility_type_id  = f.facility_type_id
+                                      where fti.facility_type = 'processor' and fc.io = 'o'
+                                      order by facility_name
+                                      ;"""
+        
+
         # initialize SQL queries that vary based on whether routes are used or not
         logger.debug("start: summarize optimal supply and demand")
         if the_scenario.ndrOn:
@@ -1144,7 +1183,7 @@ def make_optimal_scenario_results_db(the_scenario, logger):
                                       where ov.o_facility > 0 and ov.edge_type = 'transport' and fti.facility_type = 'processor'
                                       group by ov.commodity_name, ov.o_facility, ne.mode_source, ov.units
                                       ;"""
-
+            
             sql_processor_input = """insert into optimal_scenario_results
                                      select
                                      "facility_summary",
@@ -1323,26 +1362,65 @@ def make_optimal_scenario_results_db(the_scenario, logger):
         db_con.execute(sql_destination_demand_optimal)
         db_con.execute(sql_destination_demand_optimal_frac)
         db_con.execute(sql_rmp_supply)
+        db_con.execute(sql_proc_input_capacity)
+        db_con.execute(sql_proc_output_capacity)
+
         db_con.execute(sql_rmp_supply_optimal)
         db_con.execute(sql_rmp_supply_optimal_frac)
         db_con.execute(sql_processor_output)
         db_con.execute(sql_processor_input)
         db_con.execute(sql_processor_build_cost)
 
-        # measure totals
-        sql_total = """insert into optimal_scenario_results
+        # totals across modes
+        sql_allmodes = """insert into optimal_scenario_results
                        select table_name, commodity, facility_name, measure, "allmodes" as mode, sum(value), units, notes
                        from optimal_scenario_results
+                       where mode != "total"
                        group by table_name, commodity, facility_name, measure
                        ;"""
-        db_con.execute(sql_total)
+        db_con.execute(sql_allmodes)
 
-        # delete redundant total rows
-        sql_delete = """delete from optimal_scenario_results
-                        where (measure in ("rmp_supply_potential","destination_demand_potential","processor_amortized_build_cost"))
-                        and mode = "allmodes"
-                        ;"""
-        db_con.execute(sql_delete)
+        # processor optimal fractions, as share of input and output capacity
+        sql_proc_input_optimal_frac = """insert into optimal_scenario_results
+                                             select
+                                             "facility_summary",
+                                             ov.commodity_name,
+                                             ov.d_facility,
+                                             "processor_input_optimal_frac",
+                                             "allmodes",
+                                             IFNULL((sum(ov.variable_value) / fc.scaled_quantity), "Unconstrained"),
+                                             "fraction",
+                                             ''
+                                             from optimal_variables ov
+                                             join facilities f on f.facility_name = ov.d_facility
+                                             join facility_commodities fc on fc.facility_id = f.facility_id and fc.commodity_id = ov.commodity_id
+                                             join commodities c on c.commodity_id = fc.commodity_id
+                                             join facility_type_id fti on fti.facility_type_id  = f.facility_type_id
+                                             where ov.d_facility > 0 and ov.edge_type = 'transport' and fti.facility_type = 'processor' and fc.io = 'i'
+                                             group by ov.commodity_name, ov.d_facility
+                                             ;"""
+        db_con.execute(sql_proc_input_optimal_frac)
+
+        sql_proc_output_optimal_frac = """insert into optimal_scenario_results
+                                             select
+                                             "facility_summary",
+                                             ov.commodity_name,
+                                             ov.o_facility,
+                                             "processor_output_optimal_frac",
+                                             "allmodes",
+                                             IFNULL((sum(ov.variable_value) / fc.scaled_quantity), "Unconstrained"),
+                                             "fraction",
+                                             ''
+                                             from optimal_variables ov
+                                             join facilities f on f.facility_name = ov.o_facility
+                                             join facility_commodities fc on fc.facility_id = f.facility_id and fc.commodity_id = ov.commodity_id
+                                             join commodities c on c.commodity_id = fc.commodity_id
+                                             join facility_type_id fti on fti.facility_type_id  = f.facility_type_id
+                                             where ov.o_facility > 0 and ov.edge_type = 'transport' and fti.facility_type = 'processor' and fc.io = 'o'
+                                             group by ov.commodity_name, ov.o_facility
+                                             ;"""
+        
+        db_con.execute(sql_proc_output_optimal_frac)
 
         # manually set unit for "vehicles" for "allmodes"
         sql_vehicles = """update optimal_scenario_results
@@ -1417,12 +1495,14 @@ def generate_scenario_summary(the_scenario, logger):
             logger.result('{}_{}_{}_{}: \t {:,.2f} : \t {}'.format(row[0].upper(), row[1].upper(), row[3].upper(), row[4].upper(), row[5], row[6]))
 
         # Log facility results
-        sql = "select * from optimal_scenario_results where table_name = 'facility_summary' order by commodity, measure, mode;"
+        sql = "select * from optimal_scenario_results where table_name = 'facility_summary' order by facility_name, measure, commodity, mode;"
         db_cur = db_con.execute(sql)
         data = db_cur.fetchall()
         for row in data:
             if row[3] == 'processor_amortized_build_cost':
                 logger.result('{}_{}_{}_{}: \t {:,.2f} : \t {}'.format(row[0].upper(), row[2].upper(), row[3].upper(), row[4].upper(), row[5], row[6]))
+            elif type(row[5]) == str:
+                logger.result('{}_{}_{}_{}_{}: \t {} : \t {}'.format(row[0].upper(), row[2].upper(), row[3].upper(), row[1].upper(), row[4].upper(), row[5], row[6]))
             else:
                 logger.result('{}_{}_{}_{}_{}: \t {:,.2f} : \t {}'.format(row[0].upper(), row[2].upper(), row[3].upper(), row[1].upper(), row[4].upper(), row[5], row[6]))
 
@@ -1614,7 +1694,7 @@ def db_report_commodity_utilization(the_scenario, logger):
                                from optimal_scenario_results osr
                                where osr.commodity = c.commodity_name
                                and osr.measure = 'total_flow'
-                               and osr.mode = 'allmodes'), 2), 'No Flow') optimal_flow,
+                               and osr.mode = 'allmodes'), 2), "No Flow") as optimal_flow,
                  fc.units
                  from facility_commodities fc
                  join commodities c on fc.commodity_id = c.commodity_id
@@ -1636,7 +1716,7 @@ def db_report_commodity_utilization(the_scenario, logger):
             if (type(row[3]) == str):
                 logger.result("{:15.15} {:15.15} {:4.1} {:15.15} {:15.15}".format(row[0], row[1], row[2], row[3], row[4]))
             else:
-                logger.result("{:15.15} {:15.15} {:4.1} {:15,.2f} {:15.15}".format(row[0], row[1], row[2], row[3], row[4]))
+                logger.result("{:15.15} {:15.15} {:4.1} {:15,.2f} {:15.15}".format(row[0], row[1], row[2], row[3], row[4]))                
         logger.result("-------------------------------------------------------------------")
 
     # This query compares the optimal flow to the net available quantity of material in the scenario.
@@ -1645,11 +1725,16 @@ def db_report_commodity_utilization(the_scenario, logger):
     # -----------------------------------
     with sqlite3.connect(the_scenario.main_db) as db_con:
         sql = """select c.commodity_name, fti.facility_type, fc.io,
-                 IFNULL(round((select osr.value
-                               from optimal_scenario_results osr
-                               where osr.commodity = c.commodity_name
-                               and osr.measure = 'total_flow'
-                               and osr.mode = 'allmodes') / sum(fc.scaled_quantity), 2), 'Unconstrained') as Utilization,
+                 (CASE WHEN (select osr.value
+                            from optimal_scenario_results osr
+                            where osr.commodity = c.commodity_name and osr.measure = 'total_flow'
+                            and osr.mode = 'allmodes') IS NULL THEN "No Flow"
+					   WHEN sum(CASE WHEN fc.scaled_quantity IS NULL THEN 1 ELSE 0 END) <> 0 THEN "Unconstrained"
+                       ELSE round((select osr.value
+                            from optimal_scenario_results osr
+                            where osr.commodity = c.commodity_name
+                            and osr.measure = 'total_flow'
+                            and osr.mode = 'allmodes') / sum(fc.scaled_quantity), 2) END) as Utilization,
                  'fraction'
                  from facility_commodities fc
                  join commodities c on fc.commodity_id = c.commodity_id
@@ -1660,6 +1745,8 @@ def db_report_commodity_utilization(the_scenario, logger):
                  order by c.commodity_name, fc.io asc;"""
         db_cur = db_con.execute(sql)
 
+        # TODO: Change sum of unconstrained and constrained to unconstrained
+
         db_data = db_cur.fetchall()
         logger.result("-------------------------------------------------------------------")
         logger.result("Scenario Total Utilization of Supply and Demand")
@@ -1678,8 +1765,7 @@ def db_report_commodity_utilization(the_scenario, logger):
 # ===================================================================================================
 
 
-def make_optimal_route_segments_featureclass_from_db(the_scenario,
-                                                     logger):
+def make_optimal_route_segments_featureclass_from_db(the_scenario, logger):
 
     logger.info("starting make_optimal_route_segments_featureclass_from_db")
 
diff --git a/program/ftot_processor.py b/program/ftot_processor.py
index c3daa7c..b313485 100644
--- a/program/ftot_processor.py
+++ b/program/ftot_processor.py
@@ -208,6 +208,8 @@ def populate_candidate_process_list_table(the_scenario, candidate_process_list,
                 logger.warning("the units for the max_size and min_size candidate process do not match!")
             if max_size_units != min_aggregation_units:
                 logger.warning("the units for the max_size and min_aggregation candidate process do not match!")
+            if max_size_units != cost_formula_units.split(" / ")[-1]:
+                logger.warning("the units for the max_size and cost formula candidate process do not match!")
             if min_size == '':
                 logger.warning("the min_size is set to Null")
             if max_size == '':
@@ -434,8 +436,8 @@ def processor_candidates(the_scenario, logger):
                     cpl.maxsize max_input,
                     cpl.minsize min_input
                     from candidate_nodes cn
-                    join candidate_process_commodities cpc on cpc.commodity_id = cn.commodity_id and cpc.process_id = cn.process_id
-                    join candidate_process_list cpl on cpl.process_id = cn.process_id
+                    join candidate_process_commodities cpc on cpc.commodity_id = cn.o_commodity_id and cpc.process_id = cn.o_process_id
+                    join candidate_process_list cpl on cpl.process_id = cn.o_process_id
                     join networkx_nodes xy on cn.node_id = xy.node_id
                     join commodities c on c.commodity_id = cpc.commodity_id
                     join schedule_names sn on sn.schedule_id = cpl.schedule_id
diff --git a/program/ftot_pulp.py b/program/ftot_pulp.py
index b95695b..8e05945 100644
--- a/program/ftot_pulp.py
+++ b/program/ftot_pulp.py
@@ -1,16 +1,15 @@
 # ---------------------------------------------------------------------------------------------------
 # Name: ftot_pulp
 #
-# Purpose: PulP optimization - create and run a modified facility location problem.
+# Purpose: PuLP optimization - create and run a modified facility location problem.
 # Take NetworkX and GIS scenario data as recorded in main.db and convert to a structure of edges, nodes, vertices.
 # Create variables for flow over edges, unmet demand, processor use, and candidate processors to build if present
 # Solve cost minimization for unmet demand, transportation, and facility build costs
-# Constraints ensure compliance with scenario requirements (e.g. max_route_capacity)
-# as well as general problem structure (e.g. conservation_of_flow)
+# Constraints ensure compliance with scenario requirements (e.g., max_route_capacity)
+# as well as general problem structure (e.g., conservation_of_flow)
 # ---------------------------------------------------------------------------------------------------
 
 import datetime
-import pdb
 import re
 import sqlite3
 from collections import defaultdict
@@ -39,7 +38,7 @@
 default_max_capacity = 10000000000
 default_min_capacity = 0
 
-zero_threshold=0.00001
+zero_threshold = 0.00001
 
 
 def o1(the_scenario, logger):
@@ -1113,7 +1112,6 @@ def generate_first_edges_from_source_facilities(the_scenario, schedule_length, l
                 for tariff_row in db_cur.execute(sql):
                     tariff_id = tariff_row[0]
 
-
             if mode in the_scenario.permittedModes and (mode, commodity_id) in commodity_mode_dict.keys() \
                     and commodity_mode_dict[mode, commodity_id] == 'Y':
 
@@ -1138,7 +1136,7 @@ def generate_first_edges_from_source_facilities(the_scenario, schedule_length, l
                                 # than checking if facility type is 'ultimate destination'
                                 # only connect to vertices with matching source_facility_id
                                 # source_facility_id is zero for commodities without source tracking
-
+                                
                                 main_db_con.execute("""insert or ignore into edges (from_node_id, to_node_id,
                                     start_day, end_day, commodity_id,
                                     o_vertex_id,
@@ -1197,6 +1195,7 @@ def generate_first_edges_from_source_facilities(the_scenario, schedule_length, l
         for row in db_cur.execute("select count(distinct edge_id) from edges where edge_type = 'transport';"):
             transport_edges_created = row[0]
             logger.info('{} transport edges created'.format(transport_edges_created))
+
         for row in db_cur.execute("""select count(distinct edge_id), children_created from edges
          where edge_type = 'transport'
          group by children_created;"""):
@@ -1243,7 +1242,6 @@ def generate_all_edges_from_source_facilities(the_scenario, schedule_length, log
             transport_edges_created = row_d[0]
             nx_edge_count = row_d[1]
 
-
         commodity_mode_data = main_db_con.execute("select * from commodity_mode;")
         commodity_mode_data = commodity_mode_data.fetchall()
         commodity_mode_dict = {}
@@ -1629,7 +1627,7 @@ def generate_all_edges_from_source_facilities(the_scenario, schedule_length, log
                     '{} transport edges on {} nx edges,  created in {} loops, {} edges_requiring_children'.format(
                         transport_edges_created, nx_edge_count, while_count, edges_requiring_children))
 
-        #  edges going in to the facility by re-running "generate first edges
+        # edges going in to the facility by re-running "generate first edges
         # then re-run this method
 
         logger.info('{} transport edges on {} nx edges,  created in {} loops, {} edges_requiring_children'.format(
@@ -1952,7 +1950,7 @@ def generate_edges_from_routes(the_scenario, schedule_length, logger):
         join route_edges re2 on r1.scenario_rt_id = re2.scenario_rt_id and r1.num_edges = re2.rt_order_ind) r2
         where r2.scenario_rt_id = odp.scenario_rt_id and r2.phase_of_matter = odp.phase_of_matter
         ;""")
-
+        
         route_data = main_db_con.execute("select * from route_reference where route_type = 'transport';")
 
         # Add an edge for each route, (applicable) vertex, day, commodity
@@ -1966,39 +1964,78 @@ def generate_edges_from_routes(the_scenario, schedule_length, logger):
             phase_of_matter = row_a[11]
             cost = row_a[12]
             miles = row_a[13]
-            source_facility_id = 0
+            source_facility_id = 0  # reset to default - will be updated if a source facility is relevant
+
+            # get vertices for from facility, if applicable
+            db_cur4 = main_db_con.cursor()
+            if from_location != None:
+                from_vertices = db_cur4.execute("""select vertex_id, schedule_day, source_facility_id
+                            from vertices v, facility_commodities fc
+                            where v.location_id = {}
+                            and v.commodity_id = {} 
+                            and v.storage_vertex = 1
+                            and v.facility_id = fc.facility_id
+                            and v.commodity_id = fc.commodity_id
+                            and fc.io = 'o'""".format(from_location, commodity_id))
+            else:
+                from_vertices = {}
+            
+            # save vertices for from_location in a dict
+            from_location_vertices = {}
+            from_count = 0
+            for row in from_vertices:
+                from_count = from_count + 1
+                vertex = row[0]
+                day = row[1]
+                source_facility_id = row[2]
+                from_location_vertices[day] = vertex
+
+            # get vertices for to facility, if applicable
+            db_cur5 = main_db_con.cursor()
+            if to_location != None:
+                to_vertices = db_cur5.execute("""select vertex_id, schedule_day, source_facility_id
+                            from vertices v, facility_commodities fc
+                            where v.location_id = {}
+                            and v.commodity_id = {}
+                            and v.storage_vertex = 1
+                            and v.facility_id = fc.facility_id
+                            and v.commodity_id = fc.commodity_id
+                            and fc.io = 'i'""".format(to_location, commodity_id))
+            else:
+                to_vertices = {}
+
+            # save vertices for to_location in a dict  
+            to_location_vertices = {}
+            to_count = 0
+            for row in to_vertices:
+                to_count = to_count + 1
+                vertex = row[0]
+                day = row[1]
+                # source_facility_id = row[2]
+                to_location_vertices[day] = vertex          
 
             for day in range(1, schedule_length+1):
                 if day + fixed_route_duration <= schedule_length:
                     # add edge from o_vertex to d_vertex
                     # for each day and commodity, get the corresponding origin and destination vertex
                     # ids to include with the edge info
-                    db_cur4 = main_db_con.cursor()
-                    # TODO: are these DB calls necessary for vertices?
-                    for row_d in db_cur4.execute("""select vertex_id
-                        from vertices v, facility_commodities fc
-                        where v.location_id = {} and v.schedule_day = {}
-                        and v.commodity_id = {} and v.source_facility_id = {}
-                        and v.storage_vertex = 1
-                        and v.facility_id = fc.facility_id
-                        and v.commodity_id = fc.commodity_id
-                        and fc.io = 'o'""".format(from_location, day, commodity_id, source_facility_id)):
-                        from_vertex_id = row_d[0]
-                        db_cur5 = main_db_con.cursor()
-                        for row_e in db_cur5.execute("""select vertex_id
-                        from vertices v, facility_commodities fc
-                        where v.location_id = {} and v.schedule_day = {}
-                        and v.commodity_id = {} and v.source_facility_id = {}
-                        and v.storage_vertex = 1
-                        and v.facility_id = fc.facility_id
-                        and v.commodity_id = fc.commodity_id
-                        and fc.io = 'i'""".format(to_location, day, commodity_id, source_facility_id)):
-                            to_vertex_id = row_e[0]
-                            main_db_con.execute("""insert or ignore into edges (route_id, from_node_id, 
+                    if from_count > 0:
+                        from_vertex_id = from_location_vertices[day]
+                    else:
+                        from_vertex_id = 'Null'
+                    
+                    if to_count > 0:
+                        to_vertex_id = to_location_vertices[day]
+                    else:
+                        to_vertex_id = 'Null'
+
+                    # set vertices for from and to location = null
+                    # get vertex for to_facility and from_facility if available
+                    main_db_con.execute("""insert or ignore into edges (route_id, from_node_id, 
                                 to_node_id, start_day, end_day, commodity_id, o_vertex_id, d_vertex_id, 
                                 edge_flow_cost, edge_type, 
-                                miles,phase_of_matter) VALUES ({}, {}, {}, {}, {}, {}, {}, {}, {}, 
-                                '{}', {},'{}'
+                                miles,phase_of_matter, source_facility_id) VALUES ({}, {}, {}, {}, {}, {}, {}, {}, {}, 
+                                '{}', {},'{}', {}
                                 )""".format(route_id,
                                             from_node_id,
                                             to_node_id,
@@ -2011,9 +2048,8 @@ def generate_edges_from_routes(the_scenario, schedule_length, logger):
                                             # nx_edge_id, mode, mode_oid,
                                             miles,
                                             # simple_mode, tariff_id,
-                                            phase_of_matter))
-                                            #source_facility_id))
-
+                                            phase_of_matter,
+                                            source_facility_id))
     return
 
 
@@ -2034,7 +2070,6 @@ def set_edges_volume_capacity(the_scenario, logger):
         logger.debug("volume and capacity recorded for non-pipeline edges")
 
         logger.debug("starting to record volume and capacity for pipeline edges")
-        ##
         main_db_con.executescript("""update edges set volume =
         (select l.background_flow
          from pipeline_mapping pm,
@@ -2073,6 +2108,7 @@ def set_edges_volume_capacity(the_scenario, logger):
         where simple_mode = 'pipeline'
         ;""")
         logger.debug("volume and capacity recorded for pipeline edges")
+
         logger.debug("starting to record units and conversion multiplier")
         main_db_con.execute("""update edges
         set capacity_units =
@@ -2100,6 +2136,7 @@ def set_edges_volume_capacity(the_scenario, logger):
                     the_scenario.barge_load_liquid.magnitude,
                     the_scenario.barge_load_solid.magnitude,
                     ))
+
         logger.debug("units and conversion multiplier recorded for all edges; starting capacity minus volume")
         main_db_con.execute("""update edges
         set capac_minus_volume_zero_floor =
@@ -2176,11 +2213,10 @@ def create_flow_vars(the_scenario, logger):
             counter += 1
             # create an edge for each commodity allowed on this link - this construction may change
             # as specific commodity restrictions are added
-            # TODO4-18 add days, but have no scheduel for links currently
+            # TODO add days, but have no schedule for links currently
             # running just with nodes for now, will add proper facility info and storage back soon
             edge_list.append((row[0]))
 
-
     flow_var = LpVariable.dicts("Edge", edge_list, 0, None)
     return flow_var
 
@@ -2196,7 +2232,7 @@ def create_unmet_demand_vars(the_scenario, logger):
     with sqlite3.connect(the_scenario.main_db) as main_db_con:
         db_cur = main_db_con.cursor()
         for row in db_cur.execute("""select v.facility_id, v.schedule_day, 
-        ifnull(c.supertype, c.commodity_name) top_level_commodity_name, v.udp
+         ifnull(c.supertype, c.commodity_name) top_level_commodity_name, v.udp
          from vertices v, commodities c, facility_type_id ft, facilities f
          where v.commodity_id = c.commodity_id
          and ft.facility_type = "ultimate_destination"
@@ -2343,11 +2379,11 @@ def create_opt_problem(logger, the_scenario, unmet_demand_vars, flow_vars, proce
         logger.debug("number of candidate processors from proc.csv = {}".format(input_candidates))
 
         if input_candidates > 0:
-            #force ignore_facility = 'false' for processors input from file; it should always be set to false anyway
+            # force ignore_facility = 'false' for processors input from file; it should always be set to false anyway
             input_processor_build_cost = db_cur.execute("""
             select f.facility_id,  f.build_cost
             from facilities f
-                INNER JOIN facility_type_id ft ON f.facility_type_id = ft.facility_type_id
+            inner join facility_type_id ft ON f.facility_type_id = ft.facility_type_id
             where candidate = 1 and build_cost>0
             and facility_type = 'processor'
             and ifnull(ignore_facility, 'false') = 'false'
@@ -2437,8 +2473,6 @@ def create_constraint_unmet_demand(logger, the_scenario, prob, flow_var, unmet_d
                                                                                                               key[1],
                                                                                                               key[2])
             else:
-                if key not in actual_demand_dict:
-                    pdb.set_trace()
                 # no edges in, so unmet demand equals full demand
                 prob += actual_demand_dict[key] == unmet_demand_var[
                     key], "constraint set unmet demand variable for facility {}, day {}, " \
@@ -2494,7 +2528,7 @@ def create_constraint_max_flow_out_of_supply_vertex(logger, the_scenario, prob,
 
 def create_constraint_daily_processor_capacity(logger, the_scenario, prob, flow_var, processor_build_vars,
                                                processor_daily_flow_vars):
-    logger.debug("STARTING:  create_constraint_daily_processor_capacity")
+    logger.debug("STARTING: create_constraint_daily_processor_capacity")
     # primary vertices only
     # flow into vertex is capped at facility max_capacity per day
     # sum over all input commodities, grouped by day and facility
@@ -2551,7 +2585,6 @@ def create_constraint_daily_processor_capacity(logger, the_scenario, prob, flow_
                 logger.debug(
                     "flow in for capacity constraint on processor facility {} day {}: {}".format(facility_id, day, flow_in))
 
-
                 # capacity is set to -1 if there is no restriction, so should be no constraint
                 if min_capacity >= 0:
                     daily_inflow_min_capacity = float(min_capacity) * float(daily_activity_level)
@@ -2570,8 +2603,6 @@ def create_constraint_daily_processor_capacity(logger, the_scenario, prob, flow_
 
                 prob += lpSum(flow_in) >= daily_inflow_min_capacity * processor_daily_flow_vars[
                     (facility_id, day)], "constraint min flow into processor {}, day {}".format(facility_id, day)
-            # else:
-            #     pdb.set_trace()
 
             if is_candidate == 1:
                 # forces processor build var to be correct
@@ -2580,7 +2611,7 @@ def create_constraint_daily_processor_capacity(logger, the_scenario, prob, flow_
                     (facility_id, day)], "constraint forces processor build var to be correct {}, {}".format(
                     facility_id, processor_build_vars[facility_id])
 
-    logger.debug("FINISHED:  create_constraint_daily_processor_capacity")
+    logger.debug("FINISHED: create_constraint_daily_processor_capacity")
     return prob
 
 
@@ -2588,7 +2619,7 @@ def create_constraint_daily_processor_capacity(logger, the_scenario, prob, flow_
 
 
 def create_primary_processor_vertex_constraints(logger, the_scenario, prob, flow_var):
-    logger.debug("STARTING:  create_primary_processor_vertex_constraints - conservation of flow")
+    logger.debug("STARTING: create_primary_processor_vertex_constraints - conservation of flow")
     # for all of these vertices, flow in always  == flow out
     # node_counter = 0
     # node_constraint_counter = 0
@@ -2599,7 +2630,7 @@ def create_primary_processor_vertex_constraints(logger, the_scenario, prob, flow
         # total flow in == total flow out, subject to conversion;
         # dividing by "required quantity" functionally converts all commodities to the same "processor-specific units"
 
-        # processor primary vertices with input commodity and  quantity needed to produce specified output quantities
+        # processor primary vertices with input commodity and quantity needed to produce specified output quantities
         # 2 sets of constraints; one for the primary processor vertex to cover total flow in and out
         # one for each input and output commodity (sum over sources) to ensure its ratio matches facility_commodities
 
@@ -2728,13 +2759,13 @@ def create_primary_processor_vertex_constraints(logger, the_scenario, prob, flow
         constrained_input_flow_vars = set([])
 
         for key, value in iteritems(flow_out_lists):
-            #value is a dictionary with commodity & source as keys
+            # value is a dictionary with commodity & source as keys
             # set up a dictionary that will be filled with input lists to check ratio against
             compare_input_dict = {}
             compare_input_dict_commod = {}
             vertex_id = key
             zero_in = False
-            #value is a dictionary keyed on output commodity, quantity required, edge source
+            # value is a dictionary keyed on output commodity, quantity required, edge source
             if vertex_id in flow_in_lists:
                 in_quantity = 0
                 in_commodity_id = 0
@@ -2754,7 +2785,6 @@ def create_primary_processor_vertex_constraints(logger, the_scenario, prob, flow
             else:
                 zero_in = True
 
-
             # value is a dict - we loop once here for each output commodity and source at the vertex
             for key2, value2 in iteritems(value):
                 out_commodity_id = key2[0]
@@ -2827,7 +2857,7 @@ def create_primary_processor_vertex_constraints(logger, the_scenario, prob, flow
                                                "commodity {} with source {}".format(
                             vertex_id, commodity_id, source)
 
-    logger.debug("FINISHED:  create_primary_processor_conservation_of_flow_constraints")
+    logger.debug("FINISHED: create_primary_processor_conservation_of_flow_constraints")
     return prob
 
 
@@ -2835,7 +2865,7 @@ def create_primary_processor_vertex_constraints(logger, the_scenario, prob, flow
 
 
 def create_constraint_conservation_of_flow(logger, the_scenario, prob, flow_var, processor_excess_vars):
-    logger.debug("STARTING:  create_constraint_conservation_of_flow")
+    logger.debug("STARTING: create_constraint_conservation_of_flow")
     # node_counter = 0
     node_constraint_counter = 0
     storage_vertex_constraint_counter = 0
@@ -3014,7 +3044,7 @@ def create_constraint_conservation_of_flow(logger, the_scenario, prob, flow_var,
             source_facility_id = row_a[9]
             commodity_id = row_a[10]
 
-            # node_counter = node_counter +1
+            # node_counter = node_counter + 1
             # if node is not intermodal, conservation of flow holds per mode;
             # if intermodal, then across modes
             if intermodal == 'N':
@@ -3078,7 +3108,7 @@ def create_constraint_conservation_of_flow(logger, the_scenario, prob, flow_var,
 
         # Note: no consesrvation of flow for primary vertices for supply & demand - they have unique constraints
 
-    logger.debug("FINISHED:  create_constraint_conservation_of_flow")
+    logger.debug("FINISHED: create_constraint_conservation_of_flow")
 
     return prob
 
@@ -3087,7 +3117,7 @@ def create_constraint_conservation_of_flow(logger, the_scenario, prob, flow_var,
 
 
 def create_constraint_max_route_capacity(logger, the_scenario, prob, flow_var):
-    logger.info("STARTING:  create_constraint_max_route_capacity")
+    logger.info("STARTING: create_constraint_max_route_capacity")
     logger.info("modes with background flow turned on: {}".format(the_scenario.backgroundFlowModes))
     # min_capacity_level must be a number from 0 to 1, inclusive
     # min_capacity_level is only relevant when background flows are turned on
@@ -3142,7 +3172,7 @@ def create_constraint_max_route_capacity(logger, the_scenario, prob, flow_var):
         logger.debug("route_capacity constraints created for all storage routes")
 
         # capacity for transport routes
-        # Assumption - all flowing material is in kgal, all flow is summed on a single non-pipeline nx edge
+        # Assumption - all flowing material is in thousand_gallon, all flow is summed on a single non-pipeline nx edge
         sql = """select e.edge_id, e.nx_edge_id, e.max_edge_capacity, e.start_day, e.simple_mode, e.phase_of_matter,
          e.capac_minus_volume_zero_floor
         from edges e
@@ -3181,11 +3211,11 @@ def create_constraint_max_route_capacity(logger, the_scenario, prob, flow_var):
             else:
                 use_capacity = nx_edge_capacity
 
-            # flow is in thousand gallons (kgal), for liquid, or metric tons, for solid
+            # flow is in thousand gallons, for liquid, or metric tons, for solid
             # capacity is in truckload, rail car, barge, or pipeline movement per day
-            # if mode is road and phase is liquid, capacity is in truckloads per day, we want it in kgal
-            # ftot_supporting_gis tells us that there are 8 kgal per truckload,
-            # so capacity * 8 gives us correct units or kgal per day
+            # if mode is road and phase is liquid, capacity is in truckloads per day, we want it in thousand_gallon
+            # ftot_supporting_gis tells us that there are 8 thousand_gallon per truckload,
+            # so capacity * 8 gives us correct units or thousand_gallon per day
             # => use capacity * ftot_supporting_gis multiplier to get capacity in correct flow units
 
             multiplier = 1  # if units match, otherwise specified here
@@ -3214,7 +3244,7 @@ def create_constraint_max_route_capacity(logger, the_scenario, prob, flow_var):
 
         logger.debug("route_capacity constraints created for all non-pipeline  transport routes")
 
-    logger.debug("FINISHED:  create_constraint_max_route_capacity")
+    logger.debug("FINISHED: create_constraint_max_route_capacity")
     return prob
 
 
@@ -3222,7 +3252,7 @@ def create_constraint_max_route_capacity(logger, the_scenario, prob, flow_var):
 
 
 def create_constraint_pipeline_capacity(logger, the_scenario, prob, flow_var):
-    logger.debug("STARTING:  create_constraint_pipeline_capacity")
+    logger.debug("STARTING: create_constraint_pipeline_capacity")
     logger.debug("Length of flow_var: {}".format(len(list(flow_var.items()))))
     logger.info("modes with background flow turned on: {}".format(the_scenario.backgroundFlowModes))
     logger.info("minimum available capacity floor set at: {}".format(the_scenario.minCapacityLevel))
@@ -3272,7 +3302,7 @@ def create_constraint_pipeline_capacity(logger, the_scenario, prob, flow_var):
             # tariff_id = row_a[1]
             link_id = row_a[2]
             # Link capacity is recorded in "thousand barrels per day"; 1 barrel = 42 gall
-            # Link capacity * 42 is now in kgal per day, to match flow in kgal
+            # Link capacity * 42 is now in thousand_gallon per day, to match flow in thousand_gallon
             link_capacity_kgal_per_day = THOUSAND_GALLONS_PER_THOUSAND_BARRELS * row_a[3]
             start_day = row_a[4]
             capac_minus_background_flow_kgal = max(THOUSAND_GALLONS_PER_THOUSAND_BARRELS * row_a[5], 0)
@@ -3296,7 +3326,7 @@ def create_constraint_pipeline_capacity(logger, the_scenario, prob, flow_var):
 
         logger.debug("pipeline capacity constraints created for all transport routes")
 
-    logger.debug("FINISHED:  create_constraint_pipeline_capacity")
+    logger.debug("FINISHED: create_constraint_pipeline_capacity")
     return prob
 
 
@@ -3421,7 +3451,7 @@ def save_pulp_solution(the_scenario, prob, logger, zero_threshold):
     with sqlite3.connect(the_scenario.main_db) as db_con:
 
         db_cur = db_con.cursor()
-        # drop  the optimal_solution table
+        # drop the optimal_solution table
         # -----------------------------
         db_cur.executescript("drop table if exists optimal_solution;")
 
@@ -3476,8 +3506,10 @@ def save_pulp_solution(the_scenario, prob, logger, zero_threshold):
     logger.info(
         "FINISH: save_pulp_solution: Runtime (HMS): \t{}".format(ftot_supporting.get_total_runtime_string(start_time)))
 
+
 # ===============================================================================
 
+
 def record_pulp_solution(the_scenario, logger):
     logger.info("START: record_pulp_solution")
     non_zero_variable_count = 0
@@ -3494,21 +3526,21 @@ def record_pulp_solution(the_scenario, logger):
             null as converted_capacity,
             null as converted_volume,
             null as converted_capac_minus_volume,
-            null as  edge_type,
-            null as  commodity_name,
+            null as edge_type,
+            null as commodity_name,
             null as o_facility,
             'placeholder' as d_facility,
-            null as  o_vertex_id,
-            null as  d_vertex_id,
-            null as  from_node_id,
-            null as  to_node_id,
-            null as  time_period,
+            null as o_vertex_id,
+            null as d_vertex_id,
+            null as from_node_id,
+            null as to_node_id,
+            null as time_period,
             null as commodity_id,
             null as source_facility_id,
             null as source_facility_name,
             null as units,
             variable_name,
-            null as  nx_edge_id,
+            null as nx_edge_id,
             null as mode,
             null as mode_oid,
             null as miles,
@@ -3563,21 +3595,21 @@ def record_pulp_solution(the_scenario, logger):
             null as converted_capacity,
             null as converted_volume,
             null as converted_capac_minus_volume,
-            null as  edge_type,
-            null as  commodity_name,
+            null as edge_type,
+            null as commodity_name,
             'placeholder' as o_facility,
             'placeholder' as d_facility,
-            null as  o_vertex_id,
-            null as  d_vertex_id,
-            null as  from_node_id,
-            null as  to_node_id,
-            null as  time_period,
+            null as o_vertex_id,
+            null as d_vertex_id,
+            null as from_node_id,
+            null as to_node_id,
+            null as time_period,
             null as commodity_id,
             null as source_facility_id,
             null as source_facility_name,
             null as units,
             variable_name,
-            null as  nx_edge_id,
+            null as nx_edge_id,
             null as mode,
             null as mode_oid,
             null as miles,
@@ -3593,6 +3625,7 @@ def record_pulp_solution(the_scenario, logger):
 
     logger.info("FINISH: record_pulp_solution")
 
+
 # ===============================================================================
 
 
@@ -3691,7 +3724,6 @@ def parse_optimal_solution_db(the_scenario, logger):
                 else:
                     optimal_unmet_demand[dest_name][commodity_flowed] += int(v_value)
 
-
     logger.info("length of optimal_processors list: {}".format(len(optimal_processors)))  # a list of optimal processors
     logger.info("length of optimal_processor_flows list: {}".format(
         len(optimal_processor_flows)))  # a list of optimal processor flows
@@ -3701,5 +3733,3 @@ def parse_optimal_solution_db(the_scenario, logger):
         len(optimal_unmet_demand)))  # a dictionary of route keys and unmet demand values
 
     return optimal_processors, optimal_route_flows, optimal_unmet_demand, optimal_storage_flows, optimal_excess_material
-
-
diff --git a/program/ftot_pulp_candidate_generation.py b/program/ftot_pulp_candidate_generation.py
index de10641..9ab4c74 100644
--- a/program/ftot_pulp_candidate_generation.py
+++ b/program/ftot_pulp_candidate_generation.py
@@ -1,13 +1,12 @@
 # ---------------------------------------------------------------------------------------------------
-# Name: pulp c2-spiderweb variant
+# Name: pulp candidate generation
 #
-# Purpose: PulP optimization  - partial source facility as subcommodity variant
+# Purpose: PuLP optimization - partial source facility as subcommodity variant
 # only creates from source edges for commodities with max transport distance
 # other commodities get regular edges to reduce problem size
 # ---------------------------------------------------------------------------------------------------
 
 import datetime
-import pdb
 import re
 import sqlite3
 from collections import defaultdict
@@ -39,6 +38,8 @@
 
 
 # ===============================================================================
+
+
 def oc1(the_scenario, logger):
     # create vertices, then edges for permitted modes, then set volume & capacity on edges
     pre_setup_pulp(logger, the_scenario)
@@ -53,7 +54,15 @@ def oc2(the_scenario, logger):
 
 
 def oc3(the_scenario, logger):
+    # new method when routes are on to break candidate generation result into constituent route segments
+    # new method matches the existing SQL table given to post_optimization
     record_pulp_candidate_gen_solution(the_scenario, logger, zero_threshold)
+
+    if the_scenario.ndrOn:
+        identify_candidate_nodes_from_routes(the_scenario, logger) 
+    else:
+        identify_candidate_nodes(the_scenario, logger) 
+    
     from ftot_supporting import post_optimization
     post_optimization(the_scenario, 'oc3', logger)
     # finalize candidate creation and report out
@@ -141,7 +150,6 @@ def source_as_subcommodity_setup(the_scenario, logger):
             -- TODO as of 2-25 this will throw an error if run repeatedly, 
             -- adding columns to the candidate_process_commodities table
             
-
             drop table if exists candidate_process_commodities_temp;
 
             create table candidate_process_commodities_temp as
@@ -172,7 +180,6 @@ def source_as_subcommodity_setup(the_scenario, logger):
             
             drop table IF EXISTS candidate_process_commodities_temp;
 
-
             """.format(multi_commodity_name, multi_commodity_name)
                                   )
 
@@ -193,6 +200,7 @@ def schedule_avg_availabilities(the_scenario, schedule_dict, schedule_length, lo
 
     return avg_availabilities
 
+
 # ===============================================================================
 
 
@@ -243,6 +251,7 @@ def generate_all_edges_from_source_facilities(the_scenario, schedule_length, log
              group by children_created
              order by children_created asc;""")
         current_edge_data = current_edge_data.fetchall()
+
         for row in current_edge_data:
             if row[1] == 'N':
                 edges_requiring_children = row[0]
@@ -263,7 +272,7 @@ def generate_all_edges_from_source_facilities(the_scenario, schedule_length, log
 
         counter = 0
 
-        #       set up a table to keep track of endcap nodes
+        # set up a table to keep track of endcap nodes
         sql = """
         drop table if exists endcap_nodes;
         create table if not exists endcap_nodes(
@@ -396,7 +405,6 @@ def generate_all_edges_from_source_facilities(the_scenario, schedule_length, log
             left outer join candidate_process_commodities proc on
             (c.commodity_id = proc.input_commodity and best_process_id is not null)
 
-
             where p.to_node_id = ch.from_node_id
             --and p.mode = ch.mode_source --build across modes, control at conservation of flow
             and ch.to_node_id = chtn.node_id
@@ -414,7 +422,6 @@ def generate_all_edges_from_source_facilities(the_scenario, schedule_length, log
 
             main_db_con.execute("update edges set children_created = 'Y' where children_created = 'N';")
 
-
             # this deliberately includes updating "parent" edges that did not get chosen because they weren't the
             # current shortest path
             # those edges are still "resolved" by this batch
@@ -478,7 +485,6 @@ def generate_all_edges_from_source_facilities(the_scenario, schedule_length, log
                             VALUES ({}, '{}', {}, {},{},'{}');
                             """.format(from_node, mode, source_facility_id, commodity_id, input_commodity_process_id, destination_yn))
 
-
                         # designate leadin edge as endcap
                         # this does, deliberately, allow endcap status to be overwritten if we've found a shorter
                         # path to a previous endcap
@@ -952,6 +958,7 @@ def generate_all_edges_without_max_commodity_constraint(the_scenario, schedule_l
                                                        default_min_capacity, route_cost, dollar_cost,
                                                        'transport', nx_edge_id, mode, mode_oid, miles, simple_mode,
                                                        tariff_id, phase_of_matter, source_facility_id))
+
                                 elif from_location == 'NULL' and to_location != 'NULL':
                                     # for each day and commodity, get the corresponding destination vertex id to
                                     # include with the edge info
@@ -980,6 +987,7 @@ def generate_all_edges_without_max_commodity_constraint(the_scenario, schedule_l
                                                        default_min_capacity, route_cost, dollar_cost,
                                                        'transport', nx_edge_id, mode, mode_oid, miles, simple_mode,
                                                        tariff_id, phase_of_matter, source_facility_id))
+
                                 elif from_location != 'NULL' and to_location != 'NULL':
                                     # for each day and commodity, get the corresponding origin and destination vertex
                                     # ids to include with the edge info
@@ -1032,6 +1040,7 @@ def generate_all_edges_without_max_commodity_constraint(the_scenario, schedule_l
         edge_type, nx_edge_id, mode, mode_oid, miles, source_facility_id);""")
         db_cur.execute(sql)
         logger.info("edge_index Total Runtime (HMS): \t{} \t ".format(get_total_runtime_string(index_start_time)))
+
     return
 
 
@@ -1068,15 +1077,20 @@ def pre_setup_pulp(logger, the_scenario):
     logger.debug("----- Using generate_connector_and_storage_edges method imported from ftot_pulp ------")
     generate_connector_and_storage_edges(the_scenario, logger)
 
-    from ftot_pulp import generate_first_edges_from_source_facilities
-    logger.debug("----- Using generate_first_edges_from_source_facilities method imported from ftot_pulp ------")
-    generate_first_edges_from_source_facilities(the_scenario, schedule_avg_length, logger)
+    # pull in generate_edges_from_routes from ftot_pulp and run following code block only if NDR is off
+    if not the_scenario.ndrOn:
+        from ftot_pulp import generate_first_edges_from_source_facilities
+        logger.debug("----- Using generate_first_edges_from_source_facilities method imported from ftot_pulp ------")
+        generate_first_edges_from_source_facilities(the_scenario, schedule_avg_length, logger)
 
-    generate_all_edges_from_source_facilities(the_scenario, schedule_avg_length, logger)
+        generate_all_edges_from_source_facilities(the_scenario, schedule_avg_length, logger)
 
-    clean_up_endcaps(the_scenario, logger)
+        clean_up_endcaps(the_scenario, logger)
 
-    generate_all_edges_without_max_commodity_constraint(the_scenario, schedule_avg_length, logger)
+        generate_all_edges_without_max_commodity_constraint(the_scenario, schedule_avg_length, logger)
+    else:
+        from ftot_pulp import generate_edges_from_routes
+        generate_edges_from_routes(the_scenario, schedule_avg_length, logger)
 
     logger.info("Edges generated for modes: {}".format(the_scenario.permittedModes))
 
@@ -1126,12 +1140,12 @@ def create_unmet_demand_vars(the_scenario, logger):
     logger.info("START: create_unmet_demand_vars")
     demand_var_list = []
     # may create vertices with zero demand, but only for commodities that the facility has demand for at some point
-    #checks material incoming from storage vertex
+    # checks material incoming from storage vertex
 
     with sqlite3.connect(the_scenario.main_db) as main_db_con:
         db_cur = main_db_con.cursor()
         for row in db_cur.execute("""select v.facility_id, v.schedule_day, ifnull(c.supertype, c.commodity_name) 
-        top_level_commodity_name, v.udp
+         top_level_commodity_name, v.udp
          from vertices v, commodities c, facility_type_id ft, facilities f
          where v.commodity_id = c.commodity_id
          and ft.facility_type = "ultimate_destination"
@@ -1285,7 +1299,7 @@ def create_constraint_unmet_demand(logger, the_scenario, prob, flow_var, unmet_d
         actual_demand_dict = {}
 
         # want to specify that all edges leading into this vertex + unmet demand = total demand
-        # demand primary (non-storage)  vertices
+        # demand primary (non-storage) vertices
 
         db_cur = main_db_con.cursor()
         # each row_a is a primary vertex whose edges in contributes to the met demand of var
@@ -1333,8 +1347,6 @@ def create_constraint_unmet_demand(logger, the_scenario, prob, flow_var, unmet_d
                                                                                                               key[1],
                                                                                                               key[2])
             else:
-                if key not in actual_demand_dict:
-                    pdb.set_trace()
                 # no edges in, so unmet demand equals full demand
                 prob += actual_demand_dict[key] == unmet_demand_var[
                     key], "constraint set unmet demand variable for facility {}, day {}, commodity {} -  " \
@@ -1389,7 +1401,7 @@ def create_constraint_max_flow_out_of_supply_vertex(logger, the_scenario, prob,
 
 def create_primary_processor_vertex_constraints(logger, the_scenario, prob, flow_var):
     logger.info("STARTING:  create_primary_processor_vertex_constraints - capacity and conservation of flow")
-    # for all of these vertices, flow in always  == flow out
+    # for all of these vertices, flow in always == flow out
 
     with sqlite3.connect(the_scenario.main_db) as main_db_con:
         db_cur = main_db_con.cursor()
@@ -1684,10 +1696,11 @@ def create_constraint_conservation_of_flow_storage_vertices(logger, the_scenario
 
 # ===============================================================================
 
+
 def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, prob, flow_var,
                                                         processor_excess_vars):
     # This creates constraints for all non-vertex nodes, with variant rules for endcaps nodes
-    logger.debug("STARTING:  create_constraint_conservation_of_flow_endcap_nodes")
+    logger.debug("STARTING: create_constraint_conservation_of_flow_endcap_nodes")
     node_constraint_counter = 0
     passthrough_constraint_counter = 0
     other_endcap_constraint_counter = 0
@@ -1724,7 +1737,6 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
             process_dict.setdefault(commodity_id, [process_id, quantity])
             process_dict[commodity_id] = [process_id, quantity]
 
-
         sql = """select o.process_id, o.commodity_id, o.quantity, i.commodity_id
         from candidate_process_commodities o, candidate_process_commodities i
         where o.io = 'o'
@@ -1759,7 +1771,7 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
             (case when en.node_id is null then 0 else en.source_facility_id end) as endcap_source_facility
             from networkx_nodes nn
             join edges e on (nn.node_id = e.from_node_id or nn.node_id = e.to_node_id)
-            left outer join endcap_nodes en on (nn.node_id = en.node_id and e.commodity_id = en.commodity_id and 
+            left join endcap_nodes en on (nn.node_id = en.node_id and e.commodity_id = en.commodity_id and 
             e.source_facility_id = en.source_facility_id)
             where nn.location_id is null
             order by nn.node_id, e.commodity_id,
@@ -1841,7 +1853,7 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
             process_id = value[0]
             if process_id > 0:
                 endcap_ref[node][3] = process_outputs_dict[process_id]
-
+                
         # if this node has at least one edge flowing out
         for key, value in iteritems(flow_in_lists):
 
@@ -1872,6 +1884,7 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
             # if this commodity and source match the inputs for the endcap process
             # or if this output commodity matches an input process for this node, enter the first loop
             # if output commodity, don't do anything - will be handled via input commodity
+            
             if (node_id in endcap_ref and
                     ((endcap_input_process == endcap_ref[node_id][0] and commodity_id == endcap_ref[node_id][
                         1] and source_facility_id in endcap_ref[node_id][2])
@@ -1914,10 +1927,15 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
                     # this is an endcap for this commodity and source, so unprocessed flow out is not allowed
                     # unless it's from a different source
 
-                    prob += lpSum(flow_out_lists[k] for k in in_key_list if k in flow_out_lists) == lpSum(
-                        0), "conservation of flow (zero allowed out, must be processed), endcap, nx node {}, " \
-                            "source facility {},  commodity {}, day {}, mode {}".format(
-                        node_id, source_facility_id, commodity_id, day, node_mode)
+                    # add logic check to only create zero allowed out, must be processed constraint only when
+                    # edges are matching
+                    for key in in_key_list:
+                        if key in flow_out_lists.keys():
+                            prob += lpSum(flow_out_lists[k] for k in in_key_list if k in flow_out_lists) == lpSum(
+                                0), "conservation of flow (zero allowed out, must be processed), endcap, nx node {}, " \
+                                    "source facility {},  commodity {}, day {}, mode {}".format(
+                                node_id, source_facility_id, commodity_id, day, node_mode)
+                            break  # only make constraint when we encounter the first key in flow_out_lists
 
                     outputs_dict = {}
                     # starting with the 3rd item in the list, the first output commodity tuple
@@ -1933,7 +1951,7 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
                         node_constraint_counter = node_constraint_counter + 1
                         output_source_facility_id = 0
                         output_process_id = 0
-
+                        
                         # setting up the outflow edges to check
                         if node_mode == 'intermodal':
                             out_key = (
@@ -1986,10 +2004,6 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
                                             node_mode)
                                 other_endcap_constraint_counter = other_endcap_constraint_counter + 1
 
-
-
-
-
             else:
                 # if this node is not an endcap, or has no viable candidate process, handle as standard
                 # if this node has at least one edge flowing in, and out, ensure flow is compliant
@@ -2020,7 +2034,9 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
                 node_mode = 'intermodal'
 
             # node has edges flowing out, but not in; restrict flow out to zero
-            if key not in flow_in_lists:
+            # in non-routes, there is an abundance of edges that can transport commodity 2
+            # so this would always be a weird exception 
+            if key not in flow_in_lists and not the_scenario.ndrOn :
                 prob += lpSum(flow_out_lists[key]) == lpSum(
                     0), "conservation of flow (zero allowed out), nx node {}, source facility {}, commodity {},  " \
                         "day {}, mode {}".format(
@@ -2034,7 +2050,7 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
 
     logger.debug("endcap constraints: passthrough {}, other {}, total no reverse {}".format(
         passthrough_constraint_counter, other_endcap_constraint_counter, endcap_no_reverse_counter))
-    logger.info("FINISHED:  create_constraint_conservation_of_flow_endcap_nodes")
+    logger.info("FINISHED: create_constraint_conservation_of_flow_endcap_nodes")
 
     return prob
 
@@ -2043,7 +2059,7 @@ def create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, pr
 
 
 def create_constraint_pipeline_capacity(logger, the_scenario, prob, flow_var):
-    logger.debug("STARTING:  create_constraint_pipeline_capacity")
+    logger.debug("STARTING: create_constraint_pipeline_capacity")
     logger.debug("Length of flow_var: {}".format(len(list(flow_var.items()))))
     logger.info("modes with background flow turned on: {}".format(the_scenario.backgroundFlowModes))
     logger.info("minimum available capacity floor set at: {}".format(the_scenario.minCapacityLevel))
@@ -2093,7 +2109,7 @@ def create_constraint_pipeline_capacity(logger, the_scenario, prob, flow_var):
             tariff_id = row_a[1]
             link_id = row_a[2]
             # Link capacity is recorded in "thousand barrels per day"; 1 barrel = 42 gall
-            # Link capacity * 42 is now in kgal per day, to match flow in kgal
+            # Link capacity * 42 is now in thousand_gallon per day, to match flow in thousand_gallon
             link_capacity_kgal_per_day = THOUSAND_GALLONS_PER_THOUSAND_BARRELS * row_a[3]
             start_day = row_a[4]
             capac_minus_background_flow_kgal = max(THOUSAND_GALLONS_PER_THOUSAND_BARRELS * row_a[5], 0)
@@ -2117,7 +2133,7 @@ def create_constraint_pipeline_capacity(logger, the_scenario, prob, flow_var):
 
         logger.debug("pipeline capacity constraints created for all transport routes")
 
-    logger.debug("FINISHED:  create_constraint_pipeline_capacity")
+    logger.debug("FINISHED: create_constraint_pipeline_capacity")
     return prob
 
 
@@ -2166,6 +2182,7 @@ def setup_pulp_problem_candidate_generation(the_scenario, logger):
     prob = create_constraint_conservation_of_flow_storage_vertices(logger, the_scenario, prob,
                                                                                    flow_vars, processor_excess_vars)
 
+    # TO DO: verify that conservation of flow method for endcap nodes functions with current layout of endcap nodes
     prob = create_constraint_conservation_of_flow_endcap_nodes(logger, the_scenario, prob, flow_vars,
                                                                processor_excess_vars)
 
@@ -2193,12 +2210,14 @@ def setup_pulp_problem_candidate_generation(the_scenario, logger):
 
 # ===============================================================================
 
+
 def record_pulp_candidate_gen_solution(the_scenario, logger, zero_threshold):
     logger.info("START: record_pulp_candidate_gen_solution")
     non_zero_variable_count = 0
 
     with sqlite3.connect(the_scenario.main_db) as db_con:
-
+        
+        # TO DO: the non-zero variable count is never set to anything other than zero...
         logger.info("number of solution variables greater than zero: {}".format(non_zero_variable_count))
         sql = """
             create table optimal_variables as
@@ -2209,26 +2228,26 @@ def record_pulp_candidate_gen_solution(the_scenario, logger, zero_threshold):
             null as converted_capacity,
             null as converted_volume,
             null as converted_capac_minus_volume,
-            null as  edge_type,
-            null as  commodity_name,
+            null as edge_type,
+            null as commodity_name,
             null as o_facility,
             'placeholder' as d_facility,
-            null as  o_vertex_id,
-            null as  d_vertex_id,
-            null as  from_node_id,
-            null as  to_node_id,
-            null as  time_period,
+            null as o_vertex_id,
+            null as d_vertex_id,
+            null as from_node_id,
+            null as to_node_id,
+            null as time_period,
             null as commodity_id,
             null as source_facility_id,
             null as source_facility_name,
             null as units,
             variable_name,
-            null as  nx_edge_id,
+            null as nx_edge_id,
             null as mode,
             null as mode_oid,
             null as miles,
             null as edge_count_from_source,
-            null as  miles_travelled
+            null as miles_travelled
             from optimal_solution
             where variable_name like 'UnmetDemand%'
             union
@@ -2274,21 +2293,21 @@ def record_pulp_candidate_gen_solution(the_scenario, logger, zero_threshold):
             null as converted_capacity,
             null as converted_volume,
             null as converted_capac_minus_volume,
-            null as  edge_type,
-            null as  commodity_name,
+            null as edge_type,
+            null as commodity_name,
             'placeholder' as o_facility,
             'placeholder' as d_facility,
-            null as  o_vertex_id,
-            null as  d_vertex_id,
-            null as  from_node_id,
-            null as  to_node_id,
-            null as  time_period,
+            null as o_vertex_id,
+            null as d_vertex_id,
+            null as from_node_id,
+            null as to_node_id,
+            null as time_period,
             null as commodity_id,
             null as source_facility_id,
             null as source_facility_name,
             null as units,
             variable_name,
-            null as  nx_edge_id,
+            null as nx_edge_id,
             null as mode,
             null as mode_oid,
             null as miles,
@@ -2300,19 +2319,29 @@ def record_pulp_candidate_gen_solution(the_scenario, logger, zero_threshold):
         db_con.execute("drop table if exists optimal_variables;")
         db_con.executescript(sql)
 
+    logger.info("FINISH: record_pulp_candidate_gen_solution")
+
+
+# ===============================================================================
+
+
+def identify_candidate_nodes(the_scenario, logger):
 
+    logger.info("START: identify_candidate_nodes")
+
+    with sqlite3.connect(the_scenario.main_db) as db_con:
         sql = """
             create table optimal_variables_c as
             select * from optimal_variables
             ;
-             drop table if exists candidate_nodes;
+            drop table if exists candidate_nodes;
 
             create table candidate_nodes as
             select *
             from
 
-            (select pl.minsize, pl.process_id, ov.to_node_id, ov.mode_oid, ov.mode, ov.commodity_id,
-            sum(variable_value) agg_value, count(variable_value) num_aggregated
+            (select pl.minsize, pl.process_id as i_process_id, ov.to_node_id, ov.mode_oid, ov.mode, ov.commodity_id as i_commodity_id,
+            sum(variable_value) as agg_value, count(variable_value) as num_aggregated
             from optimal_variables ov, candidate_process_list pl, candidate_process_commodities pc
 
             where pc.io = 'i'
@@ -2322,8 +2351,8 @@ def record_pulp_candidate_gen_solution(the_scenario, logger, zero_threshold):
             and variable_value > {}
             group by ov.mode_oid, ov.to_node_id, ov.commodity_id, pl.process_id, pl.minsize, ov.mode) i,
 
-            (select pl.min_aggregation, pl.process_id, ov.from_node_id, ov.commodity_id,
-            sum(variable_value) agg_value, count(variable_value) num_aggregated, min(edge_count_from_source) as 
+            (select pl.min_aggregation, pl.process_id as o_process_id, ov.from_node_id, ov.commodity_id as o_commodity_id,
+            sum(variable_value) as agg_value, count(variable_value) as num_aggregated, min(edge_count_from_source) as 
             edge_count,
             o_vertex_id
             from optimal_variables ov, candidate_process_list pl, candidate_process_commodities pc
@@ -2336,18 +2365,17 @@ def record_pulp_candidate_gen_solution(the_scenario, logger, zero_threshold):
             on o.from_node_id = nn.node_id
 
             where (i.to_node_id = o.from_node_id
-            and i.process_id = o.process_id
-            and i.commodity_id = o.commodity_id
+            and i_process_id = o_process_id
+            and i_commodity_id = o_commodity_id
             and o.agg_value > i.agg_value
             and o.agg_value >= o.min_aggregation)
-            or (o.o_vertex_id is not null and  o.agg_value >= o.min_aggregation)
+            or (o.o_vertex_id is not null and o.agg_value >= o.min_aggregation)
 
-            group by o.min_aggregation, o.process_id, o.from_node_id, o.commodity_id,  o.agg_value, o.num_aggregated, 
+            group by o.min_aggregation, o_process_id, o.from_node_id, o_commodity_id, o.agg_value, o.num_aggregated, 
             o.o_vertex_id
             order by o.agg_value, edge_count
             ;
 
-           
             """.format(zero_threshold)
         db_con.execute("drop table if exists optimal_variables_c;")
         db_con.executescript(sql)
@@ -2361,9 +2389,124 @@ def record_pulp_candidate_gen_solution(the_scenario, logger, zero_threshold):
             logger.warning("the candidate nodes table is empty. Hints: "
                            "(1) increase the maximum raw material transport distance to allow additional flows to "
                            "aggregate and meet the minimum candidate facility size"
-                           "(2) increase the ammount of material available to flow from raw material producers")
+                           "(2) increase the amount of material available to flow from raw material producers")
 
-    logger.info("FINISH: record_pulp_candidate_gen_solution")
+    logger.info("FINISH: identify_candidate_nodes")
+
+
+# ===============================================================================
+
+
+def identify_candidate_nodes_from_routes(the_scenario, logger):
+
+    logger.info("START: identify_candidate_nodes_from_routes")
+
+
+    with sqlite3.connect(the_scenario.main_db) as db_con:
+        sql = """
+            create table optimal_variables_c as
+            select * from optimal_variables
+            ;
+            drop table if exists candidate_nodes;
+
+            -- SQL query below maps "edges" results (on level of O-D pairs) to NetworkX edges
+
+            -- BEGIN post process query
+
+            drop table if exists temp_optimal_variables;
+
+            create table temp_optimal_variables as
+            select
+            nx_e.mode_source,
+            nx_e.mode_source_oid,
+            c.commodity_name,
+            c.commodity_id,
+            sum(ov.variable_value) as variable_value,
+            ov.converted_volume,
+            ov.converted_capacity,
+            ov.converted_capac_minus_volume,
+            ov.edge_type,
+            ov.units,
+            ov.time_period,
+            ov.edge_count_from_source, --NULL
+            CASE WHEN re.rt_order_ind = 1 THEN ov.o_vertex_id ELSE NULL END AS o_vertex_id,
+            nx_e.miles,
+            c.phase_of_matter,
+            nx_e_cost.dollar_cost,
+            nx_e_cost.route_cost,
+            nx_e.artificial,
+			re.scenario_rt_id,
+			nx_e.from_node_id,
+			nx_e.to_node_id
+            from optimal_variables ov
+            join commodities as c on c.commodity_id = ov.commodity_id
+            join edges as e on e.edge_id = ov.var_id
+            join route_reference as rr on rr.route_id = e.route_id
+            join route_edges as re on re.scenario_rt_id = rr.scenario_rt_id
+            left join networkx_edges as nx_e on nx_e.edge_id = re.edge_id
+            join networkx_edge_costs as nx_e_cost on nx_e_cost.edge_id = re.edge_id and nx_e_cost.phase_of_matter_id = c.phase_of_matter
+            group by nx_e.edge_id, c.commodity_id
+            ;
+            -- END post process query
+
+            -- use temp_optimal_variables table instead of optimal_variables table to identify candidate nodes
+
+            create table candidate_nodes as
+            select *
+            from
+
+            (select pl.minsize, pl.process_id as i_process_id, ov.to_node_id, ov.mode_source_oid, ov.mode_source, ov.commodity_id as i_commodity_id,
+            sum(variable_value) as agg_value, count(variable_value) as num_aggregated
+            from temp_optimal_variables ov, candidate_process_list pl, candidate_process_commodities pc
+
+            where pc.io = 'i'
+            and pl.process_id = pc.process_id
+            and ov.commodity_id = pc.commodity_id
+            and ov.edge_type = 'transport'
+            and variable_value > {}
+            group by ov.mode_source_oid, ov.to_node_id, ov.commodity_id, pl.process_id, pl.minsize, ov.mode_source) i,
+
+            (select pl.min_aggregation, pl.process_id as o_process_id, ov.from_node_id, ov.commodity_id as o_commodity_id,
+            sum(variable_value) as agg_value, count(variable_value) as num_aggregated, min(edge_count_from_source) as 
+            edge_count,
+            o_vertex_id
+            from temp_optimal_variables ov, candidate_process_list pl, candidate_process_commodities pc
+            where pc.io = 'i'
+            and pl.process_id = pc.process_id
+            and ov.commodity_id = pc.commodity_id
+            and ov.edge_type = 'transport'
+            group by ov.from_node_id, ov.commodity_id, pl.process_id, pl.min_aggregation, ov.o_vertex_id) o
+            left outer join networkx_nodes nn
+            on o.from_node_id = nn.node_id
+
+            where (i.to_node_id = o.from_node_id
+            and i_process_id = o_process_id
+            and i_commodity_id = o_commodity_id
+            and o.agg_value > i.agg_value
+            and o.agg_value >= o.min_aggregation)
+            or (o.o_vertex_id is not null and o.agg_value >= o.min_aggregation)
+
+            group by o.min_aggregation, o_process_id, o.from_node_id, o_commodity_id, o.agg_value, o.num_aggregated, 
+            o.o_vertex_id
+            order by o.agg_value, edge_count
+            ;
+
+            """.format(zero_threshold)
+        db_con.execute("drop table if exists optimal_variables_c;")
+        db_con.executescript(sql)
+
+        sql_check_candidate_node_count = "select count(*) from candidate_nodes"
+        db_cur = db_con.execute(sql_check_candidate_node_count)
+        candidate_node_count = db_cur.fetchone()[0]
+        if candidate_node_count > 0:
+            logger.info("number of candidate nodes identified: {}".format(candidate_node_count))
+        else:
+            logger.warning("the candidate nodes table is empty. Hints: "
+                           "(1) increase the maximum raw material transport distance to allow additional flows to "
+                           "aggregate and meet the minimum candidate facility size"
+                           "(2) increase the amount of material available to flow from raw material producers")
+
+    logger.info("FINISH: identify_candidate_nodes_from_routes")
 
 
 # ===============================================================================
@@ -2391,15 +2534,15 @@ def parse_optimal_solution_db(the_scenario, logger):
         # do the Route Edges
         sql = """select
             variable_name, variable_value,
-            cast(substr(variable_name, 6) as int) edge_id,
+            cast(substr(variable_name, 6) as int) as edge_id,
             route_ID, start_day time_period, edges.commodity_id,
             o_vertex_id, d_vertex_id,
             v1.facility_id o_facility_id,
-         v2.facility_id d_facility_id
+            v2.facility_id d_facility_id
             from optimal_solution
             join edges on edges.edge_id = cast(substr(variable_name, 6) as int)
             join vertices v1 on edges.o_vertex_id = v1.vertex_id
-        join vertices v2 on edges.d_vertex_id = v2.vertex_id
+            join vertices v2 on edges.d_vertex_id = v2.vertex_id
             where variable_name like 'Edge%_' and variable_name not like 'Edge%_storage';
             """
         data = db_con.execute(sql)
@@ -2436,7 +2579,6 @@ def parse_optimal_solution_db(the_scenario, logger):
         optimal_processor_flows_sql = data.fetchall()
         for proc in optimal_processor_flows_sql:
             optimal_processor_flows.append(proc)
-            # optimal_biorefs.append(v.name[22:(v.name.find(",")-1)])# find the name from the
 
         # do the UnmetDemand
         sql = "select variable_name, variable_value from optimal_solution where variable_name like 'UnmetDemand%';"
@@ -2461,7 +2603,6 @@ def parse_optimal_solution_db(the_scenario, logger):
                 else:
                     optimal_unmet_demand[dest_name][commodity_flowed] += int(v_value)
 
-
     logger.info("length of optimal_processors list: {}".format(len(optimal_processors)))  # a list of optimal processors
     logger.info("length of optimal_processor_flows list: {}".format(
         len(optimal_processor_flows)))  # a list of optimal processor flows
diff --git a/program/ftot_pulp_sourcing.py b/program/ftot_pulp_sourcing.py
index 4891292..1c2e07e 100644
--- a/program/ftot_pulp_sourcing.py
+++ b/program/ftot_pulp_sourcing.py
@@ -1586,7 +1586,7 @@ def create_constraint_max_route_capacity(logger, the_scenario, prob, flow_var):
 
 
         #capacity for transport routes
-        #Assumption - all flowing material is in kgal, all flow is summed on a single non-pipeline nx edge
+        #Assumption - all flowing material is in thousand_gallon, all flow is summed on a single non-pipeline nx edge
         sql = """select e.edge_id, e.nx_edge_id, e.max_edge_capacity, e.start_day, e.simple_mode, e.phase_of_matter, e.capac_minus_volume_zero_floor
         from edges e
         where e.max_edge_capacity is not null
@@ -1621,10 +1621,10 @@ def create_constraint_max_route_capacity(logger, the_scenario, prob, flow_var):
                use_capacity = min_restricted_capacity
             else: use_capacity = nx_edge_capacity
 
-            #flow is in thousand gallons (kgal), for liquid, or metric tons, for solid
-            #capacity is in truckload, rail car, barge, or pipeline movement per day
-            # if mode is road and phase is liquid, capacity is in truckloads per day, we want it in kgal
-            # ftot_supporting_gis tells us that there are 8 kgal per truckload, so capacity * 8 gives us correct units or kgal per day
+            # flow is in thousand gallons, for liquid, or metric tons, for solid
+            # capacity is in truckload, rail car, barge, or pipeline movement per day
+            # if mode is road and phase is liquid, capacity is in truckloads per day, we want it in thousand_gallon
+            # ftot_supporting_gis tells us that there are 8 thousand_gallon per truckload, so capacity * 8 gives us correct units or thousand_gallon per day
             # use capacity * ftot_supporting_gis multiplier to get capacity in correct flow units
 
             multiplier = 1 #unless otherwise specified
@@ -1706,7 +1706,7 @@ def create_constraint_pipeline_capacity(logger, the_scenario, prob, flow_var):
             tariff_id = row_a[1]
             link_id = row_a[2]
             # Link capacity is recorded in "thousand barrels per day"; 1 barrel = 42 gall
-            # Link capacity * 42 is now in kgal per day, to match flow in kgal
+            # Link capacity * 42 is now in thousand_gallon per day, to match flow in thousand_gallon
             link_capacity_kgal_per_day = THOUSAND_GALLONS_PER_THOUSAND_BARRELS*row_a[3]
             start_day = row_a[4]
             capac_minus_background_flow_kgal = max(THOUSAND_GALLONS_PER_THOUSAND_BARRELS*row_a[5],0)
diff --git a/program/ftot_report.py b/program/ftot_report.py
index 572d3e8..f47f73a 100644
--- a/program/ftot_report.py
+++ b/program/ftot_report.py
@@ -25,7 +25,7 @@
 # ==================================================================
 
 
-def prepare_tableau_assets(timestamp_directory, report_file, the_scenario, logger):
+def prepare_tableau_assets(timestamp_directory, the_scenario, logger):
     logger.info("start: prepare_tableau_assets")
 
     # make tableau_output.gdb file
@@ -140,8 +140,27 @@ def prepare_tableau_assets(timestamp_directory, report_file, the_scenario, logge
     # copy tableau report to the assets location
     latest_generic_path = os.path.join(timestamp_directory, "tableau_report.csv")
     logger.debug("copying the latest tableau report csv file to the timestamped tableau report directory")
+    report_file_name = 'report_' + TIMESTAMP.strftime("%Y_%m_%d_%H-%M-%S") + ".csv"
+    report_file_name = clean_file_name(report_file_name)
+    report_file = os.path.join(timestamp_directory, report_file_name)
     copy(report_file, latest_generic_path)
 
+    # add all_routes report to assets location
+    logger.debug("adding routes report csv file to the tableau report directory")
+    if the_scenario.ndrOn:
+        # if NDR On, copy existing routes report
+        latest_routes_path = os.path.join(timestamp_directory, "all_routes.csv")
+        routes_file_name = 'all_routes_' + TIMESTAMP.strftime("%Y_%m_%d_%H-%M-%S") + ".csv"
+        routes_file_name = clean_file_name(routes_file_name)
+        routes_file = os.path.join(timestamp_directory, routes_file_name)
+        copy(routes_file, latest_routes_path)
+    else:
+        # generate placeholder all_routes report
+        routes_file = os.path.join(timestamp_directory, "all_routes.csv")
+        with open(routes_file, 'w', newline='') as f:
+            writer = csv.writer(f)
+            writer.writerow(['Run scenario with NDR on to view this dashboard.'])
+   
     # create packaged workbook for tableau reader compatibility
     twbx_dashboard_filename = os.path.join(timestamp_directory, "tableau_dashboard.twbx")
     zipObj = zipfile.ZipFile(twbx_dashboard_filename, 'w', zipfile.ZIP_DEFLATED)
@@ -153,40 +172,55 @@ def prepare_tableau_assets(timestamp_directory, report_file, the_scenario, logge
     zipObj.write(os.path.join(timestamp_directory, "volpeTriskelion.gif"), "volpeTriskelion.gif")
     zipObj.write(os.path.join(timestamp_directory, "parameters_icon.png"), "parameters_icon.png")
     zipObj.write(os.path.join(timestamp_directory, "tableau_output.gdb.zip"), "tableau_output.gdb.zip")
+    zipObj.write(os.path.join(timestamp_directory, "all_routes.csv"), "all_routes.csv")
 
     # close the Zip File
     zipObj.close()
 
-    # delete the other four files so its nice an clean.
+    # delete the other files so its nice an clean.
     os.remove(os.path.join(timestamp_directory, "tableau_dashboard.twb"))
     os.remove(os.path.join(timestamp_directory, "tableau_report.csv"))
     os.remove(os.path.join(timestamp_directory, "volpeTriskelion.gif"))
     os.remove(os.path.join(timestamp_directory, "parameters_icon.png"))
     os.remove(os.path.join(timestamp_directory, "tableau_output.gdb.zip"))
-
+    os.remove(os.path.join(timestamp_directory, "all_routes.csv"))
 
 # ==============================================================================================
 
 def generate_edges_from_routes_summary(timestamp_directory, the_scenario, logger):
 
     logger.info("start: generate_edges_from_routes_summary")
-    report_file_name = 'optimal_routes_' + TIMESTAMP.strftime("%Y_%m_%d_%H-%M-%S") + ".csv"
+    report_file_name = 'all_routes_' + TIMESTAMP.strftime("%Y_%m_%d_%H-%M-%S") + ".csv"
     report_file_name = clean_file_name(report_file_name)
     report_file = os.path.join(timestamp_directory, report_file_name)
     
     with sqlite3.connect(the_scenario.main_db) as main_db_con:
         db_cur = main_db_con.cursor()
-        summary_route_data = main_db_con.execute("""select rr.route_id, f1.facility_name as from_facility, 
-                            f2.facility_name as to_facility, rr.phase_of_matter, rr.dollar_cost, rr.cost, rr.miles 
-                            FROM route_reference rr 
-                            join facilities f1 on rr.from_facility_id = f1.facility_id
-                            join facilities f2 on rr.to_facility_id = f2.facility_id; """)
+        summary_route_data = main_db_con.execute("""select rr.route_id, f1.facility_name as from_facility, fti1.facility_type as from_facility_type,
+            f2.facility_name as to_facility, fti2.facility_type as to_facility_type,
+            c.commodity_name, c.phase_of_matter, m.mode, rr.dollar_cost, rr.cost, rr.miles,
+            case when ors.scenario_rt_id is NULL then "N" else "Y" end as in_solution
+            from route_reference rr
+            join facilities f1 on rr.from_facility_id = f1.facility_id
+            join facility_type_id fti1 on f1.facility_type_id = fti1.facility_type_id
+            join facilities f2 on rr.to_facility_id = f2.facility_id
+            join facility_type_id fti2 on f2.facility_type_id = fti2.facility_type_id
+            join commodities c on rr.commodity_id = c.commodity_ID
+            left join (select temp.scenario_rt_id, case when temp.modes_list like "%,%" then "multimodal" else temp.modes_list end as mode
+                       from (select re.scenario_rt_id, group_concat(distinct(nx_e.mode_source)) as modes_list
+                       from route_edges re
+                       left join networkx_edges nx_e on re.edge_id = nx_e.edge_id
+                       group by re.scenario_rt_id) temp) m on rr.scenario_rt_id = m.scenario_rt_id
+            left join (select distinct scenario_rt_id from optimal_route_segments) ors on rr.scenario_rt_id = ors.scenario_rt_id;""")
         
-        # Print route data to file in debug folder
+        # print route data to file in Reports folder
         with open(report_file, 'w', newline='') as f:
             writer = csv.writer(f)
-            writer.writerow(['route_id','from_facility','to_facility','phase','dollar_cost','routing_cost','miles'])
-            writer.writerows(summary_route_data)
+            writer.writerow(['scenario_name', 'route_id', 'from_facility', 'from_facility_type', 'to_facility', 'to_facility_type',
+                             'commodity_name', 'phase', 'mode', 'dollar_cost', 'routing_cost', 'miles', 'in_solution'])
+            for row in summary_route_data:
+                writer.writerow([the_scenario.scenario_name, row[0], row[1], row[2], row[3], row[4], row[5],
+                                 row[6], row[7], row[8], row[9], row[10], row[11]])
 
 
 # ==============================================================================================
@@ -245,7 +279,7 @@ def generate_artificial_link_summary(timestamp_directory, the_scenario, logger):
             else:
                 output_table[permitted_mode].append('NA')
 
-    # print artificial link data for each facility to file in debug folder
+    # print artificial link data for each facility to file in Reports folder
     with open(report_file, 'w', newline='') as f:
         writer = csv.writer(f)
         output_fields = ['facility_name', 'facility_type'] + the_scenario.permittedModes
@@ -285,7 +319,7 @@ def generate_reports(the_scenario, logger):
     
     from ftot_networkx import update_ndr_parameter
 
-    # Check NDR conditions before reporting optimal routes
+    # Check NDR conditions before reporting routes
     update_ndr_parameter(the_scenario, logger)
 
     # make overall Reports directory
@@ -520,7 +554,8 @@ def generate_reports(the_scenario, logger):
             # note: processor input and outputs are based on facility size and reflect a processing capacity,
             # not a conversion of the scenario feedstock supply
             # -------------------------------------------------------------------------
-            sql = """   select c.commodity_name, fti.facility_type, io, sum(fc.quantity), fc.units  
+            sql = """   select c.commodity_name, fti.facility_type, io, sum(fc.quantity), fc.units,
+                              (case when fc.scaled_quantity is null then 'Unconstrained' else sum(fc.scaled_quantity) end)
                               from facility_commodities fc
                               join commodities c on fc.commodity_id = c.commodity_id
                               join facilities f on f.facility_id = fc.facility_id
@@ -534,18 +569,23 @@ def generate_reports(the_scenario, logger):
                 io = row[2]
                 if facility_type == "raw_material_producer":
                     measure = "supply"
+                    value = row[3]
                 if facility_type == "processor" and io == 'o':
-                    measure = "processing output capacity"
+                    measure = "processing_output_capacity"
+                    value = row[5]
                 if facility_type == "processor" and io == 'i':
-                    measure = "processing input capacity"
+                    measure = "processing_input_capacity"
+                    value = row[5]
                 if facility_type == "ultimate_destination":
                     measure = "demand"
-                writer.writerow([the_scenario.scenario_name, "total_supply_demand_proc", row[0], "all_{}".format(row[1]), measure, io, row[3], row[4], None])
+                    value = row[3]
+                writer.writerow([the_scenario.scenario_name, "total_supply_demand_proc", row[0], "all_{}".format(row[1]), measure, io, value, row[4], None])
 
             # Scenario Stranded Supply, Demand, and Processing Capacity
             # note: stranded supply refers to facilities that are ignored from the analysis.")
             # -------------------------------------------------------------------------
-            sql = """ select c.commodity_name, fti.facility_type, io, sum(fc.quantity), fc.units, f.ignore_facility 
+            sql = """ select c.commodity_name, fti.facility_type, io, sum(fc.quantity), fc.units, f.ignore_facility,
+                      (case when fc.scaled_quantity is null then 'Unconstrained' else sum(fc.scaled_quantity) end)
                             from facility_commodities fc
                             join commodities c on fc.commodity_id = c.commodity_id
                             join facilities f on f.facility_id = fc.facility_id
@@ -560,19 +600,24 @@ def generate_reports(the_scenario, logger):
                 facility_type = row[1]
                 io = row[2]
                 if facility_type == "raw_material_producer":
-                    measure = "stranded supply"
+                    measure = "stranded_supply"
+                    value = row[3]
                 if facility_type == "processor" and io == 'o':
-                    measure = "stranded processing output capacity"
+                    measure = "stranded_processing_output_capacity"
+                    value = row[6]
                 if facility_type == "processor" and io == 'i':
-                    measure = "stranded processing input capacity"
+                    measure = "stranded_processing_input_capacity"
+                    value = row[6]
                 if facility_type == "ultimate_destination":
-                    measure = "stranded demand"
-                writer.writerow([the_scenario.scenario_name, "stranded_supply_demand_proc", row[0], "stranded_{}".format(row[1]), measure, io, row[3], row[4], row[5]]) # ignore_facility note
+                    measure = "stranded_demand"
+                    value = row[3]
+                writer.writerow([the_scenario.scenario_name, "stranded_supply_demand_proc", row[0], "stranded_{}".format(row[1]), measure, io, value, row[4], row[5]]) # ignore_facility note
 
             # report out net quantities with ignored facilities removed from the query
             # note: net supply, demand, and processing capacity ignores facilities not connected to the network
             # -------------------------------------------------------------------------
-            sql = """   select c.commodity_name, fti.facility_type, io, sum(fc.quantity), fc.units  
+            sql = """   select c.commodity_name, fti.facility_type, io, sum(fc.quantity), fc.units,
+                        (case when fc.scaled_quantity is null then 'Unconstrained' else sum(fc.scaled_quantity) end) 
                               from facility_commodities fc
                               join commodities c on fc.commodity_id = c.commodity_id
                               join facilities f on f.facility_id = fc.facility_id
@@ -586,14 +631,18 @@ def generate_reports(the_scenario, logger):
                 facility_type = row[1]
                 io = row[2]
                 if facility_type == "raw_material_producer":
-                    measure = "net supply"
+                    measure = "net_supply"
+                    value = row[3]
                 if facility_type == "processor" and io == 'o':
-                    measure = "net processing output capacity"
+                    measure = "net_processing_output_capacity"
+                    value = row[5]
                 if facility_type == "processor" and io == 'i':
-                    measure = "net processing input capacity"
+                    measure = "net_processing_input_capacity"
+                    value = row[5]
                 if facility_type == "ultimate_destination":
-                    measure = "net demand"
-                writer.writerow([the_scenario.scenario_name, "net_supply_demand_proc", row[0], "net_{}".format(row[1]), measure, io, row[3], row[4], None])
+                    measure = "net_demand"
+                    value = row[3]
+                writer.writerow([the_scenario.scenario_name, "net_supply_demand_proc", row[0], "net_{}".format(row[1]), measure, io, value, row[4], None])
 
         # REPORT OUT CONFIG FOR O2 STEP AND RUNTIMES FOR ALL STEPS
         # Loop through the list of configurations records in the message_dict['config'] and ['runtime'].
@@ -618,8 +667,6 @@ def generate_reports(the_scenario, logger):
             writer.writerow([the_scenario.scenario_name, "runtime", '', '', step, '', '', '', runtime])
 
     logger.debug("finish: main results report operation")
-
-    prepare_tableau_assets(timestamp_directory, report_file, the_scenario, logger)
     
     # -------------------------------------------------------------
 
@@ -634,4 +681,7 @@ def generate_reports(the_scenario, logger):
     if the_scenario.detailed_emissions:
         generate_detailed_emissions_summary(timestamp_directory, the_scenario, logger)
 
+    # tableau workbook
+    prepare_tableau_assets(timestamp_directory, the_scenario, logger)
+
     logger.result("Reports located here: {}".format(timestamp_directory))
diff --git a/program/ftot_scenario.py b/program/ftot_scenario.py
index fd0bdb2..aa37f6b 100644
--- a/program/ftot_scenario.py
+++ b/program/ftot_scenario.py
@@ -181,14 +181,14 @@ def check_relative_paths(mypath):
         scenario.railroadCO2Emissions = Q_(xmlScenarioFile.getElementsByTagName('Railroad_CO2_Emissions')[0].firstChild.data).to('g/{}/mi'.format(scenario.default_units_solid_phase))
         scenario.bargeCO2Emissions = Q_(xmlScenarioFile.getElementsByTagName('Barge_CO2_Emissions')[0].firstChild.data).to('g/{}/mi'.format(scenario.default_units_solid_phase))
         scenario.pipelineCO2Emissions = Q_(xmlScenarioFile.getElementsByTagName('Pipeline_CO2_Emissions')[0].firstChild.data).to('g/{}/mi'.format(scenario.default_units_solid_phase))
-        # setting density conversion based on 'Density_Conversion_Factor' field if it exists, or otherwise default to 3.33 ton/kgal
+        # setting density conversion based on 'Density_Conversion_Factor' field if it exists, or otherwise default to 3.33 ton/thousand_gallon
         if len(xmlScenarioFile.getElementsByTagName('Density_Conversion_Factor')):
             scenario.densityFactor = Q_(xmlScenarioFile.getElementsByTagName('Density_Conversion_Factor')[0].firstChild.data).to('{}/{}'.format(scenario.default_units_solid_phase, scenario.default_units_liquid_phase))
         else:
             if scenario.commodity_density_data == "None":
                 # User didn't specify a density file OR factor
-                logger.warning("FTOT is assuming a density of 3.33 ton/kgal for emissions reporting for liquids. Use scenario XML parameter 'Density_Conversion_Factor' to adjust this value.")
-            scenario.densityFactor = Q_('3.33 ton/kgal').to('{}/{}'.format(scenario.default_units_solid_phase, scenario.default_units_liquid_phase))
+                logger.warning("FTOT is assuming a density of 3.33 ton/thousand_gallon for emissions reporting for liquids. Use scenario XML parameter 'Density_Conversion_Factor' to adjust this value.")
+            scenario.densityFactor = Q_('3.33 ton/thousand_gallon').to('{}/{}'.format(scenario.default_units_solid_phase, scenario.default_units_liquid_phase))
         logger.debug("PASS: setting the vehicle emission factors with pint passed")
     
     except Exception as e:
diff --git a/program/ftot_supporting.py b/program/ftot_supporting.py
index 4dfc476..5ee9339 100644
--- a/program/ftot_supporting.py
+++ b/program/ftot_supporting.py
@@ -428,11 +428,11 @@ def get_processor_capacity(primary_processing, logger):
     capacity = 0
 
     if primary_processing == "FTx":
-        capacity = Q_(200000, "kgal")
+        capacity = Q_(200000, "thousand_gallon")
     if primary_processing == "Petroleum_Refinery":
-        capacity = Q_(7665000, "kgal")
+        capacity = Q_(7665000, "thousand_gallon")
     else:
-        capacity = Q_(200000, "kgal")
+        capacity = Q_(200000, "thousand_gallon")
 
     return capacity
 
@@ -463,13 +463,13 @@ def get_input_and_output_commodity_quantities_from_afpat(commodity, process, the
 
     if commodity.lower().find("test_liquid_none_none") > -1:
         #        print "in the right place"
-        input_commodity_quantities = Q_(1, "kgal")
-        output_commodity_quantities['test_product_liquid_None_None'] = Q_(1, "kgal")
-        output_commodity_quantities['jet'] = Q_(0, "kgal")
-        output_commodity_quantities['diesel'] = Q_(0, "kgal")
-        output_commodity_quantities['naphtha'] = Q_(0, "kgal")
-        output_commodity_quantities['aromatics'] = Q_(0, "kgal")
-        output_commodity_quantities['total_fuel'] = Q_(1, "kgal")
+        input_commodity_quantities = Q_(1, "thousand_gallon")
+        output_commodity_quantities['test_product_liquid_None_None'] = Q_(1, "thousand_gallon")
+        output_commodity_quantities['jet'] = Q_(0, "thousand_gallon")
+        output_commodity_quantities['diesel'] = Q_(0, "thousand_gallon")
+        output_commodity_quantities['naphtha'] = Q_(0, "thousand_gallon")
+        output_commodity_quantities['aromatics'] = Q_(0, "thousand_gallon")
+        output_commodity_quantities['total_fuel'] = Q_(1, "thousand_gallon")
         # hack to skip the petroleum refinery
         commodity = "hack-hack-hack"
         process = "hack-hack-hack"

From 2892343a65bb8e067d2cccc5dfb2dcaa85fddc21 Mon Sep 17 00:00:00 2001
From: Kevin Zhang <68650998+kzhang81@users.noreply.github.com>
Date: Fri, 13 Jan 2023 12:07:46 -0500
Subject: [PATCH 3/4] Add files via upload

---
 program/lib/Master_FTOT_Schema.xsd | 28 ++++++++++++------------
 program/lib/v6_temp_Scenario.xml   | 34 +++++++++++++++---------------
 program/lib/vehicle_types.csv      |  2 +-
 3 files changed, 32 insertions(+), 32 deletions(-)

diff --git a/program/lib/Master_FTOT_Schema.xsd b/program/lib/Master_FTOT_Schema.xsd
index db3c327..4cebbd5 100644
--- a/program/lib/Master_FTOT_Schema.xsd
+++ b/program/lib/Master_FTOT_Schema.xsd
@@ -6,7 +6,7 @@ elementFormDefault="qualified">
 <xs:element name="Scenario">
 	<xs:complexType>
 		<xs:sequence>
-			<xs:element name="Scenario_Schema_Version" type="xs:string" default="6.0.3" />
+			<xs:element name="Scenario_Schema_Version" type="xs:string" default="6.0.4" />
 			<xs:element name="Scenario_Name" type="xs:string" default="USER INPUT REQUIRED: replace this string with a descriptive Scenario Name" />
 			<xs:element name="Scenario_Description" type="xs:string" default="USER INPUT REQUIRED: replace this string with a descriptive Scenario Description"/>
 
@@ -28,7 +28,7 @@ elementFormDefault="qualified">
 						<xs:element name="Commodity_Mode_Data" type="xs:string" default="None"/>
 						<xs:element name="Commodity_Density_Data" type="xs:string" default="None" minOccurs="0"/>
 						<xs:element name="Default_Units_Solid_Phase" type="xs:string" default="tonnes"/>
-						<xs:element name="Default_Units_Liquid_Phase" type="xs:string" default="kgal"/>
+						<xs:element name="Default_Units_Liquid_Phase" type="xs:string" default="thousand_gallon"/>
 					</xs:sequence>
 				</xs:complexType>
 			</xs:element>
@@ -40,11 +40,11 @@ elementFormDefault="qualified">
 						<xs:element name="Railcar_Load_Solid" type="xs:string" default="82 tonne"/>
 						<xs:element name="Barge_Load_Solid" type="xs:string" default="700 tonne"/>
 
-						<xs:element name="Truck_Load_Liquid" type="xs:string" default="8 kgal"/>
-						<xs:element name="Railcar_Load_Liquid" type="xs:string" default="28.5 kgal"/>
-						<xs:element name="Barge_Load_Liquid" type="xs:string" default="2100 kgal"/>
-						<xs:element name="Pipeline_Crude_Load_Liquid" type="xs:string" default="3150 kgal"/>
-						<xs:element name="Pipeline_Prod_Load_Liquid" type="xs:string" default="3150 kgal"/>
+						<xs:element name="Truck_Load_Liquid" type="xs:string" default="8 thousand_gallon"/>
+						<xs:element name="Railcar_Load_Liquid" type="xs:string" default="28.5 thousand_gallon"/>
+						<xs:element name="Barge_Load_Liquid" type="xs:string" default="2100 thousand_gallon"/>
+						<xs:element name="Pipeline_Crude_Load_Liquid" type="xs:string" default="3150 thousand_gallon"/>
+						<xs:element name="Pipeline_Prod_Load_Liquid" type="xs:string" default="3150 thousand_gallon"/>
 
 						<xs:element name="Truck_Fuel_Efficiency" type="xs:string" default = "7.4 mi/gal"/>
 						<xs:element name="Rail_Fuel_Efficiency" type="xs:string" default = "10.15 mi/gal"/>
@@ -67,7 +67,7 @@ elementFormDefault="qualified">
 							</xs:simpleType>
 						</xs:element>
 
-						<xs:element name="Density_Conversion_Factor" type="xs:string" default = "3.33 ton/kgal" minOccurs="0"/>
+						<xs:element name="Density_Conversion_Factor" type="xs:string" default = "3.33 ton/thousand_gallon" minOccurs="0"/>
 					</xs:sequence>
 				</xs:complexType>
 			</xs:element>
@@ -87,7 +87,7 @@ elementFormDefault="qualified">
 															<xs:element name="Railroad">
 																<xs:complexType>
 																	<xs:sequence>
-																		<xs:element name="liquid_Railroad_Class_I_Cost" type="xs:string" default = "0.14 usd/kgal/mi"/>
+																		<xs:element name="liquid_Railroad_Class_I_Cost" type="xs:string" default = "0.14 usd/thousand_gallon/mi"/>
 																		<xs:element name="solid_Railroad_Class_I_Cost" type="xs:string" default = "0.047 usd/tonne/mi"/>
 																		<xs:element name="Rail_Density_Code_7_Weight" default = "1.0">
 																			<xs:simpleType>
@@ -152,7 +152,7 @@ elementFormDefault="qualified">
 															<xs:element name="Truck">
 																<xs:complexType>
 																	<xs:sequence>
-																		<xs:element name="liquid_Truck_Base_Cost" type="xs:string"  default = "0.66 usd/kgal/mi"/>
+																		<xs:element name="liquid_Truck_Base_Cost" type="xs:string"  default = "0.66 usd/thousand_gallon/mi"/>
 																		<xs:element name="solid_Truck_Base_Cost" type="xs:string" default = "0.22 usd/tonne/mi"/>
 																		<xs:element name="Truck_Interstate_Weight" default = "1.0">
 																			<xs:simpleType>
@@ -189,7 +189,7 @@ elementFormDefault="qualified">
 															<xs:element name="Barge">
 																<xs:complexType>
 																	<xs:sequence>
-																		<xs:element name="liquid_Barge_cost" type="xs:string" default = "0.097 usd/kgal/mi"/>
+																		<xs:element name="liquid_Barge_cost" type="xs:string" default = "0.097 usd/thousand_gallon/mi"/>
 																		<xs:element name="solid_Barge_cost" type="xs:string"  default = "0.032 usd/tonne/mi"/>
 																		<xs:element name="Water_High_Volume_Weight" default = "1.00">
 																			<xs:simpleType>
@@ -229,7 +229,7 @@ elementFormDefault="qualified">
 												<xs:element name="Intermodal_Transloading_Costs">
 													<xs:complexType>
 														<xs:sequence>
-															<xs:element name="liquid_Transloading_Cost" type="xs:string" default = "40.00 usd/kgal"/>
+															<xs:element name="liquid_Transloading_Cost" type="xs:string" default = "40.00 usd/thousand_gallon"/>
 															<xs:element name="solid_Transloading_Cost" type="xs:string" default = "12.35 usd/tonne"/>
 														</xs:sequence>
 													</xs:complexType>
@@ -253,9 +253,9 @@ elementFormDefault="qualified">
 									<xs:element name="Short_Haul_Penalties">
 										<xs:complexType>
 											<xs:sequence>
-												<xs:element name="liquid_Rail_Short_Haul_Penalty" type="xs:string" default = "71.8 usd/kgal"/>
+												<xs:element name="liquid_Rail_Short_Haul_Penalty" type="xs:string" default = "71.8 usd/thousand_gallon"/>
 												<xs:element name="solid_Rail_Short_Haul_Penalty" type="xs:string" default = "23.9 usd/tonne"/>
-												<xs:element name="liquid_Water_Short_Haul_Penalty" type="xs:string" default = "76.1 usd/kgal"/>
+												<xs:element name="liquid_Water_Short_Haul_Penalty" type="xs:string" default = "76.1 usd/thousand_gallon"/>
 												<xs:element name="solid_Water_Short_Haul_Penalty" type="xs:string" default = "25.4 usd/tonne"/>
 											</xs:sequence>
 										</xs:complexType>
diff --git a/program/lib/v6_temp_Scenario.xml b/program/lib/v6_temp_Scenario.xml
index e33aa84..3e76ea2 100644
--- a/program/lib/v6_temp_Scenario.xml
+++ b/program/lib/v6_temp_Scenario.xml
@@ -1,7 +1,7 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <!-- Document created with online XML Editor http://xmlgrid.net 2018/11/09 15:39:47  -->
 <Scenario xmlns="Schema_v6.0.0">
-    <Scenario_Schema_Version>6.0.3</Scenario_Schema_Version>
+    <Scenario_Schema_Version>6.0.4</Scenario_Schema_Version>
     <Scenario_Name>USER INPUT REQUIRED: replace this string with a descriptive Scenario Name</Scenario_Name>
     <Scenario_Description>USER INPUT REQUIRED: replace this string with a Scenario Description</Scenario_Description>
     <Scenario_Inputs>
@@ -29,10 +29,10 @@
         <Commodity_Density_Data>None</Commodity_Density_Data>
         <!-- DEFAULT UNITS FOR PHASE OF MATTER -->
         <!-- FTOT will check that the units are acceptable at runtime-->
-        <!-- The default units are tonnes (solids) and kgals (liquids)-->
+        <!-- The default units are tonnes (solids) and thousand_gallon (liquids)-->
         <!-- For small-scale scenarios, units can be changed to pounds and gallons.-->
         <Default_Units_Solid_Phase>tonnes</Default_Units_Solid_Phase>
-        <Default_Units_Liquid_Phase>kgal</Default_Units_Liquid_Phase>
+        <Default_Units_Liquid_Phase>thousand_gallon</Default_Units_Liquid_Phase>
     </Scenario_Inputs>
     <Assumptions>
         <!-- VEHICLE LOAD QUANTITY AND UNITS -->
@@ -40,11 +40,11 @@
         <Truck_Load_Solid>24 tonne</Truck_Load_Solid>
         <Railcar_Load_Solid>82 tonne</Railcar_Load_Solid>
         <Barge_Load_Solid>700 tonne</Barge_Load_Solid>
-        <Truck_Load_Liquid>8 kgal</Truck_Load_Liquid>
-        <Railcar_Load_Liquid>28.5 kgal</Railcar_Load_Liquid>
-        <Barge_Load_Liquid>2100 kgal</Barge_Load_Liquid>
-        <Pipeline_Crude_Load_Liquid>3150 kgal</Pipeline_Crude_Load_Liquid>
-        <Pipeline_Prod_Load_Liquid>3150 kgal</Pipeline_Prod_Load_Liquid>
+        <Truck_Load_Liquid>8 thousand_gallon</Truck_Load_Liquid>
+        <Railcar_Load_Liquid>28.5 thousand_gallon</Railcar_Load_Liquid>
+        <Barge_Load_Liquid>2100 thousand_gallon</Barge_Load_Liquid>
+        <Pipeline_Crude_Load_Liquid>3150 thousand_gallon</Pipeline_Crude_Load_Liquid>
+        <Pipeline_Prod_Load_Liquid>3150 thousand_gallon</Pipeline_Prod_Load_Liquid>
         <Truck_Fuel_Efficiency>7.4 mi/gal</Truck_Fuel_Efficiency>
         <Rail_Fuel_Efficiency>10.15 mi/gal</Rail_Fuel_Efficiency>
         <Barge_Fuel_Efficiency>5.00 mi/gal</Barge_Fuel_Efficiency>
@@ -63,17 +63,17 @@
         <!--The density conversion factor is used to calculate emissions for liquid commodities on the rail, water, and pipeline modes based on the above emission factors.-->
         <!--The density factor is applied to any liquid commodity not specified in the optional Commodity_Density_Data input file.-->
         <!--If there is no Commodity_Density_Data input file, the density factor will be applied to all liquid commodities.-->
-        <!--Units must be specified in terms of commodity mass per volume. The default value used by FTOT is 3.33 ton/kgal (as a proxy for SAF).-->
-        <!--Example density conversion factors include: gasoline = 3.14 ton/kgal, diesel = 3.52 ton/kgal, ethanol = 0.79 kg/liter, water = 1.00 g/cm^3-->
+        <!--Units must be specified in terms of commodity mass per volume. The default value used by FTOT is 3.33 ton/thousand_gallon (as a proxy for SAF).-->
+        <!--Example density conversion factors include: gasoline = 3.14 ton/thousand_gallon, diesel = 3.52 ton/thousand_gallon, ethanol = 0.79 kg/liter, water = 1.00 g/cm^3-->
         <!--Sources: BP Statistical Review of World Energy, July 2021, and International Energy Agency AMF TCP-->
-        <Density_Conversion_Factor>3.33 ton/kgal</Density_Conversion_Factor>
+        <Density_Conversion_Factor>3.33 ton/thousand_gallon</Density_Conversion_Factor>
     </Assumptions>
     <scriptParameters>
         <Create_Network_Layer_Script>
             <Network_Costs>
                 <Modal_Costs>
                     <Railroad>
-                        <liquid_Railroad_Class_I_Cost>0.14 usd/kgal/mi</liquid_Railroad_Class_I_Cost>
+                        <liquid_Railroad_Class_I_Cost>0.14 usd/thousand_gallon/mi</liquid_Railroad_Class_I_Cost>
                         <solid_Railroad_Class_I_Cost>0.047 usd/tonne/mi</solid_Railroad_Class_I_Cost>
                         <!--In the internal FTOT network, rail impedance categories are based on FRA Rail Density Codes: sensitive, non-public data-->
                         <!--In the public FTOT network, rail impedance categories are based on publicly available rail ownership data and whether or not the rail segment is part of STRACNET-->
@@ -91,7 +91,7 @@
                         <Rail_Density_Code_0_Weight>10.0</Rail_Density_Code_0_Weight>
                     </Railroad>
                     <Truck>
-                        <liquid_Truck_Base_Cost>0.66 usd/kgal/mi</liquid_Truck_Base_Cost>
+                        <liquid_Truck_Base_Cost>0.66 usd/thousand_gallon/mi</liquid_Truck_Base_Cost>
                         <solid_Truck_Base_Cost>0.22 usd/tonne/mi</solid_Truck_Base_Cost>
                         <!--Weights are used to encourage flows on lower functional class roadways (e.g. interstates and highways)-->
                         <!--Truck_Interstate includes FAF Function Class 1-->
@@ -105,7 +105,7 @@
                         <Truck_Local_Weight>1.3</Truck_Local_Weight>
                     </Truck>
                     <Barge>
-                        <liquid_Barge_cost>0.097 usd/kgal/mi</liquid_Barge_cost>
+                        <liquid_Barge_cost>0.097 usd/thousand_gallon/mi</liquid_Barge_cost>
                         <solid_Barge_cost>0.032 usd/tonne/mi</solid_Barge_cost>
                         <!--Water impedances based on categorizing publicly available USACE freight volumes over the national waterway network -->
                         <!--These weights are used to encourage flows portions of the waterway network that see higher freight use-->
@@ -122,7 +122,7 @@
                     </Barge>
                 </Modal_Costs>
                 <Intermodal_Transloading_Costs>
-                    <liquid_Transloading_Cost>40.00 usd/kgal</liquid_Transloading_Cost>
+                    <liquid_Transloading_Cost>40.00 usd/thousand_gallon</liquid_Transloading_Cost>
                     <solid_Transloading_Cost>12.35 usd/tonne</solid_Transloading_Cost>
                 </Intermodal_Transloading_Costs>
             </Network_Costs>
@@ -138,9 +138,9 @@
             </Artificial_Links>
             <Short_Haul_Penalties>
                 <!--The following are penalties applied to pairs of artificial links to prevent short movements on the rail or water network-->
-                <liquid_Rail_Short_Haul_Penalty>71.8 usd/kgal</liquid_Rail_Short_Haul_Penalty>
+                <liquid_Rail_Short_Haul_Penalty>71.8 usd/thousand_gallon</liquid_Rail_Short_Haul_Penalty>
                 <solid_Rail_Short_Haul_Penalty>23.9 usd/tonne</solid_Rail_Short_Haul_Penalty>
-                <liquid_Water_Short_Haul_Penalty>76.1 usd/kgal</liquid_Water_Short_Haul_Penalty>
+                <liquid_Water_Short_Haul_Penalty>76.1 usd/thousand_gallon</liquid_Water_Short_Haul_Penalty>
                 <solid_Water_Short_Haul_Penalty>25.4 usd/tonne</solid_Water_Short_Haul_Penalty>
             </Short_Haul_Penalties>
         </Create_Network_Layer_Script>
diff --git a/program/lib/vehicle_types.csv b/program/lib/vehicle_types.csv
index 26b3b0b..60f9cf0 100644
--- a/program/lib/vehicle_types.csv
+++ b/program/lib/vehicle_types.csv
@@ -1,6 +1,6 @@
 vehicle_label,mode,vehicle_property,value
 small_truck,road,Truck_Load_Solid,8 tonne
-small_truck,road,Truck_Load_Liquid,2.5 kgal
+small_truck,road,Truck_Load_Liquid,2.5 thousand_gallon
 small_truck,road,Truck_Fuel_Efficiency,12.1 mi/gal
 small_truck,road,Atmos_CO2_Urban_Unrestricted,957.20 g/mi
 small_truck,road,Atmos_CO2_Urban_Restricted,780.84 g/mi

From 44241825655b6dd84fdf33bee8ae0ae81ad77d59 Mon Sep 17 00:00:00 2001
From: Kevin Zhang <68650998+kzhang81@users.noreply.github.com>
Date: Fri, 13 Jan 2023 12:08:15 -0500
Subject: [PATCH 4/4] Add files via upload

---
 program/tools/ftot_tools.py            |  6 ++---
 program/tools/scenario_compare_tool.py | 36 ++++++++++++++++++++++----
 2 files changed, 34 insertions(+), 8 deletions(-)

diff --git a/program/tools/ftot_tools.py b/program/tools/ftot_tools.py
index 9ff7a75..62fddbf 100644
--- a/program/tools/ftot_tools.py
+++ b/program/tools/ftot_tools.py
@@ -8,9 +8,9 @@
 import network_disruption_tool
 from six.moves import input
 
-FTOT_VERSION = "2022.3"
-SCHEMA_VERSION = "6.0.3"
-VERSION_DATE = "9/22/2022"
+FTOT_VERSION = "2022.4"
+SCHEMA_VERSION = "6.0.4"
+VERSION_DATE = "1/13/2023"
 
 header = "\n\
  _______  _______  _______  _______    _______  _______  _______  ___      _______ \n\
diff --git a/program/tools/scenario_compare_tool.py b/program/tools/scenario_compare_tool.py
index 9d2532b..6132e3b 100644
--- a/program/tools/scenario_compare_tool.py
+++ b/program/tools/scenario_compare_tool.py
@@ -46,15 +46,22 @@ def scenario_compare_prep():
                                    out_name=output_gdb_name, out_version="CURRENT")
     output_gdb = os.path.join(output_dir, output_gdb_name + ".gdb")
 
-    # create output csv
+    # create output report csv
     print("start: Tableau results report")
     report_file_name = 'tableau_report.csv'
     report_file = os.path.join(output_dir, report_file_name)
-    # write the first header line
     wf = open(report_file, 'w')
     header_line = 'scenario_name, table_name, commodity, facility_name, measure, mode, value, units, notes\n'
     wf.write(header_line)
 
+    # create output all_routes csv
+    routes_file_name = 'all_routes.csv'
+    routes_file = os.path.join(output_dir, routes_file_name)
+    rf = open(routes_file, 'w')
+    routes_header = 'scenario_name,route_id,from_facility,from_facility_type,to_facility,to_facility_type,commodity_name,phase,mode,dollar_cost,routing_cost,miles,in_solution\n'
+    rf.write(routes_header)
+    isEmpty=True # track if anything written to all_routes
+
     # get user directory to search
     # returns list of dirs
     scenario_dirs = get_input_dirs()
@@ -116,7 +123,20 @@ def scenario_compare_prep():
             wf.write(line)
         csv_in.close()
 
-        # unzip gdb.zip locally
+        # concat all_routes.csv
+        print("look at routes csv and import")
+        rf_in = open(os.path.join(temp_folder,"all_routes.csv"))
+        for line in rf_in:
+            if line == 'Run scenario with NDR on to view this dashboard.':
+                continue # NDR off. No content.
+            elif line.startswith(routes_header):
+                continue
+            else:
+                rf.write(line)
+                isEmpty = False
+        rf_in.close()
+
+         # unzip gdb.zip locally
         with zipfile.ZipFile(os.path.join(temp_folder, 'tableau_output.gdb.zip'),
                          'r') as zipObj:
             # extract all the contents of zip file in current directory
@@ -146,6 +166,12 @@ def scenario_compare_prep():
     # close merged csv
     wf.close()
 
+    # add text line to routes report if empty, then close
+    if isEmpty:
+        line = 'Run scenario with NDR on to view this dashboard.'
+        rf.write(line)
+    rf.close()
+
     # package up the concat files into a compare.twbx
     # create the zip file for writing compressed data
 
@@ -166,8 +192,8 @@ def scenario_compare_prep():
 
     # package the workbook
     # need to specify the archive name parameter to avoid the whole path to the file being added to the archive
-    file_list = ["tableau_dashboard.twb", "tableau_output.gdb.zip", "tableau_report.csv", "volpeTriskelion.gif",
-                 "parameters_icon.png"]
+    file_list = ["tableau_dashboard.twb", "tableau_output.gdb.zip", "tableau_report.csv", "all_routes.csv",
+                 "volpeTriskelion.gif", "parameters_icon.png"]
 
     for a_file in file_list: