require(robustbase)
## see also ./lmrob-psifns.R <<<<<<<<
source(system.file("xtraR/plot-psiFun.R", package = "robustbase", mustWork=TRUE))

EQ <- function(x,y) all.equal(x,y, tolerance = 1e-13)

## Demonstrate that  one of  tukeyChi() / tukeyPsi1() is superfluous
x <- seq(-4,4, length=201)
suppressWarnings(## as tukeyPsi1(), tukeyChi() are deprecated
for(c. in c(0.1, 1:2, pi, 100)) {
    ix <- abs(x) != c.
    stopifnot(EQ(tukeyChi(x, c.),
		 6/c.^2* tukeyPsi1(x, c., deriv=-1)),
	      EQ(tukeyChi(x, c., deriv= 1),
		 6/c.^2* tukeyPsi1(x, c., deriv= 0)),
	      EQ(tukeyChi(x, c., deriv= 2),
		 6/c.^2* tukeyPsi1(x, c., deriv= 1)),
	      ## Now show equivalence with Mpsi():
	      EQ(tukeyPsi1(x,     c.),      Mpsi(x,     c., "tukey")),
	      EQ(tukeyPsi1(x,     c., d=1), Mpsi(x,     c., "tukey", d=1)),
	      EQ(tukeyPsi1(x[ix], c., d=2), Mpsi(x[ix], c., "tukey", d=2))
	      )
}
)
## Test if default arguments are used
h2Psi <- chgDefaults(huberPsi, k = 2)

x <- 1:10
stopifnot(h2Psi@ rho(x, k=2) == h2Psi@ rho(x),
          h2Psi@ psi(x, k=2) == h2Psi@ psi(x),
          h2Psi@Dpsi(x, k=2) == h2Psi@Dpsi(x),
          h2Psi@ wgt(x, k=2) == h2Psi@ wgt(x),
          h2Psi@Dwgt(x, k=2) == h2Psi@Dwgt(x))

## Test default arguments for E... slots
stopifnot(EQ(h2Psi@Erho (), 0.49423127328548),
          EQ(h2Psi@Epsi2(), 0.920536925636323),
          EQ(h2Psi@EDpsi(), 0.954499736103642))

stopifnot(EQ(1, huberPsi@psi(1, k = 1e16)),
          huberPsi@wgt(0.1591319494080224, 0.5 + 1/13) <= 1)
## both used to fail because of numeric instability in pmin2/pmax2

f1 <- function(.) rep.int(1, length(.))
F1 <- function(x, .) rep.int(1, length(x))
## correct "classical psi":
cPs <- psiFunc(rho = function(x,.) x^2 / 2, psi = function(x, .) x,
               wgt = F1, Dpsi = F1, Erho = function(.) rep.int(1/2, length(.)),
               Epsi2 = f1, EDpsi = f1, . = Inf)
validObject(cPs); cPs
## incorrect dummy psi
cP <- psiFunc(rho = F1, psi = F1, wgt = F1, Dpsi = F1,
              Erho = f1, Epsi2 = f1, EDpsi = f1, . = Inf)
cP
## Check the autogenerated  Dwgt():
x <- seq(0,2, by=1/4)
stopifnot(## strict symmetry { including Dwgt(0) == 0 } :
	  huberPsi @Dwgt(-x) == -huberPsi @Dwgt(x),
	  hampelPsi@Dwgt(-x) == -hampelPsi@Dwgt(x),
	  huberPsi @Dwgt(x)[x < 1.345] == 0,
	  hampelPsi@Dwgt(x)[x < 1.487] == 0,
	  EQ(huberPsi @Dwgt(x[x >= 1.5]),
	     c(-0.597777777777778, -0.439183673469388, -0.33625)),
	  EQ(hampelPsi@Dwgt(x[x >= 1.5]),
	     c(-0.660883932259397, -0.485547378802822, -0.371747211895911))
	  )

.defDwgt <- robustbase:::.defDwgt
(ddd <- .defDwgt(psi  = function(u, k) pmin.int(k, pmax.int(-k, u)),
		 Dpsi = function(u, k) abs(u) <= k))
stopifnot(is.function(ddd), names(formals(ddd)) == c("u","k"),
	  EQ(ddd(x, 1.345), huberPsi@Dwgt(x)))

## TODO: Provide some functionality of this as a Plot+Check function
## ----  and then call the function for all our  psiFunc objects (with different 'k')
kk <- c(1.5, 3, 8)
psiH.38 <- chgDefaults(hampelPsi, k = kk)
c1 <- curve(psiH.38@psi(x), -10, 10, n=512, col=2)
abline(h=0, v=0, lty=3, lwd=.5, col="gray25")
c2 <- curve(x * psiH.38@wgt(x), add=TRUE, n=512, col=adjustcolor("blue", .5), lwd=2)
title("psi_Hampel_(1.5, 3, 8)  :  psi(x) =  x * wgt(x)")
axis(1, at=kk, expression(k[1], k[2], k[3]), pos=0)
axis(2, at=kk[1], quote(k[1]), pos=0, las=1)
stopifnot(all.equal(c1,c2, tolerance= 1e-15))

r1 <- curve(psiH.38@rho(x), -10, 10, col=2,
            main = quote(rho(x) == integral(phi(t) * dt, 0, x)))
axis(1, at=kk, expression(k[1], k[2], k[3]), pos=0)
curve(psiH.38@psi(x), add=TRUE, n=512, col=adjustcolor("blue", .5), lwd=2)
abline(h=0, v=0, lty=3, lwd=.5, col="gray25")
## check  rho(x) = \int_0^x psi(x) dx  {slightly *more* than  rho' = psi !}
rhoH.38.int <- function(x) integrate(function(u) psiH.38@psi(u), 0, x, rel.tol=1e-10)$value
r2 <- curve(sapply(x, rhoH.38.int), add = TRUE,
            lwd=4, col=adjustcolor("red", 1/4))
## numerical integration == "formula" :
stopifnot(all.equal(r1,r2, tolerance=1e-10))

curve(psiH.38@Dpsi(x), -10, 10, n=512, col=2,
      main = quote(psi*minute(x)))
abline(h=0, v=0, lty=3, lwd=.5, col="gray25")

## check  rho'(x) = phi(x)  etc  {TODO: for all our psiFun.}
head(xx <- seq(-10, 10, length=1024))
FrhoH.38 <- splinefun(xx, rho.x <- psiH.38@rho (xx))
FpsiH.38 <- splinefun(xx, psi.x <- psiH.38@psi (xx))
F1psH.38 <- splinefun(xx, Dps.x <- psiH.38@Dpsi(xx))

curve(FpsiH.38(x, deriv=1), -10,10, n=512)
curve(F1psH.38, add=TRUE, col=4, n=512)
stopifnot(all.equal(FpsiH.38(xx, deriv=1), Dps.x,
                    tolerance = 0.02))# not better because of discontinuities

curve(FrhoH.38(x, deriv=1), -10,10, n=512)
curve(FpsiH.38, add=TRUE, col=4, n=512)
stopifnot(all.equal(FrhoH.38(xx, deriv=1), psi.x, tolerance = 1e-4))

E.norm <- function(FUN, tol=1e-12, ...) {
    integrate(function(x) FUN(x) * dnorm(x), -Inf, Inf,
              rel.tol=tol, ...)$value
}

##' asymptotic efficiency -- both integrate + "formula"(@Epsi, @EDpsi) version
aeff.P <- function(psiF, k, ...) {
    stopifnot(is(psiF, "psi_func"))
    if(!missing(k))
	psiF <- chgDefaults(psiF, k = k)
    ## E[ psi'(X) ] ^2	/  E[ psi(X) ^ 2 ] :
    c(int = E.norm(psiF@Dpsi, ...)^2 / E.norm(function(x) psiF@psi(x)^2, ...),
      form= psiF@EDpsi()^2 / psiF@Epsi2())
}


## Breakdown Point --- for redescenders only,
## both integrate + "formula"(@Erho) version
bp.P <- function(psiF, k, ...) {
    stopifnot(is(psiF, "psi_func"))
    if(!missing(k))
	psiF <- chgDefaults(psiF, k = k)
    if(!is.finite( rhoInf <- psiF@rho(Inf) ))
	stop("rho(Inf) is not finite: ", rhoInf)
    integ <- function(x) psiF@rho(x)
    c(int = E.norm(integ, ...), form= psiF@Erho()) / rhoInf
}

## Print & Check the result of  aeff.P() or bp.P()
chkP <- function(rp, tol = 1e-9) {
    print(rp)
    ae <- all.equal(rp[[1]], rp[[2]], tolerance=tol)
    if(isTRUE(ae)) invisible(rp) else stop(ae)
}

chkP(aeff.P(huberPsi))
chkP(aeff.P(huberPsi, k = 1.5))
chkP(aeff.P(huberPsi, k = 2))
chkP(aeff.P(huberPsi, k = 2.5))

chkP(aeff.P(hampelPsi))
chkP(aeff.P(hampelPsi, k = c(1.5, 3, 8)))
chkP(aeff.P(hampelPsi, k = c(2,   4, 8), tol=1e-10),# fails with tol=1e-11
     tol = 1e-4)

## Now works too:
chkP(bp.P(hampelPsi))
chkP(bp.P(hampelPsi, k = c(1.5, 3, 8)))
chkP(bp.P(hampelPsi, k = c(2,   4, 8)))


## test derivatives (adapted from ./lmrob-psifns.R)
head(x. <- seq(-5, 10, length=1501))
## [separate lines, for interactive "play": ]
stopifnot(chkPsiDeriv(plot(huberPsi, x.)))
## ToDo: improve accuracy of derivative check
stopifnot(chkPsiDeriv(plot(hampelPsi, x.), tol=c(1e-4, 1e-1)))