################################
## Generating coverage figures #
################################

## remove (almost) everything in the working environment
rm(list = ls())

# packages
library(readxl)
library(dplyr)
library(RColorBrewer)

# Set working directory based on this.path
if (!require(this.path)) {
  install.packages("this.path")
  require(this.path)
}
setwd(dirname(this.path()))



## ---------
## 1. Load data
## ---------

## Paths (adjust if needed)
path_main   <- "siblings-extended.xlsx"
path_bal    <- "sample-for-balance.xlsx"
output_pdf  <- "genealogy_coverage_vs_rejection.pdf"

## Read main scholars file
profs <- as.data.frame(read_excel(path_main))

## ---------
## 2. Recreate profsgood (the "good genealogies" subsample)
##    This matches the filtering logic in your big script.
## ---------

# Step 1: drop those with <2 sibships AND unknown father's year of death
profsgood <- profs[
  profs$SibshipSize > 2 |
    (profs$SibshipSize < 3 & !is.na(profs$YearDeathFather)),
]

# Step 2: drop those with one sibship and one descendant
profsgood <- profsgood[
  profsgood$SibshipSize > 1 |
    profsgood$NumberDescendants > 1,
]

## ---------
## 3. Time series setup
## ---------

pts <- 70                  # number of evaluation points
sha <- 1 - pts / 100       # window width in quantile terms

## Order both datasets by refdate, like in the big script
profsgood <- profsgood[order(profsgood$refdate), ]
profs     <- profs[order(profs$refdate), ]

ns    <- nrow(profsgood)
nsall <- nrow(profs)

## We'll store:
## - sharereject: share of genealogies *rejected* by filters
## - date: cohort tick labels for the x-axis
results <- data.frame(
  sharereject = numeric(pts),
  stringsAsFactors = FALSE
)

for (i in 1:pts) {
  
  idx1 <- as.integer(1 + (i - 1) * ns * (1 - sha) / (pts - 1))
  idx2 <- as.integer((i - 1) * ns * (1 - sha) / (pts - 1) + ns * sha)
  
  t1 <- profsgood$refdate[idx1]
  t2 <- profsgood$refdate[idx2]
  
  # sharereject = 1 - (good genealogies / all genealogies)
  num_good <- nrow(profsgood[profsgood$refdate >= t1 & profsgood$refdate < t2, ])
  den_all  <- nrow(profs[profs$refdate >= t1 & profs$refdate < t2, ])
  
  results[i, "sharereject"] <- if (den_all == 0) NA_real_ else 1 - num_good / den_all
}

## Build x-axis tick labels ("date" in your main script)
date <- numeric(pts)
for (i in 1:pts) {
  date[i] <- round(quantile(profsgood$refdate, i/100 + .15)[[1]], 0)
}

## ---------
## 4. Compute sharewithgenealogy from the balance file
##    (coverage: how many professors HAVE a genealogy at all)
## ---------

allprofs <- as.data.frame(read_excel(path_bal))
allprofs <- allprofs[order(allprofs$refdate), ]
ns_balance <- nrow(allprofs)

results_balance <- data.frame(
  sharewithgenealogy = numeric(pts),
  stringsAsFactors   = FALSE
)

for (i in 1:pts) {
  
  idx1b <- as.integer(1 + (i - 1) * ns_balance * (1 - sha) / (pts - 1))
  idx2b <- as.integer((i - 1) * ns_balance * (1 - sha) / (pts - 1) + ns_balance * sha)
  
  tb1 <- allprofs$refdate[idx1b]
  tb2 <- allprofs$refdate[idx2b]
  
  winb <- allprofs$refdate >= tb1 & allprofs$refdate < tb2
  
  num_with <- nrow(allprofs[winb & allprofs$exist == 1, ])
  den_all2 <- nrow(allprofs[winb, ])
  
  # share of professors WITH genealogy in that birth window
  results_balance[i, "sharewithgenealogy"] <-
    if (den_all2 == 0) NA_real_ else num_with / den_all2
}

## ---------
## 5. Plot both lines on the same figure
## ---------

# colour palette
coul <- colorRampPalette(brewer.pal(12, "Set3"))(21)

pdf(file = output_pdf, width = 10, height = 7)

# First line: sharewithgenealogy (coverage)
plot(
  results_balance$sharewithgenealogy,
  type = "l",
  xlab = "Year of Birth",
  ylab = "Share",
  lwd  = 2,
  frame.plot = TRUE,
  ylim = c(0, 1),
  col  = coul[1],
  main = "Genealogy coverage and rejection over time",
  sub  = paste("Coverage Nobs =", nrow(allprofs),
               "| Quality Nobs =", nrow(profs)),
  axes = FALSE
)

# Second line: sharereject (rejection rate of genealogies)
lines(
  results$sharereject,
  lwd = 2,
  col = coul[4]
)

# X-axis: same labelling logic you already use in your main code
axis(
  1,
  at = c(1,10,20,30,40,50,60,pts),
  labels = c(
    as.integer(date[1]),
    date[10],
    date[20],
    date[30],
    date[40],
    date[50],
    date[60],
    date[pts]
  )
)

# Y-axis: 0 to 1 in 0.1 steps
axis(
  2,
  at = seq(0, 1, by = 0.1),
  labels = c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)
)

box()

legend(
  x = 10,
  y = 0.9,
  legend = c("Share with genealogy", "Share rejected"),
  col    = c(coul[1], coul[4]),
  lwd    = 2,
  box.lty = 0
)

dev.off()
