Epanechnikov.kernel  <- function(x){
  res <- ((3/4)*(1-x^2))*(abs(x)<=1)
  return(res)
}

Epanechnikov.kernel.2d  <- function(x,y){
  res <- (9/16)*(1-x^2)*(1-y^2)*(abs(x)<=1)*(abs(y)<=1)
  return(res)
}

mu.estimate <- function(h, xi, yi, times = sort(unique(xi))){
  mu.hat <- c()
  for(x in times){
    X <- x - xi
    W <- Epanechnikov.kernel((x - xi)/h)
    if (length(W[which(W!=0)]) == 0){
      # zvysit h
      h.new <- h
      while (length(W[which(W!=0)]) == 0) {
        h.new <- h.new*1.3
        print(paste(c("zvysuji na",h.new)))
        W <- Epanechnikov.kernel((x - xi)/h.new)
      }
    }
    model <- lm(yi ~ X, weights = W)
    mu.hat <- c(mu.hat,model$coefficients[1])
  }
  return(mu.hat)
}

cov.hat.fda <- function(tecka, ypska, munbasis = 14, covnbasis = 14, varnbasis = munbasis, lambda = 1e4){
  # vypocte odhad kovariancniho operatoru pomoci splajnoveho vyhlazovani
  # tecka, ypska ... listy pro t a y
  # munbasis .... pocet bazovych splajnu pro odhad str. hodnoty
  # covnbasis ... pocet bazovych splajnu pro odhad kovariancniho op.
  n <- length(tecka)
  BTime <- unlist(tecka)
  BAmount <- unlist(ypska)
  rknots <- range(BTime)
  bbasis <- create.bspline.basis(rknots, nbasis = munbasis, norder = 4)
  
  palmfd <- smooth.basis(BTime, BAmount, bbasis)
  #plot(palmfd$fd)
  errs <- BAmount - eval.fd(BTime, palmfd$fd)
  
  # Setting Up Bivariate Smoothing for one Auction
  AuctionErrs <- list()
  covtt <- data.frame()
  Xmat <- 0
  Ymat <- 0
  Auction <- vector("double")
  for(i in 1:n){
    kus <- rep(i, length(tecka[[i]]))
    Auction <- c(Auction, kus)
  }
  covbasis <- create.bspline.basis(rknots, nbasis = covnbasis, norder = 4)
  
  for(i in 1:n){
    AuctionErrs[[i]] <- cbind(BTime[Auction == i], errs[Auction == i])
    bvals <- eval.basis(AuctionErrs[[i]][,1], covbasis) %x%
      eval.basis(AuctionErrs[[i]][,1], covbasis)
    bigerrs <- AuctionErrs[[i]][,2]%x%AuctionErrs[[i]][,2]
    Xmat <- Xmat + t(bvals)%*%bvals
    Ymat <- Ymat + t(bvals)%*%bigerrs
    covtt <- rbind(covtt, data.frame(x = BTime[Auction == i], y = errs[Auction == i]*errs[Auction == i]))
  }
  
  R1 <- eval.penalty(covbasis,0)%x%eval.penalty(covbasis,2)
  R2 <- eval.penalty(covbasis,2)%x%eval.penalty(covbasis,0)
  
  Bigmat <- Xmat #t(Xmat)%*%Xmat
  nbasis <- covnbasis
  coefs <- solve(Bigmat + lambda*(R1+R2), Ymat)
  coefs <- matrix(coefs, nbasis, nbasis)
  covfd <- bifd(coefs, covbasis, covbasis)
  t0 <- seq(rknots[1], rknots[2], length.out = 50)
  # odhad rozptylu
  bbasis <- create.bspline.basis(rknots, nbasis = varnbasis, norder = 4)
  varfd <- smooth.basis(covtt$x, covtt$y, bbasis)
  cov.tilde <- eval.fd(t0, varfd$fd)
  temp.cov <- eval.bifd(t0, t0, covfd)
  cov.hat <- matrix(temp.cov, ncol = length(t0))
  var.hat <- cov.tilde - diag(cov.hat)
  # jen pro vystup
  bbasis <- create.bspline.basis(rknots, breaks = t0, norder = 2)
  sigmafd <- smooth.basis(t0, c(var.hat), bbasis)
  Z <- list(t = t0, mu.hat = eval.fd(t0, palmfd$fd), cov.hat = cov.hat + diag(c(var.hat)),
            mu.hat.fd = palmfd$fd, cov.hat.fd = covfd, var.hat.fd = sigmafd$fd)
  return(Z)
}

cov.hat.kok <- function(tecka, ypska, lambda1 = 0.1, lambda2 = 0.2){
  # podle knihy Kokoszka
  x <- unlist(tecka)
  y <- unlist(ypska)
  library(mgcv)
  gmf <- gam(y ~ s(x, sp = lambda1))

  N <- length(tecka)
  mu.hat <- gmf$fitted.values
  errs <- y - mu.hat
  
  covtt <- data.frame()
  Xmat <- vector("double")
  Ymat <- vector("double")
  Auction <- vector("double")
  for(i in 1:N){
    kus <- rep(i, length(tecka[[i]]))
    Auction <- c(Auction, kus)
  }
  # vytvoreni dat, na ktere se aplikuje gam odhad
  for(i in 1:N){
    tt <- x[Auction == i]
    yy <- errs[Auction == i]
    comb.times <- combn(tt, 2)
    comb.ind <- combn(1:length(tt), 2)
    comb.y <- yy[comb.ind[1,]] * yy[comb.ind[2,]]
    Xmat <- cbind(Xmat, comb.times)
    Ymat <- c(Ymat, comb.y)
    covtt <- rbind(covtt, data.frame(x = tt, y = yy*yy))
  }
  
  cov.hat.val <- Ymat
  tt.val <- t(Xmat)
  
  x1 <- tt.val[, 1]
  x2 <- tt.val[, 2]
  y <- cov.hat.val
  rozsah <- range(c(x1, x2))
  tt0 <- seq(rozsah[1], rozsah[2], length.out = 50)
  tt0.df <- merge(tt0, tt0)
  gmsig <- gam(y ~ s(x1, sp = lambda2) + s(x2, sp = lambda2))
  zz <- predict(gmsig, data.frame(x1 = tt0.df[, 1], x2 = tt0.df[, 2]))
  cov.hat <- matrix(zz, ncol = length(tt0))
  # odhad rozptylu
  x <- covtt$x
  y <- covtt$y
  gmvar <- gam(y ~ s(x, sp = lambda1))
  cov.tilde <- predict(gmvar, data.frame(x = tt0))
  var.hat <- cov.tilde - diag(cov.hat)  
  
  Z <- list(t = tt0, mu.hat = predict(gmf, data.frame(x = tt0)), cov.hat = cov.hat + diag(c(var.hat)))
  return(Z)
}

cov.hat.ks <- function(tecka, ypska, h.mu = (range(tecka)[2] - range(tecka)[1])/(sqrt(length(unlist(tecka)))), 
                       h.cov = (range(tecka)[2] - range(tecka)[1])/((length(unlist(tecka)))^(1/4)), 
                       t0 = seq(range(tecka)[1], range(tecka)[2], length.out = 50)){
  # vypocte odhad kovariancniho operatoru pomoci jadroveho vyhlazovani
  # ------------------------
  # Diagonalu odhaduje zvlast
  # ------------------------
  # tecka, ypska ... listy pro t a y
  n <- length(tecka)
  BTime <- unlist(tecka)
  BAmount <- unlist(ypska)
  mu.hat <- mu.estimate(h = h.mu, BTime, BAmount, BTime)
  errs <- BAmount - mu.hat
  
  AuctionErrs <- list()
  covtt <- data.frame()
  Xmat <- vector("double")
  Ymat <- vector("double")
  Auction <- vector("double")
  for(i in 1:n){
    kus <- rep(i, length(tecka[[i]]))
    Auction <- c(Auction, kus)
  }
  # vytvoreni dat, na ktere se aplikuje jadrovy odhad
  for(i in 1:n){
    times <- BTime[Auction == i]
    y <- errs[Auction == i]
    comb.times <- combn(times, 2)
    comb.ind <- combn(1:length(times), 2)
    comb.y <- y[comb.ind[1,]] * y[comb.ind[2,]]
    Xmat <- cbind(Xmat, comb.times)
    Ymat <- c(Ymat, comb.y)
    covtt <- rbind(covtt, data.frame(x = times, y = y*y))
  }
  n.t0 <- length(t0)
  dens <- density(BTime, from = t0[1], to = t0[n.t0], n = n.t0)
  alfa <- dens$y^(1/5)
  cov.f <- matrix(nrow = length(t0), ncol = length(t0))
  for(s in 1:length(t0)){
    for(t in s:length(t0)){
      temp <- Xmat - c(t0[s],t0[t])
      X <- data.frame(x1 = temp[1,], x2 = temp[2,], y = Ymat)
      W <- alfa[s]*alfa[t]*Epanechnikov.kernel.2d(alfa[s]*X$x1/h.cov, alfa[t]*X$x2/h.cov)
      if (length(W[which(W!=0)]) == 0){
        # zvysit h.cov
        h.new <- h.cov
        while (length(W[which(W!=0)]) == 0) {
          h.new <- h.new*1.3
          print(paste(c("zvysuji na",h.new)))
          W <- alfa[s]*alfa[t]*Epanechnikov.kernel.2d(alfa[s]*X$x1/h.new, alfa[t]*X$x2/h.new)
        }
      }
      model <- lm(y ~ x1 + x2, weights = W[which(W!=0)], data = X[which(W!=0), ])
      cov.f[s,t] <- model$coefficients[1]
      cov.f[t,s] <- cov.f[s,t]
    }
  }
  # odhad rozptylu
  cov.tilde <- mu.estimate(h = h.mu, covtt$x, covtt$y, t0)
  cov.hat <- cov.f
  var.hat <- cov.tilde - diag(cov.hat)  
  Z <- list(t = t0, mu.hat = mu.estimate(h = h.mu, BTime, BAmount, t0), cov.hat = cov.hat, var.hat = var.hat)
  return(Z)
}

data.fpca <- function(tecka, ypska, mu.hat, Cmat, var.hat,
                         t0 = seq(range(tecka)[1], range(tecka)[2], length.out = length(mu.hat)), K = 6){
  N <- length(tecka)
  if (N!=length(ypska)){stop("Lengths differ!")}

  # Setting Up fd object for mu.hat
  bbasis <- create.bspline.basis(range(t0), breaks = t0, norder = 2)
  mu.hat.fd <- smooth.basis(t0, mu.hat, bbasis)
  
  # Setting Up fd object for var.hat
  if (!hasArg(var.hat)){
    var.hat <- diag(Cmat)
  }
  if (is.fd(var.hat)){
    var.hat.fd <- var.hat
  }else{
    var.hat.fd <- smooth.basis(t0, var.hat, bbasis)$fd
  }
  
  # if (min(var.hat) < 0){
  #   print("var.hat is negative!")
  #   plot(var.hat)
  #   var.hat <- abs(var.hat)
  # }
  # Setting Up bifd object for Cmat
  if (class(Cmat) == "bifd"){
    cov.hat.fd <- Cmat
    Cmat <- eval.bifd(t0, t0, cov.hat.fd)
  }else{
    covbasis <- create.bspline.basis(range(t0), breaks = t0, norder = 2)
    bvals <- eval.basis(t0, covbasis)%x%eval.basis(t0, covbasis)
    Xmat <- t(bvals)%*%bvals
    Ymat <- t(bvals)%*%matrix(c(Cmat), ncol = 1)
    nbasis <- covbasis$nbasis
    coefs <- solve(Xmat, Ymat)
    coefs <- matrix(coefs, nbasis, nbasis)
    cov.hat.fd <- bifd(coefs, covbasis, covbasis)
  }
  ypska.hat <- matrix(nrow = length(t0), ncol = N)

  #Cmat <- (Cmat + t(Cmat))/2
  pca <- eigen(Cmat)
  U <- pca$vectors
  sumvecc <- apply(U, 2, sum)
  U[, sumvecc < 0] <- -U[, sumvecc < 0]
  Lambda <- pca$values[1:K]
  UK <- U[, 1:K]
  rknots <- range(tecka)
  bbasis <- create.bspline.basis(range(t0), breaks = t0, norder = 2)
  harmfd <- smooth.basis(t0, UK, bbasis)
  
  # FPCA
  mu <- eval.fd(t0, mu.hat.fd$fd)
  for (k in 1:N){
    t <- tecka[[k]]
    y <- ypska[[k]]
    M <- eval.bifd(t, t, cov.hat.fd) + diag(c(eval.fd(t, var.hat.fd)))
    mu.k <- eval.fd(t, mu.hat.fd$fd)
    V <- eval.fd(t, harmfd$fd)
    sc.hat <- Lambda * t(V) %*% solve(M) %*% (y - mu.k)
    ypska.hat[,k] <- mu + eval.fd(t0, harmfd$fd) %*% sc.hat
  }
  return(ypska.hat)
}

cov.hat.bks <- function(tecka, ypska, h.mu = (range(tecka)[2] - range(tecka)[1])/(sqrt(length(unlist(tecka)))), 
                       h.cov = (range(tecka)[2] - range(tecka)[1])/((length(unlist(tecka)))^(1/4)), 
                       t0 = seq(range(tecka)[1], range(tecka)[2], length.out = 50)){
  # vypocte odhad kovariancniho operatoru pomoci jadroveho vyhlazovani
  # tecka, ypska ... listy pro t a y
  n <- length(tecka)
  BTime <- unlist(tecka)
  BAmount <- unlist(ypska)
  mu.hat <- mu.estimate(h = h.mu, BTime, BAmount, BTime)
  errs <- BAmount - mu.hat
  
  AuctionErrs <- list()
  covtt <- data.frame()
  Xmat <- vector("double")
  Ymat <- vector("double")
  Auction <- vector("double")
  for(i in 1:n){
    kus <- rep(i, length(tecka[[i]]))
    Auction <- c(Auction, kus)
  }
  # vytvoreni dat, na ktere se aplikuje jadrovy odhad
  for(i in 1:n){
    times <- BTime[Auction == i]
    y <- errs[Auction == i]
    comb.times <- combn(times, 2)
    comb.ind <- combn(1:length(times), 2)
    comb.y <- y[comb.ind[1,]] * y[comb.ind[2,]]
    Xmat <- cbind(Xmat, comb.times)
    Ymat <- c(Ymat, comb.y)
    covtt <- rbind(covtt, data.frame(x = times, y = y*y))
  }
  n.t0 <- length(t0)
  dens <- density(BTime, from = t0[1], to = t0[n.t0], n = n.t0)
  alfa <- dens$y^(1/5)
  # nebo
  #alfa <- (dens$y*(t0[n.t0] - t0[1])/n.t0)^(1/5)
  cov.f <- matrix(nrow = length(t0), ncol = length(t0))
  for(s in 1:length(t0)){
    for(t in s:length(t0)){
            if(s!=t){
      temp <- Xmat - c(t0[s],t0[t])
      X <- data.frame(x1 = temp[1,], x2 = temp[2,], y = Ymat)
      W <- alfa[s]*alfa[t]*Epanechnikov.kernel.2d(alfa[s]*X$x1/h.cov, alfa[t]*X$x2/h.cov)
      if (length(W[which(W!=0)]) == 0){
        # zvysit h.cov
        h.new <- h.cov
        while (length(W[which(W!=0)]) == 0) {
          h.new <- h.new*1.3
          print(paste(c("zvysuji na",h.new)))
          W <- alfa[s]*alfa[t]*Epanechnikov.kernel.2d(alfa[s]*X$x1/h.new, alfa[t]*X$x2/h.new)
        }
      }
      model <- lm(y ~ x1 + x2, weights = W[which(W!=0)], data = X[which(W!=0), ])
      cov.f[s,t] <- model$coefficients[1]
      cov.f[t,s] <- cov.f[s,t]
      }else{
        Xmat.rotated <- matrix(c(sqrt(2)/2, -sqrt(2)/2, sqrt(2)/2, sqrt(2)/2), nrow=2, byrow=T)%*%Xmat
        temp <- c(0, t0[t]*sqrt(2)) - Xmat.rotated
        X <- data.frame(x1 = temp[1,], x2 = temp[2,], y = Ymat)
        W <- alfa[s]^2*Epanechnikov.kernel.2d(alfa[s]*X$x1/(h.cov/2), alfa[t]*X$x2/(h.cov/2))
        X$x2 <- X$x2^2
        model <- lm(y ~ x1 + x2, weights = W[which(W!=0)], data = X[which(W!=0),])
        cov.f[s,t] <- model$coefficients[1]
      }
    }
  }
  # odhad rozptylu
  
  cov.tilde <- mu.estimate(h = h.mu, covtt$x, covtt$y, t0)
  cov.hat <- cov.f
  var.hat <- cov.tilde - diag(cov.hat)  
  
  # pca <- eigen(cov.f)
  # Lambda <- diag(c(pca$values[pca$values>0],rep(0,sum(pca$values<=0))))
  # cov.f.fitted <- pca$vectors%*%Lambda%*%t(pca$vectors)
  # cov.f.fitted <- cov.f
  Z <- list(t = t0, mu.hat = mu.estimate(h = h.mu, BTime, BAmount, t0), cov.hat = cov.hat, var.hat = var.hat)
  return(Z)
}