openIO.py

"""
Class that creates symmetric Input-Output tables for Canada and its provinces
author: maxime.agez@polymtl.ca
"""

import pandas as pd
import numpy as np
import re
import pkg_resources
import os
import pymrio
import json
import country_converter as coco
import logging
import warnings
import gzip
import pickle
import math


class IOTables:
    def __init__(self, folder_path, exiobase_folder, endogenizing_capitals=True):
        """
        :param folder_path: [string] the path to the folder with the economic data (e.g. /../Detail level/)
        :param exiobase_folder: [string] path to exiobase folder for international imports
        :param endogenizing_capitals: [boolean] True if you want to endogenize capitals
        :param aggregated_ghgs: [boolean] True to work with aggregated GHG physical flow accounts. False requires you to
        have access to the disaggregated files provided by StatCan.
        """

        # ignoring some warnings
        warnings.filterwarnings(action='ignore', category=FutureWarning)
        warnings.filterwarnings(action='ignore', category=np.VisibleDeprecationWarning)
        warnings.filterwarnings(action='ignore', category=pd.errors.PerformanceWarning)

        # set up logging tool
        logger = logging.getLogger('openIO-Canada')
        logger.setLevel(logging.INFO)
        logger.handlers = []
        formatter = logging.Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s")
        ch = logging.StreamHandler()
        ch.setLevel(logging.DEBUG)
        ch.setFormatter(formatter)
        logger.addHandler(ch)
        logger.propagate = False
        # set up logging tool for country_converter
        coco_logger = coco.logging.getLogger()
        coco_logger.setLevel(logging.CRITICAL)

        logger.info('Reading all the Excel files...')

        self.level_of_detail = [i for i in os.path.normpath(folder_path).split(os.sep) if 'level' in i][-1]
        self.exiobase_folder = exiobase_folder
        self.endogenizing = endogenizing_capitals
        self.folder_path = folder_path

        # values
        self.V = pd.DataFrame()
        self.U = pd.DataFrame()
        self.A = pd.DataFrame()
        self.K = pd.DataFrame()
        self.Z = pd.DataFrame()
        self.W = pd.DataFrame()
        self.R = pd.DataFrame()
        self.Y = pd.DataFrame()
        self.WY = pd.DataFrame()
        self.g = pd.DataFrame()
        self.inv_g = pd.DataFrame()
        self.q = pd.DataFrame()
        self.inv_q = pd.DataFrame()
        self.F = pd.DataFrame()
        self.S = pd.DataFrame()
        self.FY = pd.DataFrame()
        self.C = pd.DataFrame()
        self.INT_imports = pd.DataFrame()
        self.L = pd.DataFrame()
        self.E = pd.DataFrame()
        self.D = pd.DataFrame()
        self.who_uses_int_imports_U = pd.DataFrame()
        self.who_uses_int_imports_K = pd.DataFrame()
        self.who_uses_int_imports_Y = pd.DataFrame()
        self.A_exio = pd.DataFrame()
        self.Z_exio = pd.DataFrame()
        self.x_exio = pd.DataFrame()
        self.Y_exio = pd.DataFrame()
        self.K_exio = pd.DataFrame()
        self.S_exio = pd.DataFrame()
        self.F_exio = pd.DataFrame()
        self.C_exio = pd.DataFrame()
        self.link_openio_exio_A = pd.DataFrame()
        self.link_openio_exio_K = pd.DataFrame()
        self.link_openio_exio_Y = pd.DataFrame()
        self.imports = pd.DataFrame()
        self.imports_scaled_U = pd.DataFrame()
        self.imports_scaled_K = pd.DataFrame()
        self.imports_scaled_Y = pd.DataFrame()
        self.minerals = pd.DataFrame()
        self.emission_factors_basic_price = pd.DataFrame()
        self.emission_factors_purchaser_price = pd.DataFrame()

        # metadata
        self.emission_metadata = pd.DataFrame()
        self.unit_exio = pd.DataFrame()
        self.methods_metadata = pd.DataFrame()
        self.industries = []
        self.commodities = []
        self.factors_of_production = []

        self.matching_dict = {'AB': 'Alberta',
                              'BC': 'British Columbia',
                              'CE': 'Canadian territorial enclaves abroad',
                              'MB': 'Manitoba',
                              'NB': 'New Brunswick',
                              'NL': 'Newfoundland and Labrador',
                              'NS': 'Nova Scotia',
                              'NT': 'Northwest Territories',
                              'NU': 'Nunavut',
                              'ON': 'Ontario',
                              'PE': 'Prince Edward Island',
                              'QC': 'Quebec',
                              'SK': 'Saskatchewan',
                              'YT': 'Yukon'}

        files = [i for i in os.walk(folder_path)]
        files = [i for i in files[0][2] if i[:2] in self.matching_dict.keys() and 'SUT' in i]
        files.sort()
        self.year = int(files[0].split('SUT_C')[1].split('_')[0])

        try:
            self.NPRI = pd.read_excel(pkg_resources.resource_stream(
                __name__, '/Data/Environmental_data/NPRI-INRP_DataDonnées_' + str(self.year) + '.xlsx'), None)
            self.NPRI_file_year = self.year
        # 2018 by default (for older years)
        except FileNotFoundError:
            self.NPRI = pd.read_excel(pkg_resources.resource_stream(
                __name__, '/Data/Environmental_data/NPRI-INRP_DataDonnées_2018.xlsx'), None)
            self.NPRI_file_year = 2018

        logger.info("Formatting the Supply and Use tables...")
        for province_data in files:
            supply = pd.read_excel(os.path.join(folder_path, province_data), 'Supply')
            use = pd.read_excel(os.path.join(folder_path, province_data), 'Use_Basic')
            region = province_data[:2]
            self.format_tables(supply, use, region)

        self.W = self.W.fillna(0)
        self.WY = self.WY.fillna(0)
        self.Y = self.Y.fillna(0)
        self.q = self.q.fillna(0)
        self.g = self.g.fillna(0)
        self.U = self.U.fillna(0)
        self.V = self.V.fillna(0)

        # output of the industries (g) should match between V and U+W
        assert np.isclose((self.U.sum()+self.W.sum()).sum(), self.g.sum().sum())
        # output of the commodities (q) should match between V and U+Y adjusted for international imports
        assert np.isclose((self.U.sum(1) + self.Y.sum(1) - self.INT_imports.sum(1)).sum(), self.q.sum(1).sum())

        logger.info("Modifying names of duplicated sectors...")
        self.dealing_with_duplicated_names()

        logger.info('Organizing final demand sectors...')
        self.organize_final_demand()

        logger.info('Removing IOIC codes from index...')
        self.remove_codes()

        if self.endogenizing:
            logger.info('Endogenizing capitals of OpenIO-Canada...')
            self.endogenizing_capitals()

        logger.info("Balancing inter-provincial trade...")
        self.province_import_export(
            pd.read_excel(
                os.path.join(folder_path, [i for i in [j for j in os.walk(folder_path)][0][2] if 'Provincial_trade_flow' in i][0]),
                'Data'))

        logger.info("Loading Exiobase...")
        self.load_exiobase()

        logger.info('Treatment of international trade data...')
        self.determine_sectors_importing()
        self.load_and_correct_international_trade_data()
        logger.info("Linking international trade data to openIO-Canada...")
        self.link_imports_to_openio()

        logger.info("Building the symmetric tables...")
        self.gimme_symmetric_iot()

        logger.info("Linking openIO-Canada to Exiobase...")
        self.link_international_trade_data_to_exiobase()

        self.remove_abroad_enclaves()

        self.concatenate_matrices()

        logger.info("Extracting and formatting environmental data from the NPRI file...")
        self.extract_environmental_data()

        logger.info("Matching emission data from NPRI to IOT sectors...")
        self.match_npri_data_to_iots()

        logger.info("Matching GHG accounts to IOT sectors...")
        self.match_ghg_accounts_to_iots()

        logger.info("Matching water accounts to IOT sectors...")
        self.match_water_accounts_to_iots()

        logger.info("Matching energy accounts to IOT sectors...")
        self.match_energy_accounts_to_iots()

        logger.info("Matching mineral extraction data to IOT sectors...")
        self.match_mineral_extraction_to_iots()

        logger.info("Creating the characterization matrix...")
        self.characterization_matrix()

        logger.info("Refining the GHG emissions for the agriculture sector...")
        self.better_distribution_for_agriculture_ghgs()

        logger.info("Cleaning province and country names...")
        self.differentiate_country_names_openio_exio()

        logger.info("Refining the GHG emissions for the meat sector...")
        self.refine_meat_sector()

        self.convert_F_to_commodity()

        logger.info("Adding HFP and PFC flows...")
        self.add_hfc_emissions()

        logger.info("Refining water consumption of livestock and crops...")
        self.add_water_consumption_flows_for_livestock_and_crops()

        logger.info("Adding plastic waste flows...")
        self.add_plastic_emissions()

        logger.info("Normalizing emissions...")
        self.normalize_flows()

        logger.info("Differentiating biogenic from fossil CO2 emissions...")
        self.differentiate_biogenic_carbon_emissions()

        logger.info("Done extracting openIO-Canada!")

    def format_tables(self, supply, use, region):
        """
        Extracts the relevant dataframes from the Excel files in the Stat Can folder
        :param supply: the supply table
        :param use: the use table
        :param region: the province of Canada to compile data for
        :return: self.W, self.WY, self.Y, self.g, self.q, self.V, self.U
        """

        supply_table = supply
        use_table = use

        if self.year in [2014, 2015, 2016, 2017]:
            # starting_line is the line in which the Supply table starts (the first green row)
            starting_line = 11
            # starting_line_values is the line in which the first value appears
            starting_line_values = 16

        elif self.year in [2018, 2019, 2020, 2021, 2022]:
            # starting_line is the line in which the Supply table starts (the first green row)
            starting_line = 3
            # starting_line_values is the line in which the first value appears
            starting_line_values = 7

        if not self.industries:
            for i in range(0, len(supply_table.columns)):
                if supply_table.iloc[starting_line, i] == 'Total':
                    break
                if supply_table.iloc[starting_line, i] not in [np.nan, 'Industries']:
                    # tuple with code + name (need code to deal with duplicate names in detailed levels)
                    self.industries.append((supply_table.iloc[starting_line+1, i],
                                            supply_table.iloc[starting_line, i]))
            # remove fictive sectors
            self.industries = [i for i in self.industries if not re.search(r'^F', i[0])]

        if not self.commodities:
            for i, element in enumerate(supply_table.iloc[:, 0].tolist()):
                if type(element) == str:
                    # identify by their codes
                    if re.search(r'^[M,F,N,G,I,E]\w*\d', element):
                        self.commodities.append((element, supply_table.iloc[i, 1]))
                    elif re.search(r'^P\w*\d', element) or re.search(r'^GVA', element):
                        self.factors_of_production.append((element, supply_table.iloc[i, 1]))

        final_demand = []
        for i in range(0, len(use_table.columns)):
            if use_table.iloc[starting_line, i] in self.matching_dict.values():
                break
            if use_table.iloc[starting_line, i] not in [np.nan, 'Industries']:
                final_demand.append((use_table.iloc[starting_line+1, i],
                                     use_table.iloc[starting_line, i]))

        final_demand = [i for i in final_demand if i not in self.industries and i[1] != 'Total']

        tables = [supply_table, use_table]

        for i, table in enumerate(tables):
            df = table.iloc[starting_line_values - 2:, 2:]
            df.index = list(zip(table.iloc[starting_line_values - 2:, 0].tolist(),
                                table.iloc[starting_line_values - 2:, 1].tolist()))
            df.columns = list(zip(table.iloc[starting_line + 1, 2:].tolist(),
                                  table.iloc[starting_line, 2:].tolist()))
            df = df.astype(str).replace(to_replace=['.'], value='0').astype(float)
            tables[i] = df  # Update the corresponding table in the list

        # Unpack the tables back into their original variables
        supply_table, use_table = tables

        if self.level_of_detail == 'Detail level':
            # tables from k$ to $
            supply_table *= 1000
            use_table *= 1000
        else:
            # tables from M$ to $
            supply_table *= 1000000
            use_table *= 1000000

        # check calculated totals matched displayed totals
        assert np.allclose(use_table.iloc[:, use_table.columns.get_loc(('TOTAL', 'Total'))],
                           use_table.iloc[:, :use_table.columns.get_loc(('TOTAL', 'Total'))].sum(axis=1), atol=1e-5)
        assert np.allclose(supply_table.iloc[supply_table.index.get_loc(('TOTAL', 'Total'))],
                           supply_table.iloc[:supply_table.index.get_loc(('TOTAL', 'Total'))].sum(), atol=1e-5)

        # extract the tables we need
        W = use_table.loc[self.factors_of_production, self.industries]
        W.drop(('GVA', 'Gross value-added at basic prices'), inplace=True)
        Y = use_table.loc[self.commodities, final_demand]
        WY = use_table.loc[self.factors_of_production, final_demand]
        WY.drop(('GVA', 'Gross value-added at basic prices'), inplace=True)
        g = use_table.loc[[('TOTAL', 'Total')], self.industries]
        q = supply_table.loc[self.commodities, [('TOTAL', 'Total')]]
        V = supply_table.loc[self.commodities, self.industries]
        U = use_table.loc[self.commodities, self.industries]
        INT_imports = supply_table.loc[self.commodities, [('INTIM000', 'International imports')]]

        # create multiindex with region as first level
        for matrix in [W, Y, WY, g, q, V, U, INT_imports]:
            matrix.columns = pd.MultiIndex.from_product([[region], matrix.columns]).tolist()
            matrix.index = pd.MultiIndex.from_product([[region], matrix.index]).tolist()

        # concat the region tables with the all the other tables
        self.W = pd.concat([self.W, W])
        self.WY = pd.concat([self.WY, WY])
        self.Y = pd.concat([self.Y, Y])
        self.q = pd.concat([self.q, q])
        self.g = pd.concat([self.g, g])
        self.U = pd.concat([self.U, U])
        self.V = pd.concat([self.V, V])
        self.INT_imports = pd.concat([self.INT_imports, INT_imports])

    def dealing_with_duplicated_names(self):
        """
        IOIC classification has duplicate names, so we rename when it's the case
        :return: updated dataframes
        """

        # reindexing to fix the order of the columns
        self.V = self.V.T.reindex(pd.MultiIndex.from_product([self.matching_dict.keys(), self.industries]).tolist()).T
        self.U = self.U.T.reindex(pd.MultiIndex.from_product([self.matching_dict.keys(), self.industries]).tolist()).T
        self.g = self.g.T.reindex(pd.MultiIndex.from_product([self.matching_dict.keys(), self.industries]).tolist()).T
        self.W = self.W.T.reindex(pd.MultiIndex.from_product([self.matching_dict.keys(), self.industries]).tolist()).T

        if self.level_of_detail in ['Link-1961 level', 'Link-1997 level', 'Detail level']:
            self.industries = [(i[0], i[1] + ' (private)') if re.search(r'^BS61', i[0]) else i for i in
                               self.industries]
            self.industries = [(i[0], i[1] + ' (non-profit)') if re.search(r'^NP61|^NP71', i[0]) else i for i in
                               self.industries]
            self.industries = [(i[0], i[1] + ' (public)') if re.search(r'^GS61', i[0]) else i for i in
                               self.industries]
        if self.level_of_detail in ['Link-1997 level', 'Detail level']:
            self.industries = [(i[0], i[1] + ' (private)') if re.search(r'^BS623|^BS624', i[0]) else i for i in
                               self.industries]
            self.industries = [(i[0], i[1] + ' (non-profit)') if re.search(r'^NP624', i[0]) else i for i in
                               self.industries]
            self.industries = [(i[0], i[1] + ' (public)') if re.search(r'^GS623', i[0]) else i for i in
                               self.industries]

        # applying the change of names to columns
        for df in [self.V, self.U, self.g, self.W]:
            df.columns = pd.MultiIndex.from_product([self.matching_dict.keys(), self.industries]).tolist()

    def organize_final_demand(self):
        """
        Extract the final demand sectors. These will be disaggregated. If you do not want the detail, use
        self.aggregated_final_demand()
        Provincial exports will be included in self.U and are thus excluded from self.Y
        :return: self.Y & self.WY updated
        """

        # dealing with duplicate names of disaggregated final demand sector names separately
        self.Y.columns = [(i[0], (i[1][0], i[1][1] + ' (private)')) if re.search(r'^COB61|^MEB61|^IPB61|^MEBU', i[1][0])
                          else i for i in self.Y.columns]
        self.Y.columns = [(i[0], (i[1][0], i[1][1] + ' (public)')) if re.search(r'^COG61|^MEG61|^IPG61|^MEGU', i[1][0])
                          else i for i in self.Y.columns]
        self.Y.columns = [(i[0], (i[1][0], i[1][1] + ' (non-profit)')) if re.search(r'^MENU', i[1][0])
                          else i for i in self.Y.columns]
        self.WY.columns = [(i[0], (i[1][0], i[1][1] + ' (private)')) if re.search(r'^COB61|^MEB61|^IPB61|^MEBU', i[1][0])
                          else i for i in self.WY.columns]
        self.WY.columns = [(i[0], (i[1][0], i[1][1] + ' (public)')) if re.search(r'^COG61|^MEG61|^IPG61|^MEGU', i[1][0])
                          else i for i in self.WY.columns]
        self.WY.columns = [(i[0], (i[1][0], i[1][1] + ' (non-profit)')) if re.search(r'^MENU', i[1][0])
                          else i for i in self.WY.columns]

        Y = pd.DataFrame()

        fd_households = [i for i in self.Y.columns if re.search(r'^PEC\w*\d', i[1][0])]
        df = self.Y.loc[:, fd_households].copy()
        df.columns = [(i[0], "Household final consumption expenditure", i[1][1]) for i in fd_households]
        Y = pd.concat([Y, df], axis=1)

        fd_npish = [i for i in self.Y.columns if re.search(r'^CEN\w*\d', i[1][0])]
        df = self.Y.loc[:, fd_npish].copy()
        df.columns = [(i[0], i[1][1], 'NPISH') for i in fd_npish]
        Y = pd.concat([Y, df], axis=1)

        fd_gov = [i for i in self.Y.columns if re.search(r'^CEG\w*\d', i[1][0])]
        df = self.Y.loc[:, fd_gov].copy()
        df.columns = [(i[0], "Governments final consumption expenditure", i[1][1]) for i in fd_gov]
        Y = pd.concat([Y, df], axis=1)

        fd_construction = [i for i in self.Y.columns if re.search(r'^CO\w*\d', i[1][0])]
        df = self.Y.loc[:, fd_construction].copy()
        df.columns = [(i[0], "Gross fixed capital formation, Construction", i[1][1]) for i in fd_construction]
        Y = pd.concat([Y, df], axis=1)

        fd_machinery = [i for i in self.Y.columns if re.search(r'^ME\w*\d', i[1][0])]
        df = self.Y.loc[:, fd_machinery].copy()
        df.columns = [(i[0], "Gross fixed capital formation, Machinery and equipment", i[1][1]) for i in fd_machinery]
        Y = pd.concat([Y, df], axis=1)

        fd_ip = [i for i in self.Y.columns if re.search(r'^IP\w[T]*\d', i[1][0])]
        df = self.Y.loc[:, fd_ip].copy()
        df.columns = [(i[0], "Gross fixed capital formation, Intellectual property products", i[1][1]) for i in fd_ip]
        Y = pd.concat([Y, df], axis=1)

        fd_inv = [i for i in self.Y.columns if re.search(r'^INV\w*\d', i[1][0])]
        df = self.Y.loc[:, fd_inv].copy()
        df.columns = [(i[0], "Changes in inventories", i[1][1]) for i in fd_inv]
        Y = pd.concat([Y, df], axis=1)

        fd_int = [i for i in self.Y.columns if re.search(r'^INT\w*', i[1][0])]
        df = self.Y.loc[:, fd_int].copy()
        df.columns = [(i[0], "International exports", i[1][1].split(' ')[1].capitalize()) for i in fd_int]
        Y = pd.concat([Y, df], axis=1)

        self.Y = Y

        WY = pd.DataFrame()

        fd_households = [i for i in self.WY.columns if re.search(r'^PEC\w*\d', i[1][0])]
        df = self.WY.loc[:, fd_households].copy()
        df.columns = [(i[0], "Household final consumption expenditure", i[1][1]) for i in fd_households]
        WY = pd.concat([WY, df], axis=1)

        fd_npish = [i for i in self.WY.columns if re.search(r'^CEN\w*\d', i[1][0])]
        df = self.WY.loc[:, fd_npish].copy()
        df.columns = [(i[0], i[1][1], 'NPISH') for i in fd_npish]
        WY = pd.concat([WY, df], axis=1)

        fd_gov = [i for i in self.WY.columns if re.search(r'^CEG\w*\d', i[1][0])]
        df = self.WY.loc[:, fd_gov].copy()
        df.columns = [(i[0], "Governments final consumption expenditure", i[1][1]) for i in fd_gov]
        WY = pd.concat([WY, df], axis=1)

        fd_construction = [i for i in self.WY.columns if re.search(r'^CO\w*\d', i[1][0])]
        df = self.WY.loc[:, fd_construction].copy()
        df.columns = [(i[0], "Gross fixed capital formation, Construction", i[1][1]) for i in fd_construction]
        WY = pd.concat([WY, df], axis=1)

        fd_machinery = [i for i in self.WY.columns if re.search(r'^ME\w*\d', i[1][0])]
        df = self.WY.loc[:, fd_machinery].copy()
        df.columns = [(i[0], "Gross fixed capital formation, Machinery and equipment", i[1][1]) for i in fd_machinery]
        WY = pd.concat([WY, df], axis=1)

        fd_ip = [i for i in self.WY.columns if re.search(r'^IP\w[T]*\d', i[1][0])]
        df = self.WY.loc[:, fd_ip].copy()
        df.columns = [(i[0], "Gross fixed capital formation, Intellectual property products", i[1][1]) for i in fd_ip]
        WY = pd.concat([WY, df], axis=1)

        fd_inv = [i for i in self.WY.columns if re.search(r'^INV\w*\d', i[1][0])]
        df = self.WY.loc[:, fd_inv].copy()
        df.columns = [(i[0], "Changes in inventories", i[1][1]) for i in fd_inv]
        WY = pd.concat([WY, df], axis=1)

        fd_int = [i for i in self.WY.columns if re.search(r'^INT\w*', i[1][0])]
        df = self.WY.loc[:, fd_int].copy()
        df.columns = [(i[0], "International exports", i[1][1].split(' ')[1].capitalize()) for i in fd_int]
        WY = pd.concat([WY, df], axis=1)

        self.WY = WY

        assert np.isclose((self.U.sum(1) + self.Y.sum(1) - self.INT_imports.sum(1)).sum(), self.q.sum(1).sum())

    def remove_codes(self):
        """
        Removes the IOIC codes from the index to only leave the name.
        :return: Dataframes with the code of the multi-index removed
        """
        # removing the IOIC codes
        for df in [self.W, self.g, self.V, self.U, self.INT_imports]:
            df.columns = [(i[0], i[1][1]) for i in df.columns]
        for df in [self.W, self.Y, self.WY, self.q, self.V, self.U, self.INT_imports]:
            df.index = [(i[0], i[1][1]) for i in df.index]

        # recreating MultiIndexes
        for df in [self.W, self.Y, self.WY, self.g, self.q, self.V, self.U, self.INT_imports]:
            df.index = pd.MultiIndex.from_tuples(df.index)
            df.columns = pd.MultiIndex.from_tuples(df.columns)

        # reordering columns
        reindexed_columns = pd.MultiIndex.from_product([self.matching_dict.keys(), [i[1] for i in self.industries]])
        self.W = self.W.T.reindex(reindexed_columns).T
        self.g = self.g.T.reindex(reindexed_columns).T
        self.V = self.V.T.reindex(reindexed_columns).T
        self.U = self.U.T.reindex(reindexed_columns).T

    def endogenizing_capitals(self):
        """
        Endogenize consumption of fixed capitals (CFC) of openIO-Canada. CFC is not directly available in the supply
        and use tables of StatCan. However, it is available at the total provincial level. For openIO-Canada, we want
        this data per province AND per sector AND we want the detail of what kind of capital good was purchased.
        So how do we redistribute the provincial values that we have to all sectors while keeping the inputs detail?
        We first derive historical gross fixed capital formation over the years 2014 to the latest available year from
        the SUTs. We get the average capital goods built over these last years. This distribution is then scaled to
        the provincial totals. This gives the CFC per sector per province. To get the inputs details (i.e., to know if
        it's a car or a laptop that was bought) we then apply the reference year GFCF (Gross fixed capital formation)
        distribution from the SUTs.

        Note that we excluded Used car and equipment from the GFCF and the CFC and put it in a newly created Changes in
        inventories category.
        We also replace the Gross operating surplus value added category by the Net operating surplus, which is the
        difference between the Gross operating surplus and the estimated CFC. Negative values of Net operating surplus
        simply indicate that the whole sector in that province that year operated at a loss.

        :return:
        """

        # load historical gross fixed capital formation (generated with SUTs from 2014 to 2021)
        hist_gfcf = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Capitals_endogenization/historical_gross_fixed_capital_formation.xlsx'), index_col=0)
        hist_gfcf.index = pd.MultiIndex.from_tuples([eval(i) for i in hist_gfcf.index])
        hist_gfcf.columns = pd.MultiIndex.from_tuples([eval(i) for i in hist_gfcf.columns])
        # change of classification requires some adjustments in the historical gfcf values
        if self.year == 2022:
            hist_gfcf.columns = pd.MultiIndex.from_tuples(
                [(i[0], 'Other Indigenous government services') if i[1] == 'Other aboriginal government services' else i
                 for i in hist_gfcf.columns])
            hist_gfcf.columns = pd.MultiIndex.from_tuples(
                [(i[0], 'Animal production and aquaculture') if i[1] == 'Animal production' else i for i in
                 hist_gfcf.columns])
            hist_gfcf.columns = pd.MultiIndex.from_tuples(
                [(i[0].replace('2111A', '21111'), i[1]) if '2111A' in i[0] else i for i in hist_gfcf.columns])
            hist_gfcf.columns = pd.MultiIndex.from_tuples(
                [(i[0].replace('2111B', '21114'), i[1]) if '2111B' in i[0] else i for i in hist_gfcf.columns])
            hist_gfcf.columns = pd.MultiIndex.from_tuples(
                [(i[0].replace('55100', '55000'), i[1]) if '55100' in i[0] else i for i in hist_gfcf.columns])

            df = pd.concat([hist_gfcf.loc[:, 'COG61000']] * 3, axis=1)
            df.columns = pd.MultiIndex.from_tuples([('COG61110', 'Elementary and secondary schools'),
                                                    ('COG61120', 'Community colleges and C.E.G.E.P.s'),
                                                    ('COG61130', 'Universities')])
            hist_gfcf = hist_gfcf.drop('COG61000', axis=1).join(df)

            df = pd.concat([hist_gfcf.loc[:, 'MEG61000']] * 3, axis=1)
            df.columns = pd.MultiIndex.from_tuples([('MEG61110', 'Elementary and secondary schools'),
                                                    ('MEG61120', 'Community colleges and C.E.G.E.P.s'),
                                                    ('MEG61130', 'Universities')])
            hist_gfcf = hist_gfcf.drop('MEG61000', axis=1).join(df)

            df = pd.concat([hist_gfcf.loc[:, 'IPG61000']] * 3, axis=1)
            df.columns = pd.MultiIndex.from_tuples([('IPG61110', 'Elementary and secondary schools'),
                                                    ('IPG61120', 'Community colleges and C.E.G.E.P.s'),
                                                    ('IPG61130', 'Universities')])
            hist_gfcf = hist_gfcf.drop('IPG61000', axis=1).join(df)

            df = hist_gfcf.loc[:, 'COG62200'].copy()
            df.columns = pd.MultiIndex.from_tuples(
                [('COG62A00', 'Ambulatory health care services and social assistance')])
            hist_gfcf = hist_gfcf.join(df)

            df = hist_gfcf.loc[:, 'MEG62200'].copy()
            df.columns = pd.MultiIndex.from_tuples(
                [('MEG62A00', 'Ambulatory health care services and social assistance')])
            hist_gfcf = hist_gfcf.join(df)

            df = hist_gfcf.loc[:, 'IPG62200'].copy()
            df.columns = pd.MultiIndex.from_tuples(
                [('IPG62A00', 'Ambulatory health care services and social assistance')])
            hist_gfcf = hist_gfcf.join(df)

        # differentiate between public and private capitals
        hist_gfcf.columns = [(i[0], i[1] + ' (private)') if re.search(r'^COB61|^MEB61|^IPB61|^MEBU', i[0])
                             else i for i in hist_gfcf.columns]
        hist_gfcf.columns = [(i[0], i[1] + ' (public)') if re.search(r'^COG61|^MEG61|^IPG61|^MEGU', i[0])
                             else i for i in hist_gfcf.columns]
        hist_gfcf.columns = [(i[0], i[1] + ' (non-profit)') if re.search(r'^MENU', i[0])
                             else i for i in hist_gfcf.columns]
        hist_gfcf.columns = pd.MultiIndex.from_tuples([i for i in hist_gfcf.columns])
        hist_gfcf_public = hist_gfcf.loc[:, [i for i in hist_gfcf.columns if (re.findall(r'^COG', i[0]) or
                                                                              re.findall(r'^MEG', i[0]) or
                                                                              re.findall(r'^IPG', i[0]) or
                                                                              re.findall(r'^CON', i[0]) or
                                                                              re.findall(r'^MEN', i[0]) or
                                                                              re.findall(r'^IPN', i[0]))]]

        hist_gfcf_private = hist_gfcf.loc[:, [i for i in hist_gfcf.columns if not (re.findall(r'^COG', i[0]) or
                                                                                   re.findall(r'^MEG', i[0]) or
                                                                                   re.findall(r'^IPG', i[0]) or
                                                                                   re.findall(r'^CON', i[0]) or
                                                                                   re.findall(r'^MEN', i[0]) or
                                                                                   re.findall(r'^IPN', i[0]))]]
        # load total provincial consumption of fixed capital data from StatCan
        factor_prod = pd.read_csv(pkg_resources.resource_stream(
            __name__, '/Data/Capitals_endogenization/factor_of_production_detail_StatCan.csv'))
        # rename to match with names of embassies in supply and use tables
        factor_prod.loc[factor_prod.GEO == 'Outside Canada', 'GEO'] = 'Canadian territorial enclaves abroad'
        factor_prod = factor_prod[
            (factor_prod.REF_DATE == self.year) & (factor_prod.GEO.isin(list(self.matching_dict.values())))]

        # load mapping connecting IOCC (intermediary sectors) to GFCF categories
        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            endo = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2014-2021/Endogenizing.xlsx'))
        elif self.year in [2022]:
            endo = pd.read_excel(
                pkg_resources.resource_stream(__name__, '/Data/Concordances/2022/Endogenizing.xlsx'))
        self.K = pd.DataFrame(0, self.U.index, self.U.columns, float)
        self.W = pd.concat([self.W, pd.DataFrame(0, index=[(i, 'Net mixed income') for i in self.matching_dict.keys()] +
                                                          [(i, 'Net operating surplus') for i in self.matching_dict.keys()],
                                                 columns=self.W.columns)])
        # get capitals matrix
        for province in self.matching_dict.keys():
            gfcf_public = self.Y.loc(axis=1)[province, ['Gross fixed capital formation, Construction',
                                                        'Gross fixed capital formation, Machinery and equipment',
                                                        'Gross fixed capital formation, Intellectual property products'],
                                             list(set([i[1] for i in hist_gfcf_public.columns]))].groupby(axis=1, level=[0, 2]).sum().copy()

            gfcf_private = self.Y.loc(axis=1)[province, ['Gross fixed capital formation, Construction',
                                                         'Gross fixed capital formation, Machinery and equipment',
                                                         'Gross fixed capital formation, Intellectual property products'],
                                              list(set([i[1] for i in hist_gfcf_private.columns]))].groupby(axis=1, level=[0, 2]).sum().copy()

            gfcf_public = gfcf_public.drop(
                [i for i in gfcf_public.columns if 'Used cars and equipment and scrap' in i[1]], axis=1)
            gfcf_private = gfcf_private.drop(
                [i for i in gfcf_private.columns if 'Used cars and equipment and scrap' in i[1]], axis=1)

            cfc_public = (hist_gfcf_public.loc[province].mean() /
                          hist_gfcf_public.loc[province].mean().sum() *
                          factor_prod.loc[
                              (factor_prod.GEO == self.matching_dict[province]) &
                              (factor_prod.Estimates == 'Consumption of fixed capital: general governments and non-profit institutions serving households'),
                                          'VALUE'].iloc[0] * 1000000).fillna(0)

            cfc_private_corpo = (hist_gfcf_private.loc[province].mean() /
                                 hist_gfcf_private.loc[province].mean().sum() *
                                 factor_prod.loc[(factor_prod.GEO == self.matching_dict[province]) &
                                                 (factor_prod.Estimates == 'Consumption of fixed capital: corporations'),
                                                 'VALUE'].sum() * 1000000).fillna(0)

            cfc_private_not_corpo = (hist_gfcf_private.loc[province].mean() /
                                     hist_gfcf_private.loc[province].mean().sum() *
                                     factor_prod.loc[(factor_prod.GEO == self.matching_dict[province]) &
                                                     (factor_prod.Estimates == 'Consumption of fixed capital: unincorporated businesses'),
                                                     'VALUE'].sum() * 1000000).fillna(0)

            cfc_public.index = [i[1] for i in cfc_public.index]
            cfc_public = cfc_public.groupby(cfc_public.index).sum()
            cfc_private_corpo.index = [i[1] for i in cfc_private_corpo.index]
            cfc_private_corpo = cfc_private_corpo.groupby(cfc_private_corpo.index).sum()
            cfc_private_not_corpo.index = [i[1] for i in cfc_private_not_corpo.index]
            cfc_private_not_corpo = cfc_private_not_corpo.groupby(cfc_private_not_corpo.index).sum()

            # rescale because of historical investments that are not present in year of study
            cfc_public = (cfc_public.loc[
                              [i for i in cfc_public.index if gfcf_public.sum().loc(axis=0)[province, i] != 0]] /
                          cfc_public.loc[[i for i in cfc_public.index if
                                          gfcf_public.sum().loc(axis=0)[province, i] != 0]].sum() * cfc_public.sum())

            cfc_public_rescaled = ((gfcf_public.loc[:, province] / gfcf_public.loc[:, province].sum()).fillna(0) *
                                   cfc_public.reindex([i[1] for i in gfcf_public.columns]).fillna(0))

            cfc_private_corpo = (cfc_private_corpo.loc[[i for i in cfc_private_corpo.index if
                                                        gfcf_private.sum().loc(axis=0)[province, i] != 0]] /
                                 cfc_private_corpo.loc[[i for i in cfc_private_corpo.index if
                                                        gfcf_private.sum().loc(axis=0)[province, i] != 0]].sum() *
                                 cfc_private_corpo.sum())

            cfc_private_corpo_rescaled = (
                        (gfcf_private.loc[:, province] / gfcf_private.loc[:, province].sum()).fillna(0) *
                        cfc_private_corpo.reindex([i[1] for i in gfcf_private.columns]).fillna(0))

            cfc_private_not_corpo = (
                    cfc_private_not_corpo.loc[
                        [i for i in cfc_private_not_corpo.index if gfcf_private.sum().loc(axis=0)[province, i] != 0]] /
                    cfc_private_not_corpo.loc[[i for i in cfc_private_not_corpo.index if
                                               gfcf_private.sum().loc(axis=0)[province, i] != 0]].sum() *
                    cfc_private_not_corpo.sum())

            cfc_private_not_corpo_rescaled = (
                        (gfcf_private.loc[:, province] / gfcf_private.loc[:, province].sum()).fillna(0) *
                        cfc_private_not_corpo.reindex([i[1] for i in gfcf_private.columns]).fillna(0))

            # ensure rescaling works correctly
            assert np.isclose(
                cfc_public_rescaled.sum().sum() + cfc_private_corpo_rescaled.sum().sum() + cfc_private_not_corpo_rescaled.sum().sum(),
                factor_prod.loc[(factor_prod.GEO == self.matching_dict[province]) & (
                    factor_prod.Estimates.isin(['Consumption of fixed capital: corporations',
                                                'Consumption of fixed capital: unincorporated businesses',
                                                'Consumption of fixed capital: general governments and non-profit institutions serving households'])),
                                'VALUE'].sum().sum() * 1000000)

            cfc_rescaled = pd.concat([cfc_public_rescaled, cfc_private_corpo_rescaled, cfc_private_not_corpo_rescaled],
                                     axis=1)
            cfc_rescaled = cfc_rescaled.groupby(axis=1, level=0).sum()

            for capital in cfc_rescaled.columns:
                # total K and total factor_prod don't match because IF condition not respected all the time
                if self.V.loc[:, province].loc[:, endo[endo.Capitals == capital].loc[:, 'IOIC']].sum().sum() != 0:
                    share = (self.V.loc[:, province].loc[:, endo[endo.Capitals == capital].loc[:, 'IOIC']].sum() /
                             self.V.loc[:, province].loc[:, endo[endo.Capitals == capital].loc[:, 'IOIC']].sum().sum())
                    assert np.isclose(share.sum(), 1)
                    for sector in share.index:
                        # add to capital matrix
                        self.K.loc[:, (province, sector)] += (cfc_rescaled.loc[:, capital] * share.loc[sector])

                        # remove from value-added matrix
                        if capital in cfc_private_not_corpo_rescaled.columns.tolist():
                            self.W.loc[(province, 'Net mixed income'), (province, sector)] += (
                                    self.W.loc[(province, 'Gross mixed income'), (province, sector)] -
                                    (cfc_private_not_corpo_rescaled.loc[:, capital] * share.loc[sector]).sum()
                            )

                        self.W.loc[(province, 'Net operating surplus'), (province, sector)] += (
                                self.W.loc[(province, 'Gross operating surplus'), (province, sector)] -
                                (pd.concat([cfc_private_corpo_rescaled, cfc_public_rescaled], axis=1).loc[:, capital] *
                                 share.loc[sector]).sum()
                        )
        # Gotta change Gross indicators for net ones for the residential building construction sector
        self.W.loc[[i for i in self.W.index if i[1] == 'Net mixed income'],
                   [j for j in self.W.columns if j[1] == 'Residential building construction']] = (
            self.W.loc(axis=1)[:, 'Residential building construction'].loc(axis=0)[:, 'Gross mixed income'].values)
        self.W.loc[[i for i in self.W.index if i[1] == 'Net operating surplus'],
                   [j for j in self.W.columns if j[1] == 'Residential building construction']] = (
            self.W.loc(axis=1)[:, 'Residential building construction'].loc(axis=0)[:, 'Gross operating surplus'].values)

        self.W = self.W.drop(['Gross operating surplus', 'Gross mixed income'], axis=0, level=1)

        # check balance still respected
        assert np.isclose((self.U.sum() + self.W.sum() + self.K.sum()), self.g.sum()).all()

        # move Used cars and equipment and scrap categories to changes in inventories
        df = self.Y.loc(axis=1)[:, :, ['Used cars and equipment and scrap (private)',
                                       'Used cars and equipment and scrap (public)',
                                       'Used cars and equipment and scrap (non-profit)']].copy('deep')
        df = df.T.groupby(level=0).sum().T
        df.columns = pd.MultiIndex.from_product([self.matching_dict.keys(), ['Changes in inventories'],
                                                 ['Changes in inventories, used cars and equipment and scrap']])
        self.Y = pd.concat([self.Y, df], axis=1)

        # drop the share of the GFCF that was endogenized, the redisual represents investments of that year
        # first we need to retranslate the endogenized capitals (i.e., go back to the COICOP classification)
        endogenized_capitals = self.K.copy()
        endogenized_capitals.columns = pd.MultiIndex.from_tuples(
            [(i[0], endo.set_index('IOIC').drop('IOIC code', axis=1).loc[i[1], 'Capitals']) if
             i[1] != 'Residential building construction' else i for i in self.K.columns])
        endogenized_capitals = endogenized_capitals.groupby(endogenized_capitals.columns, axis=1).sum()
        endogenized_capitals.columns = pd.MultiIndex.from_tuples(endogenized_capitals.columns)
        # isolate and format construction capitals
        buildings = list(set([i[1] for i in self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Construction'].sum(1)[
            self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Construction'].sum(1) != 0
            ].index.tolist()]))
        buildings = endogenized_capitals.loc(axis=0)[:, buildings].reindex(self.Y.index).fillna(0)
        buildings.columns = pd.MultiIndex.from_tuples(
            [(i[0], 'Gross fixed capital formation, Construction', i[1]) for i in buildings.columns])
        # isolate and format machinery and equipment capitals
        machines = list(
            set([i[1] for i in self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Machinery and equipment'].sum(1)[
                self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Machinery and equipment'].sum(1) != 0
                ].index.tolist()]))
        machines = endogenized_capitals.loc(axis=0)[:, machines].reindex(self.Y.index).fillna(0)
        machines.columns = pd.MultiIndex.from_tuples(
            [(i[0], 'Gross fixed capital formation, Machinery and equipment', i[1]) for i in machines.columns])
        # isolate and format IP capitals
        ip = list(set([i[1] for i in
                       self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Intellectual property products'].sum(1)[
                           self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Intellectual property products'].sum(
                               1) != 0
                           ].index.tolist()]))
        ip = endogenized_capitals.loc(axis=0)[:, ip].reindex(self.Y.index).fillna(0)
        ip.columns = pd.MultiIndex.from_tuples(
            [(i[0], 'Gross fixed capital formation, Intellectual property products', i[1]) for i in ip.columns])

        # do the actual removing
        self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Construction'] -= buildings
        self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Machinery and equipment'] -= machines
        self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Intellectual property products'] -= ip

        # remove negative values, they clearly do not represent an investment of this year
        self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Construction'] = (
            self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Construction'][
                self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Construction'] > 0
                ]
        ).fillna(0)
        self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Machinery and equipment'] = (
            self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Machinery and equipment'][
                self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Machinery and equipment'] > 0
                ]
        ).fillna(0)
        self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Intellectual property products'] = (
            self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Intellectual property products'][
                self.Y.loc(axis=1)[:, 'Gross fixed capital formation, Intellectual property products'] > 0
                ]
        ).fillna(0)

    def province_import_export(self, province_trade_file):
        """
        Method extracting and formatting inter province imports/exports
        :return: modified self.U, self.V, self.W, self.Y
        """

        province_trade_file = province_trade_file

        province_trade_file.Origin = [{v: k for k, v in self.matching_dict.items()}[i.split(') ')[1]] if ')' in i
                                      else i for i in province_trade_file.Origin]
        province_trade_file.Destination = [{v: k for k, v in self.matching_dict.items()}[i.split(') ')[1]] if ')' in i
                                           else i for i in province_trade_file.Destination]
        # extracting and formatting supply for each province
        province_trade = pd.pivot_table(data=province_trade_file, index='Destination', columns=['Origin', 'Product'])

        province_trade = province_trade.loc[
            [i for i in province_trade.index if i in self.matching_dict.keys()],
            [i for i in province_trade.columns if i[1] in self.matching_dict.keys()]]
        if self.level_of_detail == 'Detail level':
            province_trade *= 1000
        else:
            province_trade *= 1000000
        province_trade.columns = [(i[1], i[2].split(': ')[1]) if ':' in i[2] else i for i in
                                  province_trade.columns]
        province_trade.drop([i for i in province_trade.columns if i[1] not in [i[1] for i in self.commodities]],
                            axis=1, inplace=True)
        province_trade.columns = pd.MultiIndex.from_tuples(province_trade.columns)
        for province in province_trade.index:
            province_trade.loc[province, province] = 0

        import_markets = pd.DataFrame(0, province_trade.index, province_trade.columns, float)
        for importing_province in province_trade.index:
            for exported_product in province_trade.columns.levels[1]:
                exported_products = [i for i in import_markets.columns if i[1] == exported_product]
                import_markets.loc[importing_province, exported_products] = (
                            province_trade.loc[importing_province, exported_products] /
                            province_trade.loc[importing_province, exported_products].sum()).values

        before_distribution_U = self.U.copy('deep')
        before_distribution_K = self.K.copy('deep')
        before_distribution_Y = self.Y.copy('deep')

        if not self.endogenizing:
            for importing_province in province_trade.index:
                total_use = pd.concat([self.U.loc[importing_province, importing_province],
                                       self.Y.loc[importing_province, importing_province]], axis=1)
                total_imports = province_trade.T.groupby(level=1).sum().T.loc[importing_province]
                index_commodity = [i[1] for i in self.commodities]
                total_imports = total_imports.reindex(index_commodity).fillna(0)

                import_distribution_U = ((self.U.loc[importing_province, importing_province].T /
                                          (total_use.sum(axis=1))) *
                                         total_imports).T.fillna(0)
                import_distribution_Y = ((self.Y.loc[importing_province, importing_province].T /
                                          (total_use.sum(axis=1))) *
                                         total_imports).T.fillna(0)

                # distribution balance imports to the different exporting regions
                for exporting_province in province_trade.index:
                    if importing_province != exporting_province:
                        scaled_imports_U = ((import_distribution_U.T * import_markets.fillna(0).loc[
                            importing_province, exporting_province]).T).reindex(import_distribution_U.index).fillna(0)
                        scaled_imports_Y = ((import_distribution_Y.T * import_markets.fillna(0).loc[
                            importing_province, exporting_province]).T).reindex(import_distribution_Y.index).fillna(0)

                        self.assert_order(exporting_province, importing_province, scaled_imports_U, scaled_imports_Y)

                        # assign new values into self.U
                        self.U.loc[exporting_province, importing_province] = (
                            scaled_imports_U.loc[:, self.U.columns.levels[1]].reindex(
                                self.U.loc[exporting_province, importing_province].columns, axis=1).values)
                        # assign new values into self.Y
                        self.Y.loc[exporting_province, importing_province] = (
                            scaled_imports_Y.loc[:, self.Y.columns.levels[1]].reindex(
                                self.Y.loc[exporting_province, importing_province].columns, axis=1).values)

                # remove interprovincial from intraprovincial to not double count
                self.U.loc[importing_province, importing_province] = (
                        self.U.loc[importing_province, importing_province] - self.U.loc[
                    [i for i in self.matching_dict if i != importing_province], importing_province].groupby(
                    level=1).sum()).reindex([i[1] for i in self.commodities]).values
                self.Y.loc[importing_province, importing_province] = (
                        self.Y.loc[importing_province, importing_province] - self.Y.loc[
                    [i for i in self.matching_dict if i != importing_province], importing_province].groupby(
                    level=1).sum()).reindex([i[1] for i in self.commodities]).values

                # if some province buys more than they use, drop the value in "changes in inventories"
                df = self.U[self.U < -1].dropna(how='all', axis=1).dropna(how='all')
                if len(df):
                    self.Y.loc[:, (importing_province, 'Changes in inventories',
                                   'Changes in inventories, finished goods and goods in process')] += (
                        df.sum(1).reindex(self.Y.index).fillna(0))
                    self.U.loc[df.index, df.columns] = 0
                # removing negative values lower than 1$ (potential calculation artefacts)
                self.U = self.U[self.U > 0].fillna(0)
                # checking negative values were removed
                assert not self.U[self.U < 0].any().any()

            # check that the distribution went smoothly
            assert np.isclose(self.U.sum().sum()+self.Y.sum().sum(),
                              before_distribution_U.sum().sum()+before_distribution_Y.sum().sum())

        else:
            for importing_province in province_trade.index:
                total_use = pd.concat([self.U.loc[importing_province, importing_province] +
                                       self.K.loc[importing_province, importing_province],
                                       self.Y.loc[importing_province, importing_province]], axis=1)
                total_imports = province_trade.T.groupby(level=1).sum().T.loc[importing_province]
                index_commodity = [i[1] for i in self.commodities]
                total_imports = total_imports.reindex(index_commodity).fillna(0)

                import_distribution_U = ((self.U.loc[importing_province, importing_province].T /
                                          (total_use.sum(axis=1))) *
                                         total_imports).T.fillna(0)
                import_distribution_K = ((self.K.loc[importing_province, importing_province].T /
                                          (total_use.sum(axis=1))) *
                                         total_imports).T.fillna(0)
                import_distribution_Y = ((self.Y.loc[importing_province, importing_province].T /
                                          (total_use.sum(axis=1))) *
                                         total_imports).T.fillna(0)

                # distribution balance imports to the different exporting regions
                for exporting_province in province_trade.index:
                    if importing_province != exporting_province:
                        scaled_imports_U = ((import_distribution_U.T * import_markets.fillna(0).loc[
                            importing_province, exporting_province]).T).reindex(import_distribution_U.index).fillna(0)
                        scaled_imports_K = ((import_distribution_K.T * import_markets.fillna(0).loc[
                            importing_province, exporting_province]).T).reindex(import_distribution_K.index).fillna(0)
                        scaled_imports_Y = ((import_distribution_Y.T * import_markets.fillna(0).loc[
                            importing_province, exporting_province]).T).reindex(import_distribution_Y.index).fillna(0)

                        self.assert_order(exporting_province, importing_province, scaled_imports_U, scaled_imports_Y,
                                          scaled_imports_K)

                        # assign new values into self.U
                        self.U.loc[exporting_province, importing_province] = (
                            scaled_imports_U.loc[:, self.U.columns.levels[1]].reindex(
                                self.U.loc[exporting_province, importing_province].columns, axis=1).values)
                        # assign new values into self.K
                        self.K.loc[exporting_province, importing_province] = (
                            scaled_imports_K.loc[:, self.K.columns.levels[1]].reindex(
                                self.K.loc[exporting_province, importing_province].columns, axis=1).values)
                        # assign new values into self.Y
                        self.Y.loc[exporting_province, importing_province] = (
                            scaled_imports_Y.loc[:, self.Y.columns.levels[1]].reindex(
                                self.Y.loc[exporting_province, importing_province].columns, axis=1).values)

                # remove interprovincial from intraprovincial to not double count
                self.U.loc[importing_province, importing_province] = (
                        self.U.loc[importing_province, importing_province] - self.U.loc[
                    [i for i in self.matching_dict if i != importing_province], importing_province].groupby(
                    level=1).sum()).reindex([i[1] for i in self.commodities]).values
                self.K.loc[importing_province, importing_province] = (
                        self.K.loc[importing_province, importing_province] - self.K.loc[
                    [i for i in self.matching_dict if i != importing_province], importing_province].groupby(
                    level=1).sum()).reindex([i[1] for i in self.commodities]).values
                self.Y.loc[importing_province, importing_province] = (
                        self.Y.loc[importing_province, importing_province] - self.Y.loc[
                    [i for i in self.matching_dict if i != importing_province], importing_province].groupby(
                    level=1).sum()).reindex([i[1] for i in self.commodities]).values

                # if some province buys more than they use, drop the value in "changes in inventories"
                df = self.U[self.U < -1].dropna(how='all', axis=1).dropna(how='all')
                if len(df):
                    self.Y.loc[:, (importing_province, 'Changes in inventories',
                                   'Changes in inventories, finished goods and goods in process')] += (
                        df.sum(1).reindex(self.Y.index).fillna(0))
                    self.U.loc[df.index, df.columns] = 0
                df = self.K[self.K < -1].dropna(how='all', axis=1).dropna(how='all')
                if len(df):
                    self.Y.loc[:, (importing_province, 'Changes in inventories',
                                   'Changes in inventories, finished goods and goods in process')] += (
                        df.sum(1).reindex(self.Y.index).fillna(0))
                    self.K.loc[df.index, df.columns] = 0
                # removing negative values lower than 1$ (potential calculation artefacts)
                self.U = self.U[self.U > 0].fillna(0)
                self.K = self.K[self.K > 0].fillna(0)
                # checking negative values were removed
                assert not self.U[self.U < 0].any().any()
                assert not self.K[self.K < 0].any().any()

                # check that the distribution went smoothly
                assert np.isclose(self.U.sum().sum() + self.K.sum().sum() + self.Y.sum().sum(),
                                  before_distribution_U.sum().sum() +
                                  before_distribution_K.sum().sum() +
                                  before_distribution_Y.sum().sum())

    def load_exiobase(self):

        # loading Exiobase
        io = pymrio.parse_exiobase3(self.exiobase_folder)

        # need to calc to get access to S environmental matrices
        io.calc_all()

        if self.endogenizing:
            with gzip.open(pkg_resources.resource_stream(
                    __name__, '/Data/Capitals_endogenization/K_cfc_pxp_exio3.8.2_' +
                              str(self.year) + '.gz.pickle'), 'rb') as f:
                self.K_exio = pickle.load(f)

        # save the matrices from exiobase because we need them later
        self.A_exio = io.A.copy('deep')
        self.Z_exio = io.Z.copy('deep')
        self.x_exio = io.x.copy('deep')
        self.Y_exio = io.Y.copy('deep')

        try:
            io.satellite.S
            version_exio = '3.8'
        except AttributeError:
            version_exio = '3.9'

        if version_exio == '3.9':
            try:
                self.S_exio = pd.concat(
                    [io.air_emissions.S, io.nutrients.S, io.water.S, io.land.S, io.material.S, io.energy.S, io.employment.S,
                     io.factor_inputs.S])
                self.F_exio = pd.concat(
                    [io.air_emissions.F, io.nutrients.F, io.water.F, io.land.F, io.material.F, io.energy.F, io.employment.F,
                     io.factor_inputs.F])
                self.FY_exio = pd.concat(
                    [io.air_emissions.F_Y, io.nutrients.F_Y, io.water.F_Y, io.land.F_Y, io.material.F_Y, io.energy.F_Y,
                     io.employment.F_Y, io.factor_inputs.F_Y])
                # millions euros to euros for environmental outputs
                self.S_exio.loc[[i for i in self.S_exio.index if i not in (io.employment.S.index.tolist() +
                                                                           io.factor_inputs.S.index.tolist())]] /= 1000000
                self.unit_exio = pd.concat(
                    [io.air_emissions.unit, io.nutrients.unit, io.water.unit, io.land.unit, io.material.unit, io.energy.unit,
                     io.employment.unit, io.factor_inputs.unit])
                self.unit_exio.columns = ['Unit']
            # the 2022 reference year of EXIOBASE has no land and energy modules, unlike all other reference years
            except AttributeError:
                self.S_exio = pd.concat(
                    [io.air_emissions.S, io.nutrients.S, io.water.S, io.material.S, io.employment.S, io.factor_inputs.S])
                self.F_exio = pd.concat(
                    [io.air_emissions.F, io.nutrients.F, io.water.F, io.material.F, io.employment.F, io.factor_inputs.F])
                self.FY_exio = pd.concat(
                    [io.air_emissions.F_Y, io.nutrients.F_Y, io.water.F_Y, io.material.F_Y, io.employment.F_Y,
                     io.factor_inputs.F_Y])
                # millions euros to euros for environmental outputs
                self.S_exio.loc[[i for i in self.S_exio.index if i not in (io.employment.S.index.tolist() +
                                                                           io.factor_inputs.S.index.tolist())]] /= 1000000
                self.unit_exio = pd.concat(
                    [io.air_emissions.unit, io.nutrients.unit, io.water.unit, io.material.unit, io.employment.unit,
                     io.factor_inputs.unit])
                self.unit_exio.columns = ['Unit']

        elif version_exio == '3.8':
            self.S_exio = io.satellite.S.copy()
            self.F_exio = io.satellite.F.copy()
            self.FY_exio = io.satellite.F_Y.copy()
            # millions euros to euros for environmental outputs
            self.S_exio.iloc[9:] /= 1000000
            self.unit_exio = io.satellite.unit.copy()
            self.unit_exio.columns = ['Unit']

        del io

    def determine_sectors_importing(self):
        """
        Determine which sectors use international imports and removing international imports from use
        :return:
        """

        # aggregating international imports in 1 column
        self.INT_imports = self.INT_imports.T.groupby(level=1).sum().T
        # save total values to check if distribution was done properly later
        before_distribution_U = self.U.sum().sum()
        before_distribution_K = self.K.sum().sum()
        before_distribution_Y = self.Y.sum().sum()
        # need to flatten multiindex for the concatenation to work properly
        self.Y.columns = self.Y.columns.tolist()
        self.U.columns = self.U.columns.tolist()
        self.K.columns = self.K.columns.tolist()
        # determine total use
        if not self.endogenizing:
            total_use = pd.concat([self.U, self.Y], axis=1)
        else:
            total_use = pd.concat([self.U + self.K, self.Y], axis=1)
        # weighted average of who is requiring the international imports, based on national use
        self.who_uses_int_imports_U = (self.U.T / total_use.sum(1)).fillna(0).T * self.INT_imports.values
        self.who_uses_int_imports_Y = (self.Y.T / total_use.sum(1)).fillna(0).T * self.INT_imports.values
        if self.endogenizing:
            self.who_uses_int_imports_K = (self.K.T / total_use.sum(1)).fillna(0).T * self.INT_imports.values
        # remove international imports from national use
        self.U = self.U - self.who_uses_int_imports_U.reindex(self.U.columns, axis=1).fillna(0)
        self.Y = self.Y - self.who_uses_int_imports_Y.reindex(self.Y.columns, axis=1).fillna(0)
        if self.endogenizing:
            self.K = self.K - self.who_uses_int_imports_K.reindex(self.K.columns, axis=1).fillna(0)

        # check totals match
        assert np.isclose(self.U.sum().sum() + self.who_uses_int_imports_U.sum().sum(), before_distribution_U)
        assert np.isclose(self.K.sum().sum() + self.who_uses_int_imports_K.sum().sum(), before_distribution_K)
        assert np.isclose(self.Y.sum().sum() + self.who_uses_int_imports_Y.sum().sum(), before_distribution_Y)

        if self.endogenizing:
            # if some province buys more than they use, drop the value in "changes in inventories"
            # also remove artefacts (between -0.001$ and 0$)
            df = self.U[self.U < -1e-3].dropna(how='all', axis=1).dropna(how='all')
            if len(df):
                df.columns = pd.MultiIndex.from_tuples(df.columns)
                for province in df.columns.levels[0]:
                    self.Y.loc[:, [(province, 'Changes in inventories',
                                    'Changes in inventories, finished goods and goods in process')]] += pd.DataFrame(
                        df.loc[:, province].sum(1).reindex(self.Y.index).fillna(0),
                        columns=[(province, 'Changes in inventories',
                                  'Changes in inventories, finished goods and goods in process')]
                    )
                self.U.loc[df.index, df.columns] = 0

            df = self.K[self.K < -1e-3].dropna(how='all', axis=1).dropna(how='all')
            if len(df):
                df.columns = pd.MultiIndex.from_tuples(df.columns)
                for province in df.columns.levels[0]:
                    self.Y.loc[:, [(province, 'Changes in inventories',
                                    'Changes in inventories, finished goods and goods in process')]] += pd.DataFrame(
                        df.loc[:, province].sum(1).reindex(self.Y.index).fillna(0),
                        columns=[(province, 'Changes in inventories',
                                  'Changes in inventories, finished goods and goods in process')]
                    )
                self.K.loc[df.index, df.columns] = 0
        assert self.U[self.U < -1e-3].dropna(how='all', axis=1).dropna(how='all', axis=0).empty
        assert self.K[self.K < -1e-3].dropna(how='all', axis=1).dropna(how='all', axis=0).empty

        # remove negative artefacts
        self.Y = pd.concat([self.Y[self.Y >= 0].fillna(0), self.Y[self.Y < -1e-3].fillna(0)], axis=1)
        self.Y = self.Y.T.groupby(by=self.Y.columns).sum().T
        self.Y.columns = pd.MultiIndex.from_tuples(self.Y.columns)

        assert np.isclose(self.U.sum().sum() + self.K.sum().sum() + self.Y.sum().sum() + self.INT_imports.sum().sum(),
                          before_distribution_U + before_distribution_K + before_distribution_Y)

    def load_and_correct_international_trade_data(self):
        """
        Loading and treating the international trade merchandise database of Statistics Canada.
        Original source: https://open.canada.ca/data/en/dataset/b1126a07-fd85-4d56-8395-143aba1747a4
        :return:
        """

        # load concordance between HS classification and IOIC classification
        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            conc = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2014-2021/HS-IOIC.xlsx'))

        elif self.year in [2022]:
            conc = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2022/HS-IOIC.xlsx'))

        # load database
        self.imports = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Imports_data/Imports_' +
                                                                   str(self.year) + '_HS06_treated.xlsx'))

        self.imports = self.imports.ffill()
        self.imports.columns = ['Province', 'Country', 'HS6', 'Value']

        # apply concordance
        self.imports = self.imports.merge(conc, on='HS6', how='left')

        # only keep useful information
        self.imports = self.imports.loc[:, ['IOIC', 'Province', 'Country', 'Value']]

        # remove HS sectors that cannot be matched to IOIC
        self.imports = self.imports.dropna(subset=['IOIC'])

        # change IOIC codes to sector names
        code_to_name = {j[0]: j[1] for j in self.commodities}
        self.imports.loc[:, 'IOIC'] = [code_to_name[i] for i in self.imports.IOIC]

        # set MultiIndex with country and classification
        self.imports = self.imports.set_index(['Province', 'Country', 'IOIC'])

        # regroup purchases together (on province + country + IOIC sector)
        self.imports = self.imports.groupby(self.imports.index).sum()

        # set Multi-index
        self.imports.index = pd.MultiIndex.from_tuples(self.imports.index)

        # reset the index to apply country converter
        self.imports = self.imports.reset_index()
        # apply country converter
        self.imports.level_1 = coco.convert(self.imports.level_1, to='EXIO3')
        # restore index
        self.imports = self.imports.set_index(['level_0', 'level_1', 'level_2'])
        self.imports.index.names = None, None, None

        # groupby on country/sector (e.g., there were multiple 'WL' after applying coco)
        self.imports = self.imports.groupby(self.imports.index).sum()

        # restore multi-index
        self.imports.index = pd.MultiIndex.from_tuples(self.imports.index)

        # the absolute values do not matter, we only need the share of origins per province
        totals = self.imports.groupby(level=[0, 2])["Value"].transform("sum")
        self.imports["Share"] = self.imports["Value"] / totals
        self.imports["Share"] = self.imports["Share"].fillna(0)
        self.imports = self.imports.drop('Value', axis=1)

        # reindexing to ensure all sectors are here
        self.imports = self.imports.reindex(pd.MultiIndex.from_product([
            self.matching_dict.keys(),
            self.imports.index.levels[1],
            [i[1] for i in self.commodities]]))

        # there are some issues with the import of electricity for some provinces, so for those, hardcode that it's coming from the US
        self.imports.loc(axis=0)[:, 'US', 'Electricity'] = self.imports.loc(axis=0)[:, 'US', 'Electricity'].fillna(1)
        # also some issue with the Nickel import of Manitoba for the year 2014. Mismatch between SUTs and the merchandise trade database -> fix to US
        self.imports.loc(axis=0)['MB', 'US', 'Nickel ores and concentrates'] = (
            self.imports.loc(axis=0)['MB', 'US', 'Nickel ores and concentrates'].fillna(1))

        # get average origin of imports for each commodity for Canada according to exiobase
        canadian_imports_exio = pd.concat([self.Z_exio.loc[:, 'CA'],
                                           self.Y_exio.loc[:, 'CA']], axis=1).sum(1).drop('CA', axis=0, level=0)

        for commodity in canadian_imports_exio.index.levels[1]:
            canadian_imports_exio[[i for i in canadian_imports_exio.index if i[1] == commodity]] = (
                    canadian_imports_exio.loc(axis=0)[:, commodity] /
                    canadian_imports_exio.loc(axis=0)[:, commodity].sum()).fillna(0)

        # loading concordances between exiobase classification and IOIC
        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            ioic_exio = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2014-2021/IOIC_EXIOBASE.xlsx'),
                                      'commodities')
        elif self.year in [2022]:
            ioic_exio = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2022/IOIC_EXIOBASE.xlsx'),
                                      'commodities')
        ioic_exio = ioic_exio[2:].drop('IOIC Detail level - EXIOBASE', axis=1).set_index(
            'Unnamed: 1').infer_objects(copy=False).fillna(0)
        ioic_exio.index.name = None
        ioic_exio.index = [{j[0]: j[1] for j in self.commodities}[i] for i in ioic_exio.index]

        # for sectors that have no match, use exiobase Canadian imports
        for province in self.matching_dict.keys():
            for sector in self.imports.index.levels[2]:
                df = self.imports.loc(axis=0)[province, :, sector].copy()
                # if empty, it means we have no match with merchandise trade database -> use exiobase values
                if df.dropna().empty:
                    # 1 for 1 mapping with exiobase
                    if ioic_exio.loc[sector].sum() == 1:
                        exio_sector = ioic_exio.loc[sector][ioic_exio.loc[sector] == 1].index[0]
                        dff = canadian_imports_exio.loc(axis=0)[:, exio_sector]
                        dff /= dff.sum()
                        dff.index = pd.MultiIndex.from_tuples([(province, i[0], sector) for i in dff.index])
                        self.imports.loc[dff.index, 'Share'] = dff
                    # not mapped -> that's for the fictional sectors of the IOCC (e.g., FIC110000)
                    elif ioic_exio.loc[sector].sum() == 0:
                        pass
                    # 1 for many with exiobase
                    else:
                        exio_sectors = ioic_exio.loc[sector][ioic_exio.loc[sector] == 1].index.tolist()
                        dff = canadian_imports_exio.loc(axis=0)[:, exio_sectors].groupby('region').sum()
                        dff /= dff.sum()
                        dff.index = pd.MultiIndex.from_tuples([(province, i, sector) for i in dff.index])
                        self.imports.loc[dff.index, 'Share'] = dff

    def link_imports_to_openio(self):
        """
        Linking the international trade merchandise database of Statistics Canada to openIO-Canada.
        :return:
        """

        def scale_international_imports(who_uses_int_imports):
            imports = self.imports.copy()
            uses = who_uses_int_imports.copy()

            imports = imports.reset_index()
            uses = uses.reset_index()

            imports = imports.rename(columns={'level_0': 'Province', 'level_1': 'Origin', 'level_2': 'Commodity'})
            uses = uses.rename(columns={'level_0': 'Province', 'level_1': 'Commodity'})

            merged = pd.merge(imports, uses, on=['Province', 'Commodity'], how='inner')

            sector_cols = uses.columns.difference(['Province', 'Commodity'])
            for col in sector_cols:
                merged[col] = merged[col] * merged['Share']

            final = merged.set_index(['Province', 'Origin', 'Commodity'])[sector_cols]

            final = final.droplevel(0)

            final = final.groupby(final.index).sum()

            final.index = pd.MultiIndex.from_tuples(final.index)

            return final

        self.imports_scaled_U = scale_international_imports(self.who_uses_int_imports_U)
        self.imports_scaled_Y = scale_international_imports(self.who_uses_int_imports_Y)
        if self.endogenizing:
            self.imports_scaled_K = scale_international_imports(self.who_uses_int_imports_K)
            assert np.isclose(self.who_uses_int_imports_K.sum().sum(),
                              self.imports_scaled_K.sum().sum())

        assert np.isclose(self.who_uses_int_imports_U.sum().sum(),
                          self.imports_scaled_U.sum().sum())
        assert np.isclose(self.who_uses_int_imports_Y.sum().sum(),
                          self.imports_scaled_Y.sum().sum())

    def gimme_symmetric_iot(self):
        """
        Transforms Supply and Use tables to symmetric IO tables and transforms Y from product to industries if
        selected classification is "industry"
        :return:
        """

        # recalculate g because we introduced K
        if self.endogenizing:
            g = (self.U.sum() + self.W.sum() + self.who_uses_int_imports_U.sum() + self.K.sum() +
                 self.who_uses_int_imports_K.sum())
            g = pd.concat([g] * len(self.matching_dict), axis=1).T
            g.index = self.g.index
            g.columns = pd.MultiIndex.from_tuples(g.columns)
            for province in self.matching_dict.keys():
                g.loc[[i for i in g.index.levels[0] if i != province], province] = 0
            self.g = g.copy('deep')

        self.inv_q = pd.DataFrame(np.diag((1 / self.q.sum(axis=1)).replace(np.inf, 0)), self.q.index, self.q.index)
        self.inv_g = pd.DataFrame(np.diag((1 / self.g.sum()).replace(np.inf, 0)), self.g.columns, self.g.columns)

        self.A = self.U.dot(self.inv_g.dot(self.V.T)).dot(self.inv_q)
        self.R = self.W.dot(self.inv_g.dot(self.V.T)).dot(self.inv_q)

        self.imports_scaled_U = self.imports_scaled_U.reindex(
            self.U.columns, axis=1).dot(self.inv_g.dot(self.V.T)).dot(self.inv_q)

        if self.endogenizing:
            self.imports_scaled_K = self.imports_scaled_K.reindex(
                self.U.columns, axis=1).dot(self.inv_g.dot(self.V.T)).dot(self.inv_q)
            self.K = self.K.dot(self.inv_g.dot(self.V.T)).dot(self.inv_q)
            # assert that balance is respected after changing into symmetric tables
            assert (self.A.sum() + self.R.sum() + self.K.sum() +
                    self.imports_scaled_U.sum() + self.imports_scaled_K.sum())[
                       (self.A.sum() + self.R.sum() + self.K.sum() + self.imports_scaled_U.sum() +
                        self.imports_scaled_K.sum()) < 0.999].sum() == 0
        else:
            # assert that balance is respected after changing into symmetric tables
            assert ((self.A.sum() + self.R.sum() + self.imports_scaled_U.sum())[
                        (self.A.sum() + self.R.sum() + self.imports_scaled_U.sum()) < 0.999].sum() == 0)

    def link_international_trade_data_to_exiobase(self):
        """
        Linking the data from the international merchandise trade database, which was previously linked to openIO-Canada,
        to exiobase.
        :return:
        """

        # determine the Canadian imports according to Exiobase
        canadian_imports_exio = self.A_exio.loc[:, 'CA'].sum(1).drop('CA', axis=0, level=0)

        # loading concordances between exiobase classification and IOIC
        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            ioic_exio = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2014-2021/IOIC_EXIOBASE.xlsx'),
                                      'commodities')
        elif self.year in [2022]:
            ioic_exio = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2022/IOIC_EXIOBASE.xlsx'),
                                      'commodities')
        ioic_exio = ioic_exio[2:].drop('IOIC Detail level - EXIOBASE', axis=1).set_index(
            'Unnamed: 1').infer_objects(copy=False).fillna(0)
        ioic_exio.index.name = None
        ioic_exio.index = [{j[0]: j[1] for j in self.commodities}[i] for i in ioic_exio.index]

        def link_openIO_to_exio(merchandise_imports_scaled, who_uses_int_imports):

            link_openio_exio = pd.DataFrame()

            for merchandise in merchandise_imports_scaled.index.levels[1]:
                # check if there is trading happening for the uncovered commodity or not
                if who_uses_int_imports.groupby(level=1).sum().loc[merchandise].sum() != 0:
                    # 1 for 1 with exiobase -> easy
                    if ioic_exio.loc[merchandise].sum() == 1:
                        exio_sector = ioic_exio.loc[merchandise][ioic_exio.loc[merchandise] == 1].index[0]
                        dff = merchandise_imports_scaled.loc(axis=0)[:, merchandise]
                        dff.index = [(i[0], exio_sector) for i in dff.index]
                        link_openio_exio = pd.concat([link_openio_exio, dff])
                    # 1 for many with exiobase -> headscratcher
                    elif ioic_exio.loc[merchandise].sum() > 1:
                        exio_sector = ioic_exio.loc[merchandise][ioic_exio.loc[merchandise] == 1].index.tolist()
                        dff = merchandise_imports_scaled.loc(axis=0)[:, merchandise].copy()
                        dff = pd.concat([dff] * len(exio_sector))
                        dff = dff.sort_index()
                        dff.index = pd.MultiIndex.from_product([dff.index.levels[0], exio_sector])
                        for region in dff.index.levels[0]:
                            dfff = (dff.loc[region].T *
                                    (canadian_imports_exio.loc(axis=0)[region, exio_sector] /
                                     canadian_imports_exio.loc(axis=0)[region, exio_sector].sum()).loc[region]).T
                            # if our calculations shows imports (e.g., fertilizers from Bulgaria) for a product but there
                            # are not seen in exiobase, then we rely on self.x_exio to distribute between commodities
                            if not np.isclose(
                                    merchandise_imports_scaled.loc(axis=0)[:, merchandise].loc[region].sum().sum(),
                                    dfff.sum().sum()):
                                dfff = (dff.loc[region].T *
                                        (self.x_exio.loc(axis=0)[region, exio_sector].iloc[:, 0] /
                                         self.x_exio.loc(axis=0)[region, exio_sector].iloc[:, 0].sum()).loc[region]).T
                            dfff.index = pd.MultiIndex.from_product([[region], dfff.index])
                            link_openio_exio = pd.concat([link_openio_exio, dfff])
                            link_openio_exio.index = pd.MultiIndex.from_tuples(link_openio_exio.index)
                    else:
                        print(merchandise + ' is not linked to any Exiobase sector!')

            link_openio_exio.index = pd.MultiIndex.from_tuples(link_openio_exio.index)
            link_openio_exio = link_openio_exio.groupby(link_openio_exio.index).sum()
            link_openio_exio.index = pd.MultiIndex.from_tuples(link_openio_exio.index)
            link_openio_exio = link_openio_exio.reindex(self.A_exio.index).fillna(0)

            return link_openio_exio

        self.link_openio_exio_A = link_openIO_to_exio(self.imports_scaled_U, self.who_uses_int_imports_U)
        if self.endogenizing:
            self.link_openio_exio_K = link_openIO_to_exio(self.imports_scaled_K, self.who_uses_int_imports_K)
        self.link_openio_exio_Y = link_openIO_to_exio(self.imports_scaled_Y, self.who_uses_int_imports_Y)

        self.link_openio_exio_A = self.link_openio_exio_A.reindex(self.A.columns, axis=1)
        self.link_openio_exio_K = self.link_openio_exio_K.reindex(self.K.columns, axis=1)
        self.link_openio_exio_Y = self.link_openio_exio_Y.reindex(self.Y.columns, axis=1)

        # check financial balance is respected before converting to euros
        if self.endogenizing:
            assert (self.A.sum() + self.R.sum() + self.K.sum() +
                    self.link_openio_exio_A.sum() + self.link_openio_exio_K.sum())[
                       (self.A.sum() + self.R.sum() + self.K.sum() +
                        self.link_openio_exio_A.sum() + self.link_openio_exio_K.sum()) < 0.995].sum() == 0
        else:
            assert (self.A.sum() + self.R.sum() + self.link_openio_exio_A.sum())[
                       (self.A.sum() + self.R.sum() + self.link_openio_exio_A.sum()) < 0.995].sum() == 0

        # convert from CAD to EURO (https://www.bankofcanada.ca/rates/exchange/annual-average-exchange-rates/)
        if self.year == 2014:
            self.link_openio_exio_A /= 1.4671
            self.link_openio_exio_Y /= 1.4671
            if self.endogenizing:
                self.link_openio_exio_K /= 1.4671
        elif self.year == 2015:
            self.link_openio_exio_A /= 1.4182
            self.link_openio_exio_Y /= 1.4182
            if self.endogenizing:
                self.link_openio_exio_K /= 1.4182
        elif self.year == 2016:
            self.link_openio_exio_A /= 1.466
            self.link_openio_exio_Y /= 1.466
            if self.endogenizing:
                self.link_openio_exio_K /= 1.466
        elif self.year == 2017:
            self.link_openio_exio_A /= 1.465
            self.link_openio_exio_Y /= 1.465
            if self.endogenizing:
                self.link_openio_exio_K /= 1.465
        elif self.year == 2018:
            self.link_openio_exio_A /= 1.5302
            self.link_openio_exio_Y /= 1.5302
            if self.endogenizing:
                self.link_openio_exio_K /= 1.5302
        elif self.year == 2019:
            self.link_openio_exio_A /= 1.4856
            self.link_openio_exio_Y /= 1.4856
            if self.endogenizing:
                self.link_openio_exio_K /= 1.4856
        elif self.year == 2020:
            self.link_openio_exio_A /= 1.5298
            self.link_openio_exio_Y /= 1.5298
            if self.endogenizing:
                self.link_openio_exio_K /= 1.5298
        elif self.year == 2021:
            self.link_openio_exio_A /= 1.4828
            self.link_openio_exio_Y /= 1.4828
            if self.endogenizing:
                self.link_openio_exio_K /= 1.4828
        elif self.year == 2022:
            self.link_openio_exio_A /= 1.3696
            self.link_openio_exio_Y /= 1.3696
            if self.endogenizing:
                self.link_openio_exio_K /= 1.3696

    def remove_abroad_enclaves(self):
        """
        SUT accounts include Canadian abroad enclaves (i.e., embassies). Physical flow accounts however, do not provide
        emissions for these embassies. So we exclude them from the model.
        This slightly disturbs the economic balance as we do not re-equilibrate the model. But embassies accounts for
        an extremely low amount of trade. It is not worth the effort.
        :return:
        """

        self.V = self.V.drop('CE', axis=0, level=0)
        self.V = self.V.drop('CE', axis=1, level=0)
        self.U = self.U.drop('CE', axis=0, level=0)
        self.U = self.U.drop([i for i in self.U.columns if i[0] == 'CE'], axis=1)
        self.g = self.g.drop('CE', axis=0, level=0)
        self.g = self.g.drop('CE', axis=1, level=0)
        self.q = self.q.drop('CE', axis=0, level=0)
        self.q = self.q.drop('CE', axis=1, level=0)
        self.inv_g = self.inv_g.drop('CE', axis=0, level=0)
        self.inv_g = self.inv_g.drop('CE', axis=1, level=0)
        self.inv_q = self.inv_q.drop('CE', axis=0, level=0)
        self.inv_q = self.inv_q.drop('CE', axis=1, level=0)
        self.A = self.A.drop('CE', axis=0, level=0)
        self.A = self.A.drop('CE', axis=1, level=0)
        self.Y = self.Y.drop('CE', axis=0, level=0)
        self.Y = self.Y.drop('CE', axis=1, level=0)
        self.W = self.W.drop('CE', axis=0, level=0)
        self.W = self.W.drop('CE', axis=1, level=0)
        self.WY = self.WY.drop('CE', axis=0, level=0)
        self.WY = self.WY.drop('CE', axis=1, level=0)
        self.R = self.R.drop('CE', axis=0, level=0)
        self.R = self.R.drop('CE', axis=1, level=0)
        self.link_openio_exio_A = self.link_openio_exio_A.drop('CE', axis=1, level=0)
        self.link_openio_exio_Y = self.link_openio_exio_Y.drop('CE', axis=1, level=0)
        del self.matching_dict['CE']

        if self.endogenizing:
            self.K = self.K.drop('CE', axis=0, level=0)
            self.K = self.K.drop('CE', axis=1, level=0)
            self.link_openio_exio_K = self.link_openio_exio_K.drop('CE', axis=1, level=0)

    def concatenate_matrices(self):
        """
        Concatenate openIO-Canada matrices to Exiobase matrices and the link between them.
        :return:
        """

        # concat international trade with interprovincial trade
        self.A = pd.concat([self.A, self.link_openio_exio_A])
        # concat openIO-Canada with exiobase to get the full technology matrix
        df = pd.concat([pd.DataFrame(0, index=self.A.columns, columns=self.A_exio.columns), self.A_exio])
        self.A = pd.concat([self.A, df], axis=1)

        if self.endogenizing:
            self.K = pd.concat([self.K, self.link_openio_exio_K])
            df = pd.concat([pd.DataFrame(0, index=self.K.columns, columns=self.K_exio.columns), self.K_exio])
            self.K = pd.concat([self.K, df], axis=1)

        # concat interprovincial and international trade for final demands
        self.Y = pd.concat([self.Y, self.link_openio_exio_Y])

    def extract_environmental_data(self):
        """
        Extracts the data from the NPRI file
        :return: self.F but linked to NAICS codes
        """
        # Tab name changes with selected year, so identify it using "INRP-NPRI"
        emissions = self.NPRI[[i for i in self.NPRI.keys() if "INRP-NPRI" in i][0]]
        emissions.columns = list(zip(emissions.iloc[0].ffill().tolist(), emissions.iloc[2]))
        emissions = emissions.iloc[3:]
        # selecting the relevant columns from the file
        emissions = emissions.loc[:, [i for i in emissions.columns if
                                      (i[1] in
                                       ['NAICS 6 Code', 'CAS Number', 'Substance Name (English)', 'Units', 'Province']
                                       or 'Total' in i[1] and 'Air' in i[0]
                                       or 'Total' in i[1] and 'Water' in i[0]
                                       or 'Total' in i[1] and 'Land' in i[0])]].infer_objects(copy=False).fillna(0)
        # renaming the columns
        emissions.columns = ['Province', 'NAICS 6 Code', 'CAS Number', 'Substance Name', 'Units', 'Emissions to air',
                             'Emissions to water', 'Emissions to land']

        # somehow the NPRI manages to have entries without NAICS codes... Remove them
        no_naics_code_entries = emissions.loc[:, 'NAICS 6 Code'][emissions.loc[:, 'NAICS 6 Code'] == 0].index
        emissions.drop(no_naics_code_entries, inplace=True)
        # same for Provinces
        no_province_code_entries = emissions.loc[:, 'Province'][emissions.loc[:, 'Province'] == 0].index
        emissions.drop(no_province_code_entries, inplace=True)

        # NAICS codes as strings and not floats
        emissions.loc[:, 'NAICS 6 Code'] = emissions.loc[:, 'NAICS 6 Code'].astype('str')
        emissions.loc[:, 'NAICS 6 Code'] = [i.split('.0')[0] for i in emissions.loc[:, 'NAICS 6 Code']]

        # extracting metadata for substances
        temp_df = emissions.copy()
        temp_df.set_index('Substance Name', inplace=True)
        temp_df = temp_df.groupby(temp_df.index).head(n=1)
        # separating the metadata for emissions (CAS and units)
        self.emission_metadata = pd.DataFrame('', index=temp_df.index, columns=['CAS Number', 'Unit'])
        for emission in temp_df.index:
            self.emission_metadata.loc[emission, 'CAS Number'] = temp_df.loc[emission, 'CAS Number']
            self.emission_metadata.loc[emission, 'Unit'] = temp_df.loc[emission, 'Units']
        del temp_df

        self.F = pd.pivot_table(data=emissions, index=['Province', 'Substance Name'],
                                columns=['Province', 'NAICS 6 Code'], aggfunc="sum").fillna(0)
        try:
            self.F = self.F.drop(['CAS Number', 'Units'], axis=1)
        except KeyError:
            pass
        # renaming compartments
        self.F = self.F.rename(columns={'Emissions to air': 'Air', 'Emissions to soil': 'Soil',
                                        'Emissions to water': 'Water'}, level=0)
        # renaming the names of the columns indexes
        self.F.columns = self.F.columns.rename(['compartment', 'Province', 'NAICS'])
        # reorder multi index to have province as first level
        self.F = self.F.reorder_levels(['Province', 'compartment', 'NAICS'], axis=1)
        # match compartments with emissions and not to provinces
        self.F = self.F.T.unstack('compartment').T[self.F.T.unstack('compartment').T != 0].fillna(0)
        # identify emissions that are in tonnes
        emissions_to_rescale = [i for i in self.emission_metadata.index if
                                self.emission_metadata.loc[i, 'Unit'] == 'tonnes']
        # convert them to kg
        self.F.loc(axis=0)[:, emissions_to_rescale] *= 1000
        self.emission_metadata.loc[emissions_to_rescale, 'Unit'] = 'kg'
        # same thing for emissions in grams
        emissions_to_rescale = [i for i in self.emission_metadata.index if
                                self.emission_metadata.loc[i, 'Unit'] == 'grams']
        self.F.loc(axis=0)[:, emissions_to_rescale] /= 1000
        self.emission_metadata.loc[emissions_to_rescale, 'Unit'] = 'kg'

        # harmonizing emissions across provinces, set to zero if missing initially
        new_index = pd.MultiIndex.from_product(
            [self.matching_dict, self.emission_metadata.sort_index().index, ['Air', 'Water', 'Soil']])
        self.F = self.F.reindex(new_index).fillna(0)

        # harmonizing NAICS codes across provinces this time
        self.F = self.F.T.reindex(
            pd.MultiIndex.from_product([self.F.columns.levels[0], self.F.columns.levels[1]])).T.fillna(0)

    def match_npri_data_to_iots(self):

        total_emissions_origin = self.F.sum().sum()

        # load and format concordances file
        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            concordance = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2014-2021/NPRI_concordance.xlsx'),
                                        self.level_of_detail)
        elif self.year in [2022]:
            concordance = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2022/NPRI_concordance.xlsx'),
                                        self.level_of_detail)
        concordance.set_index('NAICS 6 Code', inplace=True)
        concordance.drop('NAICS 6 Sector Name (English)', axis=1, inplace=True)

        # splitting emissions between private and public sectors
        if self.level_of_detail == 'Summary level':
            self.split_private_public_sectors(NAICS_code=['611210', '611310', '611510'], IOIC_code='GS610')
        elif self.level_of_detail == 'Link-1961 level':
            self.split_private_public_sectors(NAICS_code=['611210', '611510'], IOIC_code='GS611B0')
            self.split_private_public_sectors(NAICS_code='611310', IOIC_code='GS61130')
        elif self.level_of_detail in ['Link-1997 level', 'Detail level']:
            self.split_private_public_sectors(NAICS_code='611210', IOIC_code='GS611200')
            self.split_private_public_sectors(NAICS_code='611310', IOIC_code='GS611300')
            self.split_private_public_sectors(NAICS_code='611510', IOIC_code='GS611A00')

        # switch NAICS codes in self.F for corresponding IOIC codes (from concordances file)
        IOIC_index = []
        for NAICS in self.F.columns:
            try:
                IOIC_index.append((NAICS[0], concordance.loc[int(NAICS[1]), 'IOIC']))
            except ValueError:
                IOIC_index.append(NAICS)
        self.F.columns = pd.MultiIndex.from_tuples(IOIC_index)

        # adding emissions from same sectors together (summary level is more aggregated than NAICS 6 Code)
        self.F = self.F.T.groupby(self.F.columns).sum().T
        # reordering columns
        self.F = self.F.T.reindex(
            pd.MultiIndex.from_product([self.matching_dict, [i[0] for i in self.industries]])).T.fillna(0)
        # changing codes for actual names of the sectors
        self.F.columns = pd.MultiIndex.from_product([self.matching_dict, [i[1] for i in self.industries]])

        # assert that nearly all emissions present in the NPRI were successfully transferred in self.F
        assert self.F.sum().sum() / total_emissions_origin > 0.98
        assert self.F.sum().sum() / total_emissions_origin < 1.02

    def match_ghg_accounts_to_iots(self):
        """
        Method matching GHG accounts to IOCC classification selected by the user
        :return: self.F and self.FY with GHG flows included
        """

        GHG = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Environmental_data/GHG_by_gas_type_2009-2023.xlsx'), 'Data_Données')
        GHG = GHG.loc[[i for i in GHG.index if GHG.loc[i, 'Reference period/Période de référence'] == self.year and
                       GHG.GEO_E[i] != 'Canada']]

        # kilotonnes to kgs
        GHG.loc[:, ['CO2', 'CH4 (CO2 eq)', 'N2O (CO2 eq)']] *= 1000000
        # data from StatCan is in kgCO2eq using IPCC AR5 values, so divide to apply our own up-to-date CFs
        GHG.loc[:, 'CH4'] = GHG.loc[:, 'CH4 (CO2 eq)'] / 28
        GHG.loc[:, 'N2O'] = GHG.loc[:, 'N2O (CO2 eq)'] / 265

        Household_GHG = GHG.loc[[i for i in GHG.index if 'PEH' in GHG.loc[i, 'IOIC_L61']]]

        Household_GHG = Household_GHG.drop(['Reference period/Période de référence', 'GEO_F', 'DescE_L61', 'DescF_L61',
                                            'Unit of measure/Unité de mesure', 'Total (CO2 eq)', 'CH4 (CO2 eq)',
                                            'N2O (CO2 eq)'], axis=1)
        # assume all direct emissions from home appliances come from "Other fuels"
        Household_GHG.IOIC_L61 = ['Other fuels' if i == 'PEH1' else 'Fuels and lubricants' for i in
                                  Household_GHG.IOIC_L61]
        Household_GHG.loc[Household_GHG.GEO_E == 'Québec', 'GEO_E'] = 'Quebec'
        Household_GHG.GEO_E = [{v: k for k, v in self.matching_dict.items()}[i] for i in Household_GHG.GEO_E]
        Household_GHG = pd.pivot_table(data=Household_GHG, values=['CO2', 'CH4', 'N2O'],
                                       columns=['GEO_E', 'IOIC_L61'])
        Household_GHG.columns = [(i[0], "Household final consumption expenditure", i[1]) for i in
                                 Household_GHG.columns]
        Household_GHG = Household_GHG.reindex(self.Y.columns, axis=1).fillna(0)

        # spatialization
        Household_GHG = pd.concat([Household_GHG] * len(self.matching_dict))
        Household_GHG.index = pd.MultiIndex.from_product([self.matching_dict,
                                                          ['Methane', 'Carbon dioxide', 'Dinitrogen monoxide'],
                                                          ['Air']]).drop_duplicates()
        for province in self.matching_dict:
            Household_GHG.loc[Household_GHG.index.get_level_values(0) != province, province] = 0

        # create FY and update it with GHG emissions from households
        self.FY = pd.DataFrame(0, index=pd.MultiIndex.from_product(
            [self.matching_dict, ['Carbon dioxide', 'Methane', 'Dinitrogen monoxide'], ['Air']]).drop_duplicates(),
                               columns=self.Y.columns, dtype=object)
        self.FY.update(Household_GHG)

        # Now the emissions from production
        GHG.loc[GHG.GEO_E == 'Québec', 'GEO_E'] = 'Quebec'
        GHG.set_index(pd.MultiIndex.from_tuples(tuple(
            list(zip([{v: k for k, v in self.matching_dict.items()}[i] for i in GHG.GEO_E], GHG.IOIC_L61.tolist())))),
            inplace=True)
        GHG.drop(['IOIC_L61', 'Reference period/Période de référence', 'GEO_E', 'GEO_F', 'DescE_L61', 'DescF_L61',
                  'Unit of measure/Unité de mesure', 'Total (CO2 eq)', 'CH4 (CO2 eq)', 'N2O (CO2 eq)'], axis=1,
                 inplace=True)
        GHG.drop([i for i in GHG.index if re.search(r'^FC', i[1])
                  or re.search(r'^PEH', i[1])
                  or re.search(r'^Total', i[1])], inplace=True)
        GHG.columns = ['Carbon dioxide', 'Methane', 'Dinitrogen monoxide']

        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            concordance = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2014-2021/GHG_concordance.xlsx'),
                                        self.level_of_detail)
        elif self.year in [2022]:
            concordance = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2022/GHG_concordance.xlsx'),
                                        self.level_of_detail)
        concordance.set_index('GHG codes', inplace=True)

        # dropping empty sectors (mostly Cannabis related)
        to_drop = concordance.loc[concordance.loc[:, 'IOIC'].isna()].index
        concordance.drop(to_drop, inplace=True)
        ghgs = pd.DataFrame()
        for code in concordance.index:
            # L97 and D levels are more precise than GHG accounts, we use market share to distribute GHGs
            sectors_to_split = [i[1] for i in self.industries if
                                i[0] in concordance.loc[code].dropna().values.tolist()]
            output_sectors_to_split = self.V.loc[:,
                                      [i for i in self.V.columns if i[1] in sectors_to_split]].sum()
            share_sectors_to_split = pd.DataFrame(0, output_sectors_to_split.index,
                                                  ['Carbon dioxide', 'Methane', 'Dinitrogen monoxide'], dtype=object)
            for province in self.matching_dict:
                df = ((output_sectors_to_split.loc[province] /
                       output_sectors_to_split.loc[province].sum()).fillna(0).values)
                # hardcoded 3 because 3 GHGs: CO2, CH4, N2O
                share_sectors_to_split.loc[province] = (pd.DataFrame(
                    [df] * 3,  index=['Carbon dioxide', 'Methane', 'Dinitrogen monoxide'],
                    columns=sectors_to_split).T * GHG.loc(axis=0)[:, code].loc[province].values).values
            ghgs = pd.concat([ghgs, share_sectors_to_split])

        # spatializing GHG emissions
        list_ghgs = ghgs.columns.tolist()
        ghgs = pd.concat([ghgs] * len(self.matching_dict), axis=1)
        ghgs.columns = pd.MultiIndex.from_product(
            [list(self.matching_dict.keys()), list_ghgs, ['Air']])
        for province in self.matching_dict:
            ghgs.loc[ghgs.index.get_level_values(0) != province, province] = 0

        # adding GHG accounts to pollutants
        self.F = pd.concat([self.F, ghgs.T])

        # reindexing
        self.F = self.F.reindex(self.U.columns, axis=1)

        self.emission_metadata.loc['Carbon dioxide', 'CAS Number'] = '124-38-9'
        self.emission_metadata.loc['Methane', 'CAS Number'] = '74-82-8'
        self.emission_metadata.loc['Dinitrogen monoxide', 'CAS Number'] = '10024-97-2'
        self.emission_metadata.loc[list_ghgs, 'Unit'] = 'kg'

    def match_water_accounts_to_iots(self):
        """
        Method matching water accounts to IOIC classification selected by the user
        :return: self.F and self.FY with GHG flows included
        """

        # load the water use data from STATCAN
        water_use = pd.read_csv(pkg_resources.resource_stream(__name__, '/Data/Environmental_data/water_use.csv'))
        water_conso = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Environmental_data/Water_consumption_values.xlsx'), None)

        # Only odd years from 2009 to 2017
        match_year_data = {2014: 2015, 2015: 2015, 2016: 2015, 2017: 2017, 2018: 2017, 2019: 2019, 2020: 2019,
                           2021: 2021, 2022: 2021}
        year_for_water = match_year_data[self.year]

        # select the year of the data
        water_use = water_use.loc[
            [i for i in water_use.index if water_use.REF_DATE[i] == int(year_for_water)], ['Sector', 'VALUE']].fillna(0)
        # convert into cubic meters
        water_use.VALUE *= 1000

        # huge mess to deal with the change in classification for year 2022
        if self.year == 2022:
            replacements = {'BS312200': 'BS312A00',
                            'BS332100': 'BS332B00',
                            'BS334100': 'BS334C00',
                            'BS335100': 'BS335B00',
                            'BS336310': 'BS336300',
                            'BS442000': 'BS449000',
                            'BS446000': 'BS456000',
                            'BS447000': 'BS457000',
                            'BS448000': 'BS458000',
                            'BS451000': 'BS459A00',
                            'BS452000': 'BS455000',
                            'BS453BL0': 'BS459BL0',
                            'BS453BU0': 'BS459BU0',
                            'BS511110': 'BS513110',
                            'BS5111A0': 'BS5131A0',
                            'BS511200': 'BS513200',
                            'BS515100': 'BS516100',
                            'BS515200': 'BS516211',
                            'BS522130': 'BS522B00',
                            'BS522300': 'BS52C000',
                            'BS52A000': 'BS52D000',
                            'BS524200': 'BS526A00'}

            for replace in replacements:
                water_conso['Manufacturing Industry'].loc[:, 'Code sector industry Open IO'] = (
                    water_conso['Manufacturing Industry'].loc[:, 'Code sector industry Open IO'].replace(replace,
                                                                                                         replacements[
                                                                                                             replace]))
                water_conso['Commercial and Institutionnal'].loc[:, 'Code sector industry Open IO'] = (
                    water_conso['Commercial and Institutionnal'].loc[:, 'Code sector industry Open IO'].replace(replace,
                                                                                                                replacements[
                                                                                                                    replace]))

            water_conso['Manufacturing Industry'] = water_conso['Manufacturing Industry'][
                ~water_conso['Manufacturing Industry'].loc[:, 'Code sector industry Open IO'].isin(
                    ['BS332A00', 'BS332500', 'BS332600', 'BS334A00', 'BS335200', 'BS336320', 'BS336330', 'BS336340',
                     'BS336350', 'BS336360',
                     'BS336370', 'BS336380', 'BS336390'])]

            water_conso['Commercial and Institutionnal'] = water_conso['Commercial and Institutionnal'][
                ~water_conso['Commercial and Institutionnal'].loc[:, 'Code sector industry Open IO'].isin(
                    ['BS443000', 'BS453A00', 'BS454000', 'GS611A00'])]

            df = water_conso['Commercial and Institutionnal'][
                water_conso['Commercial and Institutionnal'].loc[:,
                'Code sector industry Open IO'] == 'GS622000'].copy()
            df.loc[:, 'Code sector industry Open IO'] = 'GS62A000'
            water_conso['Commercial and Institutionnal'] = pd.concat([water_conso['Commercial and Institutionnal'], df])

            df = water_conso['Commercial and Institutionnal'][
                water_conso['Commercial and Institutionnal'].loc[:,
                'Code sector industry Open IO'] == 'BS516100'].copy()
            df.loc[:, 'Code sector industry Open IO'] = 'BS51621A'
            water_conso['Commercial and Institutionnal'] = pd.concat([water_conso['Commercial and Institutionnal'], df])

            water_conso['Commercial and Institutionnal'] = water_conso[
                'Commercial and Institutionnal'].reset_index().drop('index', axis=1)

        # -------------------------------------- Final demand ------------------------------------------------
        fd_water = water_use.loc[water_use.Sector == 'Households']
        fd_water = fd_water.drop(['Sector'], axis=1)
        # distribute national water use to all provinces
        provincial_household_water_use = (self.Y.loc[[i for i in self.Y.index if
                                                      i[1] == 'Water delivered by water works and irrigation systems'],
                                                     [i for i in self.Y.columns if
                                                      i[2] == 'Water supply and sanitation services']].sum() /
                                          self.Y.loc[[i for i in self.Y.index if
                                                      i[1] == 'Water delivered by water works and irrigation systems'],
                                                     [i for i in self.Y.columns if
                                                      i[2] == 'Water supply and sanitation services']].sum().sum())
        fd_water = provincial_household_water_use * fd_water.iloc[0, 0]

        # spatializing
        fd_water = pd.concat([fd_water] * len(fd_water.index), axis=1)
        fd_water.columns = pd.MultiIndex.from_product([self.matching_dict.keys(), ['Water'], ['Water']])

        for province in fd_water.columns.levels[0]:
            fd_water.loc[province, fd_water.columns.get_level_values(0) != province] = 0

        fd_water = fd_water.reindex(self.Y.columns).T.fillna(0)
        # transform water use into water consumption
        fd_water *= water_conso['Households'].loc[:, 'Water consumption (%) 2013'].iloc[0] / 100
        self.FY = pd.concat([self.FY, fd_water], axis=0)

        # -----------------------------------Intermediary demand ------------------------------------------
        water_use = water_use.loc[[i for i in water_use.index if '[' in water_use.Sector[i]]]
        water_use.Sector = [i.split('[')[1].split(']')[0] for i in water_use.Sector]
        water_use.drop([i for i in water_use.index if re.search(r'^FC', water_use.Sector.loc[i])], inplace=True)
        water_use.set_index('Sector', inplace=True)

        # load concordances matching water use data classification to the different classifications used in OpenIO
        concordance = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2014-2021/water_concordance.xlsx'),
                                    self.level_of_detail)
        concordance.set_index('Sector', inplace=True)
        # dropping potential empty sectors (mostly Cannabis related)
        to_drop = concordance.loc[concordance.loc[:, 'IOIC'].isna()].index
        concordance.drop(to_drop, inplace=True)

        water_flows = pd.DataFrame()
        for code in concordance.index:
            # Detail level is more precise than water accounts, we use market share to distribute water flows
            sectors_to_split = [i[1] for i in self.industries if
                                i[0] in concordance.loc[code].dropna().values.tolist()]
            output_sectors_to_split = self.V.loc[:,
                                      [i for i in self.V.columns if i[1] in sectors_to_split]].sum()

            share_sectors_to_split = output_sectors_to_split / output_sectors_to_split.sum() * water_use.loc[
                code, 'VALUE']
            water_flows = pd.concat([water_flows, share_sectors_to_split])

        water_flows = water_flows.groupby(water_flows.index).sum().fillna(0)
        # multi index for the win
        water_flows.index = pd.MultiIndex.from_tuples(water_flows.index)
        water_flows.columns = ['Water']
        # apply water consumption % for specific industries
        relative_water_consumption_data = pd.concat([water_conso['Manufacturing Industry'],
                                                     water_conso['Commercial and Institutionnal'],
                                                     water_conso['Oil and Gas'],
                                                     water_conso['Mining']])
        relative_water_consumption_data.index = [i for i in range(0, len(relative_water_consumption_data.index))]
        for i in relative_water_consumption_data.index:
            relative_water_consumption_data.loc[i, 'Industry name'] = {_[0]: _[1] for _ in self.industries}[
                relative_water_consumption_data.loc[i, 'Code sector industry Open IO']]

        water_flows = water_flows.reset_index().merge(relative_water_consumption_data, left_on='level_1',
                                                      right_on='Industry name', how='left')
        water_flows.loc[:, 'Average water consumption'] = water_flows.loc[:, [i for i in water_flows.columns if
                                                                              'Water consumption' in i]].mean(axis=1)
        water_flows = water_flows.set_index(['level_0', 'level_1']).loc[:, ['Water', 'Average water consumption']]
        water_flows.index.names = None, None
        water_flows.loc[:, 'Water'] *= water_flows.loc[:, 'Average water consumption'] / 100
        water_flows.drop('Average water consumption', axis=1, inplace=True)

        # For electricity we have absolute water consumption data (not applying a %)
        water_conso['Electricity'] = water_conso['Electricity'][water_conso['Electricity'].Year == self.year]
        water_conso['Electricity'].loc[:, 'Code product Open IO'] = [{i[0]: i[1] for i in self.industries}[i] for i in
                                                                     water_conso['Electricity'].loc[:,
                                                                     'Code product Open IO']]
        water_flows.loc(axis=0)[:, 'Electric power generation, transmission and distribution'] = (
        water_conso['Electricity'].set_index(
            ['Canada Province', 'Code product Open IO']).drop(['Year'], axis=1).loc[:, 'Water Consumption (m3)'])
        # crop and livestock water consumption data operate at the product level. For now set these values to zero.
        water_flows = water_flows.fillna(0)
        # spatializing
        water_flows = pd.concat([water_flows.T] * len(water_flows.index.levels[0]))
        water_flows.index = pd.MultiIndex.from_product([self.matching_dict.keys(), ['Water'], ['Water']])
        water_flows = water_flows.T.reindex(self.F.columns).T
        for province in water_flows.index.levels[0]:
            water_flows.loc[province, water_flows.columns.get_level_values(0) != province] = 0

        self.F = pd.concat([self.F, water_flows]).infer_objects(copy=False).fillna(0)

        self.emission_metadata.loc['Water', 'Unit'] = 'm3'

    def match_energy_accounts_to_iots(self):
        """
        Method matching energy accounts to IOIC classification selected by the user
        :return: self.F and self.FY with GHG flows included
        """
        NRG = pd.read_csv(pkg_resources.resource_stream(__name__, '/Data/Environmental_data/Energy_use.csv'))
        # select year of study
        NRG = NRG.loc[[i for i in NRG.index if NRG.REF_DATE[i] == self.year]]
        # keep households energy consumption in a specific dataframe
        NRG_FD = NRG.loc[[i for i in NRG.index if 'Households' in NRG.Sector[i]]]
        # keep industry energy consumption
        NRG = NRG.loc[[i for i in NRG.index if '[' in NRG.Sector[i]]]
        # extract sector codes
        NRG.Sector = [i.split('[')[1].split(']')[0] for i in NRG.Sector]
        # pivot into a dataframe
        NRG = NRG.pivot_table(values='VALUE', index=['Sector'], dropna=False).fillna(0)
        # remove fictive sectors
        NRG.drop([i for i in NRG.index if re.search(r'^FC', i)], inplace=True)

        # ------------ Industries ----------------

        # load concordance file
        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            concordance = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2014-2021/Energy_concordance.xlsx'),
                                        self.level_of_detail)
        elif self.year in [2022]:
            concordance = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2022/Energy_concordance.xlsx'),
                                        self.level_of_detail)
        concordance.set_index('NRG codes', inplace=True)
        # dropping empty sectors (mostly Cannabis related)
        to_drop = concordance.loc[concordance.loc[:, 'IOIC'].isna()].index
        concordance.drop(to_drop, inplace=True)

        # distributing energy use to more precise classifications based on market shares
        nrg = pd.DataFrame()
        for code in concordance.index:
            sectors_to_split = [i[1] for i in self.industries if
                                i[0] in concordance.loc[code].dropna().values.tolist()]
            output_sectors_to_split = self.V.loc[:,
                                      [i for i in self.V.columns if i[1] in sectors_to_split]].sum()
            output_sectors_to_split = output_sectors_to_split.groupby(axis=0, level=1).sum()
            share_sectors_to_split = output_sectors_to_split / output_sectors_to_split.sum()
            share_sectors_to_split *= NRG.loc[code, 'VALUE']
            nrg = pd.concat([nrg, share_sectors_to_split])

        # distributing national energy use to provinces based on market shares
        nrg_provincial = pd.DataFrame(0, index=pd.MultiIndex.from_product([self.matching_dict, nrg.index]),
                                      columns=['Energy'], dtype=object)

        for sector in nrg.index:
            share_province = self.g.loc(axis=1)[:, sector].sum(0) / self.g.loc(axis=1)[:, sector].sum(1).sum() * \
                             nrg.loc[sector].iloc[0]
            nrg_provincial.loc[share_province.index] = pd.DataFrame(share_province, columns=['Energy'])

        # adding to self.F
        self.F = pd.concat([self.F, nrg_provincial.reindex(self.F.columns).T])
        # cannabis stores are NaN values, we change that to zero values
        self.F = self.F.infer_objects(copy=False).fillna(0)

        # ------------- Final demand -------------

        # pivot into a dataframe
        NRG_FD = NRG_FD.pivot_table(values='VALUE', index=['Sector'], dropna=False).fillna(0)
        # rename index to IOIC FD classification
        NRG_FD.index = ['Other fuels', 'Fuels and lubricants']


        # distributing national final demand energy use to provinces based on market shares
        nrg_fd_provincial = pd.DataFrame(0, index=pd.MultiIndex.from_product(
            [self.matching_dict, ['Household final consumption expenditure'], NRG_FD.index]), columns=['Energy'],
                                         dtype=object)

        for fd_sector in NRG_FD.index:
            share_province = self.Y.loc(axis=1)[:, :, fd_sector].sum() / self.Y.loc(axis=1)[:, :,
                                                                         fd_sector].sum().sum() * NRG_FD.loc[
                                 fd_sector, 'VALUE']
            nrg_fd_provincial.loc[share_province.index] = pd.DataFrame(share_province, columns=['Energy'])

        # adding to self.FY
        self.FY = pd.concat([self.FY, nrg_fd_provincial.reindex(self.Y.columns).infer_objects(copy=False).fillna(0).T])
        # cannabis stores are NaN values, we change that to zero values
        self.FY = self.FY.infer_objects(copy=False).fillna(0)

        self.emission_metadata.loc['Energy', 'Unit'] = 'TJ'

    def match_mineral_extraction_to_iots(self):
        """
        Method matching mineral extraction data from USGS to IOIC classification selected by the user
        :return: self.F with mineral flows included
        """
        xl = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Environmental_data/Minerals_extracted_in_Canada.xlsx'), index_col=0).drop(
            'Comments (data in metric tons)', axis=1)

        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            with open(pkg_resources.resource_filename(__name__, '/Data/Concordances/2014-2021/concordance_metals.json'), 'r') as f:
                dict_data = json.load(f)
        elif self.year in [2022]:
            with open(pkg_resources.resource_filename(__name__, '/Data/Concordances/2022/concordance_metals.json'), 'r') as f:
                dict_data = json.load(f)

        distrib_minerals = pd.DataFrame()
        for mineral_sector in list(set(dict_data.values())):
            df = self.q.sum(1).loc(axis=0)[:, mineral_sector].copy()
            df /= df.sum()
            distrib_minerals = pd.concat([distrib_minerals, df])

        distrib_minerals.index = pd.MultiIndex.from_tuples(distrib_minerals.index)

        self.minerals = pd.DataFrame(0, index=dict_data, columns=self.q.index, dtype=float)

        for mineral in dict_data:
            # check if data for year is available
            if not math.isnan(xl.loc[mineral, self.year]):
                df = xl.loc[mineral, self.year] * distrib_minerals.loc(axis=0)[:, dict_data[mineral]]
            # if it's not, take latest available year value
            else:
                df = xl.loc[mineral].dropna().iloc[-1] * distrib_minerals.loc(axis=0)[:, dict_data[mineral]]
            df.columns = [mineral]
            df = df.T
            self.minerals.loc[mineral, df.columns] = df.loc[mineral]

        # convert from thousand carats to metric tons
        self.minerals.loc['Diamond'] *= 0.0002
        # from metric tons to kgs
        self.minerals *= 1000

        self.emission_metadata = pd.concat([self.emission_metadata, pd.DataFrame('kg', index=self.minerals.index,
                                                                                 columns=['Unit'])])

    def characterization_matrix(self):
        """
        Produces a characterization matrix from IMPACT World+ file
        :return: self.C, self.methods_metadata
        """

        IW = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Characterization_factors/impact_world_plus_2.1_dev.xlsx'))

        IW = pd.pivot_table(IW, values='CF value', index=('Impact category', 'CF unit'),
                            columns=['Elem flow name', 'Compartment', 'Sub-compartment']).fillna(0)

        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            concordance = pd.read_excel(pkg_resources.resource_stream(
                __name__, '/Data/Concordances/2014-2021/openIO_IW_concordance.xlsx'), 'Total')
        elif self.year in [2022]:
            concordance = pd.read_excel(pkg_resources.resource_stream(
                __name__, '/Data/Concordances/2022/openIO_IW_concordance.xlsx'), 'Total')
        concordance.set_index('OpenIO flows', inplace=True)

        # applying concordance
        hfcs = ['CF4', 'C2F6', 'SF6', 'NF3', 'c-C4F8', 'C3F8', 'HFC-125', 'HFC-134a', 'HFC-143', 'HFC-143a', 'HFC-152a',
                'HFC-227ea', 'HFC-23', 'HFC-32', 'HFC-41', 'HFC-134', 'HFC-245fa', 'HFC-43-10mee', 'HFC-365mfc', 'HFC-236fa']

        hfcs_idx = pd.MultiIndex.from_product(
            [hfcs, [i for i in self.matching_dict.keys()], ['Air']]).swaplevel(0, 1)

        self.C = pd.DataFrame(0, IW.index, self.F.index.tolist() + self.minerals.index.tolist() +
                              hfcs_idx.tolist(), dtype=float)
        for flow in self.C.columns:
            if type(flow) == tuple:
                try:
                    if concordance.loc[flow[1], 'IMPACT World+ flows'] is not None:
                        self.C.loc[:, [flow]] = IW.loc[:, [(concordance.loc[flow[1], 'IMPACT World+ flows'], flow[2],
                                                            '(unspecified)')]].values
                except KeyError:
                    pass
            # if type == str -> we are looking at the minerals extension -> hardcode Raw/in ground comp/subcomp
            elif type(flow) == str:
                try:
                    if concordance.loc[flow, 'IMPACT World+ flows'] is not None:
                        self.C.loc[:, [flow]] = IW.loc[:, [(concordance.loc[flow, 'IMPACT World+ flows'],
                                                            'Raw', 'in ground')]].values
                except KeyError:
                    pass

        # Non IW characterization factors
        self.C.loc[('Energy use', 'TJ'), [i for i in self.C.columns if i == 'Energy']] = 1
        self.C = self.C.infer_objects(copy=False).fillna(0)

        self.emission_metadata.loc['CF4', 'CAS Number'] = '75-73-0'
        self.emission_metadata.loc['C2F6', 'CAS Number'] = '76-16-4'
        self.emission_metadata.loc['C3F8', 'CAS Number'] = '76-19-7'
        self.emission_metadata.loc['HFC-125', 'CAS Number'] = '354-33-6'
        self.emission_metadata.loc['HFC-134', 'CAS Number'] = '359-35-3'
        self.emission_metadata.loc['HFC-134a', 'CAS Number'] = '811-97-2'
        self.emission_metadata.loc['HFC-143', 'CAS Number'] = '430-66-0'
        self.emission_metadata.loc['HFC-143a', 'CAS Number'] = '420-46-2'
        self.emission_metadata.loc['HFC-152a', 'CAS Number'] = '75-37-6'
        self.emission_metadata.loc['HFC-227ea', 'CAS Number'] = '431-89-0'
        self.emission_metadata.loc['HFC-23', 'CAS Number'] = '75-46-7'
        self.emission_metadata.loc['HFC-236fa', 'CAS Number'] = '690-39-1'
        self.emission_metadata.loc['HFC-245fa', 'CAS Number'] = '460-73-1'
        self.emission_metadata.loc['HFC-32', 'CAS Number'] = '75-10-5'
        self.emission_metadata.loc['HFC-365mfc', 'CAS Number'] = '406-58-6'
        self.emission_metadata.loc['HFC-41', 'CAS Number'] = '593-53-3'
        self.emission_metadata.loc['HFC-43-10mee', 'CAS Number'] = '138495-42-8'
        self.emission_metadata.loc['NF3', 'CAS Number'] = '7783-54-2'
        self.emission_metadata.loc['SF6', 'CAS Number'] = '2551-62-4'
        self.emission_metadata.loc['c-C4F8', 'CAS Number'] = '115-25-3'

        self.emission_metadata.loc[hfcs, 'Unit'] = 'kg'

        # some methods of IMPACT World+ do not make sense in our context, remove them
        self.C.drop(['Fossil and nuclear energy use',
                     'Ionizing radiations',
                     'Ionizing radiations, ecosystem quality',
                     'Ionizing radiations, human health',
                     'Land occupation, biodiversity',
                     'Land transformation, biodiversity',
                     'Thermally polluted water'], axis=0, level=0, inplace=True)

        # importing characterization matrix IMPACT World+/exiobase
        self.C_exio = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Characterization_factors/impact_world_plus_2.1_expert_version_exiobase.xlsx'),
            index_col='Unnamed: 0')
        self.C_exio.index = pd.MultiIndex.from_tuples(list(zip(
            [i.split(' (')[0] for i in self.C_exio.index],
            [i.split(' (')[1].split(')')[0] for i in self.C_exio.index])))

        self.C_exio.drop(['Fossil and nuclear energy use',
                          'Ionizing radiations',
                          'Ionizing radiations, ecosystem quality',
                          'Ionizing radiations, human health',
                          'Land occupation, biodiversity',
                          'Land transformation, biodiversity',
                          'Thermally polluted water'], axis=0, level=0, inplace=True)

        # dealing with water characterization factors
        def regionalize_water_extension_exiobase(extension_matrix):

            for country in extension_matrix.columns.levels[0]:
                df_agri = pd.DataFrame(
                    extension_matrix.loc[[i for i in extension_matrix.index if
                                          'Water Consumption Blue - Agriculture' in i], country].sum())
                df_agri.columns = pd.MultiIndex.from_product(
                    [[country], ['Total Water Consumption Blue - Agriculture']])
                df_agri = pd.concat([df_agri], keys=[country])

                df_non_agri = pd.DataFrame(extension_matrix.loc[[i for i in extension_matrix.index if (
                            'Water Consumption Blue' in i and 'Agriculture' not in i)], country].sum())
                df_non_agri.columns = pd.MultiIndex.from_product(
                    [[country], ['Total Water Consumption Blue - Non-agriculture']])
                df_non_agri = pd.concat([df_non_agri], keys=[country])

                extension_matrix = pd.concat([extension_matrix, df_agri.T, df_non_agri.T])

            return extension_matrix.fillna(0)

        self.F_exio = regionalize_water_extension_exiobase(self.F_exio)
        self.S_exio = regionalize_water_extension_exiobase(self.S_exio)
        self.C_exio = self.C_exio.T.reindex(self.S_exio.index).T.fillna(0)

        openio_water = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Characterization_factors/Regionalized_water_flows_CF_openio.xlsx'))
        openio_water.columns = list(zip(openio_water.columns, openio_water.iloc[0]))
        openio_water.index = openio_water.iloc[:, 0]
        openio_water = openio_water.iloc[1:, 1:]
        openio_water.index.name = None
        openio_water.columns = pd.MultiIndex.from_tuples(openio_water.columns)
        openio_water.index = [eval(i) for i in openio_water.index]
        self.C.loc[openio_water.columns, openio_water.index] = openio_water.T.values

        exio_water = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Characterization_factors/Regionalized_water_flows_CF_exiobase.xlsx'))
        exio_water.columns = list(zip(exio_water.columns, exio_water.iloc[0]))
        exio_water.index = exio_water.iloc[:, 0]
        exio_water = exio_water.iloc[2:, 1:]
        exio_water.index.name = None
        exio_water.columns = pd.MultiIndex.from_tuples(exio_water.columns)
        exio_water.index = [eval(i) for i in exio_water.index]
        self.C_exio.loc[exio_water.columns, exio_water.index] = exio_water.T.values
        # need to set to zero the other water consumption flows to not double count
        self.C_exio.loc[('Water scarcity', 'm3 world-eq'),
                        [i for i in self.C_exio.columns if 'Water Consumption Blue' in i]] = 0

        # adding energy use to exiobase flows to match with energy use from STATCAN physical accounts
        # energy use in exiobase is identified through "Energy Carrier Use: Total"
        # note that STATCAN only covers energy use, thus energy supply, loss, etc. flows from exiobase are excluded
        adding_energy_use = pd.DataFrame(0, index=pd.MultiIndex.from_product([['Energy'], ['TJ']]),
                                         columns=self.S_exio.index)
        adding_energy_use.loc[:, [i for i in self.S_exio.index if 'Energy Carrier Use: Total' in i]] = 1
        self.C_exio = pd.concat([self.C_exio, adding_energy_use])
        # forcing the match with self.C (annoying parentheses for climate change long and short term)
        self.C_exio.index = self.C.index
        self.C_exio = self.C_exio.fillna(0)

        self.methods_metadata = pd.DataFrame(self.C.index.tolist(), columns=['Impact category', 'unit'])
        self.methods_metadata = self.methods_metadata.set_index('Impact category')

        self.balance_flows(concordance)

        self.C = self.C.join(self.C_exio)

    def better_distribution_for_agriculture_ghgs(self):
        """
        GHG physical flow accounts from StatCan only provide the GHG emissions for Crop and animal production aggregated.
        By default, an economic allocation is applied to distribute these emissions to the corresponding sectors.
        However, this is too constraining, as direct emissions vary vastly between animal and crop productions.
        We rely on Exiobase to get a better distribution for GHGs in these sectors.
        :return:
        """

        # separate crops from animal breeding in Exiobase sectors
        crops_exio = ['Paddy rice', 'Wheat', 'Cereal grains nec', 'Vegetables, fruit, nuts', 'Oil seeds',
                      'Sugar cane, sugar beet', 'Plant-based fibers', 'Crops nec']
        animals_exio = ['Cattle', 'Pigs', 'Poultry', 'Meat animals nec', 'Raw milk', 'Wool, silk-worm cocoons']

        # identify the three GHGs that are covered by openIO
        CO2 = [i for i in self.F_exio.index if 'CO2' in i]
        CH4 = [i for i in self.F_exio.index if 'CH4' in i]
        N2O = [i for i in self.F_exio.index if 'N2O' in i]
        # isolate GHG emissions from crop production in Exiobase
        crops_emissions = pd.concat(
            [self.F_exio.loc(axis=1)[:, crops_exio].T.groupby(level=0).sum().T.loc[CO2].sum(),
             self.F_exio.loc(axis=1)[:, crops_exio].T.groupby(level=0).sum().T.loc[CH4].sum(),
             self.F_exio.loc(axis=1)[:, crops_exio].T.groupby(level=0).sum().T.loc[N2O].sum()],
            axis=1)
        crops_emissions.columns = ['Carbon dioxide', 'Methane', 'Dinitrogen monoxide']
        crops_emissions = crops_emissions.loc['CA']
        # isolate GHG emissions from meat production in Exiobase
        meat_emissions = pd.concat(
            [self.F_exio.loc(axis=1)[:, animals_exio].T.groupby(level=0).sum().T.loc[CO2].sum(),
             self.F_exio.loc(axis=1)[:, animals_exio].T.groupby(level=0).sum().T.loc[CH4].sum(),
             self.F_exio.loc(axis=1)[:, animals_exio].T.groupby(level=0).sum().T.loc[N2O].sum()], axis=1)
        meat_emissions.columns = ['Carbon dioxide', 'Methane', 'Dinitrogen monoxide']
        meat_emissions = meat_emissions.loc['CA']
        # get the totals per GHG
        tot_co2 = crops_emissions.loc['Carbon dioxide'] + meat_emissions.loc['Carbon dioxide']
        tot_ch4 = crops_emissions.loc['Methane'] + meat_emissions.loc['Methane']
        tot_n2o = crops_emissions.loc['Dinitrogen monoxide'] + meat_emissions.loc['Dinitrogen monoxide']
        # calculate the distribution, according to Exiobase
        crops_emissions.loc['Carbon dioxide'] /= tot_co2
        crops_emissions.loc['Methane'] /= tot_ch4
        crops_emissions.loc['Dinitrogen monoxide'] /= tot_n2o
        # calculate the distribution, according to Exiobase
        meat_emissions.loc['Carbon dioxide'] /= tot_co2
        meat_emissions.loc['Methane'] /= tot_ch4
        meat_emissions.loc['Dinitrogen monoxide'] /= tot_n2o
        # store it in a single dataframe
        ghgs_exio_distribution = pd.concat([crops_emissions, meat_emissions], axis=1)
        ghgs_exio_distribution.columns = ['Crops', 'Meat']

        # now that we have the distribution of GHG per big sector, we apply this distribution to openIO data
        for ghg in ['Carbon dioxide', 'Methane', 'Dinitrogen monoxide']:

            tot = self.F.loc[[i for i in self.F.index if i[1] == ghg], [i for i in self.F.columns if i[1] in [
                'Crop production (except cannabis, greenhouse, nursery and floriculture production)',
                'Greenhouse, nursery and floriculture production (except cannabis)',
                'Animal production (except aquaculture)', 'Aquaculture']]].T.groupby(level=0).sum().T

            crops = self.F.loc[[i for i in self.F.index if i[1] == ghg], [i for i in self.F.columns if i[1] in [
                'Crop production (except cannabis, greenhouse, nursery and floriculture production)',
                'Greenhouse, nursery and floriculture production (except cannabis)']]]

            animals = self.F.loc[[i for i in self.F.index if i[1] == ghg], [i for i in self.F.columns if i[1] in [
                'Animal production (except aquaculture)', 'Aquaculture']]]

            for province in tot.columns:
                tot_prod_crop_and_meat_province = tot.loc[[(province, ghg, 'Air')], province].iloc[0]

                exio_crop_distrib = ghgs_exio_distribution.loc[ghg, 'Crops']
                crops.loc[[(province, ghg, 'Air')]] = (crops.loc[[(province, ghg, 'Air')]] /
                                                       crops.loc[[(province, ghg, 'Air')]].sum().sum() *
                                                       exio_crop_distrib * tot_prod_crop_and_meat_province)
                self.F.loc[[i for i in self.F.index if i[1] == ghg and i[0] == province], [
                    i for i in self.F.columns if i[1] in [
                        'Crop production (except cannabis, greenhouse, nursery and floriculture production)',
                        'Greenhouse, nursery and floriculture production (except cannabis)']]] = crops.loc[[(province, ghg, 'Air')]]

                exio_animal_distrib = ghgs_exio_distribution.loc[ghg, 'Meat']
                animals.loc[[(province, ghg, 'Air')]] = (animals.loc[[(province, ghg, 'Air')]] /
                                                         animals.loc[[(province, ghg, 'Air')]].sum().sum() *
                                                         exio_animal_distrib * tot_prod_crop_and_meat_province)

                self.F.loc[[i for i in self.F.index if i[1] == ghg and i[0] == province], [
                    i for i in self.F.columns if i[1] in [
                        'Animal production (except aquaculture)','Aquaculture']]] = animals.loc[[(province, ghg, 'Air')]]

    def differentiate_country_names_openio_exio(self):
        """
        Some province names are identical to country names in exiobase (e.g., 'SK' and 'NL'). So we changed province
        names to, e.g., 'CA-SK'.
        :return:
        """

        self.A.index = (pd.MultiIndex.from_product([['CA-' + i for i in self.matching_dict.keys()],
                                                    [i[1] for i in self.commodities]]).tolist() +
                        self.A_exio.index.tolist())
        self.A.index = pd.MultiIndex.from_tuples(self.A.index)
        self.A.columns = self.A.index
        self.K.index = self.A.index
        if self.endogenizing:
            self.K.columns = self.A.columns
        self.Y.index = self.A.index
        self.Y.columns = [('CA-' + i[0], i[1], i[2]) for i in self.Y.columns]
        self.Y.columns = pd.MultiIndex.from_tuples(self.Y.columns)
        self.R.columns = pd.MultiIndex.from_product([['CA-' + i for i in self.matching_dict.keys()],
                                                     [i[1] for i in self.commodities]]).tolist()
        self.R.columns = pd.MultiIndex.from_tuples(self.R.columns)
        self.R.index = [('CA-' + i[0], i[1]) for i in self.R.index]
        self.R.index = pd.MultiIndex.from_tuples(self.R.index)
        self.W.columns = [('CA-' + i[0], i[1]) for i in self.W.columns]
        self.W.columns = pd.MultiIndex.from_tuples(self.W.columns)
        self.W.index = [('CA-' + i[0], i[1]) for i in self.W.index]
        self.W.index = pd.MultiIndex.from_tuples(self.W.index)
        self.WY.columns = [('CA-' + i[0], i[1], i[2]) for i in self.WY.columns]
        self.WY.columns = pd.MultiIndex.from_tuples(self.WY.columns)
        self.WY.index = [('CA-' + i[0], i[1]) for i in self.WY.index]
        self.WY.index = pd.MultiIndex.from_tuples(self.WY.index)
        self.U.columns = [('CA-' + i[0], i[1]) for i in self.U.columns]
        self.U.columns = pd.MultiIndex.from_tuples(self.U.columns)
        self.U.index = [('CA-' + i[0], i[1]) for i in self.U.index]
        self.U.index = pd.MultiIndex.from_tuples(self.U.index)
        self.V.columns = [('CA-' + i[0], i[1]) for i in self.V.columns]
        self.V.columns = pd.MultiIndex.from_tuples(self.V.columns)
        self.V.index = [('CA-' + i[0], i[1]) for i in self.V.index]
        self.V.index = pd.MultiIndex.from_tuples(self.V.index)
        self.C.columns = [('CA-' + i[0], i[1], i[2]) if (type(i) == tuple and 'Total Water Consumption' not in i[1])
                          else i for i in self.C.columns]
        self.F.columns = [('CA-' + i[0], i[1]) for i in self.F.columns]
        self.F.columns = pd.MultiIndex.from_tuples(self.F.columns)
        self.F.index = [('CA-' + i[0], i[1], i[2]) if type(i) == tuple else i for i in self.F.index]
        self.minerals.columns = [('CA-' + i[0], i[1]) for i in self.minerals.columns]
        self.minerals.columns = pd.MultiIndex.from_tuples(self.minerals.columns)
        self.FY.columns = [('CA-' + i[0], i[1], i[2]) for i in self.FY.columns]
        self.FY.columns = pd.MultiIndex.from_tuples(self.FY.columns)
        self.FY.index = [('CA-' + i[0], i[1], i[2]) if type(i) == tuple else i for i in self.FY.index]
        self.g.index = [('CA-' + i[0], i[1]) for i in self.g.index]
        self.g.columns = [('CA-' + i[0], i[1]) for i in self.g.columns]
        self.g.index = pd.MultiIndex.from_tuples(self.g.index)
        self.g.columns = pd.MultiIndex.from_tuples(self.g.columns)
        self.inv_g.columns = self.g.columns
        self.inv_g.index = self.g.columns
        self.q.index = [('CA-' + i[0], i[1]) for i in self.q.index]
        self.q.columns = [('CA-' + i[0], i[1]) for i in self.q.columns]
        self.q.index = pd.MultiIndex.from_tuples(self.q.index)
        self.q.columns = pd.MultiIndex.from_tuples(self.q.columns)
        self.inv_q.columns = self.q.index
        self.inv_q.index = self.q.index

    def refine_meat_sector(self):
        """
        Because the meat sector is aggregated into one sector, the economic allocation from the technology-industry
        construct creates some issues. For instance, the Quebec products of cattle sector is composed of 85% purchases
        of Hogs, because Quebec mainly produces Hogs and not cattle. So, we refine the definition of these meat sectors
        by forcing the products of cattle sectors to only buy Cattle and not Hogs.
        :return:
        """

        meat_transfo = ['Fresh and frozen beef and veal', 'Fresh and frozen pork',
                        'Fresh and frozen poultry of all types',
                        'Products of meat cattle', 'Products of meat pigs', 'Products of meat poultry']
        meat_breeding = ['Cattle and calves', 'Hogs', 'Poultry', 'Pigs', 'Cattle']

        for province in ['CA-AB', 'CA-BC', 'CA-MB', 'CA-NB', 'CA-NL', 'CA-NS', 'CA-NT',
                         'CA-NU', 'CA-ON', 'CA-PE', 'CA-QC', 'CA-SK', 'CA-YT']:
            meat_sector = 'Fresh and frozen beef and veal'
            total_meat_breeding = self.A.loc(axis=0)[:, meat_breeding].loc(axis=1)[province, meat_sector].sum()
            total_meat_transfo = self.A.loc(axis=0)[:, meat_transfo].loc(axis=1)[province, meat_sector].sum()
            # rescale meat_sector
            if self.A.loc[[i for i in self.A.index if i[1] in ['Cattle and calves', 'Cattle']],
                          [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum().sum() != 0:
                self.A.loc[[i for i in self.A.index if i[1] in ['Cattle and calves', 'Cattle']],
                           [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]] /= (
                        self.A.loc[[i for i in self.A.index if i[1] in ['Cattle and calves', 'Cattle']],
                                   [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum() /
                        total_meat_breeding)
            if self.A.loc[[i for i in self.A.index if i[1] in ['Fresh and frozen beef and veal',
                                                               'Products of meat cattle']],
                          [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum().sum() != 0:
                self.A.loc[[i for i in self.A.index if i[1] in ['Fresh and frozen beef and veal',
                                                                'Products of meat cattle']],
                           [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]] /= (
                        self.A.loc[[i for i in self.A.index if i[1] in ['Fresh and frozen beef and veal',
                                                                        'Products of meat cattle']],
                                   [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum() /
                        total_meat_transfo)
            # remove other meats
            self.A.loc[[i for i in self.A.index if i[1] in ['Pigs',
                                                            'Hogs',
                                                            'Poultry',
                                                            'Fresh and frozen pork',
                                                            'Fresh and frozen poultry of all types',
                                                            'Products of meat pigs',
                                                            'Products of meat poultry']],
                       [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]] = 0

            meat_sector = 'Fresh and frozen pork'
            total_meat_breeding = self.A.loc(axis=0)[:, meat_breeding].loc(axis=1)[province, meat_sector].sum()
            total_meat_transfo = self.A.loc(axis=0)[:, meat_transfo].loc(axis=1)[province, meat_sector].sum()
            # rescale meat_sector
            if self.A.loc[[i for i in self.A.index if i[1] in ['Pigs', 'Hogs']],
                          [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum().sum() != 0:
                self.A.loc[[i for i in self.A.index if i[1] in ['Pigs', 'Hogs']],
                           [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]] /= (
                        self.A.loc[[i for i in self.A.index if i[1] in ['Pigs', 'Hogs']],
                                   [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum() /
                        total_meat_breeding)
            if self.A.loc[[i for i in self.A.index if i[1] in ['Fresh and frozen pork',
                                                               'Products of meat pigs']],
                          [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum().sum() != 0:
                self.A.loc[[i for i in self.A.index if i[1] in ['Fresh and frozen pork',
                                                                'Products of meat pigs']],
                           [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]] /= (
                        self.A.loc[[i for i in self.A.index if i[1] in ['Fresh and frozen pork',
                                                                        'Products of meat pigs']],
                                   [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum() /
                        total_meat_transfo)
            # remove other meats
            self.A.loc[[i for i in self.A.index if i[1] in ['Cattle and calves',
                                                            'Cattle',
                                                            'Poultry',
                                                            'Fresh and frozen beef and veal',
                                                            'Fresh and frozen poultry of all types',
                                                            'Products of meat cattle',
                                                            'Products of meat poultry']],
                       [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]] = 0

            meat_sector = 'Fresh and frozen poultry of all types'
            total_meat_breeding = self.A.loc(axis=0)[:, meat_breeding].loc(axis=1)[province, meat_sector].sum()
            total_meat_transfo = self.A.loc(axis=0)[:, meat_transfo].loc(axis=1)[province, meat_sector].sum()
            # rescale meat_sector
            if self.A.loc[[i for i in self.A.index if i[1] in ['Poultry']],
                          [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum().sum() != 0:
                self.A.loc[[i for i in self.A.index if i[1] in ['Poultry']],
                           [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]] /= (
                        self.A.loc[[i for i in self.A.index if i[1] in ['Poultry']],
                                   [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum() /
                        total_meat_breeding)
            if self.A.loc[[i for i in self.A.index if i[1] in ['Fresh and frozen poultry of all types',
                                                               'Products of meat poultry']],
                          [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum().sum() != 0:
                self.A.loc[[i for i in self.A.index if i[1] in ['Fresh and frozen poultry of all types',
                                                                'Products of meat poultry']],
                           [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]] /= (
                        self.A.loc[[i for i in self.A.index if i[1] in ['Fresh and frozen poultry of all types',
                                                                        'Products of meat poultry']],
                                   [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]].sum() /
                        total_meat_transfo)

            # remove other meats
            self.A.loc[[i for i in self.A.index if i[1] in ['Cattle and calves',
                                                            'Cattle',
                                                            'Pigs',
                                                            'Hogs',
                                                            'Fresh and frozen beef and veal',
                                                            'Fresh and frozen pork',
                                                            'Products of meat cattle',
                                                            'Products of meat pigs']],
                       [i for i in self.A.columns if (i[0] == province and i[1] == meat_sector)]] = 0

    def convert_F_to_commodity(self):
        """
        Converts F to emissions x commodity format.
        :return:
        """
        self.F = self.F.dot(self.V.dot(self.inv_g).T)
        self.F = pd.concat([self.F, self.minerals])

    def normalize_flows(self):
        """
        Produce normalized environmental extensions
        :return: self.S and self.F with product classification if it's been selected
        """
        self.F = self.F.astype(float)
        self.S = self.F.dot(self.inv_q)
        self.S = pd.concat([self.S, self.S_exio]).fillna(0)
        self.S = self.S.groupby(self.S.index).sum()
        self.S = self.S.reindex(self.A.columns, axis=1)
        # change provinces metadata for S here
        self.S.columns = self.A.columns

        # adding empty flows to FY to allow multiplication with self.C
        self.FY = pd.concat([pd.DataFrame(0, self.F.index, self.Y.columns), self.FY])
        self.FY = self.FY.groupby(self.FY.index).sum()
        self.FY = self.FY.reindex(self.C.columns).fillna(0)

        self.emission_metadata = pd.concat([self.emission_metadata, self.unit_exio])

    def add_hfc_emissions(self):
        """
        Method matching HFCs accounts to IOIC product classification.
        Source of the data: https://doi.org/10.5281/zenodo.7432088
        :return: self.F with HFCs flows included
        """

        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            with open(pkg_resources.resource_filename(__name__, '/Data/Concordances/2014-2021/concordance_HFCs.json'), 'r') as f:
                industry_mapping = json.load(f)
        elif self.year in [2022]:
            with open(pkg_resources.resource_filename(__name__, '/Data/Concordances/2022/concordance_HFCs.json'), 'r') as f:
                industry_mapping = json.load(f)

        # environmental values
        all_GES = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Environmental_data/SF6_HFC_PFC_emissions.xlsx'), 'Canada_UNFCCC_emissions', index_col=0)
        all_GES = all_GES.dropna(axis=0, subset=['numberValue'])

        # keep hfcs only
        hfcs_list = ['C10F18', 'C2F6', 'C3F8', 'C4F10', 'C5F12', 'C6F14', 'c-C3F6', 'c-C4F8', 'CF4', 'CO', 'HFC-125',
                     'HFC-134', 'HFC-134a', 'HFC-143', 'HFC-143a', 'HFC-152', 'HFC-152a', 'HFC-161', 'HFC-227ea', 'HFC-23',
                     'HFC-236cb', 'HFC-236ea', 'HFC-236fa', 'HFC-245ca', 'HFC-245fa', 'HFC-32', 'HFC-365mfc', 'HFC-41',
                     'HFC-43-10mee', 'NF3', 'SF6']
        hfcs = all_GES.loc[all_GES['gas'].isin(hfcs_list)].copy()
        hfcs['category'] = hfcs['category'].replace("\s+", " ", regex=True).str.strip()
        hfcs = hfcs.loc[hfcs['category'].isin([k for k in industry_mapping.keys()])].reset_index(drop=True)

        # convert to kg
        mass_conversion_factors = {"kg": 1, "t": 1000, "kt": 1000000}
        hfcs['numberValue'] = hfcs['unit'].map(mass_conversion_factors).mul(hfcs['numberValue'])

        if self.year != 2022:
            hfcs = hfcs.loc[(hfcs['measure'] == 'Net emissions/removals') & (hfcs['year'] == self.year),
                            ['category', 'gas', 'numberValue']]
        else:
            hfcs = hfcs.loc[(hfcs['measure'] == 'Net emissions/removals') & (hfcs['year'] == 2021),
                            ['category', 'gas', 'numberValue']]

        hfcs = hfcs.set_index(['gas', 'category'])
        hfcs.rename_axis(['gas', 'category'])

        # economic values
        V = self.V.T.groupby(level=0).sum().T

        R = pd.DataFrame(columns=V.index)

        for k, io_sectors in industry_mapping.items():
            temp = V.iloc[V.index.get_level_values(1).isin([i for i in io_sectors])]
            # economic allocation
            temp = temp.div(temp.values.sum())

            temp = temp.reindex(V.index).fillna(0)
            temp = temp.T
            temp = temp.assign(category=k).set_index('category', append=True)
            R = pd.concat([R, temp])

        idx = pd.MultiIndex.from_tuples(R.index, names=['province', 'category'])
        R = R.reset_index(drop=True).set_index(idx)
        R.columns = R.columns.values
        # combine economic and environmental values (left outer join)
        W = hfcs.join(R)
        W.loc[:, W.columns != 'numberValue'] = W.loc[:, W.columns != 'numberValue'].mul(W['numberValue'], 0)
        W = W.drop(columns='numberValue').groupby(level=['gas', 'province']).sum()
        columns = pd.MultiIndex.from_tuples(W.columns)
        W = W.reindex(columns, axis='columns')
        W = W.assign(compartiment='Air').set_index('compartiment', append=True).reorder_levels(['province', 'gas',
                                                                                                'compartiment'])
        W.index.names = len(W.index.names)*[None]

        if W.columns.equals(self.F.columns):
            self.F = pd.concat([self.F, W])

    def add_water_consumption_flows_for_livestock_and_crops(self):
        """
        Method adding water consumption data from livestocks and crops, to openIO-Canada's own livestock and crop
        commodities.
        :return: self.F
        """

        water_conso = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Environmental_data/Water_consumption_values.xlsx'), None)

        # ------------------------------------------ Livestocks ------------------------------------------------------
        # split water consumption between provinces
        gotta_split = water_conso['Livestock'].loc[[i for i in water_conso['Livestock'].index
                                                    if (',' in water_conso['Livestock'].loc[i, 'Canada Province'] or
                                                        water_conso['Livestock'].loc[
                                                            i, 'Canada Province'] == 'Canada')]].copy('deep')

        for ix in gotta_split.index:
            if gotta_split.loc[ix, 'Canada Province'] == 'Canada':
                # territories are excluded from water consumption data
                provinces = ['AB', 'BC', 'MB', 'NB', 'NL', 'NS', 'PE', 'ON', 'QC', 'SK']
            else:
                provinces = gotta_split.loc[ix, 'Canada Province'].split(', ')
            # select the supply of the commodity by the provinces to be split
            df = self.V.loc[[i for i in self.V.index if (i[0].split('CA-')[1] in provinces and
                                                         i[1] == {i[0]: i[1] for i in self.commodities}[
                                                             gotta_split.loc[ix, 'Code product Open IO']])]].sum(1)
            # determine share
            df /= df.sum()
            # apply share to water consumption
            df *= gotta_split.loc[ix, 'Water Consumption (m3)']
            # reformat data
            df = df.reset_index()
            df.columns = ['Canada Province', 'Code product Open IO', 'Water Consumption (m3)']
            df.loc[:, 'Canada Province'] = [i.split('CA-')[1] for i in df.loc[:, 'Canada Province']]
            df.loc[:, 'Code product Open IO'] = [{j[1]: j[0] for j in self.commodities}[i] for i in
                                                 df.loc[:, 'Code product Open IO']]
            df.loc[:, 'Livestock Category'] = gotta_split.loc[ix, 'Livestock Category']
            df.loc[:, 'Year'] = gotta_split.loc[ix, 'Year']
            # concat with original water data
            water_conso['Livestock'] = pd.concat([water_conso['Livestock'], df])

        water_conso['Livestock'].index = [i for i in range(0, len(water_conso['Livestock'].index))]
        water_conso['Livestock'] = water_conso['Livestock'].drop(
            [i for i in water_conso['Livestock'].index if (',' in water_conso['Livestock'].loc[i, 'Canada Province'] or
                                                           water_conso['Livestock'].loc[i, 'Canada Province'] == 'Canada')])
        water_conso['Livestock'].index = [i for i in range(0, len(water_conso['Livestock'].index))]
        # for livestocks where only one year is available, copy paste for other years
        unique_livestocks = set(water_conso['Livestock'].loc[:, 'Livestock Category'])

        for livestock in unique_livestocks:
            if len(set(water_conso['Livestock'].loc[
                           water_conso['Livestock'].loc[:, 'Livestock Category'] == livestock, 'Year'])) == 1:
                for year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2022]:
                    df = water_conso['Livestock'].loc[
                        water_conso['Livestock'].loc[:, 'Livestock Category'] == livestock].loc[
                        water_conso['Livestock'].loc[:, 'Year'] == 2021].copy('deep')
                    df.Year = year
                    water_conso['Livestock'] = pd.concat([water_conso['Livestock'], df])
                    water_conso['Livestock'].index = [i for i in range(0, len(water_conso['Livestock'].index))]

        # select the year of study
        water_conso['Livestock'] = water_conso['Livestock'].loc[water_conso['Livestock'].Year == self.year]
        # groupby openIO commodity
        water_conso['Livestock'] = water_conso['Livestock'].drop(['Livestock Category', 'Year'], axis=1)
        water_conso['Livestock'] = water_conso['Livestock'].groupby(
            by=['Canada Province', 'Code product Open IO']).sum().reset_index()
        # reformat data
        water_conso['Livestock'].loc[:, 'Canada Province'] = [
            'CA-' + i for i in water_conso['Livestock'].loc[:, 'Canada Province']]
        water_conso['Livestock'].index = [
            (i, 'Water', 'Water') for i in water_conso['Livestock'].loc[:, 'Canada Province']]
        water_conso['Livestock'].loc[:, 'Code product Open IO'] = [{i[0]: i[1] for i in self.commodities}[i] for i in
                                                                   water_conso['Livestock'].loc[
                                                                   :, 'Code product Open IO']]
        water_conso['Livestock'] = water_conso['Livestock'].reset_index().pivot(values='Water Consumption (m3)',
                                                                                columns=['Canada Province',
                                                                                         'Code product Open IO'],
                                                                                index=['index']).fillna(0)
        water_conso['Livestock'].index.name = None
        water_conso['Livestock'].columns.names = [None, None]

        self.F.loc[[i for i in self.F.index if i[1] == 'Water'], water_conso['Livestock'].columns] = water_conso[
            'Livestock']
        # fillna(0) because there is no data for water consumption of livestocks in the territories of Canada
        self.F = self.F.fillna(0)

        # ----------------------------------------------- Crops ------------------------------------------------------
        water_conso['Crop'] = water_conso['Crop'].drop(
            water_conso['Crop'][water_conso['Crop'].loc[:, 'Crop Category'] == 'Total'].index)
        water_conso['Crop'] = water_conso['Crop'].drop(['Crop Category'], axis=1)
        # select year
        if self.year == 2022:
            water_conso['Crop'] = water_conso['Crop'].loc[:, ['Code product Open IO'] +
                                                             [i for i in water_conso['Crop'].columns if str(2021) in i]]
        else:
            water_conso['Crop'] = water_conso['Crop'].loc[:, ['Code product Open IO'] +
                                                             [i for i in water_conso['Crop'].columns if
                                                              str(self.year) in i]]
        water_consumption_crop = pd.DataFrame()

        for comm in water_conso['Crop'].loc[:, 'Code product Open IO']:
            df = self.V.loc(axis=0)[:, {i[0]: i[1] for i in self.commodities}[comm]].sum(1)
            df /= df.sum()
            if self.year != 2022:
                df *= water_conso['Crop'].loc[
                    water_conso['Crop'].loc[:, 'Code product Open IO'] == comm, 'Water consumption (m3) ' + str(
                        self.year)].iloc[0]
            elif self.year == 2022:
                df *= water_conso['Crop'].loc[
                    water_conso['Crop'].loc[:, 'Code product Open IO'] == comm, 'Water consumption (m3) ' + str(
                        2021)].iloc[0]
            water_consumption_crop = pd.concat([water_consumption_crop, df])

        water_consumption_crop.index = pd.MultiIndex.from_tuples(water_consumption_crop.index)
        # spatializing
        water_consumption_crop = pd.concat([water_consumption_crop] * len(self.matching_dict), axis=1)
        water_consumption_crop.columns = [('CA-' + i, 'Water', 'Water') for i in self.matching_dict.keys()]
        for province in water_consumption_crop.index.levels[0]:
            water_consumption_crop.loc[province, water_consumption_crop.columns.get_level_values(0) != province] = 0
        water_consumption_crop = water_consumption_crop.T
        self.F.loc[
            [i for i in self.F.index if i[1] == 'Water'], water_consumption_crop.columns] = water_consumption_crop

    def add_plastic_emissions(self):
        """
        Method adding plastic waste generated from producing certain commodities.
        :return: self.F
        """

        # load data
        plastics_data = pd.read_csv(pkg_resources.resource_stream(
            __name__, '/Data/Environmental_data/plastic_account_by_product_category.csv'), low_memory=False)
        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            with open(pkg_resources.resource_filename(__name__, "Data/Concordances/2014-2021/plastic_kpis.json"), 'r') as f:
                kpis = json.load(f)
            with open(pkg_resources.resource_filename(__name__, "Data/Concordances/2014-2021/plastic_mapping_ppfa_openio.json"), 'r') as f:
                map_plastic_data_to_io = json.load(f)
            with open(pkg_resources.resource_filename(__name__, "Data/Concordances/2014-2021/plastic_mapping_exio_openio.json"), 'r') as f:
                plastic_openio_to_exio = json.load(f)
            oecd_data = pd.read_excel(pkg_resources.resource_stream(
                __name__, '/Data/Environmental_data/OECD_plastic_waste_management.xlsx'), None)
            with open(pkg_resources.resource_filename(__name__, "Data/Concordances/2014-2021/country_mapping_exio_oecd.json"), 'r') as f:
                map_exio_to_oecd = json.load(f)
        elif self.year in [2022]:
            with open(pkg_resources.resource_filename(__name__, "Data/Concordances/2022/plastic_kpis.json"), 'r') as f:
                kpis = json.load(f)
            with open(pkg_resources.resource_filename(__name__, "Data/Concordances/2022/plastic_mapping_ppfa_openio.json"), 'r') as f:
                map_plastic_data_to_io = json.load(f)
            with open(pkg_resources.resource_filename(__name__, "Data/Concordances/2022/plastic_mapping_exio_openio.json"), 'r') as f:
                plastic_openio_to_exio = json.load(f)
            oecd_data = pd.read_excel(pkg_resources.resource_stream(
                __name__, '/Data/Environmental_data/OECD_plastic_waste_management.xlsx'), None)
            with open(pkg_resources.resource_filename(__name__, "Data/Concordances/2022/country_mapping_exio_oecd.json"), 'r') as f:
                map_exio_to_oecd = json.load(f)

        # change IOCC codes to names within mapping
        mapping_commodities = {i[0]: i[1] for i in self.commodities}
        map_plastic_data_to_io = {k: [mapping_commodities[i] for i in v] for k, v in map_plastic_data_to_io.items()}

        # ---------------------------------- Data pre-treatment ---------------------------------------------
        # select year of study
        if self.year != 2022:
            plastics_data = plastics_data.loc[plastics_data.REF_DATE ==
                                              min(plastics_data.REF_DATE, key=lambda x:abs(x-self.year))].copy('deep')
        else:
            plastics_data = plastics_data.loc[plastics_data.REF_DATE ==
                                              min(plastics_data.REF_DATE, key=lambda x:abs(x-2021))].copy('deep')

        # format province names
        plastics_data.GEO = [
            'CA-' + {v: k for k, v in self.matching_dict.items()}[i] if i in list(self.matching_dict.values()) else i
            for i in plastics_data.GEO]

        # pivot table to manipulate data easily
        plastics_data = plastics_data.pivot_table(values='VALUE', index=['Product category', 'Variable'],
                                                  columns=['GEO'], dropna=False, fill_value=0)

        # Some packaging items (bottles, film, ...) are only available for Canada. We estimate them for provinces.
        packaging_distribution = (
                    plastics_data.loc[['Bottles', 'Film', 'Non-bottle rigid', 'Other packaging products'], 'Canada'] /
                    plastics_data.loc['Packaging', 'Canada'])

        for province in ['CA-' + i for i in self.matching_dict]:
            for packaging in ['Bottles', 'Film', 'Non-bottle rigid', 'Other packaging products']:
                plastics_data.loc[packaging, province] = (
                            plastics_data.loc['Packaging', province] * packaging_distribution.loc[packaging]).values

        landfill_rate = (plastics_data.loc(axis=0)['Total, all product categories',
             'Disposed plastic waste and scrap sent to landfill or incinerated without energy recovery'].loc['Canada'] /
                         plastics_data.loc(axis=0)['Total, all product categories',
                                                   'Total disposed plastic waste and scrap'].loc['Canada'])

        incineration_rate = (plastics_data.loc(axis=0)['Total, all product categories',
           'Disposed plastic waste and scrap sent for incineration or gasification with energy recovery'].loc['Canada'] /
                             plastics_data.loc(axis=0)['Total, all product categories',
                                                       'Total disposed plastic waste and scrap'].loc['Canada'])

        plastics_data.loc(axis=0)[:,
        'Disposed plastic waste and scrap sent to landfill or incinerated without energy recovery'] = (
                plastics_data.loc(axis=0)[:, 'Total disposed plastic waste and scrap'] * landfill_rate).values

        plastics_data.loc(axis=0)[:,
        'Disposed plastic waste and scrap sent for incineration or gasification with energy recovery'] = (
                plastics_data.loc(axis=0)[:, 'Total disposed plastic waste and scrap'] * incineration_rate).values
        plastics_data = plastics_data.drop(['Canada', 'Canadian territorial enclaves abroad'], axis=1)

        self.F = pd.concat([self.F, pd.DataFrame(0, index=kpis, columns=self.F.columns)])
        self.F_exio = pd.concat([self.F_exio, pd.DataFrame(0, index=kpis, columns=self.F_exio.columns)])

        # ----------------------------- Plastic physical flow account ------------------------------------------------
        Z = self.A.loc[:, ['CA-' + i for i in self.matching_dict]] * self.q.sum(1)

        for category in map_plastic_data_to_io.keys():
            for province in self.matching_dict:

                # determine in which products of the category is plastic embedded (1.5 for broad € to CAD conversion)
                where_is_plastic = Z.loc(axis=0)[:, ['Plastic resins']].sum() + Z.loc(axis=0)[:,
                                                                                ['Plastics, basic']].sum() * 1.5
                if where_is_plastic.loc(axis=0)['CA-' + province, map_plastic_data_to_io[category]].sum() != 0:
                    where_is_plastic = (
                                where_is_plastic.loc(axis=0)['CA-' + province, map_plastic_data_to_io[category]] /
                                where_is_plastic.loc(axis=0)['CA-' + province, map_plastic_data_to_io[category]].sum())
                # if the province does not purchase any of the products of the category, use the national distribution
                else:
                    where_is_plastic = (where_is_plastic.loc(axis=0)[:, map_plastic_data_to_io[category]] /
                                        where_is_plastic.loc(axis=0)[:, map_plastic_data_to_io[category]].sum())

                for product in map_plastic_data_to_io[category]:
                    # determine production within Canada of studied product
                    if self.endogenizing:
                        product_from_canada = (self.U.loc(axis=0)[:, product].loc[:, 'CA-' + province].sum(1) +
                                               self.Y.drop([i for i in self.Y.columns if "Changes in inventories" in i[1]],
                                                           axis=1).loc(axis=0)[:, product].loc[:, 'CA-' + province].sum(1) +
                                               self.K.loc(axis=0)[:, product].loc[:, 'CA-' + province].sum(1))

                        # determine imports of studied product
                        product_imports = (
                                    self.imports_scaled_U.loc(axis=0)[:, product].loc[:, province].sum(1) +
                                    self.imports_scaled_K.loc(axis=0)[:, product].loc[:, province].sum(1) +
                                    self.imports_scaled_Y.loc[:,
                                    [i for i in self.imports_scaled_Y.columns if (
                                            "Changes in inventories" not in i[1] and i[0] == province)]].loc(axis=0)[:,
                                    product].sum(1))
                    else:
                        product_from_canada = (self.U.loc(axis=0)[:, product].loc[:, 'CA-' + province].sum(1) +
                                               self.Y.drop([i for i in self.Y.columns if "Changes in inventories" in i[1]],
                                                   axis=1).loc(axis=0)[:, product].loc[:, 'CA-' + province].sum(1))

                        # determine imports of studied product
                        product_imports = (
                                self.imports_scaled_U.loc(axis=0)[:, product].loc[:, province].sum(1) +
                                self.imports_scaled_Y.loc[:,
                                [i for i in self.imports_scaled_Y.columns if (
                                        "Changes in inventories" not in i[1] and i[0] == province)]].loc(axis=0)[:,
                                product].sum(1))

                    # translate imports of products of openIO into products of exiobase
                    distrib_market_imported_category = (self.A.loc[
                        self.A_exio.index, ['CA-' + i for i in self.matching_dict]].dot(
                        self.q.sum(1))).loc(axis=0)[:, plastic_openio_to_exio[category]]
                    distrib_market_imported_category.index = pd.MultiIndex.from_tuples(
                        distrib_market_imported_category.index)
                    for ix in distrib_market_imported_category.index.levels[0]:
                        distrib_market_imported_category.loc[ix] = (distrib_market_imported_category.loc[ix] /
                                                                    distrib_market_imported_category.loc[
                                                                        ix].sum()).values
                    for ix in product_imports.index.levels[0]:
                        distrib_market_imported_category.loc[(ix, plastic_openio_to_exio[category])] = (
                                distrib_market_imported_category.loc[(ix, plastic_openio_to_exio[category])] *
                                product_imports.loc[ix, product])
                    distrib_market_imported_category = distrib_market_imported_category.fillna(0)

                    # get final market distribution of studied product (local + imports)
                    product_market_in_province = pd.concat([product_from_canada, distrib_market_imported_category])
                    product_market_in_province /= product_market_in_province.sum()
                    # if 0/0 (provinces not using any of the plastic product) fill NaN to 0s
                    product_market_in_province = product_market_in_province.fillna(0)
                    product_market_in_province = pd.concat([product_market_in_province] * len(kpis), axis=1)
                    product_market_in_province.columns = pd.MultiIndex.from_product([[province], kpis])

                    if where_is_plastic.loc(axis=0)['CA-' + province, map_plastic_data_to_io[category]].sum() != 0:
                        # apply distribution to plastic emissions
                        product_market_in_province = (product_market_in_province.loc[:, province] *
                                                      plastics_data.loc[category].loc[kpis].loc[:, 'CA-' + province] *
                                                      where_is_plastic.loc(axis=0)['CA-' + province, product]).fillna(0).T
                    else:
                        product_market_in_province = (product_market_in_province.loc[:, province] *
                                                      plastics_data.loc[category].loc[kpis].loc[:, 'CA-' + province] *
                                                      where_is_plastic.groupby(axis=0, level=1).sum().loc[
                                                          product]).fillna(0).T

                    # add plastic emissions to emitting production processes
                    self.F.loc[product_market_in_province.index,
                               [i for i in product_market_in_province.columns if 'CA-' in i[0]]] += (
                        product_market_in_province.loc[:, [i for i in product_market_in_province.columns if 'CA-' in i[0]]])

                    self.F_exio.loc[product_market_in_province.index,
                                    [i for i in product_market_in_province.columns if  'CA-' not in i[0]]] += (
                        product_market_in_province.loc[:, [i for i in product_market_in_province.columns if 'CA-' not in i[0]]])

        del Z

        # ---------------------------------------- OECD data --------------------------------------
        def format_oecd_data(dataframe):
            oecd_data_index = list(zip(dataframe.iloc[2:, 0].ffill(),
                                       dataframe.iloc[2:, 1].ffill(),
                                       dataframe.iloc[2:, 2].ffill()))
            oecd_data_index = [(i[0], i[1].lstrip(), i[2].lstrip()) for i in oecd_data_index]
            oecd_data_index = pd.MultiIndex.from_tuples(oecd_data_index)

            oecd_data_cols = dataframe.iloc[0, 4:].astype(int).values

            dataframe = dataframe.iloc[2:, 4:]
            dataframe.index = oecd_data_index
            dataframe.columns = oecd_data_cols

            return dataframe

        for tab in oecd_data:
            oecd_data[tab] = format_oecd_data(oecd_data[tab])
            # OECD data in millions of tonnes of plastic waste -> tonnes of plastic waste
            oecd_data[tab] *= 1000000
            oecd_data[tab] = oecd_data[tab].loc[:, min(oecd_data[tab].columns, key=lambda x: abs(x - self.year))]

        plastic_waste_oecd = pd.DataFrame()
        for tab in oecd_data:
            plastic_waste_oecd = pd.concat([plastic_waste_oecd, oecd_data[tab]], axis=1)
        plastic_waste_oecd.columns = list(oecd_data.keys())

        # add mismanaged plastic waste category for waste outside of Canada (through OECD data)
        self.F = pd.concat([self.F, pd.DataFrame(0, index=[
            'Disposed plastic waste and scrap mismanaged (open dump, open pits, unsanitary landfills)'],
                                                 columns=self.F.columns)])
        self.F_exio = pd.concat([self.F_exio, pd.DataFrame(0, index=[
            'Disposed plastic waste and scrap mismanaged (open dump, open pits, unsanitary landfills)'],
                                                           columns=self.F_exio.columns)])

        for country in self.Z_exio.columns.levels[0]:
            if country != 'CA':
                # determine the distribution for conversion from region (e.g., OECD EU) to countries (AT, BE, etc.)
                plastic_waste_to_incineration_in_region = self.Z_exio.loc(axis=0)[:,
                                                          'Plastic waste for treatment: incineration'].loc[:,
                                                         list({k for k, v in map_exio_to_oecd.items() if
                                                           v == map_exio_to_oecd[country]})].sum().groupby(level=0).sum()
                plastic_waste_to_incineration_in_region /= plastic_waste_to_incineration_in_region.sum()
                plastic_waste_to_landfill_in_region = self.Z_exio.loc(axis=0)[:,
                                                      'Plastic waste for treatment: landfill'].loc[:,
                                                      list({k for k, v in map_exio_to_oecd.items() if
                                                       v == map_exio_to_oecd[country]})].sum().groupby(level=0).sum()
                plastic_waste_to_landfill_in_region /= plastic_waste_to_landfill_in_region.sum()
                plastic_waste_total_in_region = self.Z_exio.loc(axis=0)[:, ['Plastic waste for treatment: incineration',
                                                                       'Plastic waste for treatment: landfill']].loc[:,
                                                list({k for k, v in map_exio_to_oecd.items() if
                                                 v == map_exio_to_oecd[country]})].sum().groupby(level=0).sum()
                plastic_waste_total_in_region /= plastic_waste_total_in_region.sum()

                # determine how much of the products were bought by the country
                diag = pd.DataFrame(np.diag(self.Z_exio.loc[:, country].sum(1)), self.Z_exio.index, self.Z_exio.index)

                # translate these transactions into amount of "Plastics, basic" bought
                distrib_country = (self.Z_exio.loc(axis=0)[:, 'Plastics, basic'].dot(diag).sum() /
                                   self.Z_exio.loc(axis=0)[:, 'Plastics, basic'].dot(diag).sum().sum())

                indicator = 'Recycled plastic pellets and flakes ready for use in production of new products'
                df = pd.DataFrame(self.F_exio.loc[indicator] + (
                            distrib_country * plastic_waste_oecd.loc[[map_exio_to_oecd[country]], 'Recycled'].iloc[0] *
                            plastic_waste_total_in_region.loc[country]), columns=[indicator])
                self.F_exio = self.F_exio.drop(indicator)
                self.F_exio = pd.concat([self.F_exio, df.T])

                indicator = 'Disposed plastic waste and scrap sent to landfill or incinerated without energy recovery'
                df = pd.DataFrame(self.F_exio.loc[indicator] + (
                            distrib_country * plastic_waste_oecd.loc[[map_exio_to_oecd[country]], 'Landfilled'].iloc[0] *
                            plastic_waste_to_landfill_in_region.loc[country]), columns=[indicator])
                self.F_exio = self.F_exio.drop(indicator)
                self.F_exio = pd.concat([self.F_exio, df.T])

                # incineration in OECD is not specifying energy recovery or not -> assume 50/50 split in first version
                split_incineration_with_or_without_energy_recovery = 0.5

                indicator = 'Disposed plastic waste and scrap sent to landfill or incinerated without energy recovery'
                df = pd.DataFrame(self.F_exio.loc[indicator] + (
                            distrib_country * plastic_waste_oecd.loc[[map_exio_to_oecd[country]], 'Incinerated'].iloc[0] *
                            plastic_waste_to_incineration_in_region.loc[country] *
                            split_incineration_with_or_without_energy_recovery), columns=[indicator])
                self.F_exio = self.F_exio.drop(indicator)
                self.F_exio = pd.concat([self.F_exio, df.T])
                indicator = 'Disposed plastic waste and scrap sent for incineration or gasification with energy recovery'
                df = pd.DataFrame(self.F_exio.loc[indicator] + (
                            distrib_country * plastic_waste_oecd.loc[[map_exio_to_oecd[country]], 'Incinerated'].iloc[0] *
                            plastic_waste_to_incineration_in_region.loc[country] *
                            split_incineration_with_or_without_energy_recovery), columns=[indicator])
                self.F_exio = self.F_exio.drop(indicator)
                self.F_exio = pd.concat([self.F_exio, df.T])

                indicator = 'Plastic leaked permanently into the environment'
                df = pd.DataFrame(self.F_exio.loc[indicator] + (
                            distrib_country * plastic_waste_oecd.loc[[map_exio_to_oecd[country]], 'Littered'].iloc[0] *
                            plastic_waste_total_in_region.loc[country]), columns=[indicator])
                self.F_exio = self.F_exio.drop(indicator)
                self.F_exio = pd.concat([self.F_exio, df.T])

                indicator = 'Disposed plastic waste and scrap mismanaged (open dump, open pits, unsanitary landfills)'
                df = pd.DataFrame(self.F_exio.loc[indicator] + (
                            distrib_country * plastic_waste_oecd.loc[[map_exio_to_oecd[country]], 'Mismanaged'].iloc[0] *
                            plastic_waste_total_in_region.loc[country]), columns=[indicator])
                self.F_exio = self.F_exio.drop(indicator)
                self.F_exio = pd.concat([self.F_exio, df.T])

        # ------------------------------- Specify plastic resins ---------------------------------

        resin_compo = pd.read_excel(pkg_resources.resource_stream(
            __name__, '/Data/Environmental_data/plastic_resin_compositions.xlsx'), None)

        def formatting(dff):
            dff = dff.set_index('Unnamed: 0')
            dff.index.name = None
            return dff.fillna(0)

        for sheet in resin_compo:
            resin_compo[sheet] = formatting(resin_compo[sheet])

        # create new indicators including information about resin (e.g., recycled - PET)
        resin_indicators = []
        for indicator in kpis + [
            'Disposed plastic waste and scrap mismanaged (open dump, open pits, unsanitary landfills)']:
            resin_indicators.append([indicator + ' - ' + i for i in resin_compo['Europe'].index])
        resin_indicators = [item for sublist in resin_indicators for item in sublist]

        # plastic resin composition for product categories of plastic physical flow accounts
        with_resins_for_Canada = pd.DataFrame(0, index=resin_indicators, columns=self.F.columns, dtype=object)
        for plastic_cat in map_plastic_data_to_io:
            for indicator in kpis:
                df = pd.concat(
                    [self.F.loc[indicator, [i for i in self.F.columns if i[1] in map_plastic_data_to_io[plastic_cat]]]] *
                    len(resin_compo['Canada'].index), axis=1)
                df *= resin_compo['Canada'].loc[:, plastic_cat].values
                df.columns = [indicator + ' - ' + i for i in resin_compo['Canada'].index]
                with_resins_for_Canada.loc[df.columns, df.index] = df.T

        # plastic resin composition overall for OECD data. Different compositions depending on regions.
        with_resins_for_international = pd.DataFrame(0, index=resin_indicators, columns=self.F_exio.columns, dtype=object)
        eu_membership = list(zip(self.F_exio.columns.levels[0], coco.convert(self.F_exio.columns.levels[0], to='EU28')))
        for indicator in kpis + [
            'Disposed plastic waste and scrap mismanaged (open dump, open pits, unsanitary landfills)']:
            # For Europe
            df = pd.concat([self.F_exio.loc[indicator, [i[0] for i in eu_membership if i[1] == 'EU28']]] *
                           len(resin_compo['Europe'].index), axis=1)
            df *= resin_compo['Europe'].loc[:, 'Composition'].values
            df.columns = [indicator + ' - ' + i for i in resin_compo['Europe'].index]
            with_resins_for_international.loc[df.columns, df.index] = df.T
            # For the US
            df = pd.concat([self.F_exio.loc[indicator, 'US']] * len(resin_compo['US'].index), axis=1)
            df *= resin_compo['US'].loc[:, 'Composition'].values
            df.columns = [indicator + ' - ' + i for i in resin_compo['US'].index]
            df = pd.concat([df], keys=['US'])
            with_resins_for_international.loc[df.columns, df.index] = df.T
            # For China
            df = pd.concat([self.F_exio.loc[indicator, 'CN']] * len(resin_compo['China'].index), axis=1)
            df *= resin_compo['China'].loc[:, 'Composition'].values
            df.columns = [indicator + ' - ' + i for i in resin_compo['China'].index]
            df = pd.concat([df], keys=['CN'])
            with_resins_for_international.loc[df.columns, df.index] = df.T
            # For India
            df = pd.concat([self.F_exio.loc[indicator, 'IN']] * len(resin_compo['India'].index), axis=1)
            df *= resin_compo['India'].loc[:, 'Composition'].values
            df.columns = [indicator + ' - ' + i for i in resin_compo['India'].index]
            df = pd.concat([df], keys=['IN'])
            with_resins_for_international.loc[df.columns, df.index] = df.T
            # For other countries
            df = pd.concat([self.F_exio.loc[indicator, [i[0] for i in eu_membership if
                                                   (i[1] != 'EU28' and i[0] not in ['US', 'CN', 'IN'])]]] *
                           len(resin_compo['Global'].index), axis=1)
            df *= resin_compo['Global'].loc[:, 'Composition'].values
            df.columns = [indicator + ' - ' + i for i in resin_compo['Global'].index]
            with_resins_for_international.loc[df.columns, df.index] = df.T

        # replace plastic waste information without resin, by information with resin
        self.F = pd.concat([self.F, with_resins_for_Canada]).drop(
            kpis + ['Disposed plastic waste and scrap mismanaged (open dump, open pits, unsanitary landfills)'])
        self.F_exio = pd.concat([self.F_exio, with_resins_for_international]).drop(
            kpis + ['Disposed plastic waste and scrap mismanaged (open dump, open pits, unsanitary landfills)'])

        # ----------------------------- Adapting the rest of the system -------------------------------
        self.FY = self.FY.reindex(self.F.index).fillna(0)

        # introduce plastic emissions to S matrix of exiobase
        self.S_exio = pd.concat([self.S_exio, self.F_exio.loc[resin_indicators] *
                                 (1/self.x_exio).replace([np.inf, -np.inf], 0).loc[:, 'indout']])
        # millions of euros (exiobase) to euros
        self.S_exio.loc[resin_indicators] /= 1000000

        # add metadata for plastic emissions
        self.emission_metadata = pd.concat([self.emission_metadata,
                                            pd.DataFrame('t', columns=['Unit'], index=resin_indicators)])

        # add plastic emissions to characterization matrix
        self.C = pd.concat([self.C, pd.DataFrame(np.eye(len(resin_indicators)),
                                                 index=pd.MultiIndex.from_product([resin_indicators, ['tonnes of plastics']]),
                                                 columns=resin_indicators)]).fillna(0)
        self.C_exio = pd.concat([self.C_exio, pd.DataFrame(np.eye(len(resin_indicators)),
                                                 index=pd.MultiIndex.from_product([resin_indicators, ['tonnes of plastics']]),
                                                 columns=resin_indicators)]).fillna(0)

        plastic_cfs = pd.DataFrame(0, index=pd.MultiIndex.from_tuples([('Physical effects on biota', 'PDF.m2.yr')]),
                                   columns=[i for i in self.C.columns if 'leak' in i], dtype=float)
        # note all emitted plastic particles will end up in the ocean!
        # source = Macroplastic: Leakage from waste export - https://www.plasticfootprint.earth/assessment-methodology/
        # from released rates to ocean and land, take medium size plastic particle x medium residual value
        ratio_ocean = 0.15
        # assume a microfibers shape (most common) and a 10um size (most common), also convert kg to tonnes
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Polyethylene terephthalate (PET) resins'] = 0.0117891847201787 * ratio_ocean * 1000
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Other thermoplastic polyester resins'] = 0.0117891847201787 * ratio_ocean * 1000
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Low-density polyethylene (LDPE) resins'] = 14.2015621983309 * ratio_ocean * 1000
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Linear low-density polyethylene (LLDPE) resins'] = 14.2015621983309 * ratio_ocean * 1000
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - High-density polyethylene (HDPE) resins'] = 14.5208111050616 * ratio_ocean * 1000
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Other polyethylene resins'] = (14.5208111050616 + 14.2015621983309) / 2 * ratio_ocean * 1000
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Polystyrene (PS) resins'] = 14.7359173605865 * ratio_ocean * 1000
        # proxy = PS, based on density
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Acrylonitrile-butadiene-styrene (ABS) resins'] = 14.7359173605865 * ratio_ocean * 1000
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Polyvinyl chloride (PVC) resins'] = 0.011600652123176* ratio_ocean * 1000
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Polypropylene (PP) resins'] = 13.5454167946045 * ratio_ocean * 1000
        # proxy = PLA, based on density
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Thermoplastic polyurethane (PU) resins'] = 0.0118266839558158 * ratio_ocean * 1000
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Polyamide (nylon) resins'] = 0.0118081147968369 * ratio_ocean * 1000
        # proxy = PLA
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - All other thermoplastic resins'] = 0.0118266839558158 * ratio_ocean * 1000
        # proxy = PLA, based on density
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Phenolic (PF) resins'] = 0.0118266839558158 * ratio_ocean * 1000
        # proxy = PLA, based on density
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Urea formaldehyde (UF) resins'] = 0.0118266839558158 * ratio_ocean * 1000
        # copy urea formaldehyde
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - All other formaldehyde based resins'] = 0.0118266839558158 * ratio_ocean * 1000
        # proxy = PLA, based on density
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Unsaturated polyester (thermosetting) resins'] = 0.0118266839558158 * ratio_ocean * 1000
        # proxy = PLA, based on density
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Thermosetting polyurethane resins'] = 0.0118266839558158 * ratio_ocean * 1000
        # copy polyurethane/polyester thermoset
        plastic_cfs.loc[('Physical effects on biota', 'PDF.m2.yr'),
                        'Plastic leaked permanently into the environment - Other thermosetting resins'] = 0.0118266839558158 * ratio_ocean * 1000

        self.C = pd.concat([self.C.drop([i for i in self.C.index if 'leak' in i[0]]), plastic_cfs]).fillna(0)
        self.C_exio = pd.concat(
            [self.C_exio.drop([i for i in self.C_exio.index if 'leak' in i[0]]), plastic_cfs]).fillna(0)

    def differentiate_biogenic_carbon_emissions(self):
        """
        The physical flow GHG accounts from StatCan do not differentiate between CO2 fossil and biogenic. We thus use
        exiobase biogenic vs fossil CO2 distribution per sector to determine the amount of CO2 biogenic in StatCan
        data.
        :return:
        """

        # loading concordances between exiobase classification and IOIC
        if self.year in [2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021]:
            ioic_exio = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2014-2021/IOIC_EXIOBASE.xlsx'),
                                      'commodities')
        elif self.year in [2022]:
            ioic_exio = pd.read_excel(pkg_resources.resource_stream(__name__, '/Data/Concordances/2022/IOIC_EXIOBASE.xlsx'),
                                      'commodities')
        ioic_exio = ioic_exio[2:].drop('IOIC Detail level - EXIOBASE', axis=1).set_index('Unnamed: 1').fillna(0)
        ioic_exio.index.name = None
        ioic_exio.index = [{j[0]: j[1] for j in self.commodities}[i] for i in ioic_exio.index]
        ioic_exio /= ioic_exio.sum()
        ioic_exio = ioic_exio.fillna(0)

        # identify biogenic and fossil CO2 emissions in Exiobase
        CO2_fossil = [i for i in self.F_exio.index if 'CO2' in i and 'bio' not in i]
        CO2_bio = [i for i in self.F_exio.index if 'CO2' in i and 'bio' in i]
        CO2 = [i for i in self.F_exio.index if 'CO2' in i]

        # apply the distribution of biogenic CO2 from Exiobase to openIO sectors
        total = self.F_exio.loc[CO2, 'CA'].dot(ioic_exio.T).sum()[
            self.F_exio.loc[CO2, 'CA'].dot(ioic_exio.T).sum() != 0]
        bio = self.F_exio.loc[CO2_bio, 'CA'].dot(ioic_exio.T).sum().loc[total.index] / total
        bio = bio.reindex([i[1] for i in self.commodities]).fillna(0)
        bio = pd.DataFrame(pd.concat([bio] * len([i for i in self.S.columns.levels[0] if 'CA-' in i])), columns=[
            'Carbon dioxide - biogenic'])
        bio.index = [i for i in self.S.columns if 'CA-' in i[0]]
        bio_openio = self.S.loc[[i for i in self.S.index if 'Carbon dioxide' == i[1]],
                                [i for i in self.S.columns if 'CA-' in i[0]]].copy()
        bio_openio = np.multiply(bio_openio, bio.iloc[:, 0])
        bio_openio.index = [(i[0], 'Carbon dioxide - biogenic', i[2]) for i in bio_openio.index]

        # apply the distribution of fossil CO2 from Exiobase to openIO sectors
        total = self.F_exio.loc[CO2, 'CA'].dot(ioic_exio.T).sum()[
            self.F_exio.loc[CO2, 'CA'].dot(ioic_exio.T).sum() != 0]
        fossil = self.F_exio.loc[CO2_fossil, 'CA'].dot(ioic_exio.T).sum().loc[total.index] / total
        fossil = fossil.reindex([i[1] for i in self.commodities]).fillna(0)
        fossil = pd.DataFrame(pd.concat([fossil] * len([i for i in self.S.columns.levels[0] if 'CA-' in i])), columns=[
            'Carbon dioxide - fossil'])
        fossil.index = [i for i in self.S.columns if 'CA-' in i[0]]
        fossil_openio = self.S.loc[[i for i in self.S.index if 'Carbon dioxide' == i[1]],
                                   [i for i in self.S.columns if 'CA-' in i[0]]].copy()
        fossil_openio = np.multiply(fossil_openio, fossil.iloc[:, 0])
        fossil_openio.index = [(i[0], 'Carbon dioxide - fossil', i[2]) for i in fossil_openio.index]

        # drop total CO2 emissions
        self.S.drop([i for i in self.S.index if 'Carbon dioxide' == i[1]], inplace=True)
        # add fossil and biogenic CO2 emissions
        self.S = pd.concat([self.S, fossil_openio.reindex(self.S.columns, axis=1).fillna(0),
                            bio_openio.reindex(self.S.columns, axis=1).fillna(0)])

        # same story for self.F
        bio_openio_scaled = self.F.loc[[i for i in self.F.index if 'Carbon dioxide' == i[1]],
                                       [i for i in self.F.columns if 'CA-' in i[0]]].copy()
        bio_openio_scaled = np.multiply(bio_openio_scaled, bio.iloc[:, 0])
        bio_openio_scaled.index = [(i[0], 'Carbon dioxide - biogenic', i[2]) for i in bio_openio_scaled.index]
        bio_openio_scaled = bio_openio_scaled.fillna(0)
        fossil_openio_scaled = self.F.loc[[i for i in self.F.index if 'Carbon dioxide' == i[1]],
                                          [i for i in self.F.columns if 'CA-' in i[0]]].copy()
        fossil_openio_scaled = np.multiply(fossil_openio_scaled, fossil.iloc[:, 0])
        fossil_openio_scaled.index = [(i[0], 'Carbon dioxide - fossil', i[2]) for i in fossil_openio_scaled.index]
        fossil_openio_scaled = fossil_openio_scaled.fillna(0)

        self.F.drop([i for i in self.F.index if 'Carbon dioxide' == i[1]], inplace=True)
        self.F = pd.concat([self.F, fossil_openio_scaled.reindex(self.F.columns, axis=1).fillna(0),
                            bio_openio_scaled.reindex(self.F.columns, axis=1).fillna(0)])

        # same story for self.FY
        bio_final_demand = (self.FY_exio.loc[CO2_bio, 'CA'].sum() / self.FY_exio.loc[CO2, 'CA'].sum()).fillna(0)
        fossil_final_demand = (self.FY_exio.loc[CO2_fossil, 'CA'].sum() / self.FY_exio.loc[CO2, 'CA'].sum()).fillna(0)

        bio_openio_final_demand = self.FY.loc[[i for i in self.FY.index if 'Carbon dioxide' == i[1]]].copy()
        bio_openio_final_demand = np.multiply(
            bio_openio_final_demand.loc(axis=1)[:, 'Household final consumption expenditure'],
            bio_final_demand.loc['Final consumption expenditure by households'])
        bio_openio_final_demand.index = [(i[0], 'Carbon dioxide - biogenic', i[2]) for i in
                                         bio_openio_final_demand.index]
        bio_openio_final_demand = bio_openio_final_demand.fillna(0)

        fossil_openio_final_demand = self.FY.loc[[i for i in self.FY.index if 'Carbon dioxide' == i[1]]].copy()
        fossil_openio_final_demand = np.multiply(
            fossil_openio_final_demand.loc(axis=1)[:, 'Household final consumption expenditure'],
            fossil_final_demand.loc['Final consumption expenditure by households'])
        fossil_openio_final_demand.index = [(i[0], 'Carbon dioxide - fossil', i[2]) for i in
                                            fossil_openio_final_demand.index]
        fossil_openio_final_demand = fossil_openio_final_demand.fillna(0)

        self.FY.drop([i for i in self.FY.index if 'Carbon dioxide' == i[1]], inplace=True)
        self.FY = pd.concat([self.FY, fossil_openio_final_demand.reindex(self.FY.columns, axis=1).fillna(0),
                             bio_openio_final_demand.reindex(self.FY.columns, axis=1).fillna(0)])

        # add "fossil" to the elementary flow name in characterization matrix
        self.C.columns = [(i[0], 'Carbon dioxide - fossil', i[2]) if i[1] == 'Carbon dioxide' else i for i in
                          self.C.columns]

        # also add an entry for biogenic carbon in characterization matrix
        provinces = [i for i in self.A.columns.levels[0] if 'CA-' in i]
        for province in provinces:
            self.C.loc[:, [(province, 'Carbon dioxide - biogenic', 'Air')]] = 0

        # reindex stuff around
        self.F = self.F.reindex(self.C.columns).fillna(0)
        self.F = self.F.reindex(self.A.index, axis=1).fillna(0)
        self.FY = self.FY.reindex(self.F.index).fillna(0)
        self.FY = self.FY.reindex(self.Y.columns, axis=1).fillna(0)

    def calc(self):
        """
        Method to calculate the Leontief inverse and get total impacts
        :return: self.L (total requirements), self.E (total emissions), self.D (total impacts)
        """
        I = pd.DataFrame(np.eye(len(self.A)), self.A.index, self.A.columns, dtype=float)

        # in case there are some unwanted NaNs along the way
        self.A = self.A.fillna(0)

        if self.endogenizing:
            self.L = pd.DataFrame(np.linalg.solve(I - (self.A + self.K), I), self.A.index, I.columns)
        else:
            self.L = pd.DataFrame(np.linalg.solve(I - self.A, I), self.A.index, I.columns)

        self.E = self.S.dot(self.L).dot(self.Y) + self.FY

        self.D = self.C.dot(self.E)

    def get_emission_factors_basic_price(self):
        """
        Function returns the emission factors at basic price, i.e., purchaser price without trade margins/downstream
        transportation margins/taxes paid by purchaser
        :return: self.emission_factors_basic_price
        """

        if self.L.empty:
            self.calc()

        self.emission_factors_basic_price = self.C.dot(self.S).dot(self.L)

        return self.emission_factors_basic_price

    def get_emission_factors_purchaser_price(self):
        """
        Function returns the emission factors at purchaser price, i.e., the price directly paid by the purchaser
        (including trade margins, downstream transportation margins and taxes)
        :return: self.emission_factors_purchaser_price
        """

        if self.L.empty:
            self.calc()

        # if user didn't get emission factors with basic price yet
        if self.emission_factors_basic_price.empty:
            self.emission_factors_basic_price = self.C.dot(self.S).dot(self.L)

        # quick function to format the trade, transportation and taxes tables
        def formatting(table):
            if self.year in [2014, 2015, 2016, 2017]:
                starting_line = 11
                starting_line_values = 16

            elif self.year in [2018, 2019, 2020, 2021, 2022]:
                starting_line = 3
                starting_line_values = 7

            df = table.iloc[starting_line_values - 2:, 2:]
            df.index = list(zip(table.iloc[starting_line_values - 2:, 0].tolist(),
                                table.iloc[starting_line_values - 2:, 1].tolist()))
            df.columns = list(zip(table.iloc[starting_line + 1, 2:].tolist(),
                                  table.iloc[starting_line, 2:].tolist()))

            # some file versions register zeros as '.' instead of '0'
            with pd.option_context('future.no_silent_downcasting', True):
                df = df.replace('.', 0)

            return df

        # go through the files
        files = [i for i in os.walk(self.folder_path)]
        files = [i for i in files[0][2] if i[:2] in self.matching_dict.keys() and 'SUT' in i]
        files.sort()

        # iterate through each province files
        for province_data in files:
            sut = pd.read_excel(os.path.join(self.folder_path, province_data), None)
            province = province_data[:2]

            with pd.option_context('future.no_silent_downcasting', True):
                total = formatting(sut['Use_Purchaser']).loc[:, [('TOTUSE', 'Total use')]][
                    formatting(sut['Use_Purchaser']).loc[:, [('TOTUSE', 'Total use')]] != 0].fillna(0)
            # calculate share of trade, transport and tax
            trade = formatting(sut['Use_Trade']).loc[total.index, [('TOTUSE', 'Total use')]]
            # drop negative values (from retail and wholesale)
            with pd.option_context('future.no_silent_downcasting', True):
                trade = trade[trade >= 0].fillna(0)
            trade = trade.where(total != 0, 0) / total.where(total != 0, 1)
            transport = formatting(sut['Use_Transport']).loc[total.index, [('TOTUSE', 'Total use')]]
            # drop negative values (from distribution margins)
            with pd.option_context('future.no_silent_downcasting', True):
                transport = transport[transport >= 0].fillna(0)
            transport = transport.where(total != 0, 0) / total.where(total != 0, 1)
            tax = formatting(sut['Use_Tax']).loc[total.index, [('TOTUSE', 'Total use')]]
            # drop negative values (from value added)
            with pd.option_context('future.no_silent_downcasting', True):
                tax = tax[tax >= 0].fillna(0)
            tax = tax.where(total != 0, 0) / total.where(total != 0, 1)
            basic = formatting(sut['Use_Basic']).loc[:, [('TOTUSE', 'Total use')]]

            # format
            trade.index = pd.MultiIndex.from_tuples(trade.index)
            transport.index = pd.MultiIndex.from_tuples(transport.index)
            tax.index = pd.MultiIndex.from_tuples(tax.index)
            basic.index = pd.MultiIndex.from_tuples(basic.index)

            # drop useless index level
            trade = trade.droplevel(0)
            transport = transport.droplevel(0)
            tax = tax.droplevel(0)
            basic = basic.droplevel(0)

            # calculate the average impact of trade
            average_trade_impact = (
                    self.emission_factors_basic_price.loc(axis=1)[
                        'CA-' + province, [i for i in self.emission_factors_basic_price.columns.levels[1] if
                                           'Retail margins' in i or 'Wholesale margins' in i]].droplevel(0, axis=1) *
                    (basic.loc[[i for i in self.emission_factors_basic_price.columns.levels[1] if
                                'Retail margins' in i or 'Wholesale margins' in i]] /
                     basic.loc[[i for i in self.emission_factors_basic_price.columns.levels[1] if
                                'Retail margins' in i or 'Wholesale margins' in i]].sum()).iloc[:, 0]
            ).sum(axis=1)
            # calculate the average impact of downstream transport
            transport_sectors_2017 = ['Air freight transportation services', 'Rail freight transportation services',
                                 'Water freight transportation services',
                                 'Road transportation services for general freight',
                                 'Road transportation services for specialized freight',
                                 'Transportation of natural gas by pipeline',
                                 'Transportation of crude oil and other commodities by pipeline',
                                 'Water transportation support, maintenance and repair services',
                                 'Road transportation support services',
                                 'Freight transportation arrangement and customs brokering services',
                                 'Other transportation support services', 'Grain storage', 'Natural gas distribution']
            transport_sectors_2014 = ['Air freight transportation services', 'Rail freight transportation services',
                                 'Water freight transportation services',
                                 'Truck transportation services for general freight',
                                 'Truck transportation services for specialized freight',
                                 'Transportation of natural gas by pipeline',
                                 'Transportation of crude oil and other commodities by pipeline',
                                 'Water transportation support, maintenance and repair services',
                                 'Road transportation support services',
                                 'Freight transportation arrangement and customs brokering services',
                                 'Other transportation support services', 'Grain storage', 'Natural gas distribution']
            if self.year > 2017:
                average_transport_impact = (
                        self.emission_factors_basic_price.loc(axis=1)['CA-' + province,
                                                                      transport_sectors_2017].droplevel(0, axis=1) *
                        (basic.loc[transport_sectors_2017] / basic.loc[transport_sectors_2017].sum()).iloc[:, 0]).sum(axis=1)
            elif self.year in [2014, 2015, 2016]:
                average_transport_impact = (
                        self.emission_factors_basic_price.loc(axis=1)['CA-' + province,
                                                                      transport_sectors_2014].droplevel(0, axis=1) *
                        (basic.loc[transport_sectors_2014] / basic.loc[transport_sectors_2014].sum()).iloc[:, 0]).sum(axis=1)

            # reformat tables
            trade_reformatted = pd.concat(
                [trade.reindex(self.emission_factors_basic_price.loc[:, 'CA-' + province].columns).iloc[:, 0]] * len(
                    average_trade_impact), axis=1).T
            trade_reformatted.index = average_trade_impact.index

            average_trade_impact_reformatted = pd.concat([average_trade_impact] * len(trade_reformatted.columns), axis=1)
            average_trade_impact_reformatted.columns = trade_reformatted.columns

            transport_reformatted = pd.concat(
                [transport.reindex(self.emission_factors_basic_price.loc[:, 'CA-' + province].columns).iloc[:, 0]] * len(
                    average_transport_impact), axis=1).T
            transport_reformatted.index = average_transport_impact.index

            average_transport_impact_reformatted = pd.concat(
                [average_transport_impact] * len(transport_reformatted.columns), axis=1)
            average_transport_impact_reformatted.columns = transport_reformatted.columns

            tax_reformatted = pd.concat(
                [tax.reindex(self.emission_factors_basic_price.loc[:, 'CA-' + province].columns).iloc[:, 0]] * len(
                    average_transport_impact), axis=1).T
            tax_reformatted.index = average_transport_impact.index

            # calculate the purchaser price emission factors
            df = ((self.emission_factors_basic_price.loc[:, 'CA-' + province] +
                   trade_reformatted * average_trade_impact_reformatted +
                   transport_reformatted * average_transport_impact_reformatted) /
                  (1 + trade_reformatted + transport_reformatted + tax_reformatted))

            # concat with master df to egt all provinces in same dataframe
            self.emission_factors_purchaser_price = pd.concat(
                [self.emission_factors_purchaser_price, pd.concat([df], keys=['CA-' + province], axis=1)], axis=1)

        return self.emission_factors_purchaser_price

# -------------------------------------------------- SUPPORT ----------------------------------------------------------

    def balance_flows(self, concordance):
        """
        Some flows from the NPRI trigger some double counting if left unattended. This method deals with these flows
        :return: balanced self.F
        """

        # we want to use handy multi-index features so we remove flows without multi-index and plug them back at the end
        F_multiindex = self.F.loc[[i for i in self.F.index if type(i) == tuple]].copy()
        F_multiindex.index = pd.MultiIndex.from_tuples(F_multiindex.index)

        # VOCs
        rest_of_voc = [i for i in concordance.index if 'Speciated VOC' in i and concordance.loc[i].isna().iloc[0]]
        df = F_multiindex.loc[[i for i in F_multiindex.index if i[1] in rest_of_voc]]
        if 'Volatile Organic Compounds (VOCs)' in F_multiindex.index.levels[1]:
            F_multiindex.loc[:, 'Volatile Organic Compounds (VOCs)', 'Air'] += df.groupby(level=0).sum().values
        # 2022 NPRI file has no VOCs entry
        else:
            dff = df.groupby(level=0).sum()
            dff = pd.concat([dff], keys=['Volatile Organic Compounds (VOCs)'])
            dff = pd.concat([dff], keys=['Air'])
            dff = dff.reorder_levels([2, 1, 0])
            F_multiindex = pd.concat([F_multiindex, dff])

            # because VOCs are not in the NPRI, the CF is not created in C matrix. Need to manually create it here...
            foo = pd.DataFrame(0, index=pd.MultiIndex.from_product(
                [list(self.matching_dict.keys()), ['Volatile Organic Compounds (VOCs)'], ['Air']]),
                              columns=[('Photochemical ozone formation', 'kgNOxeq'),
                                       ('Photochemical ozone formation, human health', 'DALY'),
                                       ('Photochemical ozone formation, ecosystem quality', 'PDF.m2.yr')])
            foo.index = foo.index.values

            foo.loc[:, [('Photochemical ozone formation', 'kgNOxeq')]] = 0.18
            foo.loc[:, [('Photochemical ozone formation, human health', 'DALY')]] = 0.0000001638
            foo.loc[:, [('Photochemical ozone formation, ecosystem quality', 'PDF.m2.yr')]] = 2.5277027027027

            foo.columns = pd.MultiIndex.from_tuples(foo.columns)

            self.C = self.C.join(foo.T.reindex(self.C.index).fillna(0)).fillna(0)
            self.FY = self.FY.reindex(self.C.columns).fillna(0)

        F_multiindex = F_multiindex.drop([i for i in rest_of_voc if i in F_multiindex.index.levels[1]], axis=0, level=1)

        # adjust characterization matrix too
        self.C = self.C.drop([i for i in self.C.columns if i[1] in rest_of_voc], axis=1)

        if self.NPRI_file_year >= 2018:
            # PMs, only take highest value flow as suggested by the NPRI team:
            # [https://www.canada.ca/en/environment-climate-change/services/national-pollutant-release-inventory/using-interpreting-data.html]
            for sector in F_multiindex.columns:
                little_pm = F_multiindex.loc[
                    (sector[0], 'PM2.5 - Particulate Matter <= 2.5 Micrometers', 'Air'), sector]
                big_pm = F_multiindex.loc[(sector[0], 'PM10 - Particulate Matter <= 10 Micrometers', 'Air'), sector]
                unknown_size = F_multiindex.loc[(sector[0], 'Total particulate matter', 'Air'), sector]
                if little_pm >= big_pm:
                    if little_pm >= unknown_size:
                        F_multiindex.loc[(sector[0], 'PM10 - Particulate Matter <= 10 Micrometers', 'Air'), sector] = 0
                        F_multiindex.loc[(sector[0], 'Total particulate matter', 'Air'), sector] = 0
                    else:
                        F_multiindex.loc[(sector[0], 'PM10 - Particulate Matter <= 10 Micrometers', 'Air'), sector] = 0
                        F_multiindex.loc[
                            (sector[0], 'PM2.5 - Particulate Matter <= 2.5 Micrometers', 'Air'), sector] = 0
                else:
                    if big_pm > unknown_size:
                        F_multiindex.loc[
                            (sector[0], 'PM2.5 - Particulate Matter <= 2.5 Micrometers', 'Air'), sector] = 0
                        F_multiindex.loc[(sector[0], 'Total particulate matter', 'Air'), sector] = 0
                    else:
                        F_multiindex.loc[(sector[0], 'PM10 - Particulate Matter <= 10 Micrometers', 'Air'), sector] = 0
                        F_multiindex.loc[
                            (sector[0], 'PM2.5 - Particulate Matter <= 2.5 Micrometers', 'Air'), sector] = 0
        else:
            # PMs, only take highest value flow as suggested by the NPRI team:
            # [https://www.canada.ca/en/environment-climate-change/services/national-pollutant-release-inventory/using-interpreting-data.html]
            for sector in F_multiindex.columns:
                little_pm = F_multiindex.loc[(sector[0], 'PM2.5', 'Air'), sector]
                big_pm = F_multiindex.loc[(sector[0], 'PM10', 'Air'), sector]
                unknown_size = F_multiindex.loc[(sector[0], 'Total particulate matter', 'Air'), sector]
                if little_pm >= big_pm:
                    if little_pm >= unknown_size:
                        F_multiindex.loc[(sector[0], 'PM10', 'Air'), sector] = 0
                        F_multiindex.loc[(sector[0], 'Total particulate matter', 'Air'), sector] = 0
                    else:
                        F_multiindex.loc[(sector[0], 'PM10', 'Air'), sector] = 0
                        F_multiindex.loc[(sector[0], 'PM2.5', 'Air'), sector] = 0
                else:
                    if big_pm > unknown_size:
                        F_multiindex.loc[(sector[0], 'PM2.5', 'Air'), sector] = 0
                        F_multiindex.loc[(sector[0], 'Total particulate matter', 'Air'), sector] = 0
                    else:
                        F_multiindex.loc[(sector[0], 'PM10', 'Air'), sector] = 0
                        F_multiindex.loc[(sector[0], 'PM2.5', 'Air'), sector] = 0

        # plug back the non multi-index flows
        self.F = pd.concat([F_multiindex, self.F.loc[[i for i in self.F.index if type(i) != tuple]].copy()])

    def split_private_public_sectors(self, NAICS_code, IOIC_code):
        """
        Support method to split equally emissions from private and public sectors
        :param NAICS_code: [string or list] the NAICS code(s) whose emissions will be split
        :param IOIC_code: [string] the IOIC_code inhereting the split emissions (will be private or public sector)
        :return: updated self.F
        """
        try:
            df = self.F.loc(axis=1)[:, NAICS_code].copy()
            df.columns = pd.MultiIndex.from_product([self.matching_dict, [IOIC_code]])
            self.F = pd.concat([self.F, df / 2], axis=1)
            self.F.loc(axis=1)[:, NAICS_code] /= 2
        except KeyError:
            pass

    def assert_order(self, exporting_province, importing_province, scaled_imports_U, scaled_imports_Y,
                     scaled_imports_K=None):
        """
        Support method that checks that the order of index is the same so that .values can be used in the interprovincial
        balance function.
        :return:
        """
        assert all(self.U.loc[exporting_province, importing_province].index ==
                   scaled_imports_U.loc[:, self.U.columns.levels[1]].reindex(
                       self.U.loc[exporting_province, importing_province].columns, axis=1).index)
        assert all(self.U.loc[exporting_province, importing_province].columns ==
                   scaled_imports_U.loc[:, self.U.columns.levels[1]].reindex(
                       self.U.loc[exporting_province, importing_province].columns, axis=1).columns)

        if isinstance(scaled_imports_K, pd.DataFrame):
            assert all(self.K.loc[exporting_province, importing_province].index ==
                       scaled_imports_K.loc[:, self.K.columns.levels[1]].reindex(
                           self.K.loc[exporting_province, importing_province].columns, axis=1).index)
            assert all(self.K.loc[exporting_province, importing_province].columns ==
                       scaled_imports_K.loc[:, self.K.columns.levels[1]].reindex(
                           self.K.loc[exporting_province, importing_province].columns, axis=1).columns)

        assert all(self.Y.loc[exporting_province, importing_province].index ==
                   scaled_imports_Y.loc[:, self.Y.columns.levels[1]].reindex(
                       self.Y.loc[exporting_province, importing_province].columns, axis=1).index)
        assert all(self.Y.loc[exporting_province, importing_province].columns ==
                   scaled_imports_Y.loc[:, self.Y.columns.levels[1]].reindex(
                       self.Y.loc[exporting_province, importing_province].columns, axis=1).columns)

    def export(self, filepath='', format=''):
        """
        Function to export in the chosen format.
        :param filepath: the path where to store the export file
        :param format: available formats 'csv', 'excel', 'pickle', 'json'
        :return: nothing
        """

        if not filepath:
            print("Please provide a filepath")
            return
        if not format:
            print("Please enter a format")
            return

        def flat_multiindex(df):
            df.index = df.index.tolist()
            df.columns = df.columns.tolist()
        flat_multiindex(self.A)
        flat_multiindex(self.Y)
        flat_multiindex(self.R)
        flat_multiindex(self.S)
        flat_multiindex(self.FY)
        flat_multiindex(self.C)

        def remove_zeros(df):
            return df.replace({0: np.nan})
        self.A = remove_zeros(self.A)
        self.Y = remove_zeros(self.Y)
        self.R = remove_zeros(self.R)
        self.S = remove_zeros(self.S)
        self.FY = remove_zeros(self.FY)
        self.C = remove_zeros(self.C)

        if format == 'excel':
            writer = pd.ExcelWriter(filepath, engine='xlsxwriter')

            self.A.to_excel(writer, 'A')
            self.Y.to_excel(writer, 'Y')
            self.R.to_excel(writer, 'R')
            self.S.to_excel(writer, 'S')
            self.FY.to_excel(writer, 'FY')
            self.C.to_excel(writer, 'C')

            writer.save()

        else:
            print('Format requested not implemented yet.')


def treatment_import_data(original_file_path):
    """Function used to treat the merchandise imports trade database file. FIle is way too big to be provided to
    users through Github, so we treat the data to only keep what is relevant."""

    # load database
    merchandise_database = pd.read_csv(original_file_path)
    # drop useless columns

    merchandise_database = merchandise_database.drop(['YearMonth/AnnéeMois', 'Province', 'State/État',
                                                      'Quantity/Quantité', 'Unit of Measure/Unité de Mesure'],
                                                     axis=1)

    # drop international imports coming from Canada
    merchandise_database = merchandise_database[merchandise_database['Country/Pays'] != 'CA']

    # also drop nan countries for obvious reasons
    merchandise_database = merchandise_database.dropna(subset=['Country/Pays'])

    # set the index as country/code multi-index
    merchandise_database = merchandise_database.set_index(['Country/Pays', 'HS6'])

    # regroup data from several months into a single yearly data
    merchandise_database = merchandise_database.groupby(merchandise_database.index).sum()

    # multi-index is cleaner
    merchandise_database.index = pd.MultiIndex.from_tuples(merchandise_database.index)

    return merchandise_database