test_that("lpcde default output", {
  set.seed(42)
  n <- 100
  x_data <- matrix(rnorm(1 * n, mean = 0, sd = 1), ncol = 1)
  y_data <- matrix(rnorm(n, mean = 0, sd = 1))
  y_grid <- stats::quantile(y_data, seq(from = 0.1, to = 0.9, by = 0.1))

  # density estimation
  model2 <- lpcde(x_data = x_data, y_data = y_data, y_grid = y_grid, x = 0, bw = 1.8)
  print(model2)
  confint(model2)
  coef(model2)
  vcov(model2)
  plot(model2)
  summary(model2, CIuniform = TRUE)
  confint(model2, CIuniform = TRUE)
  expect_equal(model2$opt$bw_type, "user provided")
  expect_equal(model2$opt$p, 2)
  expect_equal(model2$opt$q, 1)
  expect_equal(model2$opt$mu, 1)
  expect_equal(model2$opt$nu, 0)
  expect_equal(model2$opt$kernel, "epanechnikov")
  expect_equal(model2$opt$p_RBC, 3)
  expect_equal(model2$opt$q_RBC, 2)
})

test_that("lpcde modified output", {
  set.seed(42)
  n <- 100
  x_data <- matrix(rnorm(1 * n, mean = 0, sd = 1), ncol = 1)
  y_data <- matrix(rnorm(n, mean = 0, sd = 1))
  y_grid <- stats::quantile(y_data, seq(from = 0.1, to = 0.9, by = 0.1))

  # density estimation
  model2 <- lpcde(x_data = x_data, y_data = y_data, y_grid = y_grid, x = 0, mu = 0, p = 3, bw_type = "imse-rot")
  summary(model2)
  expect_equal(model2$opt$bw_type, "imse-rot")
  expect_equal(model2$opt$p, 3)
  expect_equal(model2$opt$q, 1)
  expect_equal(model2$opt$mu, 0)
  expect_equal(model2$opt$nu, 0)
  expect_equal(model2$opt$kernel, "epanechnikov")
  expect_equal(model2$opt$p_RBC, 4)
  expect_equal(model2$opt$q_RBC, 2)
})


test_that("error checking", {
  set.seed(42)
  n <- 100
  x_data <- matrix(rnorm(1 * n, mean = 0, sd = 1), ncol = 1)
  y_data <- matrix(rnorm(n, mean = 0, sd = 1))
  y_grid <- stats::quantile(y_data, seq(from = 0.1, to = 0.9, by = 0.1))

  # density estimation
  #  model2 = lpcde(x_data=x_data, y_data=y_data, y_grid=y_grid, x=0, p=2, q=1, p_RBC=2, q_RBC=1, bw=1.8)
  #  expect_equal(model2$opt$p, model2$opt$p_RBC)
  #  expect_equal(model2$opt$q, model2$opt$q_RBC)

  model2 <- lpcde(x_data = x_data, y_data = y_data, y_grid = y_grid, x = 0, kernel_type = "triangular", bw = 1.8, cov_flag = "diag")
  expect_equal(model2$opt$kernel, "triangular")
  expect_equal(model2$opt$cov_flag, "diag")

  model2 <- lpcde(x_data = x_data, y_data = y_data, y_grid = y_grid, x = 0, kernel_type = "triangular", bw = 1.8, cov_flag = "off")
  expect_equal(model2$CovMat$CovMat, NA)
  expect_equal(model2$CovMat$CovMat_RBC, NA)
})

test_that("lpcde multivariate output", {
  set.seed(42)
  n <- 100
  x_data <- matrix(rnorm(2 * n, mean = 0, sd = 1), ncol = 2)
  y_data <- matrix(rnorm(n, mean = 0, sd = 1))
  y_grid <- stats::quantile(y_data, seq(from = 0.1, to = 0.9, by = 0.1))

  # density estimation
  model2 <- lpcde(x_data = x_data, y_data = y_data, y_grid = y_grid, x = matrix(c(0, 0), ncol = 2), bw_type = "imse-rot")
  summary(model2, CIuniform = TRUE)
  expect_equal(model2$opt$p, 2)
  expect_equal(model2$opt$q, 1)
  expect_equal(model2$opt$mu, 1)
  expect_equal(model2$opt$nu, 0)
  expect_equal(model2$opt$kernel, "epanechnikov")
  expect_equal(model2$opt$p_RBC, 3)
  expect_equal(model2$opt$q_RBC, 2)
})

test_that("Properties of pdf estimator", {
  # Setting up simulation
  set.seed(30)
  n <- 100
  x_data <- matrix(rnorm(n, mean = 0, sd = 1))
  y_data <- matrix(rnorm(n, mean = x_data, sd = 1))
  y_grid <- stats::quantile(y_data, seq(from = 0.1, to = 0.9, by = 0.1))
  # Regularized density estimation
  model_reg <- lpcde::lpcde(x_data = x_data, y_data = y_data, y_grid = y_grid, x = 0, bw = 1, nonneg = TRUE, normalize = TRUE)

  # check integration to 1
  grid_diff <- c(diff(y_grid), diff(utils::tail(y_grid, 2)))
  c <- sum(model_reg$Estimate[, 3] * grid_diff)
  expect_equal(1, c)

  # check nonegativity
  a_nng <- any(model_reg$Estimate[, 3] < 0)
  expect_equal(a_nng, FALSE)

  # check all probabilities are less than 1
  expect_equal(all(model_reg$Estimate[, 3] < 1), TRUE)
})