R Under development (unstable) (2024-10-16 r87241 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > if(!require("GNE"))stop("this test requires package GNE.") Loading required package: GNE Loading required package: alabama Loading required package: numDeriv Loading required package: nleqslv Loading required package: BB Loading required package: SQUAREM > > > # (1) associated objective functions > # > > dimx <- c(2, 2, 3) > > #Gr_x_j O_i(x) > grfullob <- function(x, i, j) + { + x <- x[1:7] + if(i == 1) + { + grad <- 3*(x - 1:7)^2 + } + if(i == 2) + { + grad <- 1:7*(x - 1:7)^(0:6) + } + if(i == 3) + { + s <- x[5]^2 + x[6]^2 + x[7]^2 - 5 + grad <- c(1, 0, 1, 0, 4*x[5]*s, 4*x[6]*s, 4*x[7]*s) + + } + grad[j] + } > > > #Gr_x_k Gr_x_j O_i(x) > hefullob <- function(x, i, j, k) + { + x <- x[1:7] + if(i == 1) + { + he <- diag( 6*(x - 1:7) ) + } + if(i == 2) + { + he <- diag( c(0, 2, 6, 12, 20, 30, 42)*(x - 1:7)^c(0, 0:5) ) + } + if(i == 3) + { + s <- x[5]^2 + x[6]^2 + x[7]^2 + + he <- rbind(rep(0, 7), rep(0, 7), rep(0, 7), rep(0, 7), + c(0, 0, 0, 0, 4*s+8*x[5]^2, 8*x[5]*x[6], 8*x[5]*x[7]), + c(0, 0, 0, 0, 8*x[5]*x[6], 4*s+8*x[6]^2, 8*x[6]*x[7]), + c(0, 0, 0, 0, 8*x[5]*x[7], 8*x[6]*x[7], 4*s+8*x[7]^2) ) + } + he[j,k] + } > > > dimmu <- 3 > > #constraint function h(x) > h <- function(x) + { + x <- x[1:7] + c(prod(x) - 1, sum(x^2) -2, sum(x^3) -3) + } > grh <- function(x, j) + { + x <- x[1:7] + c(prod(x[-j]), 2*x[j], 3*x[j]^2) + } > heh <- function(x, j, k) + { + x <- x[1:7] + c(prod(x[-c(j,k)]), 2*(j==k), 6*x[j]*(j==k)) + } > > > > # (3) compute H > # > > z <- rexp(sum(dimx) + 2*sum(dimmu)) > > n <- sum(dimx) > m <- 0 > x <- z[1:n] > lam <- NULL > mu <- z[(n+m+1):(n+m+dimmu)] > w <- z[(n+m+dimmu+1):(n+2*m+2*dimmu)] > > resphi <- funCER(z, dimx, dimlam, grobj=grfullob, joint=h, grjoint=grh, dimmu=dimmu, echo=TRUE) [1] 2 2 3 [1] 0 0 0 [1] 13 > > > check <- c(grfullob(x, 1, 1) + mu %*% grh(x, 1), + grfullob(x, 1, 2) + mu %*% grh(x, 2), + grfullob(x, 2, 3) + mu %*% grh(x, 3), + grfullob(x, 2, 4) + mu %*% grh(x, 4), + grfullob(x, 3, 5) + mu %*% grh(x, 5), + grfullob(x, 3, 6) + mu %*% grh(x, 6), + grfullob(x, 3, 7) + mu %*% grh(x, 7), + h(x) + w, + mu * w) > > #check > cat("\n\n________________________________________\n\n") ________________________________________ > > #part A > print(cbind(check, res=as.numeric(resphi))[1:length(x), ]) check res [1,] 2.5261566 2.5261566 [2,] 12.6623197 12.6623197 [3,] 20.1379826 20.1379826 [4,] -125.5999769 -125.5999769 [5,] 17.7127742 17.7127742 [6,] -0.4412469 -0.4412469 [7,] -0.4052538 -0.4052538 > #part B > print(cbind(check, res=as.numeric(resphi))[(length(x)+1):length(z), ]) check res [1,] -0.5465746 -0.5465746 [2,] 5.1922036 5.1922036 [3,] 6.6003540 6.6003540 [4,] 0.5847315 0.5847315 [5,] 0.5530717 0.5530717 [6,] 1.1472332 1.1472332 > > > if(sum(abs(check - resphi)) > .Machine$double.eps^(2/3)) + stop("wrong result") > > > > # (4) compute Jac H > # > > resjacphi <- jacCER(z, dimx, dimlam, heobj=hefullob, + joint=h, grjoint=grh, hejoint=heh, dimmu=dimmu) > > > #check > cat("\n\n________________________________________\n\n") ________________________________________ > > > cat("\n\npart A\n\n") part A > > > cat("\n\npart A\n\n") part A > > lam <- rep(0, 5) > > checkA <- + rbind( + c(hefullob(x, 1, 1, 1) + mu%*%heh(x, 1, 1), + hefullob(x, 1, 1, 2) + mu%*%heh(x, 1, 2), + hefullob(x, 1, 1, 3) + mu%*%heh(x, 1, 3), + hefullob(x, 1, 1, 4) + mu%*%heh(x, 1, 4), + hefullob(x, 1, 1, 5) + mu%*%heh(x, 1, 5), + hefullob(x, 1, 1, 6) + mu%*%heh(x, 1, 6), + hefullob(x, 1, 1, 7) + mu%*%heh(x, 1, 7) + ), + c(hefullob(x, 1, 2, 1) + mu%*%heh(x, 2, 1), + hefullob(x, 1, 2, 2) + mu%*%heh(x, 2, 2), + hefullob(x, 1, 2, 3) + mu%*%heh(x, 2, 3), + hefullob(x, 1, 2, 4) + mu%*%heh(x, 2, 4), + hefullob(x, 1, 2, 5) + mu%*%heh(x, 2, 5), + hefullob(x, 1, 2, 6) + mu%*%heh(x, 2, 6), + hefullob(x, 1, 2, 7) + mu%*%heh(x, 2, 7) + ), + c(hefullob(x, 2, 3, 1) + mu%*%heh(x, 3, 1), + hefullob(x, 2, 3, 2) + mu%*%heh(x, 3, 2), + hefullob(x, 2, 3, 3) + mu%*%heh(x, 3, 3), + hefullob(x, 2, 3, 4) + mu%*%heh(x, 3, 4), + hefullob(x, 2, 3, 5) + mu%*%heh(x, 3, 5), + hefullob(x, 2, 3, 6) + mu%*%heh(x, 3, 6), + hefullob(x, 2, 3, 7) + mu%*%heh(x, 3, 7) + ), + c(hefullob(x, 2, 4, 1) + mu%*%heh(x, 4, 1), + hefullob(x, 2, 4, 2) + mu%*%heh(x, 4, 2), + hefullob(x, 2, 4, 3) + mu%*%heh(x, 4, 3), + hefullob(x, 2, 4, 4) + mu%*%heh(x, 4, 4), + hefullob(x, 2, 4, 5) + mu%*%heh(x, 4, 5), + hefullob(x, 2, 4, 6) + mu%*%heh(x, 4, 6), + hefullob(x, 2, 4, 7) + mu%*%heh(x, 4, 7) + ), + c(hefullob(x, 3, 5, 1) + mu%*%heh(x, 5, 1), + hefullob(x, 3, 5, 2) + mu%*%heh(x, 5, 2), + hefullob(x, 3, 5, 3) + mu%*%heh(x, 5, 3), + hefullob(x, 3, 5, 4) + mu%*%heh(x, 5, 4), + hefullob(x, 3, 5, 5) + mu%*%heh(x, 5, 5), + hefullob(x, 3, 5, 6) + mu%*%heh(x, 5, 6), + hefullob(x, 3, 5, 7) + mu%*%heh(x, 5, 7) + ), + c(hefullob(x, 3, 6, 1) + mu%*%heh(x, 6, 1), + hefullob(x, 3, 6, 2) + mu%*%heh(x, 6, 2), + hefullob(x, 3, 6, 3) + mu%*%heh(x, 6, 3), + hefullob(x, 3, 6, 4) + mu%*%heh(x, 6, 4), + hefullob(x, 3, 6, 5) + mu%*%heh(x, 6, 5), + hefullob(x, 3, 6, 6) + mu%*%heh(x, 6, 6), + hefullob(x, 3, 6, 7) + mu%*%heh(x, 6, 7) + ), + c(hefullob(x, 3, 7, 1) + mu%*%heh(x, 7, 1), + hefullob(x, 3, 7, 2) + mu%*%heh(x, 7, 2), + hefullob(x, 3, 7, 3) + mu%*%heh(x, 7, 3), + hefullob(x, 3, 7, 4) + mu%*%heh(x, 7, 4), + hefullob(x, 3, 7, 5) + mu%*%heh(x, 7, 5), + hefullob(x, 3, 7, 6) + mu%*%heh(x, 7, 6), + hefullob(x, 3, 7, 7) + mu%*%heh(x, 7, 7) + ) + ) > > > print(resjacphi[1:n, 1:n] - checkA) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] -4.440892e-16 0 0.000000e+00 0 0 0 0 [2,] 0.000000e+00 0 0.000000e+00 0 0 0 0 [3,] 0.000000e+00 0 -1.776357e-15 0 0 0 0 [4,] 0.000000e+00 0 0.000000e+00 0 0 0 0 [5,] 0.000000e+00 0 0.000000e+00 0 0 0 0 [6,] 0.000000e+00 0 0.000000e+00 0 0 0 0 [7,] 0.000000e+00 0 0.000000e+00 0 0 0 0 > > > cat("\n\n________________________________________\n\n") ________________________________________ > cat("\n\npart B\n\n") part B > > > checkB <- rbind(grh(x, 1), grh(x, 2), grh(x, 3), grh(x, 4), grh(x, 5), grh(x, 6), grh(x, 7)) > > print(resjacphi[1:n, (n+1):(n+m+dimmu)] - checkB) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 [3,] 0 0 0 [4,] 0 0 0 [5,] 0 0 0 [6,] 0 0 0 [7,] 0 0 0 > > cat("\n\n________________________________________\n\n") ________________________________________ > cat("\n\npart C\n\n") part C > > checkC <- matrix(0, n, m+dimmu) > > print(resjacphi[1:n, (n+m+dimmu+1):(n+2*m+2*dimmu)] - checkC) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 [3,] 0 0 0 [4,] 0 0 0 [5,] 0 0 0 [6,] 0 0 0 [7,] 0 0 0 > > cat("\n\n________________________________________\n\n") ________________________________________ > cat("\n\npart D\n\n") part D > > checkD <- cbind(grh(x, 1), grh(x, 2), grh(x, 3), grh(x, 4), grh(x, 5), grh(x, 6), grh(x, 7)) > > > print(resjacphi[(n+1):(n+m+dimmu), 1:n] - checkD) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0 0 0 0 0 0 0 [2,] 0 0 0 0 0 0 0 [3,] 0 0 0 0 0 0 0 > > > cat("\n\n________________________________________\n\n") ________________________________________ > cat("\n\npart E\n\n") part E > > > checkE <- matrix(0, m+dimmu, m+dimmu) > > print(resjacphi[(n+1):(n+m+dimmu), (n+1):(n+m+dimmu)] - checkE) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 [3,] 0 0 0 > > > cat("\n\n________________________________________\n\n") ________________________________________ > cat("\n\npart F\n\n") part F > > > checkF <- diag(m+dimmu) > > print(resjacphi[(n+1):(n+m+dimmu), (n+m+dimmu+1):(n+2*m+2*dimmu)] - checkF) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 [3,] 0 0 0 > > > cat("\n\n________________________________________\n\n") ________________________________________ > cat("\n\npart G\n\n") part G > > > checkG <- matrix(0, m+dimmu, n) > > print(resjacphi[(n+m+dimmu+1):(n+2*m+2*dimmu), 1:n] - checkG) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0 0 0 0 0 0 0 [2,] 0 0 0 0 0 0 0 [3,] 0 0 0 0 0 0 0 > > cat("\n\n________________________________________\n\n") ________________________________________ > cat("\n\npart H\n\n") part H > > > checkH <- diag(w) > > print(resjacphi[(n+m+dimmu+1):(n+2*m+2*dimmu), (n+1):(n+m+dimmu)] - checkH) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 [3,] 0 0 0 > > > cat("\n\n________________________________________\n\n") ________________________________________ > cat("\n\npart I\n\n") part I > > > checkI <- diag(mu) > > print(resjacphi[(n+m+dimmu+1):(n+2*m+2*dimmu), (n+m+dimmu+1):(n+2*m+2*dimmu)] - checkI) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 [3,] 0 0 0 > > checkjack <- rbind(cbind(checkA, checkB, checkC), cbind(checkD, checkE, checkF), cbind(checkG, checkH, checkI)) > > if(sum(abs(checkjack - resjacphi)) > .Machine$double.eps^(2/3)) + stop("wrong result") > > proc.time() user system elapsed 0.43 0.17 0.53