<jupyter_start><jupyter_text>$A_\mathrm{ITC-cxc} = U \cdot \hat{g}^{-1} \cdot V' \cdot \hat{q}^{-1}$<jupyter_code>U.dot(np.diag(1/g)).dot((V.T).dot(np.diag(1/q)))
SUT = pysut.SupplyUseTable(V=V, U=U)
SUT.build_mr_Gamma()
SUT.build_E_bar()
A, __, __, __, Z, __ = SUT.aac_agg(keep_size=False)
print(Z)
SUT.E_bar
SUT.Gamma
V_A = np.array([[1,0,0],[0.6,0,0],[0,0,1]])
V_B = np.array([[1,0,0],[0,1,0],[0,0,1]])
U_A = np.array([[0,0,0.5],[0,0,0.5],[0.5,0,0]])
U_B = np.array([[0,0.3,0],[1,0,0],[0.5,0,0]])

V_AB = np.array([[0,0,0],[0,0,0],[0,0,0]])
V_BA = np.array([[0,0,0],[0,0,0],[0,0,0]])
U_AB = np.array([[0,0,0],[0,0,0],[0.5,0,0]])
U_BA = np.array([[0.1,0,0],[0,0,0],[0,0,0]])
V_total = np.concatenate((np.concatenate((V_A,V_AB), axis=1),np.concatenate((V_BA,V_B),axis=1)),axis=0)
U_total = np.concatenate((np.concatenate((U_A,U_AB), axis=1),np.concatenate((U_BA,U_B),axis=1)),axis=0)
V_total
U_total
SUT_MR = pysut.SupplyUseTable(V = V_total, U=U_total, regions=2)
SUT_MR.build_mr_Gamma()
SUT_MR.build_E_bar()
SUT_MR.build_mr_Xi()
#A_MR,_,_,_,_,_,Z_MR,_ = SUT_MR.psc_agg()
print(SUT_MR.Xi)
SUT_MR.aac_agg(keep_size=False)
print("V")
print(V_total)
print("U")
print(U_total)
SUT_MR.Gamma
SUT_MR.E_bar
SupplyDiag_Eval = SUT_MR.supply_diag_check()
print(SupplyDiag_Eval)

#A_MR,_,_,_,_,_,Z_MR,_ = 
SUT_MR.psc_agg()
print("V")
print(V_total)
print("U")
print(U_total)

SUT_MR.btc()
np.dot(SUT_MR.Xi,V_total-np.diag(np.diag(V_total)))
np.dot(SUT_MR.Xi,U_total)
print("V")
print(V_total)
print("U")
print(U_total)


np.dot(SUT_MR.Xi,U_total)-np.dot(SUT_MR.Xi,V_total-np.diag(np.diag(V_total)))
SUT_MR.psc_agg()[0]-np.dot(SUT_MR.Xi,U_total)-np.dot(SUT_MR.Xi,V_total-np.diag(np.diag(V_total)))
exio_v = pd.read_csv('/home/jakobs/Documents/IndEcol/BONSAI/data/exiobase_hsut_3317/MR_HSUP_2011_v3_3_17_bonsai.csv', header=[0,1,2,3], index_col=[0,1,2,3,4])
exio_u = pd.read_csv('/home/jakobs/Documents/IndEcol/BONSAI/data/exiobase_hsut_3317/MR_HUSE_2011_v3_3_17_bonsai.csv', header=[0,1,2,3], index_col=[0,1,2,3,4])
exio_u.values.shape
exio_v.values.shape
exio_u.loc[exio_u.index.get_level_values(2) == 'p40.2.b']
np.where(exio_u.loc[exio_u.index.get_level_values(2) == 'p40.2.b'].values >0)
(exio_v.loc[exio_v.index.get_level_values(1) == 'Biogas'].values >0)
exioSUT = pysut.SupplyUseTable(V=exio_v.values, U=exio_u.values)
calval = pd.read_csv('/home/jakobs/Documents/IndEcol/BONSAI/data/exiobase_hsut_3317/Calorific_values.csv', header=[4], index_col=[0,1,2,3])
calval.fillna(1, inplace=True)
aggregation_matrix = pd.read_csv('/home/jakobs/Documents/IndEcol/BONSAI/data/exiobase_hsut_3317/aggregation_matrix_exiobase.csv', header=[0,1,2], index_col=[0,1,2,3])
aggregation_matrix.ix[aggregation_matrix.iloc[:,109]==1,109]
calval.iloc[141:146]
calval.loc[calval.index.get_level_values(3) == 'C_GASE']
Natural_gas_calval = calval.loc[calval.index.get_level_values(3) == 'C_GASE'].values[0] #natural gas production has the same caloric value for now but we keep the a specific value for now
print(Natural_gas_calval)
calval.iloc[141:146]/Natural_gas_calval

new_aggregation_matrix = np.zeros((48*200,48*164))
for i in range(48):
    new_aggregation_matrix[i*200:i*200+200,i*164:i*164+164] = aggregation_matrix
new_aggregation_matrix[200,164]
for i in range(48):
    #print(calval.iloc[141:146,i]/Natural_gas_calval[i])
    new_aggregation_matrix[i*200+141:i*200+146,i*164+109] = calval.iloc[141:146,i]/Natural_gas_calval[i]
j=10
new_aggregation_matrix[141+j*200,109+164*j]
exio_vagg = new_aggregation_matrix.T.dot(exio_v.values)
j=1
print(exio_vagg[24+j*164,24+j*164])
exio_vagg.shape
exio_v.loc[exio_v.index.get_level_values(3)=='C_IRON','AU'].values
aggregation_matrix
exio_uagg = new_aggregation_matrix.T.dot(exio_u.values)
exio_SUTagg_Obj = pysut.SupplyUseTable(V=exio_vagg, U=exio_uagg, regions=48)
iot_names = pd.read_csv('/home/jakobs/data/EXIOBASE/EXIOBASE_3.3.17_hsut_2011/Classifications_v_3_3_17_IOTnames.csv')

country_list = iot_names['Country code'].unique()
country_dic = {}
for i,c in enumerate(country_list):
    country_dic[c] = i


products = iot_names[['Product name', 'Product code 1']]
products.drop_duplicates(inplace=True)

prod_dic = {}
for i,prod in enumerate(products['Product name'].values):
    prod_dic[prod] = i
products[products['Product name'].str.contains('Distribution')]
supply_diag_check_eval = exio_SUTagg_Obj.supply_diag_check()
supply_diag_check_eval.shape
prod_names = iot_names.values[supply_diag_check_eval[:,2]==1,:]
print(prod_names)
market_name_indices = np.unique(prod_names[:,1], return_index=True)[-1]

market_names = prod_names[market_name_indices]
market_names[:,0] = 'GLO'
for i in range(len(market_names)):
    market_names[i,1] = 'Market for ' + market_names[i,1]
    market_names[i,2] = market_names[i,2].replace('p','m')
    market_names[i,3] = market_names[i,3].replace('C_','M_')
    if 'm40.11' in market_names[i,2]:
        market_names[i,2] = 'm40.11'
        market_names[i,1] = 'Market for Electricity'
        market_names[i,3] = 'M_ELEC'
market_name_indices = np.unique(market_names[:,1], return_index=True)[-1]
market_names = market_names[market_name_indices]
elec_markets = np.array(list(market_names[2])*48).reshape(48,5)
elec_markets[:,0] = iot_names['Country code'].unique()
market_names = np.delete(market_names, np.where(market_names[:,2]=='m40.11')[0], axis=0)
#np.append(market_names, elec_markets, axis=0)
print(elec_markets)
market_names
U_without_elec = exio_uagg.copy()
V_without_elec = exio_vagg.copy()
diag_dummy_indices = np.arange(V_without_elec.shape[0])
V_without_elec[diag_dummy_indices, diag_dummy_indices] = 0 #set to zero for principle production so we can sum off diagonal electricity

V_elecmarkets = np.zeros(len(elec_markets)) #These are the totals of the electricity used in a country, as this is the total that a national grid will provide. Format is (N_countries x 1), will be diagonalized in final V' table
U_elecmarkets = np.zeros((exio_uagg.shape[0],len(elec_markets))) #This defines the electricity mix in a national grid. This will be updated with entso data. Format (N_products*N_countries x N_countries)
elec_indices = iot_names['Product code 1'].str.contains('p40.11').values #boolean array with the electricity commodities
elec_martket_product_use = np.zeros((len(elec_markets), exio_uagg.shape[0])) #The input vector for grid electricity for the different activities, as they now draw from the grid instead of directly from producers. Format (N_countries, N_countries*N_products+N_countries)
elec_martket_product_supply = np.zeros((len(elec_markets), exio_uagg.shape[0]))

for i in range(len(elec_markets)): #loop over countries
    U_elecmarkets[elec_indices,i] = exio_uagg[elec_indices,i*164:i*164+164].sum(axis=1)
    V_elecmarkets[i] = U_elecmarkets[elec_indices,i].sum()
    elec_martket_product_use[i,i*164:i*164+164] = exio_uagg[elec_indices,i*164:i*164+164].sum(axis=0)
    elec_martket_product_supply[i,i*164:i*164+164] = V_without_elec[elec_indices,i*164:i*164+164].sum(axis=0)
    V_without_elec[elec_indices,i*164:i*164+164] = 0

V_without_elec[diag_dummy_indices, diag_dummy_indices] = exio_vagg[diag_dummy_indices, diag_dummy_indices] 
U_without_elec[elec_indices,:] = 0 #is the new partial Use table where electricity use from producers has been set to 0
elec_martket_product_supply[1,96+15+28+164]

unique_prod_indices = np.unique(prod_names[:,1], return_index=True)[-1]
excl_byproducts = prod_names[unique_prod_indices]
rm_index = []
for i,prod in enumerate(excl_byproducts):
    if 'p40.11.' in prod[2]:
        rm_index.append(i)
excl_byproducts = np.delete(excl_byproducts, np.array(rm_index), axis=0)
print(market_names)
excl_byproducts
v_market_excl_byproduct = np.zeros(len(excl_byproducts)) #create a vector for the market supply
u_market_excl_byproduct = np.zeros((U_without_elec.shape[0], len(excl_byproducts)))
excl_market_products_use = np.zeros((len(excl_byproducts),U_without_elec.shape[0]))
excl_market_products_supply = np.zeros((len(excl_byproducts),U_without_elec.shape[0]))

U_markets = U_without_elec.copy() #this now will become the final use table (without markets)
V_markets = V_without_elec.copy() #this now will become the final supply table without markets

for i,excl_prod in enumerate(excl_byproducts):
        prod_ind = prod_dic[excl_prod[1]]
        excl_prod_indices = np.where(iot_names.as_matrix()[:,2]==excl_byproducts[i,2])[0]
        u_market_excl_byproduct[excl_prod_indices,i] = exio_vagg[excl_prod_indices,excl_prod_indices]
        v_market_excl_byproduct[i] = u_market_excl_byproduct[:,i].sum()
        
        excl_byprod_countries = prod_names[np.where(prod_names[:,2] == excl_prod[2]),0][0]
        c_index = np.array([country_dic[x] for x in excl_byprod_countries])
        excl_market_products_use[i,:] = exio_uagg[c_index*164+prod_ind,:].sum(axis=0) #every acticity that normally buys the product from a country where this is only a by product now needs to buy this from the market.
        U_markets[c_index*164+prod_ind,:] = 0 #We then set the use from the particular countries to 0 as they only buy from the market.
        
        excl_market_products_supply[i,:] = exio_vagg[c_index*164+prod_ind,:].sum(axis=0) #We now need to move the byproduction in these countries to the market. We can simpy sum over the columns as the production of the other countries in this country will always be 0 and won't affect the sum
        V_markets[c_index*164+prod_ind,:] = 0 #as above for the use but now for the byproduct supply
#here we assemble everything into one Use and one Supply table
exio_dim = exio_uagg.shape[0]
n_elecmarket = len(elec_markets)
n_excl_bp = len(excl_byproducts)
full_dimension =exio_dim + n_elecmarket + n_excl_bp

U_full = np.zeros((full_dimension,full_dimension))
V_full = np.zeros((full_dimension,full_dimension))
#first insert the main use and supply tables
U_full[:exio_dim,:exio_dim] = U_markets
V_full[:exio_dim,:exio_dim] = V_markets

#now insert_electricity markets
U_full[exio_dim:exio_dim+n_elecmarket,:exio_dim] = elec_martket_product_use
U_full[:exio_dim,exio_dim:exio_dim+n_elecmarket] = U_elecmarkets
V_full[exio_dim:exio_dim+n_elecmarket,:exio_dim] = elec_martket_product_supply
V_full[exio_dim:exio_dim+n_elecmarket,exio_dim:exio_dim+n_elecmarket] = np.diag(V_elecmarkets)

#now insert exclusive byproduct markets
U_full[exio_dim+n_elecmarket:,:exio_dim] = excl_market_products_use
U_full[:exio_dim,exio_dim+n_elecmarket:] = u_market_excl_byproduct

V_full[exio_dim+n_elecmarket:,:exio_dim] = excl_market_products_supply
V_full[exio_dim+n_elecmarket:,exio_dim+n_elecmarket:] = np.diag(v_market_excl_byproduct)



print(elec_markets)
print(market_names)
#create market electricity product names<jupyter_output>[['AU' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['AT' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['BE' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['BR' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['BG' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['CA' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['CN' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['CY' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['CZ' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['DK' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['EE' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['FI' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['FR' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['DE' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['GR' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['HU' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['HR' 'Market for Electricity' 'm40.11' 'M_ELEC' 'TJ']
 ['IN' 'Market for Electricity' 'm40.11' 'M_ELEC[...]<jupyter_text>$Z = U = $<jupyter_code>def make_IOT(U_full, V_full):
    Z_full = U_full - V_full + np.diag(np.diag(V_full))
    x_dummy = np.diag(V_full).copy()
    x_dummy[x_dummy == 0] = 1
    A_full = Z_full/x_dummy
    return Z_full, A_full
Z_full, A_full = make_IOT(U_full, V_full)<jupyter_output><empty_output><jupyter_text>* check for every non zero entry in entsoe if it has principle production exiobase* update information<jupyter_code>#Z_e_update, A_e_update = update_electricity_mix()

def update_electricity_mix(entso_u, entso_v, entso_act, entso_prod, U,V):
    
    country_ind = country_dic[country_code]
    return U_update, V_update



bentso_prods = json.load(open('/home/jakobs/Documents/IndEcol/BONSAI/bentso-data/entsoe-beebee-2016/entsoe-products.json'))
bentso_acts = json.load(open('/home/jakobs/Documents/IndEcol/BONSAI/bentso-data/entsoe-beebee-2016/entsoe-activities.json'))
bentso_use = np.load('/home/jakobs/Documents/IndEcol/BONSAI/bentso-data/entsoe-beebee-2016/entsoe-use.npy')
bentso_supply = np.load('/home/jakobs/Documents/IndEcol/BONSAI/bentso-data/entsoe-beebee-2016/entsoe-supply.npy')
bentso_acts
print(32*12)
print(672/32.)
pd.DataFrame(columns=np.array(bentso_acts).T.tolist(), index=np.array(bentso_prods).T.tolist(), data=bentso_use)
set(country_concord['bentso_Country_code'].unique()).difference(set(europe_exio))
country_concord
prod_cor = pd.read_csv('/home/jakobs/Documents/IndEcol/BONSAI/Correspondence-tables/final_tables/tables/exiobase_to_bentso_activities.csv')
len(prod_cor['Bentso technology'].dropna())
ACTIVITY_MAPPING = {
    'Fossil Hard coal': 'Production of electricity by coal',
    'Fossil Brown coal/Lignite': 'Production of electricity by coal',
    'Fossil Gas': 'Production of electricity by gas',
    'Fossil Coal-derived gas': 'Production of electricity by gas',
    'Nuclear': 'Production of electricity by nuclear',
    'Hydro Pumped Storage': 'Production of electricity by hydro',
    'Hydro Run-of-river and poundage': 'Production of electricity by hydro',
    'Hydro Water Reservoir': 'Production of electricity by hydro',
    'Wind Offshore': 'Production of electricity by wind',
    'Wind Onshore': 'Production of electricity by wind',
    'Fossil Oil': 'Production of electricity by petroleum and other oil derivatives',
    'Fossil Oil shale': 'Production of electricity by petroleum and other oil derivatives',
    'Biomass': 'Production of electricity by biomass and waste',
    'Fossil Peat': 'Production of electricity by coal',
    'Waste': 'Production of electricity by biomass and waste',
    'Solar': 'Production of electricity by solar thermal',
    'Other renewable': 'Production of electricity by tide, wave, ocean',
    'Geothermal': 'Production of electricity by Geothermal',
    'Other': 'Production of electricity nec',
    'Marine': 'Production of electricity nec',
    'Grid': "Market for Electricity",
}
FLOW_MAPPING = {
    'Fossil Hard coal': 'Electricity by coal',
    'Fossil Brown coal/Lignite': 'Electricity by coal',
    'Fossil Gas': 'Electricity by gas',
    'Fossil Coal-derived gas': 'Electricity by gas',
    'Nuclear': 'Electricity by nuclear',
    'Hydro Pumped Storage': 'Electricity by hydro',
    'Hydro Run-of-river and poundage': 'Electricity by hydro',
    'Hydro Water Reservoir': 'Electricity by hydro',
    'Wind Offshore': 'Electricity by wind',
    'Wind Onshore': 'Electricity by wind',
    'Fossil Oil': 'Electricity by petroleum and other oil derivatives',
    'Fossil Oil shale': 'Electricity by petroleum and other oil derivatives',
    'Biomass': 'Electricity by biomass and waste',
    'Fossil Peat': 'Electricity by coal',
    'Waste': 'Electricity by biomass and waste',
    'Solar': 'Electricity by solar thermal',
    'Other renewable': 'Electricity by tide; wave; ocean',
    'Geothermal': 'Electricity by Geothermal',
    'Other': 'Electricity nec',
    'Marine': 'Electricity nec',
    'Grid': "Electricity",
}


np.unique([ACTIVITY_MAPPING[x] for x in ACTIVITY_MAPPING.keys()])
np.unique([x for x in ACTIVITY_MAPPING.keys()])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1470067080116441, 'fraction_numerical': 0.035122136438425516, 'mean_word_length': 6.4824810606060606, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '267609', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6788, 'n_tokens_neox': 6360, 'n_words': 1147, 'path': 'notebooks/.ipynb_checkpoints/BONSAI system model development-checkpoint.ipynb', 'repo_name': 'BONSAMURAIS/system-model', 'size': '15802'}