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.021006, coef.cov = c(3.0000000, 0.2972), logLik = -3.670613, multistart = 3, bounds = TRUE) ## without parameters m3 <- gpcm(f, Xf, seed = seed, multistart = multistart) test_that.gpcm(m3, coef.m = 1.021004, coef.cov = c(3.0000000, 0.2972), logLik = -3.670613, 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.6122591,2.428811,1.997571,2.442421,0.336785,0.1060837,1.327671,1.331993,0.7110898,0.1769581,1.107768,2.074128,1.165799,1.229597,0.6580723,1.856563,1.032118,1.150544,1.531873,0.5616142, -0.1246298,-0.2853614,-0.1866487,-0.0385878,-0.4608072,-0.6607498,-1.340368,-0.6508122,-1.352587,-0.4046904,-0.5033067,-0.4191615,-0.4765397,-0.1175244,-0.469454,-0.4788665,-0.8011347,-0.9972914,-0.5902873,-0.8942917, 0.5738097,1.283116,1.788634,1.810229,1.066771,1.54102,1.188013,1.864328,0.9547219,1.124077,1.018869,0.3699232,0.6534465,1.90713,0.2789438,0.6997045,1.383371,2.148386,2.292424,3.649971, -0.07846662,-0.3732359,-0.005633736,-0.1526163,-0.6118102,-1.065158,-1.311049,-2.244952,-0.5557593,-0.2382527,-0.05559122,-0.05603574,-0.7521991,-0.206962,-0.8261989,-0.9331043,-0.8474916,-0.4877425,-0.8651346,-0.5046836, 1.752005,0.3721733,2.093942,0.0881462,0.300172,1.437612,1.747223,1.268184,0.2084241,0.5297287,1.106177,1.388482,0.8906913,1.536042,2.061202,0.7800325,0.6617,1.250413,1.496399,3.019615, 0.149084,2.269732,1.981262,1.306786,0.3918181,0.7286072,1.142088,0.810268,0.9624176,0.5072006,2.191569,0.335042,2.849596,1.881522,0.4828249,0.830806,1.333814,1.848544,1.55988,2.197801, 0.580544,1.065849,0.8243611,1.568932,0.547366,0.3911982,1.727271,1.722911,2.67665,2.877148,2.247637,1.766568,1.592245,1.431793,1.192576,0.9996787,1.184742,0.6401481,0.7008235,1.012382, 1.500143,0.128217,0.9797243,0.7868665,0.6408012,0.9300409,1.817421,1.231923,2.893815,1.300211,1.955581,2.580743,0.4040945,1.689755,0.2111144,0.4080143,0.5274577,0.7876883,0.7958807,1.126261, 1.183406,0.8658004,0.9141417,1.012117,0.7292691,2.449395,1.925404,1.115166,2.205335,1.156689,0.9662961,1.196432,0.3380529,3.260172,0.525666,0.9945199,1.092513,2.520948,2.778238,1.072873, 1.319323,2.046838,0.8204979,0.3991332,1.496501,0.9785637,1.461782,0.4903173,0.4091533,1.209295,0.5846927,1.218553,1.480441,1.286601,2.181852,0.5204864,0.5829928,1.298688,2.252989,2.283047 ), nrow = length(f), byrow = T) 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.8818225,0.9006519,0.9200455,0.9406197,0.9617743,0.9796086,0.9955049,1,0.9999559,0.9989786,0.9970662,0.9952128,0.9937844,0.9927826,0.9921668,0.9919432,0.9921848,0.9931207,0.9954744,0.9784429, 0.897293,0.8412018,0.7955558,0.7517794,0.7048585,0.6558185,0.6059192,0.5564689,0.5083929,0.4620727,0.4175114,0.3745306,0.3328958,0.2923681,0.2527158,0.2137291,0.1753048,0.1376718,0.1017101,0.06898273, 0.04075796,0.01767854,0.03123219,0.0900382,0.1437594,0.2082203,0.2851535,0.3740603,0.4668158,0.5586763,0.6497255,0.7384461,0.8193349,0.8877473,0.9419587,0.9861229,0.9845294,0.9609721,0.9355078,0.9064373, 0.8757535,0.8439045,0.8103883,0.774849,0.7373871,0.6983698,0.6582317,0.6173904,0.5762339,0.5351331,0.4944409,0.4544461,0.4152566,0.3766114,0.3377174,0.2972001,0.2528754,0.2034154,0.1533127,0.08785299, 0.007553127,0.0003375928,0.001277776,0.001798843,0.001806656,0.001433082,0.0008380293,0.0002746644,3.794402e-05,0,0.02048948,0.03251266,0.05088442,0.06962142,0.08779956,0.1060452,0.1242481,0.1421343,0.159552,0.1764993 ))) < 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, 2.2213) coef.m = 0.8025925 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)))