## ----echo = FALSE, results = "hide", message = FALSE------------------------------------------------------------------
oldoptions <- options()
oldpar <- par()
options(width = 120)
require(salad)

## ---------------------------------------------------------------------------------------------------------------------
f1 <- function(x) sin(x**2)

## ---------------------------------------------------------------------------------------------------------------------
x <- dual(pi)
y <- f1(x)
y 
d(y)

## ---------------------------------------------------------------------------------------------------------------------
x <- dual(c(0, 1, sqrt(pi)))
y <- f1(x)
y
# get value and derivative
value(y)
d(y)

## ---------------------------------------------------------------------------------------------------------------------
x <- dual( matrix( c(1, 2, 4, 7), 2, 2))
x

## ---------------------------------------------------------------------------------------------------------------------
varnames(x)
# derivative along x1.1
d(x, "x1.1")

## ---------------------------------------------------------------------------------------------------------------------
y <- x %*% x
y
d(y, "x1.1")

## ---------------------------------------------------------------------------------------------------------------------
det(x)
d(det(x))

## ---------------------------------------------------------------------------------------------------------------------
z <- solve(x)
z
d(z, "x1.1")

## ---------------------------------------------------------------------------------------------------------------------
f2 <- function(x) {
  a <- x**(1:2)
  b <- ifelse(a > 1, sin(a), 1 - a) 
  C <- crossprod( cbind(a,b) )
  apply(C, 2, function(x) sum(x^2))
}

## ---------------------------------------------------------------------------------------------------------------------
# creating a dual number for x = 0.2
x <- dual(0.2)
y <- f2(x)
y
# get value and the derivative 
value(y)
d(y)

## ---------------------------------------------------------------------------------------------------------------------
a <- dual(c(1,2), dx = list("x" = c(1,1)))
b <- dual(c(2,1), dx = list("y" = c(2,1)))

## ---------------------------------------------------------------------------------------------------------------------
a + b
d(a + b)

## ---------------------------------------------------------------------------------------------------------------------
a <-  dual(c(1,2), dx = list("x" = c(1,1), "y" = c(0,0)))
b <-  dual(c(2,1), dx = list("x" = c(0,0), "y" = c(2,1)))

## ---------------------------------------------------------------------------------------------------------------------
a + b
d(a + b)

## ---------------------------------------------------------------------------------------------------------------------
v <- dual( c(1,1), varnames = c("x", "y"))
v
d(v)

## ---------------------------------------------------------------------------------------------------------------------
x <- v[1]
y <- v[2]

## ---------------------------------------------------------------------------------------------------------------------
a <- c(x, x+1)
d(a)

## ---------------------------------------------------------------------------------------------------------------------
b <- c(2*y, y)
d(b)

## ---------------------------------------------------------------------------------------------------------------------
a + b
d(a + b)

## ---------------------------------------------------------------------------------------------------------------------
x <- dual( matrix( c(1, 2, 4, 7), 2, 2))
as.vector(x)

## ---------------------------------------------------------------------------------------------------------------------
dim(x) <- NULL
x

## ---------------------------------------------------------------------------------------------------------------------
x <- dual(0) + c(-1,0,1)
x
d(x)
abs(x)
d(abs(x))

## ---------------------------------------------------------------------------------------------------------------------
y <- max( dual(c(1, 2)) )
y
d(y)

## ---------------------------------------------------------------------------------------------------------------------
y <- max( dual(c(2, 2)) )
y
d(y)

## ---------------------------------------------------------------------------------------------------------------------
f <- function(x) x**2 + x + 1
x <- dual(4)
f(x)
d(f(x))

## ---------------------------------------------------------------------------------------------------------------------
f1 <- dualFun1(f, \(x) 2*x + 1)
f1(x)
d(f1(x))

## ---------------------------------------------------------------------------------------------------------------------
system.time( for(i in 1:500) f(x) )
system.time( for(i in 1:500) f1(x) )