#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Aug  8 12:36:18 2025

@author: jdaniel
"""
import numpy as np
from matplotlib import pyplot as plt
import pickle

with open('emResults_revised.pkl', 'rb') as f:   
    resultsMcEmiss = pickle.load(f)
    
print(np.shape(resultsMcEmiss))

emCat = ['Installation','Use','Landfill','Decomissioning','Weibull Scale','Lifetime']
marCat = ['Domestic Refrig','Commercial Refrig','Refrigerated Containers',\
          'Discontinous Panels','Spray Foam','Boardstock',\
          'Block and Pipe','Pipe in Pipe','Block - Pipe','Block Foam Slab','Integral Skin',\
              'Continuous Panels']
names=[]
namesbar = []
for i in range(76):
    if (i == 0):
        names.append('Baseline')
        namesbar.append('Baseline')
    elif (i == 1):
        names.append('Market Segmentation')
        namesbar.append('Market Segmentation')
    elif (i == 2):
        names.append('Production Emission')
        namesbar.append('Production Emission')
    elif (i == 3):
        names.append('Solvent Use')
        namesbar.append('Solvent Use')
    elif (i >= 4):
        iuse = int((i-4)/12)
        juse = i-4-iuse*12
        print(iuse, juse)
        names.append(str(emCat[iuse] + '\n' + marCat[juse]))
        namesbar.append(str(emCat[iuse] + ' - ' + marCat[juse]))
nyrmax = 35
year = np.zeros(nyrmax)
for i in range(nyrmax):
    year[i] = 1989+i
a = np.shape(resultsMcEmiss)
std = np.zeros(a[0])   #number of parameters perturbed
val = np.zeros(a[1])   #number of monte carlo simulations
for i in range(a[0]):
    for j in range(a[1]):
        if (np.max(resultsMcEmiss[i,j,:]) >= 1):
            tmp = resultsMcEmiss[i,j,0:nyrmax]-resultsMcEmiss[0,j,0:nyrmax]
        else: 
            tmp = [0,0]
        val[j] = np.std(tmp)
  #      plt.plot(tmp)
    std[i] = np.mean(val)
  #  plt.title(i)
  #  plt.show() 
index = np.argsort(std)
index = index[::-1]

fig, ax = plt.subplots(3,3,figsize=(5,6),dpi=300,)
plt.subplots_adjust(wspace=.4,hspace=.6)
for i in range(9):
    plt.subplot(3,3,i+1)
    for j in range(a[1]):
   # for j in range(9):
        tmp = resultsMcEmiss[index[i],j,0:nyrmax]-resultsMcEmiss[0,j,0:nyrmax]#/1000
        plt.plot(year,tmp/1000,color='black',linewidth=.5)
    print(index[i])
#    plt.title(index[i]) 
    plt.plot([1989,2025],[0,0])
    plt.xticks(fontsize=7)
    plt.yticks(fontsize=7)
    plt.title(names[index[i]],fontsize=7)
    if (i == 0 or i == 3 or i == 6):
        plt.ylabel('Gg/yr',fontsize=7)
    if (i == 6 or i == 7 or i == 8):
        plt.xlabel('Year',fontsize=7)
    plt.ylim(-11,15)
plt.suptitle('Emissions Variability from Individual Input Paramters')
plt.savefig('sensitivity.png')
plt.show()

std = std/1000
fig, ax = plt.subplots(figsize=(5,6),dpi=300)
num=30
for i in range(num):
    plt.barh(num-1-i,std[index[i]],color='black')
    plt.text(std[index[i]]+.050,num-1-i-.2,namesbar[index[i]],fontsize=6)
plt.xlim(0,2)
ax.set_yticks([])
plt.title('Average Emissions Variability')
plt.xlabel('Standard Deviation (Gg/yr)')
plt.savefig('variabilityBar.png')
plt.show()