<jupyter_start><jupyter_text>Monthly TOD MEF Comparisons<jupyter_code>import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
%matplotlib inline
from IPython.display import display
from IPython.core.debugger import set_trace
import calendar
LABELS = ['so2_kg', 'nox_kg', 'pm25_kg', 'co2_kg', 'so2_dam_ap2', 'nox_dam_ap2',
       'pm25_dam_ap2', 'so2_dam_eas', 'nox_dam_eas', 'pm25_dam_eas',
       'co2_dam']
LABELS.sort()
LABELS_FULL = LABELS + ['dam_eas', 'dam_ap2']
dam_cols_ap2 = ['co2_dam', 'so2_dam_ap2', 'nox_dam_ap2', 'pm25_dam_ap2']
dam_cols_eas = ['co2_dam', 'so2_dam_eas', 'nox_dam_eas', 'pm25_dam_eas']
GROUPING_NAMES = ['SeasonalTOD', 'MonthTOD', 'TOD', 'YearOnly', 'Month']
GROUPING_COLS = [['year', 'season', 'hour'], ['year', 'month', 'hour'], 
        ['year', 'hour'], ['year'], ['year', 'month']]
GROUPING_NAMES_COLS = dict(zip(GROUPING_NAMES, GROUPING_COLS))<jupyter_output><empty_output><jupyter_text>Read in data<jupyter_code>def get_factor_df(kind='MEF', time='MonthTOD', region='PJM', fuel_type='FossilOnly'):
    kind_folder = 'mefs' if kind=='MEF' else 'aefs'
    
    # Read in file
    if fuel_type == 'FossilOnly':
        region_breakdown = 'isorto' if region == 'PJM' else 'nerc'
        df = pd.read_csv(os.path.join(os.pardir, 'calculated_factors', kind_folder, time, 
                                      '{}_{}.csv'.format(region_breakdown, kind_folder)),
                         index_col=GROUPING_NAMES_COLS[time])
        df = df[df[region_breakdown] == region].drop(region_breakdown, axis=1)
    else:
        if region != 'PJM':
            raise NotImplementedError('fossil-plus factors are only available for PJM')
        df = pd.read_csv(os.path.join(os.pardir, 'calculated_factors', kind_folder, time, 
                                      'pjm_fplus_{}.csv'.format(kind_folder)),
                         index_col=GROUPING_NAMES_COLS[time])
        
    # Filter MEF columns
    if kind == 'MEF':
        df = df.drop([x for x in df.columns if '-r' in x or '-int' in x], axis=1)
        df.columns = [x.replace('-est', '') for x in df.columns]
        
    # Ensure columns have numeric type
    df = df.apply(pd.to_numeric, axis=1)
    
    return df
all_dfs = {}
for kind in ['MEF', 'AEF']:
    for region in ['PJM', 'RFC']:
        for fuel_type in ['FossilOnly', 'FossilPlus']:
            for time in ['MonthTOD']:
                if region == 'RFC' and fuel_type == 'FossilPlus': 
                    continue
                print(kind, region, fuel_type, time)
                df = get_factor_df(kind=kind, time=time, region=region, fuel_type=fuel_type)
                all_dfs[(kind, region, fuel_type, time)] = df
                print(df.index.get_level_values('year').unique().values)<jupyter_output>MEF PJM FossilOnly MonthTOD
[2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017]
MEF PJM FossilPlus MonthTOD
[2016 2017]
MEF RFC FossilOnly MonthTOD
[2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017]
AEF PJM FossilOnly MonthTOD
[2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017]
AEF PJM FossilPlus MonthTOD
[2016 2017]
AEF RFC FossilOnly MonthTOD
[2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017]<jupyter_text>Full plots for SI<jupyter_code>FULL_DAMS = ['dam_ap2', 'dam_eas']
def format_title(label):
    l = label.split('_')
    if label in FULL_DAMS:
        t = 'Total damages ({})'.format('AP2' if l[1] == 'ap2' else 'EASIUR')
    else:
        t = '{0}$_{{{1}}}$ {2}'.format(l[0][:2].upper(), l[0][2:], 'emissions' if l[1] == 'kg' else 'damages')
        if len(l) > 2: t += ' ({})'.format('AP2' if l[2] == 'ap2' else 'EASIUR')
    return t
def month_tod_plot_helper(label, df_dict, ylim, year,
                     figsize, outer_dims, inner_dims, months_to_plot,
                     hspace, subtitle_break,
                     lgd_ncol, lgd_bbox_to_anchor, use_suptitle):
    fontsize=16
    fig = plt.figure(figsize=figsize)
    outer = gridspec.GridSpec(*outer_dims, wspace=-0.45, hspace=hspace)
    
    def get_month_only(df, month, year):
        df = df.reset_index().query('month=={}'.format(month)).query('year=={}'.format(year))
        df = df.drop(['year', 'month'], axis=1).set_index('hour')
        return df
    
    for i, month in enumerate(months_to_plot):
        inner = gridspec.GridSpecFromSubplotSpec(*inner_dims, subplot_spec=outer[i], wspace=0.1, hspace=0.1)
        for j, kind in enumerate([('PJM', 'FossilPlus'), ('PJM', 'FossilOnly'), ('RFC', 'FossilOnly')]):
            ax = plt.Subplot(fig, inner[j])
            fig.add_subplot(ax)
            ncol = outer_dims[1]
            
            # Plot average factors for this month
            avg = get_month_only(df_dict[('AEF',) + kind + ('MonthTOD',)], month, year)    
            avg[label].plot(ax=ax, marker='*', legend=False, label='Average',
                                  ls='--', markersize=5.5, color='green')
            
            # Plot marginal factors for this month
            marg = get_month_only(df_dict[('MEF',) + kind + ('MonthTOD',)], month, year)
            marg[label].plot(ax=ax, marker='o', legend=False, label='Marginal',
                        ls='-', markersize=3, color='blue',
                        yerr=marg['{}-se'.format(label)])
            
            # Set y limits
            ax.set_ylim(*ylim)
            
            # Add titles to plots (for both month and region/gen mix)
            #   For plots at top, indicate which region/gen mix
            if month in months_to_plot[:ncol]:
                ax.set_title('{}{}({}-{})\n\n'.format(
                    kind[0], subtitle_break, kind[1][:-4].lower(), kind[1][-4:].lower()), fontweight='bold')
            #  For plots in the middle, indicate which month (and deal with cases where
            #   we've already used the title attribute for region/gen mix labels)
            if j == 1:
                if month in months_to_plot[:ncol]:
                    ax.text(0.5,1.09, calendar.month_abbr[month],
                            horizontalalignment='center', verticalalignment='center', 
                            transform=ax.transAxes, fontsize=fontsize)
                else:
                    ax.set_title(calendar.month_abbr[month])
        
            # Add x-axis label to middle plots at bottom
            if month in months_to_plot[-ncol:] and j == 1:
                ax.set_xlabel('Hour of day (UTC-5)')
            else:
                ax.set_xlabel('')
        
            # Add y-axis label on left column only, and y-axis ticks on left plots within each month
            if (i % ncol) == 0 and j == 0:
                ylabel = 'Damage factor\n(\$/MWh)' if 'dam' in label else 'Emissions factor\n(kg/MWh)'
                ax.set_ylabel(ylabel)
            if j > 0:
                ax.set_yticklabels([])
            
            # Set larger font size
            for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
                         ax.get_xticklabels() + ax.get_yticklabels()):
                item.set_fontsize(fontsize)

            # Add legend at bottom (to last plot)
            if i == len(months_to_plot)-1 and j == 2:
                lgd = ax.legend(ncol=lgd_ncol, bbox_to_anchor=lgd_bbox_to_anchor, 
                                      loc='upper left', fontsize=fontsize)
                
    # Overall title for plot 
    if use_suptitle:
        sup = plt.suptitle('{} ({})'.format(format_title(label), year),  x=0.375, y=0.97,
                           fontsize=fontsize, fontweight='bold', fontstyle='italic')
    
#     if not os.path.exists(os.path.join('figures', 'month_tod')):
#         os.makedirs(os.path.join('figures', 'month_tod'))
#     fig.savefig(os.path.join('figures', 'month_tod','si-{}.pdf'.format(label)), 
#                 bbox_extra_artists=(lgd,sup), bbox_inches='tight')
                    
    return fig
def month_tod_si_plot(label, df_dict, ylim, year=2017):
    figsize=(22,16)
    outer_dims=(6,2)
    inner_dims=(1,6)
    months_to_plot=range(1,13)
    lgd_ncol = 2
    lgd_bbox_to_anchor=(-3.35, -0.35)
    use_suptitle=True
    hspace=0.4
    subtitle_break='\n'
    return month_tod_plot_helper(label, df_dict, ylim, year, 
                     figsize, outer_dims, inner_dims, months_to_plot,
                     hspace, subtitle_break, 
                     lgd_ncol, lgd_bbox_to_anchor, use_suptitle)
ylims_dict = dict([
    ['co2_kg', (0,900)], ['so2_kg', (0,3)], ['nox_kg', (0,1.3)], ['pm25_kg', (0,0.15)],
    ['co2_dam', (0,40)], ['so2_dam_eas', (0,100)], ['so2_dam_ap2', (0,120)],
    ['nox_dam_eas', (0,15)], ['nox_dam_ap2', (0,7)], 
    ['pm25_dam_eas', (0,15)], ['pm25_dam_ap2', (0,20)],
    ['dam_eas', (0,150)], ['dam_ap2', (0,150)]
])
month_tod_si_plot('co2_kg', all_dfs, (0,900));
# for label, ylim in ylims_dict.items():
#     display(month_tod_si_plot(label, all_dfs, ylim))
def month_tod_paper_plot(label, df_dict, ylim, year=2017):
    figsize=(12,4)
    outer_dims=(2,1)
    inner_dims=(1,3)
    months_to_plot=[5,11]
    lgd_ncol=1
    lgd_bbox_to_anchor=(1.05, 1.6)
    use_suptitle=False
    hspace=0.5
    subtitle_break=' '
    return month_tod_plot_helper(label, df_dict, ylim, year, 
                     figsize, outer_dims, inner_dims, months_to_plot,
                     hspace, subtitle_break,
                    lgd_ncol, lgd_bbox_to_anchor, use_suptitle)
month_tod_paper_plot('co2_kg', all_dfs, (0,900));
month_tod_paper_plot('so2_kg', all_dfs, (0,3));
month_tod_paper_plot('dam_ap2', all_dfs, (0,150));<jupyter_output><empty_output><jupyter_text>Presentations<jupyter_code>def pres_plot(df_dict, kind):
    fig, axes = plt.subplots(1, 2, figsize=(10, 1.5), sharey=True)
    fontsize=16
    label = 'dam_ap2'

    for i, month in enumerate([5,11]):
        ax = axes[i]

        def get_month_only(df, month, year):
                df = df.reset_index().query('month=={}'.format(month)).query('year=={}'.format(year))
                df = df.drop(['year', 'month'], axis=1).set_index('hour')
                return df
        
        # Plot average factors for this month
        avg = get_month_only(df_dict[('AEF',) + kind + ('MonthTOD',)], month, 2017)    
        avg[label].plot(ax=ax, marker='*', legend=False, label='Average',
                              ls='--', markersize=5.5, color='green')

        # Plot marginal factors for this month
        marg = get_month_only(df_dict[('MEF',) + kind + ('MonthTOD',)], month, 2017)
        marg[label].plot(ax=ax, marker='o', legend=False, label='Marginal',
                    ls='-', markersize=3, color='blue',
                    yerr=marg['{}-se'.format(label)])

        ax.set_ylim(0, 150)
        ax.set_xlabel('Hour of day (UTC-5)')
        ax.set_ylabel('Damage factor\n(\$/MWh)')

        ax.set_title(calendar.month_name[month])

        for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
                     ax.get_xticklabels() + ax.get_yticklabels()):
            item.set_fontsize(fontsize)
    
#     fig.savefig(os.path.join('figures', 'month_tod','presentation', '{}.pdf'.format(ftype)), bbox_inches='tight')
pres_plot(all_dfs, ('PJM', 'FossilOnly'))
pres_plot(all_dfs, ('PJM', 'FossilPlus'))
pres_plot(all_dfs, ('RFC', 'FossilOnly'))<jupyter_output><empty_output><jupyter_text>Stats for 2017  Helpers<jupyter_code># Mean percent error
def get_mpe(df1, df2):
    return 100 * ((df1-df2)/df2).mean()

# Mean absolute percent error
def get_mape(df1, df2):
    return 100 * (abs((df1-df2)/df2)).mean()<jupyter_output><empty_output><jupyter_text>Restrict to 2017<jupyter_code>all_dfs_2017 = {}
for key in all_dfs.keys():
    df = all_dfs[key]
    idx = df.index.names
    df = df.reset_index().query('year==2017').set_index(idx)
    all_dfs_2017[key] = df<jupyter_output><empty_output><jupyter_text>Average vs. marginal<jupyter_code>avg_mar_stats = {}
for region in ['PJM', 'RFC']:
    for fuel_type in ['FossilOnly', 'FossilPlus']:
        if region == 'RFC' and fuel_type == 'FossilPlus': 
            continue
        marg = all_dfs_2017[('MEF', region, fuel_type, 'MonthTOD')][LABELS_FULL]
        avg = all_dfs_2017[('AEF', region, fuel_type, 'MonthTOD')][LABELS_FULL]
        avg_mar_stats[('MPE', region, fuel_type)] = get_mpe(avg, marg)
        avg_mar_stats[('MAPE', region, fuel_type)] = get_mape(avg, marg)
pd.DataFrame(avg_mar_stats).round(0).astype(int)<jupyter_output><empty_output><jupyter_text>RFC vs. PJM<jupyter_code>pjm_rfc_stats = {}
for kind in ['AEF', 'MEF']:
    pjm = all_dfs_2017[(kind, 'PJM', 'FossilOnly', 'MonthTOD')][LABELS_FULL]
    rfc = all_dfs_2017[(kind, 'RFC', 'FossilOnly', 'MonthTOD')][LABELS_FULL]
    pjm_rfc_stats[(kind, 'MPE')] = get_mpe(rfc, pjm)
    pjm_rfc_stats[(kind, 'MAPE')] = get_mape(rfc, pjm)
pd.DataFrame(pjm_rfc_stats).round(0).astype(int)
# Highlighted in main text
pd.DataFrame(pjm_rfc_stats).round(0).astype(int).loc['dam_ap2', 'MEF']<jupyter_output><empty_output><jupyter_text>PJM vs. PJM+RE<jupyter_code>fossil_fplus_stats = {}
for kind in ['AEF', 'MEF']:
    fossil = all_dfs_2017[(kind, 'PJM', 'FossilOnly', 'MonthTOD')][LABELS_FULL]
    fplus = all_dfs_2017[(kind, 'PJM', 'FossilPlus', 'MonthTOD')][LABELS_FULL]
    fossil_fplus_stats[(kind, 'MPE')] = get_mpe(fplus, fossil)
    fossil_fplus_stats[(kind, 'MAPE')] = get_mape(fplus, fossil)
pd.DataFrame(fossil_fplus_stats).round(0).astype(int)
# Highlighted in main text
pd.DataFrame(fossil_fplus_stats).round(0).astype(int).loc['dam_ap2', 'MEF']<jupyter_output><empty_output><jupyter_text>Within-month variation Variation from monthly means of month TOD factors<jupyter_code>mtod_var_stats = {}
for kind in ['AEF', 'MEF']:
    for region in ['PJM', 'RFC']:
        for fuel_type in ['FossilOnly', 'FossilPlus']:
            if region == 'RFC' and fuel_type == 'FossilPlus': 
                    continue 
            df = all_dfs_2017[(kind, region, fuel_type, 'MonthTOD')][LABELS_FULL]

            # Get monthly means and put in same format as monthly TOD factors
            month_means = df.groupby([df.index.get_level_values('year'), df.index.get_level_values('month')]).mean()
            month_means['hour'] = 0
            month_means.set_index([month_means.index, 'hour'], inplace=True)
            month_means = month_means.reindex(df.index).ffill()

            # Calculate MAPE. Note that MPE is always 0 by definition.
            pes = 100 * (df - month_means) / month_means
            monthly_mape = abs(pes).groupby(
                [pes.index.get_level_values('year'), pes.index.get_level_values('month')]).mean()
            
            get_min_max = lambda d: d.apply(lambda x: '{} to {}'.format(
                int(x.min().round(0)), int(x.max().round(0))))
            mtod_var_stats[(kind, region, fuel_type, 'MAPE')] = get_min_max(monthly_mape)
pd.DataFrame(mtod_var_stats)
# Highlighted in main text
mtod_var_stats[('MEF', 'PJM', 'FossilOnly', 'MAPE')]<jupyter_output><empty_output><jupyter_text>Maximum variation<jupyter_code>mtod_var_stats = {}
for kind in ['AEF', 'MEF']:
    for region in ['PJM', 'RFC']:
        for fuel_type in ['FossilOnly', 'FossilPlus']:
            if region == 'RFC' and fuel_type == 'FossilPlus': 
                    continue 
            df = all_dfs_2017[(kind, region, fuel_type, 'MonthTOD')][LABELS_FULL]

            # Get monthly means and put in same format as monthly TOD factors
            month_means = df.groupby([df.index.get_level_values('year'), df.index.get_level_values('month')]).mean()
            month_means['hour'] = 0
            month_means.set_index([month_means.index, 'hour'], inplace=True)
            month_means = month_means.reindex(df.index).ffill()

            # Calculate MAPE. Note that MPE is always 0 by definition.
            pes = 100 * (df - month_means) / month_means
            monthly_mape = abs(pes).groupby(
                [pes.index.get_level_values('year'), pes.index.get_level_values('month')]).max()
            
            get_min_max = lambda d: d.apply(lambda x: '{} to {}'.format(
                int(x.min().round(0)), int(x.max().round(0))))
            mtod_var_stats[(kind, region, fuel_type, 'MAPE')] = get_min_max(monthly_mape)
pd.DataFrame(mtod_var_stats)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15054621333006016, 'fraction_numerical': 0.03271142751933227, 'mean_word_length': 2.848606518658479, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 46, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '574416', 'license': 'Apache-2.0', 'n_tokens_mistral': 6337, 'n_tokens_neox': 5760, 'n_words': 1154, 'path': 'factor_estimates/notebooks/MonthlyTODMEFComparisons.ipynb', 'repo_name': 'priyald17/emissions-assumptions', 'size': '16294'}