tests/testthat/test-GPCmodel.R

In GPCsign: Gaussian Process Classification as Described in Bachoc et al. (2020)

require("GPCsign")
require("testthat")

seed = 123


# -------------------------------------------------------------
# -------------------------------------------------------------
# ------------------------Test gpcm function ------------------
# -------------------------------------------------------------
# -------------------------------------------------------------


test_that.gpcm <- function(model, coef.m = NULL, coef.cov = NULL,logLik = NULL, precision = 1e-4, multistart = 1, bounds = FALSE, Z_obs = NULL) {
  
  if(multistart==1){
    if(!is.null(coef.m) && !is.null(coef.cov))
      test_that(desc = "Check GPC coef.m and coef.cov",
                expect_true((abs((model@coef.m - coef.m)) < precision) &&  (max(abs((model@coef.cov - coef.cov)/coef.cov)) < precision)))
    
    if(is.null(coef.m) || is.null(coef.cov)){
      test_that(desc = "Optimization works when at least one parameter is missed",
                expect_true(model@param.estim == T))
    }
    else{
      test_that(desc = "No optimization when both parameters are provided and multistart=1 ",
                   expect_true(model@param.estim == F))}
  }
  if(multistart>1){
    test_that(desc = "Optimization works when multistart>1",
              expect_true(model@param.estim == T))
    test_that(desc = "Check GPC coef.m and coef.cov",
            expect_true((abs((model@coef.m - coef.m)) < precision) &&  (max(abs((model@coef.cov - coef.cov)/coef.cov)) < precision)))}

  
  if(bounds)
    test_that(desc = "Check GPC bounds",
              expect_true(unique(((model@y == 1) == (model@u == Inf)) & ((model@y == -1) == (model@l == -Inf)))))
  if (!is.null(Z_obs))
    test_that(desc = "Check GPC Z_obs",
              expect_true(max(abs((model@Z_obs - Z_obs)/Z_obs)) < precision))
  
  if (!is.null(logLik))
    test_that(desc = "Check GPC logLik",
              expect_true(max(abs((model@logLik - logLik)/logLik)) < precision))
}


# ----------------------------------
# A 2D example - Branin-Hoo function
# ----------------------------------

# 10-points DoE, and the corresponding response
d <- 2
nb_PX <- 10
x <- matrix(c(0.205953271570615,0.0427366598043591,
            0.610758095560595,0.16720792807173,
            0.195444350061007,0.875391226564534,
            0.921169486455619,0.204316665465012,
            0.411698259599507,0.604549635085277,
            0.871242247987539,0.795794046646915,
            0.05910230781883,0.410017502959818,
            0.354338526469655,0.354666584380902,
            0.5471894511953,0.989707531733438,
            0.744856498553418,0.532242936454713), ncol = d, byrow = T)

require(DiceKriging)
fx <- apply(x, 1, branin)
f <- ifelse(fx < 14, -1, 1)
Xf <- as.matrix(x)

# Test optimisation, multistart = 1
coef.cov <- c(0.3,0.2)
coef.m <- NULL
m1 <- gpcm(f, Xf, seed = seed, coef.cov = coef.cov, coef.m = coef.m)
test_that.gpcm(m1,
               coef.m = coef.m,
               coef.cov = coef.cov,
               bounds = TRUE)


# Test optimisation, multistart > 1
multistart <- 3
coef.cov <- c(0.3,0.2) 
coef.m <- 0
m2 <- gpcm(f, Xf, seed = seed, coef.cov = coef.cov , coef.m = coef.m, multistart = multistart)
m2
test_that.gpcm(m2,
               coef.m = 1.024028,
               coef.cov = c(3.0000000, 0.3046948),
               logLik = -3.66642,
               multistart = 3,
               bounds = TRUE)

## without parameters
m3 <- gpcm(f, Xf, seed = seed, multistart = multistart)
test_that.gpcm(m3,
               coef.m = 1.024029,
               coef.cov = c(3.0000000, 0.3047),
               logLik = -3.66642,
               multistart = 3,
               bounds = TRUE)


# Test Z_obs with multistart > 1
m4 <- gpcm(f, Xf, seed = seed, multistart = multistart, nsimu = 20)

Z_obs <- matrix(c(0.6065117,-0.1232498,0.564445,-0.07775007,1.725963,0.1490353,0.5739855,1.469259,1.173223,1.312926,2.407662,-0.2785747,1.273053,-0.3711358,0.3702328,2.258483,
1.05638,0.1263969,0.8672171,2.019331,1.983481,-0.1832863,1.783672,-0.005671975,2.094133,1.989232,0.8155662,0.9623516,0.9207395,0.8321107,2.420045,-0.03761402,1.814366,-0.1525295,
0.09056849,1.307406,1.553658,0.7797381,1.023494,0.3931547,0.3340786,-0.4567696,1.069838,-0.6076337,0.2979702,0.3906739,0.5461513,0.6279915,0.7312548,1.477448,0.1047055,-0.6591164,
1.527485,-1.057407,1.426642,0.730582,0.3907899,0.9093211,2.444633,0.9755318,1.314244,-1.325511,1.197695,-1.306888,1.731924,1.150123,1.714716,1.786977,1.925169,1.462034,1.313841,
-0.6507607,1.87047,-2.228603,1.271865,0.8200389,1.724701,1.217291,1.126575,0.4967319,0.7133302,-1.349531,0.9590339,-0.5680959,0.2080481,0.9577018,2.654053,2.864128,2.200684,0.4092046,
0.1712899,-0.4099546,1.131511,-0.2369335,0.5227997,0.507785,2.86789,1.311164,1.159168,1.201136,1.090928,-0.4969704,1.01599,-0.05570711,1.103215,2.182993,2.24725,1.953898,0.9653986,
0.5908197,2.04882,-0.4106678,0.3748812,-0.05590933,1.384168,0.334814,1.778675,2.566496,1.188518,1.229792,1.150435,-0.4689835,0.6499797,-0.7425705,0.890075,2.834854,1.611013,0.4144046,
0.3367303,1.473362,1.219667,-0.116078,1.898891,-0.2066555,1.547525,1.886384,1.42774,1.672934,3.258192,1.298007,0.6519265,-0.464371,0.30107,-0.8174835,2.06121,0.4876322,1.203768,
0.2109134,0.5241103,2.185415,1.842686,-0.471029,0.6937389,-0.9308983,0.8006155,0.8254235,0.9964199,0.4014273,0.9900444,0.5247607,1.02709,-0.7963877,1.375042,-0.8488883,0.658907,
1.330128,1.174522,0.5160265,1.093934,0.5773668,1.136948,-0.9915276,2.140794,-0.4937453,1.241677,1.853937,0.6342596,0.7686075,2.526545,1.289942,1.508144,-0.5864517,2.308951,-0.8620192,
1.494225,1.57945,0.6966728,0.7779959,2.793423,2.240587,0.5557361,-0.8913785,3.663125,-0.505999,3.016084,2.240008,1.008555,1.106413,1.110652,2.292958
), 
               nrow = length(f), byrow = F)

test_that.gpcm(m4,Z_obs = Z_obs)




# -------------------------------------------------------------
# -------------------------------------------------------------
# ------------------------Test predict function ---------------
# -------------------------------------------------------------
# -------------------------------------------------------------

# ----------------------------------
# A 1D example - sinusoidal function
# ----------------------------------

sinusoidal_function <- function(x) sin(4 * pi * x)

Xf <- as.matrix(c(0.07, 0.19, 0.42, 0.56, 0.81, 0.90))
f <- rep(1, length(Xf))
f[(sinusoidal_function(Xf) < 0)] <- -1

## parameters not given & without normalization
m5 <- gpcm(f, Xf, seed = seed, normalize = F)
newdata <- as.matrix(seq(0, 1, length.out = 100))
pred <- predict(m5, newdata = newdata)$prob

test_that(desc="Test predict",
          expect_true((max(abs(pred - c(0.8818044,0.9006364,0.9200324,0.9406092,0.9617672,0.9796049,0.9955032,1,0.9999558,0.9989784,0.9970659,0.9952124,0.9937841,
          0.9927823,0.9921665,0.9919429,0.9921846,0.9931205,0.9954743,0.9784406,0.8972897,0.8411962,0.7955495,0.7517711,0.7048484,0.6558069,0.6059068,0.5564561,0.50838,
          0.4620602,0.4174992,0.3745191,0.3328849,0.292358,0.2527066,0.2137208,0.1752975,0.1376657,0.1017053,0.06897952,0.04075618,0.01767788,0.03123165,0.09003728,0.1437574,
          0.2082169,0.2851483,0.3740537,0.4668087,0.5586693,0.6497188,0.7384403,0.8193306,0.8877447,0.9419577,0.9861225,0.9845287,0.9609706,0.9355048,0.9064329,0.8757478,0.8438976,
          0.8103801,0.7748394,0.7373762,0.6983578,0.6582188,0.6173768,0.5762199,0.5351191,0.494427,0.4544326,0.4152436,0.3765991,0.337706,0.2971896,0.252866,0.203407,0.1533051,0.0878467,
          0.007552753,0.0003375535,0.001277704,0.001798767,0.00180658,0.00143300,0.0008379434,0.0002746014,3.792301e-05,0,0.02048777,0.03250841,0.05087609,0.06961104,0.08778701,0.1060301,
          0.1242306,0.1421142,0.1595297,0.1764747
                                        
))) < 1e-6)))



# -------------------------------------------------------------
# -------------------------------------------------------------
# ------------------------Test update function ---------------
# -------------------------------------------------------------
# -------------------------------------------------------------

# ----------------------------------
# A 2D example - Branin-Hoo function
# ----------------------------------

# 20-points DoE, and the corresponding response
d <- 2
nb_PX <- 20

x <- matrix(c(0.402976635785308,0.422796698764432,
            0.116121468227357,0.948769315774552,
            0.27136832990218,0.585542013135273,
            0.577333292190451,0.250286511168815,
            0.735621123993769,0.0761102014454082,
            0.0550087514799088,0.605523034196813,
            0.227169263234828,0.789179603208322,
            0.983603964035865,0.714573476591613,
            0.797722175030503,0.312077023123857,
            0.305849129799753,0.0179746593115851,
            0.0448537658667192,0.192644317634404,
            0.92359472559765,0.134090949618258,
            0.46058474322781,0.86535982969217,
            0.629551153909415,0.470292898977641,
            0.852158332732506,0.538418710732367,
            0.522428249276709,0.660226629115641,
            0.352274817542639,0.217214710044209,
            0.687695613282267,0.801848788373172,
            0.805379047780298,0.965464736078866,
            0.197897023323458,0.354885047744028), ncol = 2, byrow = T)
fx <- apply(x, 1, branin)
f <- ifelse(fx < 14, -1, 1)
Xf <- as.matrix(x)

## build the first model
m7 <- gpcm(f, Xf, seed = seed, multistart = 5)

## add new points
newX <- matrix(c(0.87344913368579,0.87940254251007,
                 0.462787610036321,0.0823443247471005,
                 0.94271466366481,0.331297715334222,
                 0.00917922100052238,0.261589628178626,
                 0.679831806896254,0.523015858000144,
                 0.71016103150323,0.150230075791478,
                 0.505532473139465,0.428238149173558,
                 0.133920220751315,0.700809390516952,
                 0.237088595610112,0.961996482056566,
                 0.393821372953244,0.62225547786802), ncol = 2, byrow = T)
newfx <- apply(newX, 1, branin)
newf <- ifelse(newfx < 14, -1, 1)
newXf <- as.matrix(newX)

## build the updated model
m8 <- update(object = m7, newf = newf, newXf = newXf, seed = seed, nsimu  = 50)

coef.cov = c(0.2000, 1.57391)
coef.m = 0.7720781
test_that(desc = "Check GPC coef.m and coef.cov",
          expect_true((abs((m8@coef.m - coef.m)) < 1e-4) &&  (max(abs((m8@coef.cov - coef.cov)/coef.cov)) < 1e-4)))



# test update with covandmean.reestim = F
m9 <- update(object = m7, newf = newf,covandmean.reestim = F, newXf = newXf, seed = seed, nsimu  = 50)

test_that(desc = "Test re-estimation of cov and mean ",
          expect_equal(c(m9@coef.cov,m9@coef.m),c(m7@coef.cov, m7@coef.m)))