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)))