test_that("gmix matches the hardcoded reference from the original R ECM implementation", {
  # Reference values were derived once from:
  #   tests/originalR/constraints-dev.R
  #   tests/originalR/ecm-dev.R
  # using ECM_v2(dat, init, erc = 2, iter.max = 200, tol = 1e-8).
  # The package reports the converged code as 1 instead of the original 2,
  # maps the original ERC flag "4" to flag 2, and stores loglik scaled by N.
  dat <- matrix(
    c(
      -3, 0,
      -2, 0.2,
      -1, -0.1,
      4, 4,
      8, 4.1,
      12, 3.9
    ),
    ncol = 2,
    byrow = TRUE
  )
  init <- matrix(
    c(
      1, 0,
      1, 0,
      1, 0,
      0, 1,
      0, 1,
      0, 1
    ),
    ncol = 2,
    byrow = TRUE
  )

  fit <- gmix(dat, K = 2, erc = 2, init = init, save_taus = TRUE, iter.max = 200, tol = 1e-8)

  expected_proportion <- c(0.500000115012709, 0.499999884987291)
  expected_mean <- matrix(
    c(
      -1.99999861156069, 0.0333342524433696,
      8.00000091181476, 3.99999999332391
    ),
    nrow = 2
  )
  expected_cov <- array(
    c(
      1.50534963958242, 1.39669229830353e-18,
      -7.06402492440929e-19, 1.50534963958242,
      3.01046388452237, -0.0188227476626236,
      -0.0188227476626236, 1.5055850342249
    ),
    dim = c(2, 2, 2)
  )
  expected_posterior <- matrix(
    c(
      0.999999999992446, 7.55430430289607e-12,
      0.999999999692042, 3.07957878004184e-10,
      0.999999995407711, 4.59228892524982e-09,
      6.95007646597057e-07, 0.999999304992353,
      2.19454872673204e-17, 1,
      7.51017418651596e-30, 1
    ),
    ncol = 2,
    byrow = TRUE
  )

  expect_equal(fit$info$code, 1)
  expect_equal(fit$info$flag, 2)
  expect_equal(fit$iter, 2L)
  expect_equal(fit$loglik, -24.68000666757756, tolerance = 1e-12)
  expect_equal(unname(fit$params$proportion), expected_proportion, tolerance = 1e-10)
  expect_equal(unname(fit$params$mean), expected_mean, tolerance = 1e-10)
  expect_equal(unname(fit$params$cov), expected_cov, tolerance = 1e-7)
  expect_equal(unname(fit$posterior), expected_posterior, tolerance = 1e-10)
  expect_equal(unname(fit$cluster), c(1, 1, 1, 2, 2, 2))
  expect_equal(unname(fit$size), c(3, 3))
})