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
> # ?GNE
> 
> xmax <- 10
> 
> x <- seq(-xmax, xmax, length=31)
> 
> 
> a3Dplot <- function(x, phi, main, ...)
+ {
+ 	if(!is.function(phi))
+ 		stop("wrong argument phi.")
+ z <- outer(x, x, phi, ...)
+ nrz <- nrow(z)
+ ncz <- ncol(z)
+ # Create a function interpolating colors in the range of specified colors
+ jet.colors <- colorRampPalette( c("red", "orange", "yellow", "green","blue", "blueviolet", "magenta", "pink") ) 
+ # Generate the desired number of colors from this palette
+ nbcol <- 100
+ color <- jet.colors(nbcol)
+ # Compute the z-value at the facet centres
+ zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
+ # Recode facet z-values into color indices
+ facetcol <- cut(zfacet, nbcol)
+ par(bg = "white", mar=c(2, 2, 2, 2))
+ persp(x, x, z, col=color[facetcol], phi=30, theta=65, ticktype="detail",
+ 	main=main, xlab="a", ylab="b", zlab="phi(a,b)")
+ 
+ }
> 
> par(mfrow=c(1, 3))
> a3Dplot(x, phiFB, main="phiFB, lam->2")
> a3Dplot(x, phiKK, lambda=.5, main="phiKK, lam=.5")
> a3Dplot(x, phiKK, lambda=.1, main="phiKK, lam=.1")
> 
> par(mfrow=c(1, 3))
> a3Dplot(x, phiFB, main="phiFB(a,b)")
> a3Dplot(x, GrAphiFB, main="Grad a phiFB(a,b)")
> a3Dplot(x, GrBphiFB, main="Grad b phiFB(a,b)")
> 
> sum( outer(x, x, GrAphiFB) <= 0 )/(31^2)
[1] 1
> 
> max(abs(outer(x, x, GrAphiFB)))
[1] 2
> 
> 
> par(mfrow=c(1, 3))
> a3Dplot(x, phiFB, main="phiFB(a,b)")
> a3Dplot(x, phi=phipFB, main="phipFB(a,b), p=1/4", p=1/4)
> a3Dplot(x, phi=phipFB, main="phipFB(a,b), p=1", p=1)
> 
> 
> par(mfrow=c(1, 3))
> a3Dplot(x, GrAphiFB, main="phiFB(a,b)")
> a3Dplot(x, phi=GrAphipFB, main="phipFB(a,b), p=1/4", p=1/4)
> a3Dplot(x, phi=GrAphipFB, main="phipFB(a,b), p=1", p=1)
> 
> 
> 
> proc.time()
   user  system elapsed 
   0.39    0.23    0.62