From 553f449a079609f96f153dfd014afbe024b08ccb Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 18:02:57 +0100
Subject: [PATCH 01/15] llm I choose you

---
 R/fitdistdoublecens.R     | 31 +++++++++----
 R/primary_censored_dist.R | 91 +++++++++++++++++++++++++++++++++++----
 man/dot-dpcens.Rd         | 25 ++++++++++-
 man/dot-ppcens.Rd         | 24 ++++++++++-
 man/fitdistdoublecens.Rd  |  7 +++
 5 files changed, 159 insertions(+), 19 deletions(-)

diff --git a/R/fitdistdoublecens.R b/R/fitdistdoublecens.R
index 72631d7..16d1f8d 100644
--- a/R/fitdistdoublecens.R
+++ b/R/fitdistdoublecens.R
@@ -55,6 +55,7 @@
 fitdistdoublecens <- function(censdata, distr,
                               pwindow = 1, D = Inf,
                               dprimary = stats::dunif,
+                              dprimary_name = NULL,
                               dprimary_args = list(), ...) {
   # Check if fitdistrplus is available
   if (!requireNamespace("fitdistrplus", quietly = TRUE)) {
@@ -67,9 +68,13 @@ fitdistdoublecens <- function(censdata, distr,
     stop("censdata must contain 'left' and 'right' columns")
   }
 
-  # Get the distribution functions
-  pdist <- get(paste0("p", distr))
+  if (is.null(dprimary_name)) {
+    dprimary_name <- .extract_function_name(dprimary)
+  }
 
+  # Get the distribution functions
+  pdist_name <- paste0("p", distr)
+  pdist <- get(pdist_name)
   swindows <- censdata$right - censdata$left
 
   # Create the function definition with named arguments for dpcens
@@ -83,7 +88,9 @@ fitdistdoublecens <- function(censdata, distr,
         pwindow = pwindow,
         D = D,
         dprimary = dprimary,
-        dprimary_args = dprimary_args
+        dprimary_args = dprimary_args,
+        pdist_name = pdist_name,
+        dprimary_name = dprimary_name
       )
     ))
   }
@@ -99,7 +106,9 @@ fitdistdoublecens <- function(censdata, distr,
         pwindow = pwindow,
         D = D,
         dprimary = dprimary,
-        dprimary_args = dprimary_args
+        dprimary_args = dprimary_args,
+        pdist_name = pdist_name,
+        dprimary_name = dprimary_name
       )
     ))
   }
@@ -125,14 +134,15 @@ fitdistdoublecens <- function(censdata, distr,
 #' each element in x
 #' @keywords internal
 .dpcens <- function(x, swindows, pdist, pwindow, D, dprimary,
-                    dprimary_args, ...) {
+                    dprimary_args, pdist_name, dprimary_name, ...) {
   tryCatch(
     {
       if (length(unique(swindows)) == 1) {
         dprimarycensoreddist(
           x, pdist,
           pwindow = pwindow, swindow = swindows[1], D = D, dprimary = dprimary,
-          dprimary_args = dprimary_args, ...
+          dprimary_args = dprimary_args, pdist_name = pdist_name,
+          dprimary_name = dprimary_name, ...
         )
       } else {
         # Group x and swindows by unique swindow values
@@ -144,7 +154,8 @@ fitdistdoublecens <- function(censdata, distr,
           result[mask] <- dprimarycensoreddist(
             x[mask], pdist,
             pwindow = pwindow, swindow = sw, D = D,
-            dprimary = dprimary, dprimary_args = dprimary_args, ...
+            dprimary = dprimary, dprimary_args = dprimary_args,
+            pdist_name = pdist_name, dprimary_name = dprimary_name, ...
           )
         }
 
@@ -160,13 +171,15 @@ fitdistdoublecens <- function(censdata, distr,
 #' Define a fitdistrplus compatible wrapper around pprimarycensoreddist
 #' @inheritParams pprimarycensoreddist
 #' @keywords internal
-.ppcens <- function(q, pdist, pwindow, D, dprimary, dprimary_args, ...) {
+.ppcens <- function(q, pdist, pwindow, D, dprimary, dprimary_args,
+                    pdist_name, dprimary_name, ...) {
   tryCatch(
     {
       pprimarycensoreddist(
         q, pdist,
         pwindow = pwindow,
-        D = D, dprimary = dprimary, dprimary_args = dprimary_args, ...
+        D = D, dprimary = dprimary, dprimary_args = dprimary_args,
+        pdist_name = pdist_name, dprimary_name = dprimary_name, ...
       )
     },
     error = function(e) {
diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index 303cfb9..6250cac 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -124,16 +124,53 @@ primary_censored_cdf.pcens_numeric <- function(
 #' @export
 primary_censored_cdf.pcens_pgamma_dunif <- function(
     object, q, pwindow, use_numeric = FALSE) {
-  use_numeric <- TRUE
   if (isTRUE(use_numeric)) {
     return(
       primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
     )
   }
+  # Extract Gamma distribution parameters
+  shape <- object$args$shape
+  scale <- object$args$scale
+  # if we don't have scale get fromm rate
+  if (is.null(scale)) {
+    scale <- 1 / object$args$rate
+  }
+  if (is.null(shape)) {
+    stop("shape parameter is required for Gamma distribution")
+  }
+  if (is.null(scale)) {
+    stop("scale or rate parameter is required for Gamma distribution")
+  }
 
-  result <- vapply(q, function(n) {
-    # Implement analytical solution here
-  }, numeric(1))
+  partial_pgamma <- function(q) {
+    pgamma(q, shape = shape, scale = scale)
+  }
+
+  # Handle negative q values safely
+  result <- ifelse(q < 0, 0, NA)
+
+  valid_q <- q[q >= 0]
+
+  if (length(valid_q) > 0) {
+    # Compute necessary survival and distribution functions
+    pgamma_q <- partial_pgamma(valid_q)
+    pgamma_q_pwindow <- partial_pgamma(valid_q + pwindow)
+    pgamma_q_pwindow_1 <- partial_pgamma(valid_q + pwindow + 1)
+    pgamma_q_1 <- partial_pgamma(valid_q + 1)
+
+    Q_T <- 1 - pgamma_q_pwindow
+    Delta_F_T_kp1 <- pgamma_q_pwindow_1 - pgamma_q_1
+    Delta_F_T_k <- pgamma_q_pwindow - pgamma_q
+
+    # Calculate Q_Splus using the analytical formula
+    Q_Splus <- Q_T +
+      (shape * scale / pwindow) * Delta_F_T_kp1 -
+      (valid_q / pwindow) * Delta_F_T_k
+
+    # Compute the CDF as 1 - Q_Splus
+    result[q >= 0] <- 1 - Q_Splus
+  }
 
   return(result)
 }
@@ -147,16 +184,54 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
 #' @export
 primary_censored_cdf.pcens_plnorm_dunif <- function(
     object, q, pwindow, use_numeric = FALSE) {
-  use_numeric <- TRUE
   if (isTRUE(use_numeric)) {
     return(
       primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
     )
   }
 
-  result <- vapply(q, function(n) {
-    # Implement analytical solution here
-  }, numeric(1))
+  # Extract Log-Normal distribution parameters
+  mu <- object$args$meanlog
+  sigma <- object$args$sdlog
+  if (is.null(mu)) {
+    stop("meanlog parameter is required for Log-Normal distribution")
+  }
+  if (is.null(sigma)) {
+    stop("sdlog parameter is required for Log-Normal distribution")
+  }
+
+  partial_plnorm <- function(q) {
+    plnorm(q, meanlog = mu, sdlog = sigma)
+  }
+
+  # Handle negative q values safely
+  result <- ifelse(q < 0, 0, NA)
+
+  valid_q <- q[q >= 0]
+
+  if (length(valid_q) > 0) {
+    # Compute necessary survival and distribution functions
+    plnorm_q <- partial_plnorm(valid_q)
+    plnorm_q_pwindow <- partial_plnorm(valid_q + pwindow)
+    plnorm_q_pwindow_sigma2 <- partial_plnorm(
+      valid_q + pwindow + sigma^2
+    )
+    plnorm_q_sigma2 <- partial_plnorm(valid_q + sigma^2)
+
+    # Compute necessary survival and distribution functions
+    Q_T <- 1 - plnorm_q_pwindow
+    Delta_F_T_mu_sigma <- plnorm_q_pwindow_sigma2 -
+      plnorm_q_sigma2
+    Delta_F_T <- plnorm_q_pwindow - plnorm_q
+
+    # Calculate Q_Splus using the analytical formula
+    Q_Splus <- Q_T +
+      (exp(mu + 0.5 * sigma^2) / pwindow) * Delta_F_T_mu_sigma -
+      (valid_q / pwindow) * Delta_F_T
+
+    # Compute the CDF as 1 - Q_Splus
+    result[q >= 0] <- 1 - Q_Splus
+  }
 
   return(result)
 }
diff --git a/man/dot-dpcens.Rd b/man/dot-dpcens.Rd
index 7d4fc60..135baaf 100644
--- a/man/dot-dpcens.Rd
+++ b/man/dot-dpcens.Rd
@@ -4,7 +4,18 @@
 \alias{.dpcens}
 \title{Define a fitdistrplus compatible wrapper around dprimarycensoreddist}
 \usage{
-.dpcens(x, swindows, pdist, pwindow, D, dprimary, dprimary_args, ...)
+.dpcens(
+  x,
+  swindows,
+  pdist,
+  pwindow,
+  D,
+  dprimary,
+  dprimary_args,
+  pdist_name,
+  dprimary_name,
+  ...
+)
 }
 \arguments{
 \item{x}{Vector of quantiles}
@@ -32,6 +43,18 @@ dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
 and rate parameters}
 
+\item{pdist_name}{A string specifying the name of the delay distribution
+function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{pdist} is
+passed a pre-assigned variable rather than a function name.}
+
+\item{dprimary_name}{A string specifying the name of the primary event
+distribution function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{dprimary} is
+passed a pre-assigned variable rather than a function name.}
+
 \item{...}{Additional arguments to be passed to the distribution function}
 }
 \description{
diff --git a/man/dot-ppcens.Rd b/man/dot-ppcens.Rd
index 3988afc..65f5428 100644
--- a/man/dot-ppcens.Rd
+++ b/man/dot-ppcens.Rd
@@ -4,7 +4,17 @@
 \alias{.ppcens}
 \title{Define a fitdistrplus compatible wrapper around pprimarycensoreddist}
 \usage{
-.ppcens(q, pdist, pwindow, D, dprimary, dprimary_args, ...)
+.ppcens(
+  q,
+  pdist,
+  pwindow,
+  D,
+  dprimary,
+  dprimary_args,
+  pdist_name,
+  dprimary_name,
+  ...
+)
 }
 \arguments{
 \item{q}{Vector of quantiles}
@@ -29,6 +39,18 @@ dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,
 and rate parameters}
 
+\item{pdist_name}{A string specifying the name of the delay distribution
+function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{pdist} is
+passed a pre-assigned variable rather than a function name.}
+
+\item{dprimary_name}{A string specifying the name of the primary event
+distribution function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{dprimary} is
+passed a pre-assigned variable rather than a function name.}
+
 \item{...}{Additional arguments to be passed to pdist}
 }
 \description{
diff --git a/man/fitdistdoublecens.Rd b/man/fitdistdoublecens.Rd
index 7e3f63a..48dad8a 100644
--- a/man/fitdistdoublecens.Rd
+++ b/man/fitdistdoublecens.Rd
@@ -10,6 +10,7 @@ fitdistdoublecens(
   pwindow = 1,
   D = Inf,
   dprimary = stats::dunif,
+  dprimary_name = NULL,
   dprimary_args = list(),
   ...
 )
@@ -35,6 +36,12 @@ distribution over [0, pwindow]. Users can provide custom functions or use
 helper functions like \code{dexpgrowth} for an exponential growth distribution.
 See \code{primary_dists.R} for examples.}
 
+\item{dprimary_name}{A string specifying the name of the primary event
+distribution function. If NULL, the function name is extracted using
+\code{\link[=.extract_function_name]{.extract_function_name()}}. Used to determine if a analytical solution
+exists for the primary censored distribution. Must be set if \code{dprimary} is
+passed a pre-assigned variable rather than a function name.}
+
 \item{dprimary_args}{List of additional arguments to be passed to
 dprimary. For example, when using \code{dexpgrowth}, you would
 pass \code{list(min = 0, max = pwindow, r = 0.2)} to set the minimum, maximum,

From 1e89436f496d49a5f0ad74880368497f6ba07a50 Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 18:07:56 +0100
Subject: [PATCH 02/15] get rid of pcens_numeric in favour of just using the
 default

---
 NAMESPACE                                     |  1 -
 R/primary_censored_dist.R                     | 21 ++-------
 man/new_primary_censored_dist.Rd              |  1 -
 man/primary_censored_cdf.Rd                   |  1 -
 man/primary_censored_cdf.default.Rd           | 10 ++++-
 man/primary_censored_cdf.pcens_numeric.Rd     | 43 -------------------
 ...primary_censored_cdf.pcens_pgamma_dunif.Rd |  1 -
 ...primary_censored_cdf.pcens_plnorm_dunif.Rd |  1 -
 8 files changed, 12 insertions(+), 67 deletions(-)
 delete mode 100644 man/primary_censored_cdf.pcens_numeric.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 989308e..93b5927 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(primary_censored_cdf,default)
-S3method(primary_censored_cdf,pcens_numeric)
 S3method(primary_censored_cdf,pcens_pgamma_dunif)
 S3method(primary_censored_cdf,pcens_plnorm_dunif)
 export(check_dprimary)
diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index 6250cac..7c036ee 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -64,21 +64,6 @@ primary_censored_cdf <- function(
 #' the numeric integration method.
 #'
 #' @inheritParams primary_censored_cdf
-#'
-#' @family primary_censored_dist
-#'
-#' @export
-primary_censored_cdf.default <- function(
-    object, q, pwindow, use_numeric = FALSE) {
-  primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
-}
-
-#' Numeric method for computing primary event censored CDF
-#'
-#' This method uses numerical integration to compute the primary event censored
-#' CDF for any combination of delay distribution and primary event distribution.
-#'
-#' @inheritParams primary_censored_cdf
 #' @inheritParams pprimarycensoreddist
 #'
 #' @details
@@ -93,7 +78,7 @@ primary_censored_cdf.default <- function(
 #' @family primary_censored_dist
 #'
 #' @export
-primary_censored_cdf.pcens_numeric <- function(
+primary_censored_cdf.default <- function(
     object, q, pwindow, use_numeric = FALSE) {
   result <- vapply(q, function(d) {
     if (d <= 0) {
@@ -126,7 +111,7 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
     object, q, pwindow, use_numeric = FALSE) {
   if (isTRUE(use_numeric)) {
     return(
-      primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+      primary_censored_cdf.default(object, q, pwindow, use_numeric)
     )
   }
   # Extract Gamma distribution parameters
@@ -186,7 +171,7 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
     object, q, pwindow, use_numeric = FALSE) {
   if (isTRUE(use_numeric)) {
     return(
-      primary_censored_cdf.pcens_numeric(object, q, pwindow, use_numeric)
+      primary_censored_cdf.default(object, q, pwindow, use_numeric)
     )
   }
 
diff --git a/man/new_primary_censored_dist.Rd b/man/new_primary_censored_dist.Rd
index 7c68b33..b207dab 100644
--- a/man/new_primary_censored_dist.Rd
+++ b/man/new_primary_censored_dist.Rd
@@ -53,7 +53,6 @@ S3 class for primary event censored distribution computation
 Low level primary event censored distribution objects and methods
 \code{\link{primary_censored_cdf}()},
 \code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
 \code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
 }
diff --git a/man/primary_censored_cdf.Rd b/man/primary_censored_cdf.Rd
index 4dd2d7d..002f81a 100644
--- a/man/primary_censored_cdf.Rd
+++ b/man/primary_censored_cdf.Rd
@@ -29,7 +29,6 @@ Compute primary event censored CDF
 Low level primary event censored distribution objects and methods
 \code{\link{new_primary_censored_dist}()},
 \code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
 \code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
 }
diff --git a/man/primary_censored_cdf.default.Rd b/man/primary_censored_cdf.default.Rd
index 2359b20..648c05e 100644
--- a/man/primary_censored_cdf.default.Rd
+++ b/man/primary_censored_cdf.default.Rd
@@ -24,11 +24,19 @@ This method serves as a fallback for combinations of delay and primary
 event distributions that don't have specific implementations. It uses
 the numeric integration method.
 }
+\details{
+This method implements the numerical integration approach for computing
+the primary event censored CDF. It uses the same mathematical formulation
+as described in the details section of \code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}}, but
+applies numerical integration instead of analytical solutions.
+}
 \seealso{
+\code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}} for the mathematical details of the
+primary event censored CDF computation.
+
 Low level primary event censored distribution objects and methods
 \code{\link{new_primary_censored_dist}()},
 \code{\link{primary_censored_cdf}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
 \code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
 }
diff --git a/man/primary_censored_cdf.pcens_numeric.Rd b/man/primary_censored_cdf.pcens_numeric.Rd
deleted file mode 100644
index b43974f..0000000
--- a/man/primary_censored_cdf.pcens_numeric.Rd
+++ /dev/null
@@ -1,43 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/primary_censored_dist.R
-\name{primary_censored_cdf.pcens_numeric}
-\alias{primary_censored_cdf.pcens_numeric}
-\title{Numeric method for computing primary event censored CDF}
-\usage{
-\method{primary_censored_cdf}{pcens_numeric}(object, q, pwindow, use_numeric = FALSE)
-}
-\arguments{
-\item{object}{A \code{primary_censored_dist} object as created by
-\code{\link[=new_primary_censored_dist]{new_primary_censored_dist()}}.}
-
-\item{q}{Vector of quantiles}
-
-\item{pwindow}{Primary event window}
-
-\item{use_numeric}{Logical, if TRUE forces use of numeric integration
-even for distributions with analytical solutions. This is primarily
-useful for testing purposes or for settings where the analytical solution
-breaks down.}
-}
-\description{
-This method uses numerical integration to compute the primary event censored
-CDF for any combination of delay distribution and primary event distribution.
-}
-\details{
-This method implements the numerical integration approach for computing
-the primary event censored CDF. It uses the same mathematical formulation
-as described in the details section of \code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}}, but
-applies numerical integration instead of analytical solutions.
-}
-\seealso{
-\code{\link[=pprimarycensoreddist]{pprimarycensoreddist()}} for the mathematical details of the
-primary event censored CDF computation.
-
-Low level primary event censored distribution objects and methods
-\code{\link{new_primary_censored_dist}()},
-\code{\link{primary_censored_cdf}()},
-\code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_pgamma_dunif}()},
-\code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
-}
-\concept{primary_censored_dist}
diff --git a/man/primary_censored_cdf.pcens_pgamma_dunif.Rd b/man/primary_censored_cdf.pcens_pgamma_dunif.Rd
index 1ee894f..c333bd1 100644
--- a/man/primary_censored_cdf.pcens_pgamma_dunif.Rd
+++ b/man/primary_censored_cdf.pcens_pgamma_dunif.Rd
@@ -27,7 +27,6 @@ Low level primary event censored distribution objects and methods
 \code{\link{new_primary_censored_dist}()},
 \code{\link{primary_censored_cdf}()},
 \code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_plnorm_dunif}()}
 }
 \concept{primary_censored_dist}
diff --git a/man/primary_censored_cdf.pcens_plnorm_dunif.Rd b/man/primary_censored_cdf.pcens_plnorm_dunif.Rd
index 387070f..e041be4 100644
--- a/man/primary_censored_cdf.pcens_plnorm_dunif.Rd
+++ b/man/primary_censored_cdf.pcens_plnorm_dunif.Rd
@@ -27,7 +27,6 @@ Low level primary event censored distribution objects and methods
 \code{\link{new_primary_censored_dist}()},
 \code{\link{primary_censored_cdf}()},
 \code{\link{primary_censored_cdf.default}()},
-\code{\link{primary_censored_cdf.pcens_numeric}()},
 \code{\link{primary_censored_cdf.pcens_pgamma_dunif}()}
 }
 \concept{primary_censored_dist}

From e3b7e8716d6f33072fe8b728f5debe2a335e192e Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 18:14:08 +0100
Subject: [PATCH 03/15] human read llm and help

---
 R/primary_censored_dist.R | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index 7c036ee..efcf7a4 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -131,6 +131,9 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
   partial_pgamma <- function(q) {
     pgamma(q, shape = shape, scale = scale)
   }
+  partial_pgamm_k_1 <- function(q) {
+    pgamma(q, shape = shape + 1, scale = scale)
+  }
 
   # Handle negative q values safely
   result <- ifelse(q < 0, 0, NA)
@@ -141,8 +144,8 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
     # Compute necessary survival and distribution functions
     pgamma_q <- partial_pgamma(valid_q)
     pgamma_q_pwindow <- partial_pgamma(valid_q + pwindow)
-    pgamma_q_pwindow_1 <- partial_pgamma(valid_q + pwindow + 1)
-    pgamma_q_1 <- partial_pgamma(valid_q + 1)
+    pgamma_q_pwindow_1 <- partial_pgamm_k_1(valid_q + pwindow)
+    pgamma_q_1 <- partial_pgamm_k_1(valid_q + 1)
 
     Q_T <- 1 - pgamma_q_pwindow
     Delta_F_T_kp1 <- pgamma_q_pwindow_1 - pgamma_q_1
@@ -188,7 +191,9 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
   partial_plnorm <- function(q) {
     plnorm(q, meanlog = mu, sdlog = sigma)
   }
-
+  partial_plnorm_sigma2 <- function(q) {
+    plnorm(q, meanlog = mu + sigma^2, sdlog = sigma)
+  }
   # Handle negative q values safely
   result <- ifelse(q < 0, 0, NA)
 
@@ -198,10 +203,8 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
     # Compute necessary survival and distribution functions
     plnorm_q <- partial_plnorm(valid_q)
     plnorm_q_pwindow <- partial_plnorm(valid_q + pwindow)
-    plnorm_q_pwindow_sigma2 <- partial_plnorm(
-      valid_q + pwindow + sigma^2
-    )
-    plnorm_q_sigma2 <- partial_plnorm(valid_q + sigma^2)
+    plnorm_q_pwindow_sigma2 <- partial_plnorm_sigma2(valid_q + pwindow)
+    plnorm_q_sigma2 <- partial_plnorm_sigma2(valid_q)
 
     # Compute necessary survival and distribution functions
     Q_T <- 1 - plnorm_q_pwindow

From 39252aa808579780d38fd9a119ed1f0f3fbb42ad Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 21:17:54 +0100
Subject: [PATCH 04/15] catch random 1

---
 R/primary_censored_dist.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index efcf7a4..4baff9b 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -144,8 +144,8 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
     # Compute necessary survival and distribution functions
     pgamma_q <- partial_pgamma(valid_q)
     pgamma_q_pwindow <- partial_pgamma(valid_q + pwindow)
+    pgamma_q_1 <- partial_pgamm_k_1(valid_q)
     pgamma_q_pwindow_1 <- partial_pgamm_k_1(valid_q + pwindow)
-    pgamma_q_1 <- partial_pgamm_k_1(valid_q + 1)
 
     Q_T <- 1 - pgamma_q_pwindow
     Delta_F_T_kp1 <- pgamma_q_pwindow_1 - pgamma_q_1
@@ -203,8 +203,8 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
     # Compute necessary survival and distribution functions
     plnorm_q <- partial_plnorm(valid_q)
     plnorm_q_pwindow <- partial_plnorm(valid_q + pwindow)
-    plnorm_q_pwindow_sigma2 <- partial_plnorm_sigma2(valid_q + pwindow)
     plnorm_q_sigma2 <- partial_plnorm_sigma2(valid_q)
+    plnorm_q_pwindow_sigma2 <- partial_plnorm_sigma2(valid_q + pwindow)
 
     # Compute necessary survival and distribution functions
     Q_T <- 1 - plnorm_q_pwindow

From c6ae9fac260cfacf198d6b8ce8b58d075fd6b608 Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 21:41:23 +0100
Subject: [PATCH 05/15] think about inequality

---
 R/primary_censored_dist.R                   | 12 ++++++------
 tests/testthat/test-primary_censored_dist.R |  6 +++---
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index 4baff9b..d667e4f 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -136,9 +136,9 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
   }
 
   # Handle negative q values safely
-  result <- ifelse(q < 0, 0, NA)
+  result <- ifelse(q <= 0, 0, NA)
 
-  valid_q <- q[q >= 0]
+  valid_q <- q[q > 0]
 
   if (length(valid_q) > 0) {
     # Compute necessary survival and distribution functions
@@ -157,7 +157,7 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
       (valid_q / pwindow) * Delta_F_T_k
 
     # Compute the CDF as 1 - Q_Splus
-    result[q >= 0] <- 1 - Q_Splus
+    result[q > 0] <- 1 - Q_Splus
   }
 
   return(result)
@@ -195,9 +195,9 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
     plnorm(q, meanlog = mu + sigma^2, sdlog = sigma)
   }
   # Handle negative q values safely
-  result <- ifelse(q < 0, 0, NA)
+  result <- ifelse(q <= 0, 0, NA)
 
-  valid_q <- q[q >= 0]
+  valid_q <- q[q > 0]
 
   if (length(valid_q) > 0) {
     # Compute necessary survival and distribution functions
@@ -218,7 +218,7 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
       (valid_q / pwindow) * Delta_F_T
 
     # Compute the CDF as 1 - Q_Splus
-    result[q >= 0] <- 1 - Q_Splus
+    result[q > 0] <- 1 - Q_Splus
   }
 
   return(result)
diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R
index 921a14b..1006f78 100644
--- a/tests/testthat/test-primary_censored_dist.R
+++ b/tests/testthat/test-primary_censored_dist.R
@@ -40,10 +40,10 @@ test_that("primary_censored_cdf.pcens_numeric computes correct values", {
     shape = shape, rate = rate
   )
 
-  q_values <- c(0, 5, 10, 15)
-  pwindow <- 10
+  q_values <- c(0, 3, 5, 7, 10, 15)
+  pwindow <- 1
 
-  result <- primary_censored_cdf(obj, q = q_values, pwindow = pwindow)
+  result <- primary_censored_cdf(obj, q = q_values, pwindow = pwindow, use_numeric = TRUE)
 
   expect_true(all(result >= 0))
   expect_true(all(result <= 1))

From ae7a8f2b33399931de897ff6a638763d194fa408 Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 22:18:36 +0100
Subject: [PATCH 06/15] shift indexing

---
 R/primary_censored_dist.R                   | 14 ++++++++------
 tests/testthat/test-primary_censored_dist.R |  4 ++--
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index d667e4f..658ac64 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -136,9 +136,9 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
   }
 
   # Handle negative q values safely
-  result <- ifelse(q <= 0, 0, NA)
+  result <- ifelse(q < 0, 0, NA)
 
-  valid_q <- q[q > 0]
+  valid_q <- q[q >= 0]
 
   if (length(valid_q) > 0) {
     # Compute necessary survival and distribution functions
@@ -157,7 +157,7 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
       (valid_q / pwindow) * Delta_F_T_k
 
     # Compute the CDF as 1 - Q_Splus
-    result[q > 0] <- 1 - Q_Splus
+    result[q >= 0] <- 1 - Q_Splus
   }
 
   return(result)
@@ -194,10 +194,12 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
   partial_plnorm_sigma2 <- function(q) {
     plnorm(q, meanlog = mu + sigma^2, sdlog = sigma)
   }
+  # Adjust q so that we have [q-pwindow, q]
+  q <- q - pwindow
   # Handle negative q values safely
-  result <- ifelse(q <= 0, 0, NA)
+  result <- ifelse(q < 0, 0, NA)
 
-  valid_q <- q[q > 0]
+  valid_q <- q[q >= 0]
 
   if (length(valid_q) > 0) {
     # Compute necessary survival and distribution functions
@@ -218,7 +220,7 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
       (valid_q / pwindow) * Delta_F_T
 
     # Compute the CDF as 1 - Q_Splus
-    result[q > 0] <- 1 - Q_Splus
+    result[q >= 0] <- 1 - Q_Splus
   }
 
   return(result)
diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R
index 1006f78..627d120 100644
--- a/tests/testthat/test-primary_censored_dist.R
+++ b/tests/testthat/test-primary_censored_dist.R
@@ -31,7 +31,7 @@ test_that("primary_censored_cdf.pcens_numeric computes correct values", {
   dprimary_name <- "dunif"
   dprimary <- dunif
   shape <- 2
-  rate <- 1
+  rate <- 0.5
 
   obj <- new_primary_censored_dist(
     pdist,
@@ -40,7 +40,7 @@ test_that("primary_censored_cdf.pcens_numeric computes correct values", {
     shape = shape, rate = rate
   )
 
-  q_values <- c(0, 3, 5, 7, 10, 15)
+  q_values <- 0:20
   pwindow <- 1
 
   result <- primary_censored_cdf(obj, q = q_values, pwindow = pwindow, use_numeric = TRUE)

From 3b6cdda90fca09db55fe66b0fa74f69115304688 Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 22:33:23 +0100
Subject: [PATCH 07/15] complete index rejig

---
 R/primary_censored_dist.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index 658ac64..5cf3f67 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -134,7 +134,8 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
   partial_pgamm_k_1 <- function(q) {
     pgamma(q, shape = shape + 1, scale = scale)
   }
-
+  # Adjust q so that we have [q-pwindow, q]
+  q <- q - pwindow
   # Handle negative q values safely
   result <- ifelse(q < 0, 0, NA)
 

From 7445409cea342fa5cc4545348d2962b0cd0e7dd3 Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 22:47:43 +0100
Subject: [PATCH 08/15] add partial window support

---
 R/primary_censored_dist.R | 72 ++++++++++++++++++++++++++-------------
 1 file changed, 49 insertions(+), 23 deletions(-)

diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index 5cf3f67..2f20643 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -136,17 +136,27 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
   }
   # Adjust q so that we have [q-pwindow, q]
   q <- q - pwindow
-  # Handle negative q values safely
-  result <- ifelse(q < 0, 0, NA)
+  # Adjust q and pwindow for cases crossing zero
+  adjusted_q <- pmax(q, 0)
+  adjusted_pwindow <- pmin(pwindow, pmax(0, q + pwindow) - adjusted_q)
 
-  valid_q <- q[q >= 0]
+  # Initialize result vector
+  result <- numeric(length(q))
+
+  # Handle cases where q + pwindow <= 0
+  zero_cases <- q + pwindow <= 0
+  result[zero_cases] <- 0
+
+  # Process non-zero cases only if there are any
+  if (any(!zero_cases)) {
+    non_zero_q <- adjusted_q[!zero_cases]
+    non_zero_pwindow <- adjusted_pwindow[!zero_cases]
 
-  if (length(valid_q) > 0) {
     # Compute necessary survival and distribution functions
-    pgamma_q <- partial_pgamma(valid_q)
-    pgamma_q_pwindow <- partial_pgamma(valid_q + pwindow)
-    pgamma_q_1 <- partial_pgamm_k_1(valid_q)
-    pgamma_q_pwindow_1 <- partial_pgamm_k_1(valid_q + pwindow)
+    pgamma_q <- partial_pgamma(non_zero_q)
+    pgamma_q_pwindow <- partial_pgamma(non_zero_q + non_zero_pwindow)
+    pgamma_q_1 <- partial_pgamm_k_1(non_zero_q)
+    pgamma_q_pwindow_1 <- partial_pgamm_k_1(non_zero_q + non_zero_pwindow)
 
     Q_T <- 1 - pgamma_q_pwindow
     Delta_F_T_kp1 <- pgamma_q_pwindow_1 - pgamma_q_1
@@ -155,10 +165,13 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
     # Calculate Q_Splus using the analytical formula
     Q_Splus <- Q_T +
       (shape * scale / pwindow) * Delta_F_T_kp1 -
-      (valid_q / pwindow) * Delta_F_T_k
+      (non_zero_q / pwindow) * Delta_F_T_k
 
     # Compute the CDF as 1 - Q_Splus
-    result[q >= 0] <- 1 - Q_Splus
+    non_zero_result <- 1 - Q_Splus
+
+    # Assign non-zero results back to the main result vector
+    result[!zero_cases] <- non_zero_result
   }
 
   return(result)
@@ -197,31 +210,44 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
   }
   # Adjust q so that we have [q-pwindow, q]
   q <- q - pwindow
-  # Handle negative q values safely
-  result <- ifelse(q < 0, 0, NA)
+  # Adjust q and pwindow for cases crossing zero
+  adjusted_q <- pmax(q, 0)
+  adjusted_pwindow <- pmin(pwindow, pmax(0, q + pwindow) - adjusted_q)
 
-  valid_q <- q[q >= 0]
+  # Initialize result vector
+  result <- numeric(length(q))
 
-  if (length(valid_q) > 0) {
-    # Compute necessary survival and distribution functions
-    plnorm_q <- partial_plnorm(valid_q)
-    plnorm_q_pwindow <- partial_plnorm(valid_q + pwindow)
-    plnorm_q_sigma2 <- partial_plnorm_sigma2(valid_q)
-    plnorm_q_pwindow_sigma2 <- partial_plnorm_sigma2(valid_q + pwindow)
+  # Handle cases where q + pwindow <= 0
+  zero_cases <- q + pwindow <= 0
+  result[zero_cases] <- 0
+
+  # Process non-zero cases only if there are any
+  if (any(!zero_cases)) {
+    non_zero_q <- adjusted_q[!zero_cases]
+    non_zero_pwindow <- adjusted_pwindow[!zero_cases]
 
     # Compute necessary survival and distribution functions
+    plnorm_q <- partial_plnorm(non_zero_q)
+    plnorm_q_pwindow <- partial_plnorm(non_zero_q + non_zero_pwindow)
+    plnorm_q_sigma2 <- partial_plnorm_sigma2(non_zero_q)
+    plnorm_q_pwindow_sigma2 <- partial_plnorm_sigma2(
+      non_zero_q + non_zero_pwindow
+    )
+
     Q_T <- 1 - plnorm_q_pwindow
-    Delta_F_T_mu_sigma <- plnorm_q_pwindow_sigma2 -
-      plnorm_q_sigma2
+    Delta_F_T_mu_sigma <- plnorm_q_pwindow_sigma2 - plnorm_q_sigma2
     Delta_F_T <- plnorm_q_pwindow - plnorm_q
 
     # Calculate Q_Splus using the analytical formula
     Q_Splus <- Q_T +
       (exp(mu + 0.5 * sigma^2) / pwindow) * Delta_F_T_mu_sigma -
-      (valid_q / pwindow) * Delta_F_T
+      (non_zero_q / pwindow) * Delta_F_T
 
     # Compute the CDF as 1 - Q_Splus
-    result[q >= 0] <- 1 - Q_Splus
+    non_zero_result <- 1 - Q_Splus
+
+    # Assign non-zero results back to the main result vector
+    result[!zero_cases] <- non_zero_result
   }
 
   return(result)

From d2900c6e0d25e27a205794081e8b3042e64ec30f Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 22:55:54 +0100
Subject: [PATCH 09/15] catch outstanding non-adjusted window

---
 R/primary_censored_dist.R                   | 8 ++++----
 tests/testthat/test-primary_censored_dist.R | 2 +-
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index 2f20643..b42e6b3 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -164,8 +164,8 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
 
     # Calculate Q_Splus using the analytical formula
     Q_Splus <- Q_T +
-      (shape * scale / pwindow) * Delta_F_T_kp1 -
-      (non_zero_q / pwindow) * Delta_F_T_k
+      (shape * scale / non_zero_pwindow) * Delta_F_T_kp1 -
+      (non_zero_q / non_zero_pwindow) * Delta_F_T_k
 
     # Compute the CDF as 1 - Q_Splus
     non_zero_result <- 1 - Q_Splus
@@ -240,8 +240,8 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
 
     # Calculate Q_Splus using the analytical formula
     Q_Splus <- Q_T +
-      (exp(mu + 0.5 * sigma^2) / pwindow) * Delta_F_T_mu_sigma -
-      (non_zero_q / pwindow) * Delta_F_T
+      (exp(mu + 0.5 * sigma^2) / non_zero_pwindow) * Delta_F_T_mu_sigma -
+      (non_zero_q / non_zero_pwindow) * Delta_F_T
 
     # Compute the CDF as 1 - Q_Splus
     non_zero_result <- 1 - Q_Splus
diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R
index 627d120..0d5f348 100644
--- a/tests/testthat/test-primary_censored_dist.R
+++ b/tests/testthat/test-primary_censored_dist.R
@@ -41,7 +41,7 @@ test_that("primary_censored_cdf.pcens_numeric computes correct values", {
   )
 
   q_values <- 0:20
-  pwindow <- 1
+  pwindow <- 2
 
   result <- primary_censored_cdf(obj, q = q_values, pwindow = pwindow, use_numeric = TRUE)
 

From fe03c0fc14587de9b65fa0cec9daec2f630ae974 Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Thu, 19 Sep 2024 23:13:34 +0100
Subject: [PATCH 10/15] tweak tests:

---
 tests/testthat/test-primary_censored_dist.R | 24 +++++++++++++++------
 1 file changed, 17 insertions(+), 7 deletions(-)

diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R
index 0d5f348..c3d2e89 100644
--- a/tests/testthat/test-primary_censored_dist.R
+++ b/tests/testthat/test-primary_censored_dist.R
@@ -25,7 +25,7 @@ test_that("new_primary_censored_dist creates object with correct structure", {
   expect_identical(obj, new_obj)
 })
 
-test_that("primary_censored_cdf.pcens_numeric computes correct values", {
+test_that("primary_censored_cdf.default computes correct values", {
   pdist_name <- "pgamma"
   pdist <- pgamma
   dprimary_name <- "dunif"
@@ -42,14 +42,24 @@ test_that("primary_censored_cdf.pcens_numeric computes correct values", {
 
   q_values <- 0:20
   pwindow <- 2
+  result_numeric <- primary_censored_cdf(
+    obj,
+    q = q_values, pwindow = pwindow, use_numeric = TRUE
+  )
+  result_analytical <- primary_censored_cdf(
+    obj,
+    q = q_values, pwindow = pwindow, use_numeric = FALSE
+  )
 
-  result <- primary_censored_cdf(obj, q = q_values, pwindow = pwindow, use_numeric = TRUE)
+  # Check properties of numeric result
+  expect_true(all(result_numeric >= 0))
+  expect_true(all(result_numeric <= 1))
+  expect_type(result_numeric, "double")
+  expect_length(result_numeric, length(q_values))
+  expect_true(all(diff(result_numeric) >= 0)) # Ensure CDF is non-decreasing
 
-  expect_true(all(result >= 0))
-  expect_true(all(result <= 1))
-  expect_type(result, "double")
-  expect_length(result, length(q_values))
-  expect_true(all(diff(result) >= 0)) # Ensure CDF is non-decreasing
+  # Check that analytical and numeric results are the same
+  expect_equal(result_numeric, result_analytical, tolerance = 1e-6)
 })
 
 test_that("primary_censored_cdf methods dispatch correctly", {

From fa619f99e52ccd8fe5bbaaf98052974ba9c85a0b Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Fri, 20 Sep 2024 12:27:31 +0100
Subject: [PATCH 11/15] fix partial windows

---
 R/primary_censored_dist.R                   | 43 +++++++--------------
 tests/testthat/test-primary_censored_dist.R | 37 ++++++++++++++++++
 2 files changed, 51 insertions(+), 29 deletions(-)

diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index b42e6b3..97746d0 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -136,27 +136,19 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
   }
   # Adjust q so that we have [q-pwindow, q]
   q <- q - pwindow
-  # Adjust q and pwindow for cases crossing zero
-  adjusted_q <- pmax(q, 0)
-  adjusted_pwindow <- pmin(pwindow, pmax(0, q + pwindow) - adjusted_q)
-
-  # Initialize result vector
-  result <- numeric(length(q))
-
   # Handle cases where q + pwindow <= 0
   zero_cases <- q + pwindow <= 0
-  result[zero_cases] <- 0
+  result <- ifelse(zero_cases, 0, NA)
 
   # Process non-zero cases only if there are any
-  if (any(!zero_cases)) {
-    non_zero_q <- adjusted_q[!zero_cases]
-    non_zero_pwindow <- adjusted_pwindow[!zero_cases]
+  if (!all(zero_cases)) {
+    non_zero_q <- q[!zero_cases]
 
     # Compute necessary survival and distribution functions
     pgamma_q <- partial_pgamma(non_zero_q)
-    pgamma_q_pwindow <- partial_pgamma(non_zero_q + non_zero_pwindow)
+    pgamma_q_pwindow <- partial_pgamma(non_zero_q + pwindow)
     pgamma_q_1 <- partial_pgamm_k_1(non_zero_q)
-    pgamma_q_pwindow_1 <- partial_pgamm_k_1(non_zero_q + non_zero_pwindow)
+    pgamma_q_pwindow_1 <- partial_pgamm_k_1(non_zero_q + pwindow)
 
     Q_T <- 1 - pgamma_q_pwindow
     Delta_F_T_kp1 <- pgamma_q_pwindow_1 - pgamma_q_1
@@ -164,8 +156,8 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
 
     # Calculate Q_Splus using the analytical formula
     Q_Splus <- Q_T +
-      (shape * scale / non_zero_pwindow) * Delta_F_T_kp1 -
-      (non_zero_q / non_zero_pwindow) * Delta_F_T_k
+      (shape * scale / pwindow) * Delta_F_T_kp1 -
+      (non_zero_q / pwindow) * Delta_F_T_k
 
     # Compute the CDF as 1 - Q_Splus
     non_zero_result <- 1 - Q_Splus
@@ -210,28 +202,21 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
   }
   # Adjust q so that we have [q-pwindow, q]
   q <- q - pwindow
-  # Adjust q and pwindow for cases crossing zero
-  adjusted_q <- pmax(q, 0)
-  adjusted_pwindow <- pmin(pwindow, pmax(0, q + pwindow) - adjusted_q)
-
-  # Initialize result vector
-  result <- numeric(length(q))
 
   # Handle cases where q + pwindow <= 0
   zero_cases <- q + pwindow <= 0
-  result[zero_cases] <- 0
+  result <- ifelse(zero_cases, 0, NA)
 
   # Process non-zero cases only if there are any
-  if (any(!zero_cases)) {
-    non_zero_q <- adjusted_q[!zero_cases]
-    non_zero_pwindow <- adjusted_pwindow[!zero_cases]
+  if (!all(zero_cases)) {
+    non_zero_q <- q[!zero_cases]
 
     # Compute necessary survival and distribution functions
     plnorm_q <- partial_plnorm(non_zero_q)
-    plnorm_q_pwindow <- partial_plnorm(non_zero_q + non_zero_pwindow)
+    plnorm_q_pwindow <- partial_plnorm(non_zero_q + pwindow)
     plnorm_q_sigma2 <- partial_plnorm_sigma2(non_zero_q)
     plnorm_q_pwindow_sigma2 <- partial_plnorm_sigma2(
-      non_zero_q + non_zero_pwindow
+      non_zero_q + pwindow
     )
 
     Q_T <- 1 - plnorm_q_pwindow
@@ -240,8 +225,8 @@ primary_censored_cdf.pcens_plnorm_dunif <- function(
 
     # Calculate Q_Splus using the analytical formula
     Q_Splus <- Q_T +
-      (exp(mu + 0.5 * sigma^2) / non_zero_pwindow) * Delta_F_T_mu_sigma -
-      (non_zero_q / non_zero_pwindow) * Delta_F_T
+      (exp(mu + 0.5 * sigma^2) / pwindow) * Delta_F_T_mu_sigma -
+      (non_zero_q / pwindow) * Delta_F_T
 
     # Compute the CDF as 1 - Q_Splus
     non_zero_result <- 1 - Q_Splus
diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R
index c3d2e89..e6efca2 100644
--- a/tests/testthat/test-primary_censored_dist.R
+++ b/tests/testthat/test-primary_censored_dist.R
@@ -62,6 +62,43 @@ test_that("primary_censored_cdf.default computes correct values", {
   expect_equal(result_numeric, result_analytical, tolerance = 1e-6)
 })
 
+test_that("primary_censored_cdf.pcens_plnorm_dunif uses numeric method", {
+  pdist_name <- "plnorm"
+  pdist <- plnorm
+  dprimary_name <- "dunif"
+  dprimary <- dunif
+  meanlog <- 0
+  sdlog <- 1
+
+  obj <- new_primary_censored_dist(
+    pdist,
+    dprimary, list(),
+    pdist_name, dprimary_name,
+    meanlog = meanlog, sdlog = sdlog
+  )
+
+  q_values <- 0:20
+  pwindow <- 2
+  result_numeric <- primary_censored_cdf(
+    obj,
+    q = q_values, pwindow = pwindow, use_numeric = TRUE
+  )
+  result_analytical <- primary_censored_cdf(
+    obj,
+    q = q_values, pwindow = pwindow, use_numeric = FALSE
+  )
+
+  # Check properties of numeric result
+  expect_true(all(result_numeric >= 0))
+  expect_true(all(result_numeric <= 1))
+  expect_type(result_numeric, "double")
+  expect_length(result_numeric, length(q_values))
+  expect_true(all(diff(result_numeric) >= 0)) # Ensure CDF is non-decreasing
+
+  # Check that analytical and numeric results are the same
+  expect_equal(result_numeric, result_analytical, tolerance = 1e-6)
+})
+
 test_that("primary_censored_cdf methods dispatch correctly", {
   pdist_name <- "pgamma"
   pdist <- pgamma

From dfdb09a911ae95ad35696c0d15c21974522cb76c Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Fri, 20 Sep 2024 12:41:03 +0100
Subject: [PATCH 12/15] make the comparison between numeric and non-numeric a
 lot more extensive

---
 R/pprimarycensoreddist.R                    |   2 +-
 tests/testthat/test-primary_censored_dist.R | 235 +++++++++-----------
 2 files changed, 109 insertions(+), 128 deletions(-)

diff --git a/R/pprimarycensoreddist.R b/R/pprimarycensoreddist.R
index a7a9033..b28d1a2 100644
--- a/R/pprimarycensoreddist.R
+++ b/R/pprimarycensoreddist.R
@@ -121,7 +121,7 @@ pprimarycensoreddist <- function(
   )
 
   # Compute the CDF using the S3 method
-  result <- primary_censored_cdf(pcens_obj, q, pwindow)
+  result <- primary_censored_cdf(pcens_obj, q, pwindow, use_numeric = TRUE)
 
   if (!is.infinite(D)) {
     # Compute normalization factor for finite D
diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R
index e6efca2..1ea869d 100644
--- a/tests/testthat/test-primary_censored_dist.R
+++ b/tests/testthat/test-primary_censored_dist.R
@@ -25,80 +25,6 @@ test_that("new_primary_censored_dist creates object with correct structure", {
   expect_identical(obj, new_obj)
 })
 
-test_that("primary_censored_cdf.default computes correct values", {
-  pdist_name <- "pgamma"
-  pdist <- pgamma
-  dprimary_name <- "dunif"
-  dprimary <- dunif
-  shape <- 2
-  rate <- 0.5
-
-  obj <- new_primary_censored_dist(
-    pdist,
-    dprimary, list(),
-    pdist_name, dprimary_name,
-    shape = shape, rate = rate
-  )
-
-  q_values <- 0:20
-  pwindow <- 2
-  result_numeric <- primary_censored_cdf(
-    obj,
-    q = q_values, pwindow = pwindow, use_numeric = TRUE
-  )
-  result_analytical <- primary_censored_cdf(
-    obj,
-    q = q_values, pwindow = pwindow, use_numeric = FALSE
-  )
-
-  # Check properties of numeric result
-  expect_true(all(result_numeric >= 0))
-  expect_true(all(result_numeric <= 1))
-  expect_type(result_numeric, "double")
-  expect_length(result_numeric, length(q_values))
-  expect_true(all(diff(result_numeric) >= 0)) # Ensure CDF is non-decreasing
-
-  # Check that analytical and numeric results are the same
-  expect_equal(result_numeric, result_analytical, tolerance = 1e-6)
-})
-
-test_that("primary_censored_cdf.pcens_plnorm_dunif uses numeric method", {
-  pdist_name <- "plnorm"
-  pdist <- plnorm
-  dprimary_name <- "dunif"
-  dprimary <- dunif
-  meanlog <- 0
-  sdlog <- 1
-
-  obj <- new_primary_censored_dist(
-    pdist,
-    dprimary, list(),
-    pdist_name, dprimary_name,
-    meanlog = meanlog, sdlog = sdlog
-  )
-
-  q_values <- 0:20
-  pwindow <- 2
-  result_numeric <- primary_censored_cdf(
-    obj,
-    q = q_values, pwindow = pwindow, use_numeric = TRUE
-  )
-  result_analytical <- primary_censored_cdf(
-    obj,
-    q = q_values, pwindow = pwindow, use_numeric = FALSE
-  )
-
-  # Check properties of numeric result
-  expect_true(all(result_numeric >= 0))
-  expect_true(all(result_numeric <= 1))
-  expect_type(result_numeric, "double")
-  expect_length(result_numeric, length(q_values))
-  expect_true(all(diff(result_numeric) >= 0)) # Ensure CDF is non-decreasing
-
-  # Check that analytical and numeric results are the same
-  expect_equal(result_numeric, result_analytical, tolerance = 1e-6)
-})
-
 test_that("primary_censored_cdf methods dispatch correctly", {
   pdist_name <- "pgamma"
   pdist <- pgamma
@@ -136,56 +62,111 @@ test_that("primary_censored_cdf methods dispatch correctly", {
   )
 })
 
-test_that("primary_censored_cdf.pcens_pgamma_dunif uses numeric method", {
-  pdist_name <- "pgamma"
-  pdist <- pgamma
-  dprimary_name <- "dunif"
-  dprimary <- dunif
-  shape <- 2
-  rate <- 1
-
-  obj <- new_primary_censored_dist(
-    pdist,
-    dprimary, list(),
-    pdist_name, dprimary_name,
-    shape = shape, rate = rate
-  )
-
-  q_values <- c(5, 10)
-  pwindow <- 10
-
-  expect_identical(
-    primary_censored_cdf(obj, q = q_values, pwindow = pwindow),
-    primary_censored_cdf(
-      obj,
-      q = q_values, pwindow = pwindow, use_numeric = TRUE
-    )
-  )
-})
-
-test_that("primary_censored_cdf.pcens_plnorm_dunif uses numeric method", {
-  pdist_name <- "plnorm"
-  pdist <- plnorm
-  dprimary_name <- "dunif"
-  dprimary <- dunif
-  meanlog <- 0
-  sdlog <- 1
-
-  obj <- new_primary_censored_dist(
-    pdist,
-    dprimary, list(),
-    pdist_name, dprimary_name,
-    meanlog = meanlog, sdlog = sdlog
-  )
-
-  q_values <- c(5, 10)
-  pwindow <- 10
-
-  expect_identical(
-    primary_censored_cdf(obj, q = q_values, pwindow = pwindow),
-    primary_censored_cdf(
-      obj,
-      q = q_values, pwindow = pwindow, use_numeric = TRUE
-    )
-  )
-})
+test_that(
+  "primary_censored_cdf.default computes correct values for various
+   parameters",
+  {
+    pdist_name <- "pgamma"
+    pdist <- pgamma
+    dprimary_name <- "dunif"
+    dprimary <- dunif
+
+    shapes <- c(0.5, 1, 2, 5)
+    rates <- c(0.1, 0.5, 1, 2)
+    pwindows <- c(1, 2, 5, 10)
+
+    for (shape in shapes) {
+      for (rate in rates) {
+        for (pwindow in pwindows) {
+          obj <- new_primary_censored_dist(
+            pdist,
+            dprimary, list(),
+            pdist_name, dprimary_name,
+            shape = shape, rate = rate
+          )
+
+          q_values <- seq(0, 30, by = 0.1)
+          result_numeric <- primary_censored_cdf(
+            obj,
+            q = q_values, pwindow = pwindow, use_numeric = TRUE
+          )
+          result_analytical <- primary_censored_cdf(
+            obj,
+            q = q_values, pwindow = pwindow, use_numeric = FALSE
+          )
+
+          # Check properties of numeric result
+          expect_type(result_numeric, "double")
+          expect_length(result_numeric, length(q_values))
+          expect_true(
+            all(diff(result_numeric) >= 0)
+          ) # Ensure CDF is non-decreasing
+
+          # Check that analytical and numeric results are the same
+          expect_equal(
+            result_numeric, result_analytical,
+            tolerance = 1e-5,
+            info = sprintf(
+              "Mismatch for shape = %s, rate = %s, pwindow = %s",
+              shape, rate, pwindow
+            )
+          )
+        }
+      }
+    }
+  }
+)
+
+test_that(
+  "primary_censored_cdf.pcens_plnorm_dunif uses numeric method for various
+   parameters",
+  {
+    pdist_name <- "plnorm"
+    pdist <- plnorm
+    dprimary_name <- "dunif"
+    dprimary <- dunif
+
+    meanlogs <- c(-1, 0, 1, 2)
+    sdlogs <- c(0.5, 1, 1.5)
+    pwindows <- c(1, 2, 5, 8)
+
+    for (meanlog in meanlogs) {
+      for (sdlog in sdlogs) {
+        for (pwindow in pwindows) {
+          obj <- new_primary_censored_dist(
+            pdist,
+            dprimary, list(),
+            pdist_name, dprimary_name,
+            meanlog = meanlog, sdlog = sdlog
+          )
+
+          q_values <- seq(0, 30, by = 0.1)
+          result_numeric <- primary_censored_cdf(
+            obj,
+            q = q_values, pwindow = pwindow, use_numeric = TRUE
+          )
+          result_analytical <- primary_censored_cdf(
+            obj,
+            q = q_values, pwindow = pwindow, use_numeric = FALSE
+          )
+          # Check properties of numeric result
+          expect_type(result_numeric, "double")
+          expect_length(result_numeric, length(q_values))
+          expect_true(
+            all(diff(result_numeric) >= 0)
+          ) # Ensure CDF is non-decreasing
+
+          # Check that analytical and numeric results are the same
+          expect_equal(
+            result_numeric, result_analytical,
+            tolerance = 1e-5,
+            info = sprintf(
+              "Mismatch for meanlog = %s, sdlog = %s, pwindow = %s",
+              meanlog, sdlog, pwindow
+            )
+          )
+        }
+      }
+    }
+  }
+)

From 074c9e8ec67875f6bebbe4574a456a57a079e6f1 Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Fri, 20 Sep 2024 12:42:31 +0100
Subject: [PATCH 13/15] switch analytical solutions back on

---
 R/pprimarycensoreddist.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/pprimarycensoreddist.R b/R/pprimarycensoreddist.R
index b28d1a2..a7a9033 100644
--- a/R/pprimarycensoreddist.R
+++ b/R/pprimarycensoreddist.R
@@ -121,7 +121,7 @@ pprimarycensoreddist <- function(
   )
 
   # Compute the CDF using the S3 method
-  result <- primary_censored_cdf(pcens_obj, q, pwindow, use_numeric = TRUE)
+  result <- primary_censored_cdf(pcens_obj, q, pwindow)
 
   if (!is.infinite(D)) {
     # Compute normalization factor for finite D

From 8dd43a3b4e37c29c151180edd69b36ba20712a02 Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Fri, 20 Sep 2024 12:47:55 +0100
Subject: [PATCH 14/15] update news

---
 NEWS.md             | 1 +
 touchstone/script.R | 6 +++---
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 626bbc2..553aa86 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,6 +7,7 @@ This is the current development version.
 * Add `{touchstone}` based benchmarks for benchmarking R utility functions, and fitting the `stan` and `fitdistplus` models.
 * Added a "How it works" vignette.
 * Added R infrastructure for analytical solutions via the `primary_censored_dist` S3 class.
+* Added analytical solutions for the gamma and lognormal distributions with uniform primary censoring.
 
 # primarycensoreddist 0.4.0
 
diff --git a/touchstone/script.R b/touchstone/script.R
index 1b7b637..74d8804 100644
--- a/touchstone/script.R
+++ b/touchstone/script.R
@@ -22,7 +22,7 @@ touchstone::benchmark_run(
     library(primarycensoreddist)
     q <- seq(0, 10, by = 0.01)
   },
-  pprimarycensoreddist_expgrowth = {
+  pprimarycensoreddist_plnrom_expgrowth = {
     pprimarycensoreddist(
       q, plnorm,
       dprimary = dexpgrowth,
@@ -51,7 +51,7 @@ touchstone::benchmark_run(
     library(primarycensoreddist)
     x <- seq(0, 10, by = 1)
   },
-  dprimarycensoreddist_expgrowth = {
+  dprimarycensoreddist_pweibull_expgrowth = {
     dprimarycensoreddist(
       x, pweibull,
       dprimary = dexpgrowth,
@@ -118,7 +118,7 @@ touchstone::benchmark_run(
       right = samples + swindow
     )
   },
-  fitdistdoublecens_gamma = {
+  fitdistdoublecens_gamma_expgrowth = {
     fitdistdoublecens(
       delay_data,
       distr = "gamma",

From 006d118e65fd4d0790f0ca13428b93238cf64538 Mon Sep 17 00:00:00 2001
From: Sam <s.e.abbott12@gmail.com>
Date: Fri, 20 Sep 2024 13:28:50 +0100
Subject: [PATCH 15/15] add a test for mistakenly passing dist parameters

---
 R/primary_censored_dist.R                   |   5 +-
 tests/testthat/test-primary_censored_dist.R | 130 +++++++++++++++-----
 2 files changed, 99 insertions(+), 36 deletions(-)

diff --git a/R/primary_censored_dist.R b/R/primary_censored_dist.R
index 97746d0..c1df072 100644
--- a/R/primary_censored_dist.R
+++ b/R/primary_censored_dist.R
@@ -117,9 +117,10 @@ primary_censored_cdf.pcens_pgamma_dunif <- function(
   # Extract Gamma distribution parameters
   shape <- object$args$shape
   scale <- object$args$scale
+  rate <- object$args$rate
   # if we don't have scale get fromm rate
-  if (is.null(scale)) {
-    scale <- 1 / object$args$rate
+  if (is.null(scale) && !is.null(rate)) {
+    scale <- 1 / rate
   }
   if (is.null(shape)) {
     stop("shape parameter is required for Gamma distribution")
diff --git a/tests/testthat/test-primary_censored_dist.R b/tests/testthat/test-primary_censored_dist.R
index 1ea869d..fd6094d 100644
--- a/tests/testthat/test-primary_censored_dist.R
+++ b/tests/testthat/test-primary_censored_dist.R
@@ -25,46 +25,108 @@ test_that("new_primary_censored_dist creates object with correct structure", {
   expect_identical(obj, new_obj)
 })
 
-test_that("primary_censored_cdf methods dispatch correctly", {
-  pdist_name <- "pgamma"
-  pdist <- pgamma
-  dprimary_name <- "dunif"
-  dprimary <- dunif
+test_that(
+  "primary_censored_cdf methods dispatch correctly to existing
+   analytical solutions",
+  {
+    pdist_name <- "pgamma"
+    pdist <- pgamma
+    dprimary_name <- "dunif"
+    dprimary <- dunif
 
-  obj_gamma <- new_primary_censored_dist(
-    pdist, dprimary, list(),
-    pdist_name, dprimary_name,
-    shape = 2, rate = 1
-  )
+    obj_gamma <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      shape = 2, rate = 1
+    )
 
-  pdist_name <- "plnorm"
-  pdist <- plnorm
-  dprimary_name <- "dunif"
-  dprimary <- dunif
+    pdist_name <- "plnorm"
+    pdist <- plnorm
+    dprimary_name <- "dunif"
+    dprimary <- dunif
 
-  obj_lnorm <- new_primary_censored_dist(
-    pdist, dprimary, list(),
-    pdist_name, dprimary_name,
-    meanlog = 0, sdlog = 1
-  )
+    obj_lnorm <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      meanlog = 0, sdlog = 1
+    )
 
-  expect_s3_class(obj_gamma, "pcens_pgamma_dunif")
-  expect_s3_class(obj_lnorm, "pcens_plnorm_dunif")
+    expect_s3_class(obj_gamma, "pcens_pgamma_dunif")
+    expect_s3_class(obj_lnorm, "pcens_plnorm_dunif")
 
-  q_values <- c(5, 10)
-  pwindow <- 10
+    q_values <- c(5, 10)
+    pwindow <- 10
 
-  expect_no_error(
-    primary_censored_cdf(obj_gamma, q = q_values, pwindow = pwindow)
-  )
-  expect_no_error(
-    primary_censored_cdf(obj_lnorm, q = q_values, pwindow = pwindow)
-  )
-})
+    expect_no_error(
+      primary_censored_cdf(obj_gamma, q = q_values, pwindow = pwindow)
+    )
+    expect_no_error(
+      primary_censored_cdf(obj_lnorm, q = q_values, pwindow = pwindow)
+    )
+  }
+)
+
+test_that(
+  "primary_censored_cdf errors as expected when the wrong distributional
+   parameters are supplied",
+  {
+    pdist_name <- "pgamma"
+    pdist <- pgamma
+    dprimary_name <- "dunif"
+    dprimary <- dunif
+
+    obj_gamma <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      rate = 1
+    )
+
+    expect_error(
+      primary_censored_cdf(obj_gamma, q = 1, pwindow = 1),
+      "shape parameter is required for Gamma distribution"
+    )
+
+    obj_gamma_no_rate <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      shape = 2
+    )
+
+    expect_error(
+      primary_censored_cdf(obj_gamma_no_rate, q = 1, pwindow = 1),
+      "scale or rate parameter is required for Gamma distribution"
+    )
+
+    pdist_name <- "plnorm"
+    pdist <- plnorm
+
+    obj_lnorm_no_meanlog <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      sdlog = 1
+    )
+
+    expect_error(
+      primary_censored_cdf(obj_lnorm_no_meanlog, q = 1, pwindow = 1),
+      "meanlog parameter is required for Log-Normal distribution"
+    )
+
+    obj_lnorm_no_sdlog <- new_primary_censored_dist(
+      pdist, dprimary, list(),
+      pdist_name, dprimary_name,
+      meanlog = 0
+    )
+
+    expect_error(
+      primary_censored_cdf(obj_lnorm_no_sdlog, q = 1, pwindow = 1),
+      "sdlog parameter is required for Log-Normal distribution"
+    )
+  }
+)
 
 test_that(
-  "primary_censored_cdf.default computes correct values for various
-   parameters",
+  "primary_censored_cdf.default computes the same values as
+   primary_censored_cdf.pcens_pgamma_dunif",
   {
     pdist_name <- "pgamma"
     pdist <- pgamma
@@ -118,8 +180,8 @@ test_that(
 )
 
 test_that(
-  "primary_censored_cdf.pcens_plnorm_dunif uses numeric method for various
-   parameters",
+  "primary_censored_cdf.default computes the same values as
+   primary_censored_cdf.pcens_plnorm_dunif",
   {
     pdist_name <- "plnorm"
     pdist <- plnorm