# Load PET data
## Amyloid
FBP <- read_excel(paste0(dfPATH,"av45.xlsx")) %>%
    ## Change AV45 in Tracer column to FBP
    mutate(Tracer = "FBP")

fbb <- read_excel(paste0(dfPATH,"fbb.xlsx"))

pib <- read_excel(paste0(dfPATH,"pib.xlsx"))

## Tau
tau <- read_excel(paste0(dfPATH,"tau.xlsx")) %>%
    select(ID, Tracer, PET_Session, PET_Date) %>%
    ## Change AV1451 in Tracer column to FTP
    mutate(Tracer = "FTP")

## FDG
fdg <- read_excel(paste0(dfPATH,"fdg.xlsx")) %>%
    select(ID, Tracer, PET_Session, PET_Date)

## Combine ID and PET_Session for each av45, fbb, and pib
PET <- bind_rows(
  FBP %>% select(ID, Tracer, PET_Session, PET_Date),
  pib %>% select(ID, Tracer, PET_Session, PET_Date),
  fbb %>% select(ID, Tracer, PET_Session, PET_Date),
  tau %>% select(ID, Tracer, PET_Session, PET_Date),
  fdg %>% select(ID, Tracer, PET_Session, PET_Date)
)

## extract PET_year from MRdate
PET <- PET %>%
    mutate(PET_Date = as.Date(PET_Date)) %>%
    mutate(PET_year = as.numeric(format(PET$PET_Date, "%Y")))

## Count the number of Tracer per PET_year
PET_counts <- PET %>%
  #Get unique PET_year and Tracer combinations
  distinct(PET_year, Tracer) %>%
  #Create a complete grid of PET_year and Tracer
  complete(PET_year, Tracer) %>%
  #Join with counts from original data
  left_join(
    PET %>%
      group_by(PET_year, Tracer) %>%
      summarise(n = n(), .groups = "drop"),
    by = c("PET_year", "Tracer")
  ) %>%
  #Replace NA with 0
  mutate(n = replace_na(n, 0)) %>%
  #Arrange and calculate cumulative sum
  arrange(Tracer, PET_year) %>%
  group_by(Tracer) %>%
  mutate(cumulative_n = cumsum(n)) %>%
  ungroup()


## Determine stacked order
PET_counts <- PET_counts %>%
    mutate(Tracer = fct_relevel(Tracer, "FDG", "AV-1451", "Florbetaben", "FBP", "PIB"))

## Edit color scale
custom_palette <- paletteer_d("rcartocolor::Prism")[c(8,1,3,4,2)]


# Plot
p <- PET_counts %>% 
  ggplot( aes(x=PET_year, y=n, fill=Tracer, text=Tracer)) +
    geom_area( ) +
    scale_fill_manual(values = custom_palette) +
    scale_x_continuous(breaks = seq(min(PET_counts$PET_year), max(PET_counts$PET_year), by = 2)) +
    labs(y = "Total Acquisitions", x="Year", fill="Tracer") +
    theme(legend.position="none") +
    ggtitle("PET Acquisitions by Year") +
    theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5))

# Plot
p2 <- PET_counts %>% 
  ggplot(aes(x = PET_year, y = cumulative_n, fill = Tracer, text = Tracer)) +
  geom_area() +
  scale_fill_manual(values = custom_palette) +
  scale_x_continuous(breaks = seq(min(PET_counts$PET_year), max(PET_counts$PET_year), by = 2)) +
  labs(y = "Total Acquisitions", x = "Year", fill = "Tracer") +
  theme(legend.position = "none") +
  ggtitle("Accumulation of PET Acquisitions Over Time") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))