simuTechProExtreme.py

# %%
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import FUNC

from process_data import process_data

#%% Load data: Benchmark
label_T = pd.read_csv("./Eora26_2016_bp/labels_T.csv",index_col=0)
label_T['Country_Sector'] = label_T.apply(lambda row: f"{row['Country']}_{row['Sector']}", axis=1)
label_T.ISO.replace({'SUD':'SDN', "SDS":'SSD'},inplace=True)
label_T = label_T.iloc[:4914,:]

IOZ_0, IOV_0, IOX_0, IOF_0 = process_data()

RR = label_T['Country'].nunique() #189
NN = label_T['Sector'].nunique() #26
RRNN = RR*NN  # 4914
FF = 6 
UU = 1 # Factor


#%% set scenarios about losses in agriculture sectors
prep_path = './secondary_data/ag_lossExtreme/'
cty_pr = pd.read_csv(prep_path + "country_name.csv")

cty_pr.ISONAME.replace({'YUG':'SRB', "MLI":'MLT','TKM':'UGA',},inplace=True)
label_T.ISO.replace({'SUD':'SDN', "SDS":'SSD'},inplace=True)
cty_io = label_T['ISO'].unique().tolist()

# cty_io_set = set(label_T['ISO'].unique().tolist())
# cty_pr_set = set(cty_pr['ISONAME'])

#cty_only_in_pr = cty_pr_set - cty_io_set
# cty_only_in_io = cty_io_set - cty_pr_set

# cty_only_in_pr = list(cty_only_in_pr)
# cty_only_in_io = list(cty_only_in_io)

# all scenarios
models = ['ACCESS-ESM1-5', 'ACCESS-CM2', 'BCC-CSM2-MR', 'CAMS-CSM1-0', 'CanESM5', 
          'CMCC-ESM2', 'MRI-ESM2-0', 'NESM3', 'TaiESM1', 'CAS-ESM2-0', 'FIO-ESM-2-0', 'MIROC6'] # total 12
ssp = ['ssp126','ssp585']
extremes = ["pExtreme", "nExtreme", "comExtreme",
            "pSevere", "nSevere", "comSevere", 
            "pModerate", "nModerate", "comModerate"] # 
ScenarioSet = [
    f"{extreme}_{model}_{ss}" for extreme in extremes
    for model in models
    for ss in ssp]

#%% Initialization

for scenario in ScenarioSet:

    print(scenario)

    #  preprocess data
    I_NN = np.eye(NN)
    I_Sum = np.tile(I_NN, (1, RR)) #Replicate diagonal matrix (NN，RRNN)
   
    IOZ_C = I_Sum.dot(IOZ_0) #(NN,RRNN) sum of intermediate use by sector
    Z_Dis_0 = IOZ_0 / np.tile(IOZ_C, (RR, 1)) #(RRNN,RRNN)，intermediate distribution coefficients, regional contribution of unique products
    Z_Dis_0[np.isnan(Z_Dis_0)] = 0  
    Z_Dis =  Z_Dis_0

    IOF_C = I_Sum.dot(IOF_0) # (NN, RR) sum of final use by sector 
    F_Dis_0 = IOF_0 / np.tile(IOF_C, (RR, 1))  #(RRNN, RR) final demand distribution coefficients
    F_Dis_0[np.isnan(F_Dis_0)] = 0
    F_Dis =  F_Dis_0

    E_CZ_0 = np.matmul(IOZ_C, np.diag(1 / IOX_0[0,0:RRNN])) # (NN,RRNN) intermediate input coefficients to produce one unit of x
    E_VA_0 = np.matmul(IOV_0, np.diag(1 / IOX_0[0,0:RRNN]))  #(1,RRNN) value-added input coefficients to produce one unit of x
    E_CZ_0[np.isnan(E_CZ_0)] = 0
    E_VA_0[np.isnan(E_VA_0)] = 0 
    E_CZ = E_CZ_0
    E_VA = E_VA_0

    no_input = np.repeat(NN, RRNN)
    mat_key = IOZ_C < 0.00001 * np.tile(IOX_0, (NN, 1)) # (NN,RRNN) # logical value, FALSE if the condition is not met

    # Overproduction module
    OverProdFa = 1.25
    OverProd = np.ones((UU, RRNN)) #(1, RRNN)
    OverProdSign = np.zeros((UU, RRNN))
    OverProdStep = np.full((UU, RRNN), (OverProdFa -1) / 35)
    OverProdUpBd = np.full((UU, RRNN), OverProdFa)

    # Stock & Order
    StockProdTime = 2 # years of production sustained by stock
    StockMat = StockProdTime * IOZ_C  #(NN,RRNN)
    StockObj = StockMat + IOZ_C ##(NN,RRNN)

    OrderMat_0 = np.hstack([IOZ_0, IOF_0]) #(RRNN,RRNN+RR) 
    OrderMat = OrderMat_0

    # # Initialization

    TT = 35# Year
    
    Sect = list(range(2)) # Agri

    shock_TT = np.zeros((NN, RR, TT)) #(NN,RR,TT)

    filename = prep_path + "w_losses_" + scenario + "_ag.csv"
    pr_shock = pd.read_csv(filename) 
    pr_shock.ISONAME.replace({'YUG':'SRB', "MLI":'MLT','TKM':'UGA',},inplace=True)
    shock = pd.DataFrame(columns=pr_shock.columns)
    for country in cty_io:
    # Check if the country exists in pr_shock
        if country in pr_shock["ISONAME"].values:
        # Append the row with the country name and values from pr_shock
            shock = pd.concat([shock, pr_shock[pr_shock["ISONAME"] == country]], ignore_index=True)
        else:
            row = [country] + [0] * (pr_shock.shape[1] - 1)
            shock = pd.concat([shock, pd.DataFrame([row], columns=pr_shock.columns)], ignore_index=True)

    # intermediate input constraints    
    for t in range(TT):
        for i in range(len(cty_io)):
            country = cty_io[i]
            if country in shock["ISONAME"].values:
                country_shock = shock.loc[shock["ISONAME"] == country].iloc[0, t+3] # begin with year 2016
                shock_TT[Sect, i, t] = country_shock / 100  # (NN,RR,TT) negative value

    #shock_TT = np.zeros((NN, RR, TT)) 

    # final demand constraints
    IOF_TT = np.zeros((RRNN, RR, TT)) # (RRNN,RR,TT)
    for i in range(TT):
        IOF_TT[:,:,i] = 1 * (i+1) * IOF_0

    ######################## Modelling ###################################################  

    IOX_TT = np.zeros((TT, RRNN))
    IOX_TT_max = np.zeros((TT, RRNN))
    IOV_Sum_TT = np.zeros((TT, RRNN))
    IOF_t = np.zeros((RRNN, RR, TT))
    OrderMat_TT = np.zeros((RRNN, RRNN+RR, TT))

    IOX_TT[0, :] = IOX_0
    IOX_TT_max[0, :] = IOX_0
    IOV_Sum_TT[0, :] = IOV_0
    IOF_t[:,:,0] = IOF_0
    OrderMat_TT[:,:,0] = OrderMat_0

    StockUse_TT = np.zeros((TT, RRNN))
    StockAdd_TT = np.zeros((TT, RRNN))
    StockMat_TT = np.zeros((TT, RRNN))

    DirectDamage = np.zeros((TT, RRNN))
    results_save =  np.empty((TT+1, 12), dtype=object)
    results_save[0,:] = ["IOV_t","IOX_t_max","IOX_t","x_dis_inter","x_dis_final", "StockUse","StockAdd",
                         "StockMat","StockGap","OrderMat_final","OrderMat_inter","OrderMat"]

    for t in range(1, TT):

        print(t)

        # Production under constraints
        
        Shock_Cons = np.zeros((UU, RRNN))  # Primary input coefficient
        Shock_Cons[0, :] = shock_TT[:, :, t].reshape(1, -1, order='F') #(1,RRNN) only production input constraint
        IOV_t = IOV_0 * (1 + Shock_Cons) * OverProd # (1,RRNN)
  
        IOX_t_max = FUNC.Production_max(StockMat, E_CZ, IOV_t, E_VA, IOX_0, NN, mat_key, OverProd)  # maximum production capacity
        IOX_t_max = np.array(IOX_t_max).reshape(1, -1) # (1,RRNN) 
  
        IOX_t = FUNC.Production(StockMat, E_CZ, IOV_t, E_VA, OrderMat, IOX_0, NN, mat_key, OverProd)  #consider orders from its clients
        IOX_t = np.array(IOX_t).reshape(1, -1) ##(1,RRNN)
        
        IOX_t_max[0,:] = np.minimum(IOX_t_max, IOX_t_max * (1 + Shock_Cons))

        IOX_t[0,:] = np.minimum(IOX_t, IOX_t * (1 + Shock_Cons))

        # AlLocation to intermediate demand
                              
        row_sums = np.sum(OrderMat, axis=1)  
        
        IOX_t_Dis_max = (OrderMat.T / row_sums).T * np.tile(IOX_t_max, (RRNN + RR, 1)).T
        
        IOX_t_Dis = (OrderMat.T / row_sums).T * np.tile(IOX_t, (RRNN + RR, 1)).T

        x_dis_inter = IOX_t_Dis[:, :RRNN]
        
        x_dis_final = IOX_t_Dis[:, RRNN:(RRNN + RR)]  
  
        # Recover of inventory and capital stock
        
        StockUse = np.multiply(np.tile(IOX_t, (NN, 1)), E_CZ) #(NN,RRNN) element-wise multiplication
  
        StockAdd = np.matmul(I_Sum, x_dis_inter) # (NN,RRNN)=(NN,RRNN)*(RRNN,RRNN)  matrix multiplication 
  
        StockMat = StockMat - StockUse + StockAdd #(NN,RRNN)

        # Stock & Order
  
        StockGap = StockObj - StockMat #(NN,RRNN)
  
        StockGap[StockGap < 0] = 0

        StockMat[StockMat < 0] = 0

        # Production technology adjustment 1 (intermediate input coeff)

        SectOrd = np.matmul(I_Sum, OrderMat[:, :RRNN]) # (NN, RRNN) order for an intermediate product i

        SectDis = np.matmul(I_Sum, x_dis_inter) # (NN, RRNN) distrubution of an intermediate product i

        E_CZ = np.matmul(SectOrd, np.diag(1 / np.sum(OrderMat, axis=1)))

        diff_Sect = SectOrd - SectDis

        diff_Coef = E_CZ_0 - E_CZ

        mask = SectOrd > SectDis # boolean array where 'True' indicates 'Case of scarcity', 'False' indicates 'Case of NO scarcity'

        TechProFa = 0.1  # Maximum rate of technical progress, 10% improvement
        TechAdapTime = 5 # timescale for this technology adaptation is set as xx years

        E_CZ = np.where(mask, E_CZ - (diff_Sect / SectOrd) * E_CZ / TechAdapTime, E_CZ + (diff_Coef / E_CZ) * diff_Coef / TechAdapTime) 
        
        E_CZ = np.maximum(E_CZ, (1 - TechProFa) * E_CZ_0) 
        E_CZ = np.minimum(E_CZ, E_CZ_0) # maximum is E_CZ_0, i.e., no technical progress
        
        # Production technology adjustment 2 (primary input coeff)

        E_VA = np.matmul(IOV_t, np.diag(1 / IOX_t[0,:]))

        mask = IOV_0 > IOV_t # boolean array where 'True' indicates 'Case of scarcity', 'False' indicates 'Case of NO scarcity'

        E_VA = np.where(mask, E_VA - (IOV_0 - IOV_t) / IOV_0 * E_VA / TechAdapTime, E_VA + (E_VA_0 - E_VA) / E_VA_0 * (E_VA_0 - E_VA) / TechAdapTime) 
        
        E_VA = np.maximum(E_VA, (1 - TechProFa) * E_VA_0)
        E_VA = np.minimum(E_VA, E_VA_0) # minimum is E_VA_0, i.e., no technical progress

        # Demand - Orders placed with suppliers based on inventory shortfalls from the previous period

        DemAdjFa = 0.5 

        Z_Dis = Z_Dis + (x_dis_inter / np.tile(np.matmul(I_Sum, x_dis_inter), (RR, 1)) - OrderMat[:, :RRNN] / np.tile(np.matmul(I_Sum, OrderMat[:, :RRNN]), (RR, 1))) #(RRNN,RRNN)

        Z_Dis = np.maximum(Z_Dis, (1 - DemAdjFa) * Z_Dis_0)
        Z_Dis = np.minimum(Z_Dis, (1 + DemAdjFa) * Z_Dis_0)

        OrderMat[:, :RRNN] = np.multiply(np.tile(StockGap, (RR, 1)), Z_Dis)

        a = OrderMat[:, 0:(RRNN)]
        
        b = IOZ_0
        
        a[b<a] = b[b<a]
        
        OrderMat[:, 0:(RRNN)] = a

        TotFinlOrd = np.matmul(I_Sum, IOF_TT[:, :, t-1]) # (NN, RR)

        F_Dis = F_Dis + (x_dis_final / np.tile(np.matmul(I_Sum, x_dis_final), (RR, 1)) - OrderMat[:, RRNN:(RRNN + RR)] / np.tile(np.matmul(I_Sum, OrderMat[:, RRNN:(RRNN + RR)]), (RR, 1))) #(RRNN,RR)

        F_Dis = np.maximum(F_Dis, (1 - DemAdjFa) * F_Dis_0)
        F_Dis = np.minimum(F_Dis, (1 + DemAdjFa) * F_Dis_0)

        OrderMat[:, RRNN:(RRNN + RR)] = np.multiply(np.tile(TotFinlOrd, (RR, 1)), F_Dis)

        a = OrderMat[:,(RRNN):(RRNN+RR)]

        b = IOF_0 #(RRNN,RR)

        a[b < a] = b[b < a]

        OrderMat[:,RRNN:(RRNN + RR)] = a

        # Overproduction module

        OverProdSign = FUNC.OverProdSignFun(StockMat, E_CZ, IOV_t, E_VA, OrderMat, IOX_0, UU, RR, NN)

        OverProd = OverProd + OverProdSign * OverProdStep
        
        OverProd[np.where(OverProd < 1)] = 1
        
        OverProd[np.where(OverProd > OverProdUpBd)] = OverProdUpBd[np.where(OverProd > OverProdUpBd)]


        IOX_TT[t, :] = IOX_t

        IOX_TT_max[t, :] = IOX_t_max

        IOV_Sum_TT[t, :] = IOX_t - np.sum(StockUse, axis=0)

        OrderMat_TT[:, :, t] = OrderMat

        StockUse_TT[t, :] = np.sum(StockUse, axis=0)
        StockAdd_TT[t, :] = np.sum(StockAdd, axis=0)
        StockMat_TT[t, :] = np.sum(StockMat, axis=0)

        Name1 = "./Results/TechProExtreme/temp/" + scenario +"_order_" +str(t) + ".csv"
        np.savetxt(Name1, OrderMat, delimiter=",")

        print(np.sum(IOX_t) / np.sum(IOX_0))

        results_save[t,:] = [np.sum(IOV_t, axis=1),np.sum(IOX_t_max, axis=1),np.sum(IOX_t, axis=1), np.sum(x_dis_inter), np.sum(x_dis_final),np.sum(StockUse),
                             np.sum(StockAdd),np.sum(StockMat),np.sum(StockGap),np.sum(OrderMat[:,RRNN:(RRNN + RR)]),
                             np.sum(OrderMat[:, 0:(RRNN)]), np.sum(OrderMat)] 
        
        DirectDamage[t,:] = IOV_0 * shock_TT[:,:,t].reshape(1, -1, order='F')

    # Results analysis
    # SectLossRate = (I_Sum.dot(IOX_TT.T).T - I_Sum.dot(IOX_0.T).T) / (I_Sum.dot(IOX_0.T).T) * 100 # (TT, NN)

    # ReglIOX = []
    # for i in range(0, 189):
    #    cols = IOX_TT[:, i*26:(i+1)*26]
    #    col_sum = cols.sum(axis=1)
    #    ReglIOX.append(col_sum[:, np.newaxis])
    # ReglIOX = np.concatenate(ReglIOX, axis=1) # (TT, RR)

    # ReglLossRate = (ReglIOX - ReglIOX[0,:]) / ReglIOX[0,:] * 100 # (TT, RR)

    # DirectDamageSum = np.sum(DirectDamage)

    # plt.plot(np.sum(IOX_TT, axis=1))

    # sum for agri sector
    csum = np.zeros((RR, RR + RR, TT))
    bs_csum = np.load("./Results/TechProExtreme/temp/ag/order_0_cty.npy")
    csum[:,:,0] = bs_csum 

    for y in range(1, TT):

        data = OrderMat_TT[:,:,y]
        inter = data[:, 0:RRNN]
        final = data[:, (RRNN):(RRNN+RR)]

        inter_csum_col = np.zeros((RRNN, RR)) # sum along column
        for i in range(0,189):
            inter_csum_col[:, i] = inter[:, (i * 26):((i + 1) * 26)].sum(axis=1) # (RRNN, RR) sum over all sectors

        inter_csum = np.zeros((RR, RR))  # sum along row
        final_csum = np.zeros((RR, RR))
        for i in range(0,189): 
            temp1 = inter_csum_col[(i * 26):((i + 1) * 26),:]
            inter_csum[i, :] = temp1[0:2, :].sum(axis=0) # (RR,RR)
            temp2 = final[(i * 26):((i + 1) * 26),:]
            final_csum[i, :] = temp2[0:2, :].sum(axis=0) # (RR,RR)

        csum[:, :, y] = np.hstack([inter_csum, final_csum])
    
    Name1 = "./Results/TechProExtreme/temp/ag/" + scenario +"_order"
    np.save(Name1, csum)
 
    Name2 = "./Results/TechProExtreme/"  + scenario + "_IOX.csv"
    np.savetxt(Name2, IOX_TT, delimiter=",")
    name3 = "./Results/TechProExtreme/" + scenario + "_DirectDamage.csv" 
    np.savetxt(name3, DirectDamage, delimiter=",")
    
# %%