#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon May 13 08:37:19 2024

@author: jdaniel
"""

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from get_weibull import *
from get_consumption import *
from scipy.stats import qmc, norm, truncnorm, lognorm
from outputPresent import *
from AltMarket import *
import pickle
import copy
from joblib import Parallel, delayed
import os

class otherParams:
    def __init__(self):
        self.prodemiss = 1-1/1.05   #0.05
        self.prodemisssig = 1.
        self.solvent = .1#0.093
        self.solventsig = .5
        self.iyears = 200   #number of years to track values after production
        
class outputInfo:
    def __init__(self, nregions,nmarkets,iyears):
        self.prodemiss = np.zeros([nregions,nmarkets,iyears])
        self.totalemiss = np.zeros([nregions,nmarkets,iyears])
        self.cprodemiss = np.zeros([nregions,nmarkets,iyears])
        self.solventemiss = np.zeros([nregions,nmarkets,iyears])
        self.csolventemiss = np.zeros([nregions,nmarkets,iyears])
        self.installemiss = np.zeros([nregions,nmarkets,iyears])
        self.cinstallemiss = np.zeros([nregions,nmarkets,iyears])
        self.activebank = np.zeros([nregions,nmarkets,iyears])
        self.totalbank = np.zeros([nregions,nmarkets,iyears])
        self.activeemiss = np.zeros([nregions,nmarkets,iyears])
        self.cactiveemiss = np.zeros([nregions,nmarkets,iyears])
        self.activedecom = np.zeros([nregions,nmarkets,iyears])
        self.decomemiss = np.zeros([nregions,nmarkets,iyears])
        self.decomemissalt = np.zeros([nregions,nmarkets,iyears])
        self.cdecomemiss = np.zeros([nregions,nmarkets,iyears])
        self.cdecomemissalt = np.zeros([nregions,nmarkets,iyears])
        self.inactivebank = np.zeros([nregions,nmarkets,iyears])
        self.inactivebankalt = np.zeros([nregions,nmarkets,iyears])
        self.inactiveemiss = np.zeros([nregions,nmarkets,iyears])
        self.cinactiveemiss = np.zeros([nregions,nmarkets,iyears])
        self.inactiveemissalt = np.zeros([nregions,nmarkets,iyears])
        self.cinactiveemissalt = np.zeros([nregions,nmarkets,iyears])

def getReleaseParams(filename,sheet):
    fileinfo = pd.read_excel(filename, sheet_name=sheet)
    markets = fileinfo.columns[1::]
    nummarkets = len(markets)
    tmp = fileinfo.iloc[0::,1::]
    vals = np.zeros(tmp.shape)
    vals[:] = fileinfo.iloc[0::,1::]
    return(nummarkets,markets,vals)


def outputKernels(i, output, params, emissvals, emissvalsNEA):
    iyears = params.iyears
    nregions = 10
    iflag = 0

    # === Initialize arrays ===
    prod_emiss_scalar = 1. / (1. - params.prodemiss) - 1.
    solvent_scalar = params.solvent / 2.

    # Shape (nregions,)
    prodemiss = np.full((nregions,), prod_emiss_scalar)
    solventemiss_0 = np.full((nregions,), solvent_scalar)
    solventemiss_1 = np.full((nregions,), solvent_scalar)

    # Emissions from install (installemiss[0])
    installemiss_0 = emissvals[0] * (1. - params.solvent)
    installed = (1. - params.solvent) * (1 - emissvals[0])

    # === Get Weibull profiles ===
    weibullactive, weibullemiss = getweibull(i, params, emissvals)  # shape: (nregions, iyears)

    # === Fill outputs ===
    output.prodemiss[:, i, 0] = prodemiss
    output.solventemiss[:, i, 0] = solventemiss_0
    output.solventemiss[:, i, 1] = solventemiss_1
    output.installemiss[:, i, 0] = installemiss_0

    output.activebank[:, i, :] = weibullactive * installed
    output.activebank[:, i, 0] = 0  # conservation
    output.activeemiss[:, i, :] = weibullemiss * installed

    # === Compute activedecom ===
    output.activedecom[:, i, 0] = installed - output.activebank[:, i, 1] - output.activeemiss[:, i, 0]
    output.activedecom[:, i, 1:iyears-1] = output.activebank[:, i, 1:iyears-1] - \
                                           output.activebank[:, i, 2:iyears] - \
                                           output.activeemiss[:, i, 1:iyears-1]

    # === Decom emissions ===
    output.decomemiss[:, i, :] = output.activedecom[:, i, :] * emissvals[3]
    output.decomemissalt[:, i, :] = output.activedecom[:, i, :] * emissvals[3] * 0

    # === Inactive banks ===
    for j in range(1, iyears):
        output.inactivebank[:, i, j] = (
            output.inactivebank[:, i, j-1] * (1. - emissvals[2]) +
            output.activedecom[:, i, j-1] -
            output.decomemiss[:, i, j-1]
        )
        output.inactiveemiss[:, i, j-1] = output.inactivebank[:, i, j-1] * emissvals[2]

        # alt emissions set to 0
        output.inactivebankalt[:, i, j] = 0
        output.inactiveemissalt[:, i, j-1] = 0

    # === Check for negative banks ===
    if np.any(output.inactivebank[:, i, :] < 0):
        iflag = 1

    # === Cumulative sums ===
    output.cprodemiss[:, i, :] = np.cumsum(output.prodemiss[:, i, :], axis=1)
    output.csolventemiss[:, i, :] = np.cumsum(output.solventemiss[:, i, :], axis=1)
    output.cinstallemiss[:, i, :] = np.cumsum(output.installemiss[:, i, :], axis=1)
    output.cactiveemiss[:, i, :] = np.cumsum(output.activeemiss[:, i, :], axis=1)
    output.cdecomemiss[:, i, :] = np.cumsum(output.decomemiss[:, i, :], axis=1)
    output.cinactiveemiss[:, i, :] = np.cumsum(output.inactiveemiss[:, i, :], axis=1)
    output.cdecomemissalt[:, i, :] = np.cumsum(output.decomemissalt[:, i, :], axis=1)
    output.cinactiveemissalt[:, i, :] = np.cumsum(output.inactiveemissalt[:, i, :], axis=1)

    # === Repeat for NEA (region index 6) ===
    # To avoid duplication, do once with adjusted emissvals
    for region in [6]:
        emissvals_region = emissvalsNEA
        prodemiss = 1. / (1. - params.prodemiss) - 1.
        installed = (1. - params.solvent) * (1 - emissvals_region[0])
        installemiss = emissvals_region[0] * (1. - params.solvent)

        output.prodemiss[region, i, 0] = prodemiss
        output.solventemiss[region, i, 0] = solvent_scalar
        output.solventemiss[region, i, 1] = solvent_scalar
        output.installemiss[region, i, 0] = installemiss

        weibullactive, weibullemiss = getweibull(i, params, emissvals_region)
        output.activebank[region, i, :] = weibullactive * installed
        output.activebank[region, i, 0] = 0
        output.activeemiss[region, i, :] = weibullemiss * installed

        output.activedecom[region, i, 0] = installed - output.activebank[region, i, 1] - output.activeemiss[region, i, 0]
        output.activedecom[region, i, 1:iyears-1] = (
            output.activebank[region, i, 1:iyears-1] -
            output.activebank[region, i, 2:iyears] -
            output.activeemiss[region, i, 1:iyears-1]
        )

        output.decomemiss[region, i, :] = output.activedecom[region, i, :] * emissvals_region[3]
        output.decomemissalt[region, i, :] = 0

        for j in range(1, iyears):
            output.inactivebank[region, i, j] = (
                output.inactivebank[region, i, j-1] * (1. - emissvals_region[2]) +
                output.activedecom[region, i, j-1] -
                output.decomemiss[region, i, j-1]
            )
            output.inactiveemiss[region, i, j-1] = output.inactivebank[region, i, j-1] * emissvals_region[2]

        output.inactivebankalt[region, i, :] = 0
        output.inactiveemissalt[region, i, :] = 0

        # Cumulative NEA
        output.cprodemiss[region, i, :] = np.cumsum(output.prodemiss[region, i, :])
        output.csolventemiss[region, i, :] = np.cumsum(output.solventemiss[region, i, :])
        output.cinstallemiss[region, i, :] = np.cumsum(output.installemiss[region, i, :])
        output.cactiveemiss[region, i, :] = np.cumsum(output.activeemiss[region, i, :])
        output.cdecomemiss[region, i, :] = np.cumsum(output.decomemiss[region, i, :])
        output.cinactiveemiss[region, i, :] = np.cumsum(output.inactiveemiss[region, i, :])
        output.cdecomemissalt[region, i, :] = 0
        output.cinactiveemissalt[region, i, :] = 0

        if np.any(output.inactivebank[region, i, :] < 0):
            iflag = 1

    return output, weibullactive, weibullemiss, iflag

def runmodel(years,nregions,nummarkets,params,con,kernel,emissvals,emissvalsNEA):
    ieurope = [2]#[0,1,2,3,4,5,6,7,8,9,10]#2   # index for Europe for capturing decomposition
    results = outputInfo(nregions,nummarkets,params.iyears)
    year0 = years[0]
    for k in range(nregions):
        for i in range(0,params.iyears):
            year1 = year0 + float(i)
            for j in range(0,i+1):
                year = year0 + float(j)
                if (year <= years[-1]):
                    results.activebank[k,:,i] = results.activebank[k,:,i] + con[j,k,:]*kernel.activebank[k,:,i-j]
                    results.inactivebank[k,:,i] = results.inactivebank[k,:,i] + con[j,k,:]*kernel.inactivebank[k,:,i-j]
                    results.csolventemiss[k,:,i] = results.csolventemiss[k,:,i] + con[j,k,:]*kernel.csolventemiss[k,:,i-j]
                    results.cprodemiss[k,:,i] = results.cprodemiss[k,:,i] + con[j,k,:]*kernel.cprodemiss[k,:,i-j]
                    results.cinstallemiss[k,:,i] = results.cinstallemiss[k,:,i] + con[j,k,:]*kernel.cinstallemiss[k,:,i-j]
                    results.cactiveemiss[k,:,i] = results.cactiveemiss[k,:,i] + con[j,k,:]*kernel.cactiveemiss[k,:,i-j]
                    results.cdecomemiss[k,:,i] = results.cdecomemiss[k,:,i] + con[j,k,:]*kernel.cdecomemiss[k,:,i-j]
                    results.cinactiveemiss[k,:,i] = results.cinactiveemiss[k,:,i] + con[j,k,:]*kernel.cinactiveemiss[k,:,i-j]
    
                    results.prodemiss[k,:,i] = results.prodemiss[k,:,i] + con[j,k,:]*kernel.prodemiss[k,:,i-j]
                    results.solventemiss[k,:,i] = results.solventemiss[k,:,i] + con[j,k,:]*kernel.solventemiss[k,:,i-j]
                    results.installemiss[k,:,i] = results.installemiss[k,:,i] + con[j,k,:]*kernel.installemiss[k,:,i-j]
                    results.activeemiss[k,:,i] = results.activeemiss[k,:,i] + con[j,k,:]*kernel.activeemiss[k,:,i-j]
                    results.decomemiss[k,:,i] = results.decomemiss[k,:,i] + con[j,k,:]*kernel.decomemiss[k,:,i-j]
                    results.inactiveemiss[k,:,i] = results.inactiveemiss[k,:,i] + con[j,k,:]*kernel.inactiveemiss[k,:,i-j]
            if (year1 >= 2002 and k==ieurope[0]):   #year1 >= 2002
                kk = k
                jj = 0   #domestic refrigeration market
                results.decomemiss[kk,jj,i] = 0
                results.cdecomemiss[kk,jj,i] = results.cdecomemiss[kk,jj,i-1]
                results.inactiveemiss[kk,jj,i] = 0
                results.cinactiveemiss[kk,jj,i] = results.inactiveemiss[kk,jj,i-1]
                if (kk != 6):
                    results.inactivebank[kk,jj,i] = results.inactivebank[kk,jj,i-1]*(1-emissvals[2,kk])# - con[j,:] * (kernel.activedecom[:,i-j]-kernel.decomemiss[:,i-j])
                elif (kk == 6):
                    results.inactivebank[kk,jj,i] = results.inactivebank[kk,jj,i-1]*(1-emissvalsNEA[2,kk])# - con[j,:] * (kernel.activedecom[:,i-j]-kernel.decomemiss[:,i-j])

            results.totalbank[k,:,i] = results.activebank[k,:,i] + results.inactivebank[k,:,i]
            results.totalemiss[k,:,i] = results.inactiveemiss[k,:,i] + results.decomemiss[k,:,i]+\
                            results.activeemiss[k,:,i]+results.installemiss[k,:,i]+results.solventemiss[k,:,i]+\
                            results.prodemiss[k,:,i]
    return(results)

def getAltParams(sample_num,params):
    normaldistr = 0   # 1 for normal, 0 for uniform
    paramsmc = [otherParams() for i in range(sample_num)]
    cliplow, cliphi = 0.0000001, 1000
    dimension = 2
    a,b = np.zeros(dimension), np.zeros(dimension)
    mean = [params.prodemiss,params.solvent]
    std = [params.prodemiss*params.prodemisssig, params.solvent*params.solventsig]
    if (normaldistr == 1):
        for i in range(0, dimension):
            a[i] = (cliplow-mean[i])/std[i]
            b[i] = (cliphi-mean[i])/std[i]
            lhd = qmc.LatinHypercube(d=dimension, optimization="random-cd").random(n=sample_num)
            sample = truncnorm(a,b, loc=mean, scale=std).ppf(lhd)
    elif (normaldistr == 0):
        l_bounds, u_bounds = np.zeros(2), np.zeros(2)
        l_bounds[0] = mean[0] - std[0]
        l_bounds[1] = mean[1] - std[1]
        u_bounds = [2*mean[0],2*mean[1]]
        u_bounds = [mean[0]+std[0],mean[1]+std[1]]
        sampler = qmc.LatinHypercube(d=dimension)
        sample_noscale = sampler.random(n=sample_num)
        sample = qmc.scale(sample_noscale, l_bounds, u_bounds)
    for i in range(1,sample_num):
        paramsmc[i].prodemiss = sample[i,0]
        paramsmc[i].solvent = sample[i,1]
   
    return(paramsmc)

def getAltReleaseParams(sample_num,emissvals):
    isize = np.shape(emissvals)
    cliplow, cliphi = 0.0000001, 1#1000
    dimension = int(isize[0] * isize[1] / 2)
    a,b = np.zeros((dimension)), np.zeros((dimension))
    mean = np.zeros(dimension) 
    std = np.zeros(dimension)
  #  print('izize',isize,emissvals[0,1])
    for ii in range(0,int(emissvals.shape[0]/2)):
        if (ii <= 3):
            cliplow, cliphi = .0000001, 1
        else:
            cliplow, cliphi = .0000001, 1000
        for j in range(0,isize[1]):
            i = ii*isize[1] + j
            mean[i] = emissvals[ii,j]
            std[i] = emissvals[ii,j] * emissvals[ii+int(emissvals.shape[0]/2),j]
      #      print(ii,j,mean[i],std[i])
            a[i] = (cliplow-mean[i])/std[i]
            b[i] = (cliphi-mean[i])/std[i]
    lhd = qmc.LatinHypercube(d=dimension, optimization="random-cd").random(n=sample_num)
    sval = np.ndarray.flatten(emissvals[6::,:])
    svalln = sval[0:48]
    svalln = np.sqrt(np.log(svalln*svalln+1))
    std = sval[48::]
    sampleln = lognorm.ppf(lhd[:,0:48], s=svalln)
    sampleln = sampleln / np.exp(svalln**2/2)
    samplenorm = norm.ppf(lhd[:,48::], loc=mean[48::], scale=std)
    sampleln, samplenorm = np.array(sampleln), np.array(samplenorm)
    sampleout = np.concatenate((sampleln,samplenorm),axis=1)
    sampleout = np.reshape(sampleout,(sample_num,int(isize[0]/2),isize[1]))
    skeep = np.copy(sampleout)
    sampleout = sampleout * emissvals[0:6,:]
    sampleout[:,4::,:] = skeep[:,4::,:]
    sampleout[0,:,:] = emissvals[0:int(isize[0]/2),:]
    for i in range(0,sample_num):
        iover0 = np.where(sampleout[i,0,:] > 1)
        iover3 = np.where(sampleout[i,3,:] > 1)
        ish0 = np.shape(iover0)
        ish3 = np.shape(iover3)
        if (ish0[1] > 0): 
            sampleout[i,0,iover0] = 1
        if (ish3[1] > 0): 
            sampleout[i,3,iover3] = 1
    
    return(sampleout)
def outputKernels_worker(iSens,ii, paramvalsmc, emissvalsmc, emissvalsmcNEA,
                          marketSegRefrigmc2001, marketSegRefrigmc2005, marketSegRefrigmc2008,
                          marketSegNoRefrigmc2001, marketSegNoRefrigmc2005, marketSegNoRefrigmc2008,
                          years, consumptionin, nregions, nummarkets, markets, check):

    import copy
    from outputPresent import makefig1

    params = copy.deepcopy(paramvalsmc[0])
    emissvals = copy.deepcopy(emissvalsmc[0,:,:])
    emissvalsNEA = copy.deepcopy(emissvalsmcNEA[0,:,:])
    if (iSens == 1):
        iiU = ii
    else:
        iiU = 0
    if (iSens == 2):
        params.prodemiss = paramvalsmc[ii].prodemiss
    elif (iSens == 3):
        params.solvent = paramvalsmc[ii].solvent
    elif (iSens >= 4 and iSens <= 15):
        emissvals[0,iSens-4] = emissvalsmc[ii,0,iSens-4]
    elif (iSens >= 16 and iSens <= 27):
        emissvals[1,iSens-16] = emissvalsmc[ii,1,iSens-16]
    elif (iSens >= 28 and iSens <= 39):
        emissvals[2,iSens-28] = emissvalsmc[ii,2,iSens-28]
    elif (iSens >= 40 and iSens <= 51):
        emissvals[3,iSens-40] = emissvalsmc[ii,3,iSens-40]
    elif (iSens >= 52 and iSens <= 63):
        emissvals[4,iSens-52] = emissvalsmc[ii,4,iSens-52]
    elif (iSens >= 64 and iSens <= 75):
        emissvals[5,iSens-64] = emissvalsmc[ii,5,iSens-64]
    
    marketSegRefrig2001, marketSegRefrig2005, marketSegRefrig2008 = \
        marketSegRefrigmc2001[iiU,:,:], marketSegRefrigmc2005[iiU,:,:], marketSegRefrigmc2008[iiU,:,:]
    marketSegNoRefrig2001, marketSegNoRefrig2005, marketSegNoRefrig2008 = \
        marketSegNoRefrigmc2001[iiU,:,:], marketSegNoRefrigmc2005[iiU,:,:], marketSegNoRefrigmc2008[iiU,:,:]

    consump_bymarket = getSegmentedValsMult(
        years, consumptionin,
        marketSegRefrig2001, marketSegRefrig2005, marketSegRefrig2008,
        marketSegNoRefrig2001, marketSegNoRefrig2005, marketSegNoRefrig2008
    )

    outputKernel = outputInfo(nregions, nummarkets, params.iyears)
    weibull, weibullem = np.zeros((nummarkets, 200)), np.zeros((nummarkets, 200))

    for i in range(nummarkets):
        outputKernel, wa, we, iflag = outputKernels(i, outputKernel, params, emissvals[:, i], emissvalsNEA[:, i])
        weibull[i, :] = wa
        weibullem[i, :] = we
        if ii == 0 and i == 0:
            makefig1(outputKernel)
        if check:
            print(markets[i])
            getcheck(i, outputKernel)

    results = runmodel(years, nregions, nummarkets, params, consump_bymarket, outputKernel, emissvals, emissvalsNEA)
    return ii, consump_bymarket, outputKernel, weibull, weibullem, results

 
#------------------------------------------------------------------------
#          MAIN PROGRAM   #for paper submitted on 1/16/2025
#------------------------------------------------------------------------
iParamsList = [0,1,2,3,9,21,20,40,4,45,8,16,5]
iParams = 76#len(iParamsList)#75
sample = 50
resultsMcEmiss = np.zeros([iParams,sample,200])
for iSens in range(iParams):
#for iSens in iParamsList:
    paramsin = otherParams()
    check = 0
    sample_num = sample #number of Monte Carlo samples
    tmp1, tmp2 = np.zeros(sample_num), np.zeros(sample_num)
    emissFile = '../data/emiss_params_alt.xlsx'    # this file determines the markets and life cycle parameters for each one
    emissSheet = 'testall'
    emissSheetNEA = 'NEA-Wang(xdecom)'
    emissSheetNEA = 'NEA-TEAP'
    emissSheetNEA = 'testall'   #if this is uncommented, NEA and the rest of the world are calculated the same
    segfilename = '../data/segmentation_allyears.xlsx' # this file has information to determine market fraction in each region and defines regions
                                    #must be same order of regions as production files used in get_consumption
    sheetname1 = 'use-2001'    #primary segmentation information
    sheetname2 = 'use-2005'  #using fraction of market that is filled by -141b
    sheetname3 = 'use-2008'
    nummarkets, markets, emissvalsin = getReleaseParams(emissFile,emissSheet) # get number of markets and their parameters
    nummarkets, markets, emissvalsinNEA = getReleaseParams(emissFile,emissSheetNEA) # get number of markets and their parameters for NEA region
    
    marketSegRefrig2001, marketSegNoRefrig2001 = getSegmentation(segfilename,sheetname1)     #first index is region, second is market
    marketSegRefrig2005, marketSegNoRefrig2005 = getSegmentation(segfilename,sheetname2)     #first index is region, second is market
    marketSegRefrig2008, marketSegNoRefrig2008 = getSegmentation(segfilename,sheetname3)     #first index is region, second is market
    
    years, production, consumptionin, feedstock, regions = get_consumption()     #get scaled consumption from production files; determines length of production time
    s = np.shape(consumptionin)
    nregions = len(regions)
    consump_bymarket_mc = np.zeros((sample_num,s[0],s[1],nummarkets))
    consumptionin = consumptionin*1.0
    
    paramvalsmc = getAltParams(sample_num,paramsin)
    emissvalsmc = getAltReleaseParams(sample_num,emissvalsin)
    emissvalsmcNEA = getAltReleaseParams(sample_num,emissvalsinNEA)
    nummarkuse = 12
    varmarket = 0.1
    
    
    marketrefrigmc2001 = AltMarket_FractionMethod(sample_num,nummarkuse,varmarket,marketSegRefrig2001)
    marketrefrigmc2005 = AltMarket_FractionMethod(sample_num,nummarkuse,varmarket,marketSegRefrig2005)
    marketrefrigmc2008 = AltMarket_FractionMethod(sample_num,nummarkuse,varmarket,marketSegRefrig2008)
    marketSegRefrigmc2001 = marketmc(marketrefrigmc2001,marketSegRefrig2001)
    marketSegRefrigmc2005 = marketmc(marketrefrigmc2005,marketSegRefrig2005)
    marketSegRefrigmc2008 = marketmc(marketrefrigmc2008,marketSegRefrig2008)
    nummarkuse = 12
    marketnorefrigmc2001 = AltMarket_FractionMethod(sample_num,nummarkuse,varmarket,marketSegNoRefrig2001)
    marketnorefrigmc2005 = AltMarket_FractionMethod(sample_num,nummarkuse,varmarket,marketSegNoRefrig2005)
    marketnorefrigmc2008 = AltMarket_FractionMethod(sample_num,nummarkuse,varmarket,marketSegNoRefrig2008)
    marketSegNoRefrigmc2001 = marketmc(marketnorefrigmc2001,marketSegNoRefrig2001)
    marketSegNoRefrigmc2005 = marketmc(marketnorefrigmc2005,marketSegNoRefrig2005)
    marketSegNoRefrigmc2008 = marketmc(marketnorefrigmc2008,marketSegNoRefrig2008)
    
#    with open('pickledParams.pkl', 'wb') as f:
#        pickle.dump(regions, f)
#        pickle.dump(paramvalsmc, f)
#        pickle.dump(emissvalsmc, f)
#        pickle.dump(emissvalsmcNEA, f)
#        pickle.dump(marketSegNoRefrigmc2008, f)
#        pickle.dump(marketSegRefrigmc2008, f)
#    f.close()
    #print('stopping')
    #input()
    
    results = outputInfo(nregions,nummarkets,paramvalsmc[0].iyears)
    resultsmc = [results for i in range(sample_num)]
    
    weibull, weibullkeep = np.zeros([nummarkets,200]),np.zeros([nummarkets,200])
    weibullem, weibullemkeep = np.zeros([nummarkets,200]),np.zeros([nummarkets,200])
    
    if __name__ == "__main__":
        import time
        t0 = time.time()
    
        n_jobs = 4  # Number of CPU cores
        results = Parallel(n_jobs=n_jobs)(
            delayed(outputKernels_worker)(
                iSens,ii, paramvalsmc, emissvalsmc, emissvalsmcNEA,
                marketSegRefrigmc2001, marketSegRefrigmc2005, marketSegRefrigmc2008,
                marketSegNoRefrigmc2001, marketSegNoRefrigmc2005, marketSegNoRefrigmc2008,
                years, consumptionin, nregions, nummarkets, markets, check
            )
            for ii in range(sample_num)
        )
    
        t1 = time.time()
        print(f"Parallel run finished in {t1 - t0:.2f} seconds.")
        
    results.sort(key=lambda x: x[0])
    
    consump_bymarket_mc = np.zeros((sample_num, *results[0][1].shape))
    outputKernel_list = []
    weibull_all = np.zeros((sample_num, nummarkets, 200))
    weibullem_all = np.zeros((sample_num, nummarkets, 200))
    resultsmc = []
    
    for ii, consump_bymarket, outputKernel, weibull, weibullem, results in results:
        consump_bymarket_mc[ii,:,:,:] = consump_bymarket
        outputKernel_list.append(outputKernel)
        weibull_all[ii, :, :] = weibull
        weibullem_all[ii, :, :] = weibullem
        resultsmc.append(results)
    consump_bymarket_keep = consump_bymarket_mc[0]
    outputKernelKeep = outputKernel_list[0]
    emissvals = emissvalsmc[0,:,:]
    results = resultsmc[0]
    for numSample in range(sample_num):
        tmp = resultsmc[numSample].totalemiss
        resultsMcEmiss[iSens,numSample,:]=np.sum(tmp,axis=(0,1))

if (os.path.isfile("emResults.pkl")): 
    os.remove("emResults.pkl")
with open('emResults.pkl','wb') as f:
    pickle.dump(resultsMcEmiss, f)
# Save or use: resultsmc, consump_bymarket_mc, weibull_all, etc.

ibestfit = outputPresent(years,markets,regions,consump_bymarket_keep,consump_bymarket_mc,results,resultsmc,weibullkeep,weibullemkeep)

with open('pickled.pkl', 'wb') as f:
    pickle.dump(years, f)
    pickle.dump(markets, f)
    pickle.dump(regions, f)
    pickle.dump(consump_bymarket_keep, f)
    pickle.dump(consump_bymarket_mc, f)
    pickle.dump(results, f)
    pickle.dump(resultsmc, f)
    pickle.dump(weibullkeep, f)
    pickle.dump(weibullemkeep, f)
    pickle.dump(ibestfit, f)
    pickle.dump(paramvalsmc, f)
    pickle.dump(emissvals, f)
    pickle.dump(emissvalsmc, f)
    pickle.dump(marketSegRefrigmc2001, f)
    pickle.dump(marketSegRefrigmc2005, f)
    pickle.dump(marketSegRefrigmc2008, f)
    pickle.dump(marketSegNoRefrigmc2001, f)
    pickle.dump(marketSegNoRefrigmc2005, f)
    pickle.dump(marketSegNoRefrigmc2008, f)
    pickle.dump(consumptionin, f)
    pickle.dump(nummarkets, f)
    pickle.dump(weibull, f)
    pickle.dump(weibullem, f)
    pickle.dump(check, f)
    pickle.dump(sample_num, f)
f.close()