---
title: "figures"
output:
  bookdown::word_document2:
    fig_caption: yes
---

```{r setup, echo = FALSE, include = FALSE, message = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  echo = FALSE,
  comment = "#>",
  fig.path = "../figures/",
  dpi = 300
)

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse,
               janitor,
               here,
               wbstats,
               ISOcodes,
               viridis,
               hrbrthemes,
               wesanderson,
               glue,
               ggridges,
               patchwork,
               kableExtra,
               readxl,
               flextable)

pal <- wes_palette("Cavalcanti1", 5, type = "discrete")
extrafont::loadfonts()

```

```{r libraries}

library(here)
```

```{r functions}
# load data wrangling functions
source(here("analysis", "R", "wrangler_functions.R"))

```

```{r load-data, include=FALSE}

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

```

# Fig. 1

```{r ntiles-total}

p1 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Expenditure (trn€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy footprint (EJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p3 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon footprint (MtCO2eq)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 


p_top = p1 + p2 + p3

```

```{r ntiles-intensity-violin}


p1 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_co2eq_kg)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_co2eq_kg = weighted.mean(pe_co2eq_kg,total_fd_me)),
               aes(y=pe_co2eq_kg, yend=pe_co2eq_kg, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per expenditure (kgCO2eq/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_energy_use_mj)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me)),
               aes(y=pe_energy_use_mj, yend=pe_energy_use_mj, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy intensity per expenditure (MJ/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

dat3 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj)

p3 = dat3 %>%
  ggplot(aes(x=factor(eu_q_rank), y=intensity_e_c)) +
  geom_violin(aes(weight=total_energy_use_tj), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat3 %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj)),
               aes(y=intensity_e_c, yend=intensity_e_c, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per energy (gCO2eq/TJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 


p_bottom = p2 + p1 + p3

#p3 & theme(plot.margin =  margin(1, 1, 1, 1, "mm"))
```

```{r , fig.width=12, fig.height=8}

a = p_top / p_bottom + plot_annotation(tag_levels = 'a') +
  plot_layout(guides = 'collect')  & 
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        legend.position = 'bottom',
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

ggsave(here("analysis", "figures", "figure1.pdf"), device=cairo_pdf)

```

```{r figure1, out.width="98%", fig.cap="Household expenditure and environmental footprints and intensities across European expenditure deciles. Total expenditures (a), energy footprint (b), and carbon footprint (c) per decile. Energy intensity of consumption as energy footprint per expenditure (d), carbon intensity of consumption as carbon footprint per expenditure (e), and carbon intensity of energy as carbon footprint per energy footprint (f)."}
knitr::include_graphics(here::here("analysis", "figures", "figure1.pdf"))
```

# Fig.2

```{r , fig.width=8, fig.height=2}

pdat_basket = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), by="sector_agg_id") %>% 
  group_by(eu_q_rank) %>%
  mutate(value = total_fd_me/sum(total_fd_me)*100) %>%
  select(five_sectors, eu_q_rank, 
         value) %>%
  mutate(indicator = "Expenditure share (%)")

pdat_int_co2eq = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), by="sector_agg_id") %>% 
  mutate(value = (total_co2eq_kg)/(total_fd_me*1000000)) %>%
  select(five_sectors, eu_q_rank, 
         value) %>%
  mutate(indicator = "Carbon intensity (kgCO2eq/€)")

pdat_int_energy = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), by="sector_agg_id") %>% 
  mutate(value = (total_energy_use_tj)/(total_fd_me)) %>%
  select(five_sectors, eu_q_rank,
         value) %>%
  mutate(indicator = "Energy intensity (MJ/€)")

library(viridis)
pdat = pdat_int_co2eq %>%
  bind_rows(pdat_int_energy) %>%
  bind_rows(pdat_basket) 

pdat$facet = factor(pdat$indicator, levels = c("Expenditure share (%)", "Energy intensity (MJ/€)", "Carbon intensity (kgCO2eq/€)"))

a = ggplot(pdat, aes(x=value, y=five_sectors, color=eu_q_rank)) +
  #geom_point(shape=18, alpha=0.75) +
  geom_text(label="I", alpha=0.75, fontface="bold") +
  scale_color_viridis(option = "A", end = 0.8,
                      direction = -1, 
                      breaks = c(1, 5, 10), 
                      name="European\ndecile") +
  facet_wrap(~facet, scales="free_x") +
  theme_minimal() +
  labs(x="",y="")+
  theme(text=element_text(size=13, family="Liberation Sans Narrow"),
        panel.spacing = unit(2.5, "lines"))


ggsave(here("analysis", "figures", "figure2.pdf"), device=cairo_pdf)
```

```{r figure2, out.width="100%", fig.cap="Sectoral expenditure shares and carbon intensities of European expenditure deciles. Share of expenditure per final demand sector of total spending per decile in percent (a) and carbon intensity per final demand sector in kgCO2eq/€."}
knitr::include_graphics(here::here("analysis", "figures", "figure2.pdf"))
```

# Fig.3

```{r , fig.width=8, fig.height=5}

pdat = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), by="sector_agg_id") %>% 
  select(five_sectors, eu_q_rank, contains("pae_co2eq"), contains("pae_energy")) %>%
  pivot_longer(cols=-c(five_sectors, eu_q_rank), names_to="indicator", values_to="value") %>%
  mutate(indicator_type = if_else(str_detect(indicator, "co2eq"), 
                                  "tCO2eq per adult eq", 
                                  "Energy GJ per adult eq")) %>%
  filter(indicator != "pae_co2eq_t", 
         indicator != "pae_energy_use_gj")
  
  
pal <- wes_palette("Zissou1", 4, type = "discrete")
pal = pal[c(2,1,3,4)]
  
p1 = ggplot(pdat %>% filter(indicator_type == "tCO2eq per adult eq"), 
       aes(x=(eu_q_rank), y=value, fill=indicator)) +
  geom_col(position = position_stack()) +
  facet_wrap(~five_sectors, ncol = 5) +
  scale_fill_manual(values = pal, name="Location", labels=c("Direct", "Domestic", "Europe", "non-Europe")) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10), 
                     labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) +
  #theme_ipsum() +
  theme_minimal() +
  theme(text=element_text(family="Liberation Sans Narrow")) +
  labs(x="", y="Carbon footprint tCO2eq per adult eq") +
  theme(legend.position = "bottom", 
          axis.text.x = element_text(angle = 90, size = 12),
          axis.text.y = element_text(size = 12),
          strip.text = element_text(size = 12),
          axis.title.y = element_text(size = 10))


p2 = ggplot(pdat %>% filter(indicator_type == "Energy GJ per adult eq"), 
       aes(x=(eu_q_rank), y=value, fill=indicator)) +
  geom_col(position = position_stack()) +
  facet_wrap(~five_sectors, ncol = 5) +
  scale_fill_manual(values = pal, name="Location", labels=c("Direct", "Domestic", "Europe", "non-Europe")) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10), 
                     labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) +
  #theme_ipsum() +
  theme_minimal() +
  theme(text=element_text(family="Liberation Sans Narrow")) +
  labs(x="", y="Energy footprint GJ per adult eq") +
  theme(legend.position = "bottom", 
          axis.text.x = element_text(angle = 90, size = 12),
          axis.text.y = element_text(size = 12),
          strip.text = element_text(size = 12),
          axis.title.y = element_text(size = 10))

a = p1 / p2 + plot_layout(guides = "collect") & theme(legend.position = 'bottom', 
                                                  plot.margin =  margin(1, 1, 1, 1, "mm"),
                                                  panel.spacing = unit(c(1, 1, 1, 1), "mm"))


hm = pdat %>% filter(indicator_type == "Energy GJ per adult eq") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value))

ggsave(here("analysis", "figures", "figure3.pdf"))
```

```{r figure3, out.width="100%", fig.cap="Energy and carbon footprints by final consumption sector and European expenditure decile in 2015, further broken down by emission source location."}
knitr::include_graphics(here::here("analysis", "figures", "figure3.pdf"))
```

# Fig.4

```{r}
pal <- wes_palette("Zissou1", 4, type = "discrete")
pal = pal[c(2,1,3,4)]

pdat_int_energy = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>% 
  mutate(intensity_energy = (total_energy_use_tj)/(total_fd_me)) %>%
  select(five_sectors, eu_q_rank, 
         intensity_energy) %>%
  filter(eu_q_rank == 10) %>%
  select(-eu_q_rank)

pdat_final_demand = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>%
  left_join(pdat_int_energy, by="five_sectors") %>%
  mutate(total_energy_use_tj_new = (total_fd_me)*intensity_energy) %>%
  mutate(total_energy_use_tj_diff = total_energy_use_tj-total_energy_use_tj_new) %>%
  select(five_sectors, eu_q_rank, 
         total_energy_use_tj_diff, total_energy_use_tj_new) %>%
  group_by(eu_q_rank) %>%
  summarise(total_energy_use_tj_new = sum(total_energy_use_tj_new)*0.000001,
            total_energy_use_tj_diff = sum(total_energy_use_tj_diff)*0.000001) %>%
  pivot_longer(-c(eu_q_rank), names_to = "indicator", values_to = "value")


p1 = ggplot(pdat_final_demand, aes(x=eu_q_rank, y=value, fill=indicator, alpha=indicator)) +
  geom_col(position = position_stack()) +
  scale_fill_manual(values=c(pal[1], pal[2]), 
                    labels=c("2015", "Best technology"), name="Energy\nfootprint") +
  scale_alpha_manual(values=c(0.3,1),labels=c("2015", "Best technology"), name="Energy\nfootprint")+
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10), 
                     labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) +
  labs(y="Energy footprint (EJ)", x="") +
  theme_minimal() +
  coord_flip() +
  theme(legend.position="bottom")

```

```{r}

pdat_energy_country = get_sector_summary_by_country_quintile_direct() %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>%
  left_join(pdat_int_energy, by="five_sectors") %>%
  mutate(total_energy_use_tj_new = (total_fd_me)*intensity_energy) %>%
  #mutate(total_energy_use_tj_save = total_energy_use_tj_new/total_energy_use_tj*100) %>%
  select(iso2, 
         total_energy_use_tj_new, total_energy_use_tj) %>%
  group_by(iso2) %>%
  summarise(total_energy_use_tj_save = 100 - sum(total_energy_use_tj_new)/sum(total_energy_use_tj)*100)

```

```{r include=FALSE}

library(rworldmap)
library(ggthemes)



df_country_mean = pdat_energy_country %>%
  mutate(iso2 = if_else(iso2 == "EL", "GR", iso2)) %>%
  mutate(iso2 = if_else(iso2 == "UK", "GB", iso2)) %>%
  left_join(ISO_3166_1 %>% select(iso2 = Alpha_2, iso3 = Alpha_3), by="iso2")


quantile_rank_map = joinCountryData2Map(df_country_mean, 
                                        joinCode = "ISO_A3", 
                                        nameJoinColumn = "iso3")
quantile_rank_map_poly = fortify(quantile_rank_map) #extract polygons 

quantile_rank_map_poly = merge(quantile_rank_map_poly, quantile_rank_map@data, by.x="id", by.y="ADMIN", all.x=T)
quantile_rank_map_poly = quantile_rank_map_poly %>% arrange(id, order)

```

```{r dle-country, fig.width=8}
map = ggplot() + 
  geom_polygon(data = quantile_rank_map_poly, aes(long, lat, group = group), color="black") +
  geom_polygon(data = quantile_rank_map_poly, aes(long, lat, group = group, 
                                        fill=total_energy_use_tj_save)) + 
  #scale_fill_gradient(name="Mean decile \nrank", low="skyblue", high = "skyblue4", na.value = "#EEEEEE") +
  scale_fill_viridis(direction = -1, discrete = F, na.value = "white", name="Energy\nsavings (%)") +
  theme_map() +
  labs(x="",y="")+
  coord_map("bonne", lat0 = 50,xlim = c(-9, 40), ylim = c(38, 68), clip="on") +
  theme(legend.position = c(0.85, 0.6))
  
```

```{r fig.height=7, fig.width=11}

a = p1 + map + plot_annotation(tag_levels = 'a') +
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        legend.position = 'bottom',
        axis.title.y = element_text(size=13, hjust = 0.5),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

ggsave(here("analysis", "figures", "figure4.pdf"))
```

```{r figure4, out.width="100%", fig.cap="Energy savings through sectoral best current technology by expenditure deciles a) and country b)."}
knitr::include_graphics(here::here("analysis", "figures", "figure4.pdf"))
```

# Fig.5

```{r }

pdat_int_energy = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>% 
  mutate(intensity_energy = (total_energy_use_tj)/(total_fd_me)) %>%
  select(five_sectors, eu_q_rank, 
         intensity_energy) %>%
  filter(eu_q_rank == 10) %>%
  select(-eu_q_rank)

pdat_final_demand = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>%
  left_join(pdat_int_energy, by="five_sectors") %>%
  mutate(total_energy_use_tj_new = (total_fd_me)*intensity_energy) %>%
  mutate(total_energy_use_tj_diff = total_energy_use_tj-total_energy_use_tj_new) %>%
  select(eu_q_rank,total_energy_use_tj_new) %>%
  group_by(eu_q_rank) %>%
  summarise(total_energy_use_tj_new = sum(total_energy_use_tj_new)) %>%
  mutate(pae_energy_use_tj = total_energy_use_tj_new/33417583,
         pae_energy_use_gj = pae_energy_use_tj*1000)

df_energy_deciles = pdat_final_demand %>%
  select(eu_q_rank, pae_energy_use_gj)

ineq_curr = df_energy_deciles$pae_energy_use_gj[10]/df_energy_deciles$pae_energy_use_gj[1]

```

```{r }

library(readxl)

df_scenario_info = read_excel(here("analysis/data/raw/scenarios.xlsx"), sheet="overview") %>%
  select(scenario, fe_gj_aeu = final_energy_gj_per_aeu_2050,
         ccs_required = primary_energy_fossil_w_ccs2050_ej,
         description) %>%
  arrange(fe_gj_aeu) %>%
  mutate(fe_gj_aeu = round(fe_gj_aeu),
         ccs_required = round(ccs_required)) 

mea = c(min(df_scenario_info$fe_gj_aeu),max(df_scenario_info$fe_gj_aeu))
mer = c(min(df_scenario_info$fe_gj_aeu),53)

c_mean_mea = round((27+mea[2])/2)
c_mean_mer = round((mer[1]+mer[2])/2)

```

```{r , eval = FALSE}

# run once to save file. If changing scenarios included or input data, need to re-run and save 'scenarios_fine.rds'

# vectorized function that returns scaled quantiles given 
#quantile index column and column with quantile averages
#qidx: quantile index
#qavg_to_scale: column to scale
#first_target: target value of first quantile
#mean_target: target mean of scaled quantiles
scaled_quantiles <- function(.data, 
                          qidx, 
                          qavg_to_scale, 
                          first_target, 
                          mean_target) {
  
  # cumbersomely extract current quantile mean
  mean_current = .data %>%
    ungroup() %>%
    summarise(mean_cur = first(mean({{qavg_to_scale}}))) %>%
    pull(mean_cur)
    
  # cumbersomely extract current first wuantile value
  first_current = .data %>%
    ungroup() %>%
    arrange({{qidx}}) %>%
    summarise(first_cur = first({{qavg_to_scale}})) %>%
    pull(first_cur)
  
  df_tmp = .data %>%
    mutate(tmp = {{qavg_to_scale}}*mean_target/mean_current)
  
  first_tmp = df_tmp$tmp[1]
  
  df_tmp = df_tmp %>%
    mutate(scaled = mean_target-(mean_target-tmp) * (mean_target-first_target)/(mean_target-first_tmp)) %>%
    select({{qidx}}, scaled) %>%
    mutate(v_mean = mean_target, v_first = first_target)
}

## run once to save file
df_all = NULL
for (min_energy in seq(from=mer[1], to=mer[2], by=0.25)) {
  for (mean_energy in seq(from=mea[1], to=mea[2], by=0.25)) {
    if (min_energy <= mean_energy) {
      df_all = df_all %>%
      bind_rows(df_energy_deciles %>%
                  scaled_quantiles(eu_q_rank, pae_energy_use_gj, min_energy, mean_energy))
    }
  }
}

saveRDS(df_all, here("analysis/data/derived/scenarios_fine.rds"))
```

```{r , fig.width=7, fig.height=5.5}

round_by = 10

df_all = readRDS(here("analysis/data/derived/scenarios_fine.rds")) %>%
  filter(eu_q_rank %in% c(1,10)) %>%
  group_by(v_mean, v_first) %>%
  summarise(ratio = last(scaled)/first(scaled)) %>%
  mutate(bin_ratio = if_else((ratio*100)%%round_by > round_by*0.5, 
                             ratio*100+(round_by-(ratio*100)%%round_by),
                             ratio*100-(ratio*100)%%round_by)) %>%
  group_by(bin_ratio, v_first) %>%
  summarise(v_mean = mean(v_mean)) %>%
  mutate(bin_ratio = bin_ratio*0.01)

df_scenario = df_all %>%
  filter(v_mean %in% df_scenario_info$fe_gj_aeu) %>%
  left_join(df_scenario_info, by=c("v_mean"="fe_gj_aeu"))

library(wesanderson)

a = df_all %>%
  ggplot(aes(x=v_first, y=bin_ratio, fill=v_mean)) +
  geom_tile() +
  geom_hline(yintercept = ineq_curr, alpha=0.8, color="grey", linetype=2) +
  geom_line(data=df_scenario, aes(color=scenario, group=scenario)) +
  annotate(geom="text", x=max(df_all$v_first)-5.7,y=ineq_curr+0.6,label = "Current (2015) 10:10 ratio") +
  #scale_fill_viridis("Mean energy\navailable") +
  scale_fill_gradient("Mean energy\navailable (GJ/ae)",
                      low=wes_palette("Chevalier1")[3], 
                      high = wes_palette("Rushmore1")[4]) +
  scale_color_manual(values=wes_palette("Darjeeling1")) +
  theme_minimal() +
  labs(x="Minimum energy requirement (GJ/ae)", y="Maximum energy inequality (10:10 ratio)", color = "Scenario")+
  theme(text=element_text(family="Liberation Sans Narrow"),
        axis.text.x = element_text(size = 13),
        axis.text.y = element_text(size = 13)) +
  scale_y_continuous(breaks = c(2.5,5,7.5,10,12.5)) # +
  #theme_ipsum()

ggsave(here("analysis", "figures", "figure5.pdf"))
```

```{r figure5, out.width="70%", fig.align="center", fig.cap="Mean energy available for Europe in decarbonisation scenarios, positioned in option space between a range of minimum energy requirements and range of associated maximum inequality. All expenditure deciles have 'best technology' already."}
knitr::include_graphics(here::here("analysis", "figures", "figure5.pdf"))
```