punnett_rectangles.R

install.packages(c('data.table','IsoSpecR','ggplot2','ggthemes'))

library(data.table)
library(IsoSpecR)
library(ggplot2)
library(ggthemes)

isoDist = function(mol, P=.99, ...) as.data.table(IsoSpecify(molecule=mol, P, showCounts=1,...))
marginal = function(mol, cutoff, space=0) isoDist(mol, cutoff, algo=3)[order(prob, decreasing=TRUE)][,rprob:=prob/prob[1]][,R:=cumsum(rprob) + (0:(nrow(.SD)-1))*space][,L:=R-rprob ]

get_R = function(mol=c(C=1000,H=2002), P=.999, space=0){
  R = isoDist(mol, P)
  R = R[order(prob, decreasing = TRUE)]
  maxIso = R[which.max(prob)] # top isotopologue
  maxIsoP = maxIso$prob
  minCisoP = R[C12==maxIso$C12][which.min(prob)]$prob # min-prob-isotopologue on the C axis
  minCsubisoRatio = minCisoP/maxIsoP/2 # 2 just to be sure
  minHisoP = R[H1==maxIso$H1][which.min(prob)]$prob # min-prob-isotopologue on the H axis
  minHsubisoRatio = minHisoP/maxIsoP/2 # 2 just to be sure
  C = marginal(mol['C'], minCsubisoRatio, space)
  H = marginal(mol['H'], minHsubisoRatio, space)
  R = merge(R, C, by=c('C12','C13'), suffixes=c('','.C'))
  R = merge(R, H, by=c('H1','H2'),   suffixes=c('','.H'))
  setnames(R, c('rprob','R','L'), c('rprob.C','R.C','L.C'))
  R = R[order(prob, decreasing = T)]
  R
}

mol = c(C=1000, H=2002)*100
P = .999
R = get_R(mol, P)
# R = get_R(c(C=10000,H=20002), .999, .1)
R$col = c(rep('red', 500), rep('orange', 500), rep('yellow', 500), rep('darkblue', nrow(R)-500*3))

ggplot(R) +
  geom_rect(aes(xmin=L.C, xmax=R.C, ymin=L.H, ymax=R.H, fill=I(col)), color='black') +
  theme_tufte() +
  theme(axis.title=element_blank(),
        axis.text=element_blank(),
        axis.ticks=element_blank())