# Illustrative GRM comparison. No IRT fitting, Q3 analysis, or MILP selection.
# Power and required N use the exact noncentral t distribution (two-sided alpha = .05,
# two equal groups), the same calculation as pwr::pwr.t.test.
t_power <- function(d, n) { # n per group; may be fractional
  df <- 2*n-2; ncp <- d*sqrt(n/2); crit <- qt(.975, df)
  pt(crit, df, ncp, lower.tail=FALSE) + pt(-crit, df, ncp)
}
required_per_group <- function(d) { # smallest whole n reaching 80% power
  n <- max(2, ceiling(2*(2.8016/d)^2))
  while(t_power(d, n) < .8 && n < 1e7) n <- n+1
  while(n > 2 && t_power(d, n-1) >= .8) n <- n-1
  n
}
simulate_shortform_power <- function(k_full=30, k_short=6, a_mean=1.6,
                                    target_zone=c(1,2.5), d_true=.4,
                                    n_total=250, penalty=.72, boost=1.1) {
  stopifnot(length(target_zone)==2,
    all(is.finite(c(k_full,k_short,a_mean,target_zone,d_true,n_total,penalty,boost))),
    k_full==floor(k_full), k_full>=2, k_full<=200,
    k_short==floor(k_short), k_short>=1, k_short<=k_full,
    n_total==floor(n_total), n_total>=4, n_total<=1000000,
    a_mean>=.2, a_mean<=4, target_zone[1]>=-3, target_zone[2]<=3,
    target_zone[1]<target_zone[2], d_true>0, d_true<=3,
    penalty>0, penalty<=1, boost>=.5, boost<=1.5)
  spread <- function(lo,hi,n) if(n==1) (lo+hi)/2 else seq(lo,hi,length.out=n)
  theta <- seq(-3,3,length.out=61)
  grm_info <- function(theta,a,b) {
    p <- cbind(1,plogis(a*outer(theta,b+c(-.9,-.3,.3,.9),'-')),0)
    dp <- a*p*(1-p)
    rowSums((dp[,1:5]-dp[,2:6])^2/pmax(p[,1:5]-p[,2:6],1e-12))
  }
  idx <- which(theta>=target_zone[1]-1e-9 & theta<=target_zone[2]+1e-9)
  stopifnot(length(idx)>0)
  specs <- list(full=list(spread(-2.2,2.2,k_full),a_mean,1),
    naive=list(spread(-.4,.4,k_short),a_mean,penalty),
    targeted=list(spread(target_zone[1]-.3,target_zone[2]-.1,k_short),a_mean*boost,1))
  forms <- lapply(specs,function(s) {
    info <- s[[3]]*Reduce('+',lapply(s[[1]],function(b) grm_info(theta,s[[2]],b)))
    rho <- mean(info[idx]/(info[idx]+1)) # Precision proxy only.
    d <- d_true*sqrt(rho)
    list(info=info,reliability=rho,d=d,sem=mean(1/sqrt(info[idx])),
      power=t_power(d,n_total/2),
      requiredN=2*required_per_group(d))
  })
  list(theta=theta,forms=forms,burden=100*(1-k_short/k_full),
    extra=forms$naive$requiredN-forms$targeted$requiredN)
}
if(sys.nframe()==0) print(simulate_shortform_power())