test_that("bp_fit works with analytic derivatives for scalar theta", {
  set.seed(101)
  n <- 200
  x <- rnorm(n)
  y <- 0.5 * x + rnorm(n)

  U_fun <- function(theta) {
    cbind(y - theta[1], x)
  }
  dU_fun <- function(theta) {
    dU <- array(0, dim = c(n, 2, 1))
    dU[, 1, 1] <- -1
    dU
  }

  fit <- bp_fit(mean(y), U_fun, dU_fun, generator = "quadratic")
  expect_s3_class(fit, "bregproj_fit")
  expect_true(fit$success)
  expect_length(fit$theta, 1)
  expect_equal(dim(vcov(fit)), c(1L, 1L))
  expect_length(fit$std_error, 1)
  expect_true(is.finite(fit$std_error))
  expect_gt(fit$std_error, 0)
  expect_lt(fit$max_abs_moment_residual, 1e-5)
})

test_that("bp_fit works with finite differences", {
  set.seed(102)
  n <- 120
  x <- rnorm(n)
  y <- 0.3 * x + rnorm(n)
  U_fun <- function(theta) cbind(y - theta[1], x)

  fit <- bp_fit(mean(y), U_fun, generator = "quadratic")
  expect_true(fit$success)
  expect_length(fit$theta, 1)
})

test_that("bp_fit supports vector theta", {
  set.seed(103)
  n <- 150
  x <- rnorm(n)
  y1 <- 0.2 * x + rnorm(n)
  y2 <- -0.1 * x + rnorm(n)

  U_fun <- function(theta) {
    cbind(y1 - theta[1], y2 - theta[2], x)
  }
  dU_fun <- function(theta) {
    dU <- array(0, dim = c(n, 3, 2))
    dU[, 1, 1] <- -1
    dU[, 2, 2] <- -1
    dU
  }

  fit <- bp_fit(c(mean(y1), mean(y2)), U_fun, dU_fun, generator = "quadratic")
  expect_true(fit$success)
  expect_length(fit$theta, 2)
  expect_equal(fit$p, 2)
  expect_equal(dim(vcov(fit)), c(2L, 2L))
  expect_equal(length(fit$std_error), 2L)
  expect_true(all(is.finite(fit$std_error)))
})