#
#  mincontrast.R
#
#  Functions for estimation by minimum contrast
#

##################  base ################################

mincontrast <- function(observed, theoretical, startpar,
                        ...,
                        ctrl=list(q = 1/4, p = 2, rmin=NULL, rmax=NULL),
                        fvlab=list(label=NULL, desc="minimum contrast fit"),
                        explain=list(dataname=NULL, modelname=NULL, fname=NULL)) {
  verifyclass(observed, "fv")
  stopifnot(is.function(theoretical))
  if(!any("par" %in% names(formals(theoretical))))
    stop(paste("Theoretical function does not include an argument called",
               sQuote("par")))

  # enforce defaults
  ctrl <- resolve.defaults(ctrl, list(q = 1/4, p = 2, rmin=NULL, rmax=NULL))
  fvlab <- resolve.defaults(fvlab,
                            list(label=NULL, desc="minimum contrast fit"))
  explain <- resolve.defaults(explain,
                              list(dataname=NULL, modelname=NULL, fname=NULL))
  
  # determine range of r values
  rmin <- ctrl$rmin
  rmax <- ctrl$rmax
  if(!is.null(rmin) && !is.null(rmax))
    stopifnot(rmin < rmax && rmin >= 0)
  else {
    alim <- attr(observed, "alim")
    if(is.null(rmin)) rmin <- alim[1]
    if(is.null(rmax)) rmax <- alim[2]
  }
  # extract vector of r values
  argu <- attr(observed, "argu")
  rvals <- observed[[argu]]
  # extract vector of observed values of statistic
  valu <- attr(observed, "valu")
  obs <- observed[[valu]]
  # restrict to [rmin, rmax]
  if(max(rvals) < rmax)
    stop(paste("rmax=", signif(rmax,4),
               "exceeds the range of available data",
               "= [", signif(min(rvals),4), ",", signif(max(rvals),4), "]"))
  sub <- (rvals >= rmin) & (rvals <= rmax)
  rvals <- rvals[sub]
  obs <- obs[sub]
  # for efficiency
  obsq <- obs^(ctrl$q)
  # define objective function
  objective <- function(par, obsq, theoretical, rvals, qq, pp, rmin, rmax, ...) {
    theo <- theoretical(par=par, rvals, ...)
    if(!is.vector(theo) || !is.numeric(theo))
      stop("theoretical function did not return a numeric vector")
    if(length(theo) != length(obs))
      stop("theoretical function did not return the correct number of values")
    discrep <- (abs(theo^qq - obsq))^pp
    return(sum(discrep))
  }
  # go
  minimum <- optim(startpar, fn=objective,
                   obsq=obsq, theoretical=theoretical,
                   rvals=rvals,
                   qq=ctrl$q, pp=ctrl$p, rmin=rmin, rmax=rmax, ...)
  # evaluate the fitted theoretical curve
  fittheo <- theoretical(minimum$par, rvals, ...)
  # pack it up as an `fv' object
  label <- fvlab$label
  desc  <- fvlab$desc
  if(is.null(label))
    label <- paste("fit(", argu, ")", collapse="")
  fitfv <- bind.fv(observed[sub, ],
                   data.frame(fit=fittheo),
                   label, desc)
  result <- list(par=minimum$par,
                 fit=fitfv,
                 opt=minimum,
                 ctrl=list(p=ctrl$p,q=ctrl$q,rmin=rmin,rmax=rmax),
                 info=explain
                 )
  
  class(result) <- c("minconfit", class(result))
  return(result)
}

print.minconfit <- function(x, ...) {
  # explanatory
  cat(paste("Minimum contrast fit ",
            "(",
            "object of class ",
            dQuote("minconfit"),
            ")",
            "\n", sep=""))
  mo <- x$info$modelname
  fu <- x$info$fname
  da <- x$info$dataname
  if(!is.null(mo))
    cat(paste("Model:", mo, "\n"))
  if(!is.null(fu) && !is.null(da))
    cat(paste("Fitted by matching theoretical", fu, "function to", da))
  else {
    if(!is.null(fu))
      cat(paste(" based on", fu))
    if(!is.null(da))
      cat(paste(" fitted to", da))
  }
  cat("\n")
  # Values
  cat("Parameters fitted by minimum contrast ($par):\n")
  print(x$par, ...)
  mp <- x$modelpar
  if(!is.null(mp)) {
    cat(paste("Derived parameters of",
              if(!is.null(mo)) mo else "model",
              "($modelpar):\n"))
    print(mp)
  }
  # Diagnostics
  cgce <- x$opt$convergence
  switch(paste(cgce),
         "0"={
           cat(paste("Converged successfully after",
                     x$opt$counts[["function"]],
                     "iterations.\n"))
         },
         "1"={
           cat("Warning: iteration limit maxit was reached.\n")
         },
         "10"={
           cat("Warning: Nelder-Mead simplex was degenerate\n")
         },
         "51"={
           cat(paste("Warning from L-BGFS-B method:\n",
                     x$opt$message, "\n"))
         },
         "52"={
           cat(paste("Error message from L-BGFS-B method:\n",
                     x$opt$message, "\n"))
         },
         cat(paste("Unrecognised error code", cgce, "\n"))
         )
  # Algorithm parameters
  ct <- x$ctrl
  cat(paste("Domain of integration:",
            "[",
            signif(ct$rmin,4),
            ",",
            signif(ct$rmax,4),
            "]\n"))
  cat(paste("Exponents:",
            "p=", paste(signif(ct$p, 3), ",",  sep=""),
            "q=", signif(ct$q,3), "\n"))
  invisible(NULL)
}
              

plot.minconfit <- function(x, ...) {
  xname <- deparse(substitute(x))
  do.call("plot.fv",
          resolve.defaults(list(x$fit),
                           list(...),
                           list(main=xname)))
}


############### applications (specific models) ##################


thomas.estK <- function(X, startpar=c(kappa=1,sigma2=1),
                        lambda=NULL, q=1/4, p=2, rmin=NULL, rmax=NULL){

  dataname <- deparse(substitute(X))

  if(inherits(X, "fv")) {
    K <- X
  } else if(inherits(X, "ppp")) {
    K <- Kest(X)
    dataname <- paste("Kest(", dataname, ")", sep="")
    lambda <- summary(X)$intensity
  } else 
    stop("Unrecognised format for argument X")

  startpar <- check.named.vector(startpar, c("kappa","sigma2"))

  theoret <- function(par,rvals){
    if(any(par <= 0))
      return(rep(Inf, length(rvals)))
    pi*rvals^2+(1-exp(-rvals^2/(4*par[2])))/par[1]
  }
  result <- mincontrast(K, theoret, startpar,
                        ctrl=list(q=q, p=p,rmin=rmin, rmax=rmax),
                        fvlab=list(label="Kfit(r)",
                          desc="minimum contrast fit of Thomas process"),
                        explain=list(dataname=dataname,
                          fname="K", modelname="Thomas process"))
  # imbue with meaning
  par <- result$par
  names(par) <- c("kappa", "sigma2")
  result$par <- par
  # infer model parameters
  mu <- if(is.numeric(lambda) && length(lambda) == 1)
           lambda/result$par[["kappa"]] else NA
  result$modelpar <- c(kappa=par[["kappa"]],
                       sigma=sqrt(par[["sigma2"]]),
                       mu=mu)
  result$internal <- list(model="Thomas")
  return(result)
}


lgcp.estK <- function(X, startpar=list(sigma2=1,alpha=1),
                      lambda=NULL, q=1/4, p=2, rmin=NULL, rmax=NULL) {
  
  dataname <- deparse(substitute(K))
  if(inherits(X, "fv")) {
    K <- X
  } else if(inherits(X, "ppp")) {
    K <- Kest(X)
    dataname <- paste("Kest(", dataname, ")", sep="")
    lambda <- summary(X)$intensity
  } else 
    stop("Unrecognised format for argument X")

  startpar <- check.named.vector(startpar, c("sigma2","alpha"))

  integrand<-function(r,par){
    2*pi*r*exp(par[1]*exp(-r/par[2]))
  }
  theoret <- function(par, rvals) {
    if(any(par <= 0))
      return(rep(Inf, length(rvals)))
    th <- numeric(length(rvals))
    th[1] <- if(rvals[1] == 0) 0 else 
             integrate(integrand,lower=0,upper=rvals[1],par=par)$value
    for (i in 2:length(rvals))
      th[i]=th[i-1]+integrate(integrand,lower=rvals[i-1],upper=rvals[i],par=par)$value
    return(th)
  }
  result <- mincontrast(K, theoret, startpar,
                        ctrl=list(q=q, p=p, rmin=rmin, rmax=rmax),
                        fvlab=list(label="Kfit(r)",
                          desc="minimum contrast fit of LGCP"),
                        explain=list(dataname=dataname,
                          fname="K",
                          modelname="log-Gaussian Cox process"))
  # imbue with meaning
  par <- result$par
  names(par) <- c("sigma2", "alpha")
  result$par <- par
  # infer model parameters
  mu <- if(is.numeric(lambda) && length(lambda) == 1 && lambda > 0)
           log(lambda) - par[["sigma2"]]/2 else NA
  result$modelpar <- c(sigma2=par[["sigma2"]],
                       alpha=par[["alpha"]],
                       mu=mu)
  result$internal <- list(model="lcgp")
  return(result)
}


matclust.estK <- function(X, startpar=c(kappa=1,R=1),
                        lambda=NULL, q=1/4, p=2, rmin=NULL, rmax=NULL){

  dataname <- deparse(substitute(X))

  if(inherits(X, "fv")) {
    K <- X
  } else if(inherits(X, "ppp")) {
    K <- Kest(X)
    dataname <- paste("Kest(", dataname, ")", sep="")
    lambda <- summary(X)$intensity
  } else 
    stop("Unrecognised format for argument X")

  startpar <- check.named.vector(startpar, c("kappa","R"))

  hfun <- function(zz) {
    ok <- (zz < 1)
    h <- numeric(length(zz))
    h[!ok] <- 1
    z <- zz[ok]
    h[ok] <- 2 + (1/pi) * (
                           (8 * z^2 - 4) * acos(z)
                           - 2 * asin(z)
                           + 4 * z * sqrt((1 - z^2)^3)
                           - 6 * z * sqrt(1 - z^2)
                           )
    return(h)
  }
  theoret <- function(par,rvals){
    if(any(par <= 0))
      return(rep(Inf, length(rvals)))
    pi * rvals^2 + (1/par[1]) * hfun(rvals/(2 * par[2]))
  }
  result <- mincontrast(K, theoret, startpar,
                        ctrl=list(q=q, p=p,rmin=rmin, rmax=rmax),
                        fvlab=list(label="Kfit(r)",
                          desc="minimum contrast fit of Matern Cluster process"),
                        explain=list(dataname=dataname,
                          fname="K", modelname="Matern Cluster process"))
  # imbue with meaning
  par <- result$par
  names(par) <- c("kappa", "R")
  result$par <- par
  # infer model parameters
  mu <- if(is.numeric(lambda) && length(lambda) == 1)
           lambda/result$par[["kappa"]] else NA
  result$modelpar <- c(kappa=par[["kappa"]],
                       R=par[["R"]],
                       mu=mu)
  result$internal <- list(model="MatClust")
  return(result)
}


check.named.vector <- function(x, nam) {
  xname <- deparse(substitute(x))
  if(!is.numeric(x) || length(x) != 2 || !all(nam %in% names(x))) {
    whinge <- paste(xname,
                    "must be a named vector, with components",
                    commasep(nam))
    stop(whinge)
  }
  return(x[nam])
}