diff --git a/dev/timing.py b/dev/timing.py
index bfd8024..68c6386 100644
--- a/dev/timing.py
+++ b/dev/timing.py
@@ -1,51 +1,43 @@
 from pathways import Pathways
+
 p = Pathways(datapackage="/Users/romain/GitHub/premise/dev/image-SSP2/datapackage.json")
 
 for scenario in [
-    #"SSP2-Base",
+    # "SSP2-Base",
     "SSP2-RCP19"
 ]:
     p.calculate(
         methods=[
-            'EF v3.1 - acidification - accumulated exceedance (AE)',
-            'EF v3.1 - climate change - global warming potential (GWP100)',
-            'EF v3.1 - ecotoxicity: freshwater - comparative toxic unit for ecosystems (CTUe)',
-            'EF v3.1 - energy resources: non-renewable - abiotic depletion potential (ADP): fossil fuels',
-            'EF v3.1 - eutrophication: freshwater - fraction of nutrients reaching freshwater end compartment (P)',
-            'EF v3.1 - human toxicity: carcinogenic - comparative toxic unit for human (CTUh)',
-            'EF v3.1 - material resources: metals/minerals - abiotic depletion potential (ADP): elements (ultimate reserves)',
-            'EF v3.1 - particulate matter formation - impact on human health',
-            'EF v3.1 - water use - user deprivation potential (deprivation-weighted water consumption)',
-            'RELICS - metals extraction - Aluminium',
-            'RELICS - metals extraction - Cobalt',
-            'RELICS - metals extraction - Copper',
-            'RELICS - metals extraction - Graphite',
-            'RELICS - metals extraction - Lithium',
-            'RELICS - metals extraction - Molybdenum',
-            'RELICS - metals extraction - Neodymium',
-            'RELICS - metals extraction - Nickel',
-            'RELICS - metals extraction - Platinum',
-            'RELICS - metals extraction - Vanadium',
-            'RELICS - metals extraction - Zinc',
+            "EF v3.1 - acidification - accumulated exceedance (AE)",
+            "EF v3.1 - climate change - global warming potential (GWP100)",
+            "EF v3.1 - ecotoxicity: freshwater - comparative toxic unit for ecosystems (CTUe)",
+            "EF v3.1 - energy resources: non-renewable - abiotic depletion potential (ADP): fossil fuels",
+            "EF v3.1 - eutrophication: freshwater - fraction of nutrients reaching freshwater end compartment (P)",
+            "EF v3.1 - human toxicity: carcinogenic - comparative toxic unit for human (CTUh)",
+            "EF v3.1 - material resources: metals/minerals - abiotic depletion potential (ADP): elements (ultimate reserves)",
+            "EF v3.1 - particulate matter formation - impact on human health",
+            "EF v3.1 - water use - user deprivation potential (deprivation-weighted water consumption)",
+            "RELICS - metals extraction - Aluminium",
+            "RELICS - metals extraction - Cobalt",
+            "RELICS - metals extraction - Copper",
+            "RELICS - metals extraction - Graphite",
+            "RELICS - metals extraction - Lithium",
+            "RELICS - metals extraction - Molybdenum",
+            "RELICS - metals extraction - Neodymium",
+            "RELICS - metals extraction - Nickel",
+            "RELICS - metals extraction - Platinum",
+            "RELICS - metals extraction - Vanadium",
+            "RELICS - metals extraction - Zinc",
         ],
-        regions=[r for r in p.scenarios.coords["region"].values if r!="World"],
+        regions=[r for r in p.scenarios.coords["region"].values if r != "World"],
         scenarios=[scenario],
-        years=[
-            2020,
-            2030,
-            2040,
-            2050,
-            2060,
-            2070,
-            2080,
-            2090,
-            2100
-        ],
+        years=[2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100],
         variables=[
-            v for v in p.scenarios.coords["variables"].values
+            v
+            for v in p.scenarios.coords["variables"].values
             if any(i in v for i in ["Industry", "Transport", "Heating"])
         ],
-        demand_cutoff=0.01
+        demand_cutoff=0.01,
     )
     arr = p.display_results(cutoff=0.0001)
 
diff --git a/pathways/lca.py b/pathways/lca.py
index 719748f..94f62f8 100644
--- a/pathways/lca.py
+++ b/pathways/lca.py
@@ -126,7 +126,11 @@ def load_matrix_and_index(
 
 
 def get_lca_matrices(
-    datapackage: str, model: str, scenario: str, year: int, data_type: np.dtype = np.float32
+    datapackage: str,
+    model: str,
+    scenario: str,
+    year: int,
+    data_type: np.dtype = np.float32,
 ) -> Tuple[sparse.csr_matrix, sparse.csr_matrix, Dict, Dict]:
     """
     Retrieve Life Cycle Assessment (LCA) matrices from disk.
@@ -146,9 +150,15 @@ def get_lca_matrices(
     A_inds = read_indices_csv(dirpath / "A_matrix_index.csv")
     B_inds = read_indices_csv(dirpath / "B_matrix_index.csv")
 
-    A = load_matrix_and_index(dirpath / "A_matrix.csv", len(A_inds), transpose=True, data_type=data_type)
+    A = load_matrix_and_index(
+        dirpath / "A_matrix.csv", len(A_inds), transpose=True, data_type=data_type
+    )
     B = load_matrix_and_index(
-        dirpath / "B_matrix.csv", len(B_inds), sign=-1, extra_indices=len(A_inds), data_type=data_type
+        dirpath / "B_matrix.csv",
+        len(B_inds),
+        sign=-1,
+        extra_indices=len(A_inds),
+        data_type=data_type,
     )
 
     return A, B, A_inds, B_inds
diff --git a/pathways/pathways.py b/pathways/pathways.py
index 81b75b0..32d16e7 100644
--- a/pathways/pathways.py
+++ b/pathways/pathways.py
@@ -12,11 +12,11 @@
 
 import numpy as np
 import pandas as pd
+import pyprind
 import xarray as xr
 import yaml
 from datapackage import DataPackage
 from premise.geomap import Geomap
-import pyprind
 
 from . import DATA_DIR
 from .lca import (
@@ -32,9 +32,9 @@
     create_lca_results_array,
     display_results,
     get_unit_conversion_factors,
+    harmonize_units,
     load_classifications,
     load_units_conversion,
-    harmonize_units
 )
 
 # if pypardiso is installed, use it
diff --git a/pathways/utils.py b/pathways/utils.py
index 121e255..17cfa7d 100644
--- a/pathways/utils.py
+++ b/pathways/utils.py
@@ -117,17 +117,15 @@ def harmonize_units(scenario: xr.DataArray, variables: list) -> xr.DataArray:
             # create vector of conversion factors
             conversion_factors = np.array([1e-3 if u == "PJ/yr" else 1 for u in units])
             # multiply scenario by conversion factors
-            scenario.loc[
-                dict(variables=variables)
-            ] *= conversion_factors[:, np.newaxis, np.newaxis]
+            scenario.loc[dict(variables=variables)] *= conversion_factors[
+                :, np.newaxis, np.newaxis
+            ]
             # update units
             scenario.attrs["units"] = {var: "EJ/yr" for var in variables}
 
-
     return scenario
 
 
-
 def get_unit_conversion_factors(
     scenario_unit: dict, dataset_unit: list, unit_mapping: dict
 ) -> np.ndarray: