# Parameter Estimation Examples

# 1D
f <- function(x) (x-1)^2 + 1
optimize(f, c(-10,10))

# 2D
fsimp <- function(x1,y1) (1-x1)^2 + 100*(y1 - x1^2)^2
fr <- function(x) {   ## Rosenbrock Banana function
  x1 <- x[1]
  x2 <- x[2]
  100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}
x <- seq(-2,2,by=.15)
y <- seq(-1,3,by=.15)
z <- outer(x,y,fsimp)
persp(x,y,z,phi=45,theta=-45,col="yellow",shade=.00000001,ticktype="detailed")
optim(c(0,0), fr)
optim(c(10,10), fr)

# Lactate measurements
data <- read.csv('Stufentest_laktat_leistung.csv', header = FALSE, sep = ';')
subj <- levels(data$V1)
subj[1] <- "Power"
data <- data.table::transpose(data[,2:ncol(data)])
colnames(data) <- subj

Lactate <-  data$`09PB`
Lactate <- Lactate[!is.na(Lactate)]
Power = data$Power[1:length(Lactate)]

plot(Power, Lactate)
lines(Power, Lactate)

# 2(3) lines
lin3 <- function(x, power)
{
  p <- x[1]
  a1 <- x[2]
  a2 <- x[3]
  d1 <- x[4]
  # a3 <- x[4]
  # d1 <- x[5]
  # d2 <- x[6]
  la <- numeric(length(power))
  for (i in 1:length(power))
  {
    if (power[i] <= d1)
    {
      la[i] <- p + a1 * power[i];
    }
    # else if (power[i] <= d2)
    else
    {
      la[i] <- p + a1 * d1 + a2 * (power[i] - d1)
    }
    # else
    # {
    #   la[i] <- p + a1 * d1 + a2 * (d2 - d1) + a3 * (power[i] - d2)
    # }
  }
  return(la)
}

fun <- function(x, power, lactate)
{
  return(mean((lin3(x, power)-lactate)^2)) # least squares fit
}

x0 <- c(Lactate[1], 0, 0, Power[3])
lower <- c(Lactate[1]-1, -1, -1, Power[3])
upper <- c(2*Lactate[1], 1, 10, Power[length(Power-1)])
# x0 <- c(Lactate[1], 0, 0, 0, Power[2], Power[3])
# lower <- c(Lactate[1]-1, -1, -1, -1, Power[2], Power[3])
# upper <- c(2*Lactate[1], 1, 5, 10, Power[length(Power-2)], Power[length(Power-1)])

# optim: general purpose optimization (no modern methods)
res <- optim(x0, fun, power=Power, lactate = Lactate)
#res <- optim(x0, fun, power=Power, lactate = Lactate, lower = lower, upper = upper, method = 'L-BFGS-B')
lines(Power[1]:Power[length(Power)], lin3(res[[1]], Power[1]:Power[length(Power)]))

# nls: nolinear least squares
# model <- nls(Lactate ~ lin3(c(p, a1, a2, d1), Power), start = list(p = Lactate[1], a1 = 0, a2 = 0, d1 = Power[3]))

# simple linear fit: lm
# polynomial fit: pracma::polyfit (lm with polynomial function (poly))