# source("test.mle.R")

if (EXTENDED.TESTING <- file.exists("source.R")) {
print("A")  
  source("source.R")
print("B")  
  source("mle.old.R")
}
print("C")  

if (EXTENDED.TESTING) {

try(load("rs"))

options(warn=1) ## 4
## sequence in `estparam' is
## mean, variance, nugget, scale, (+ further model parameters)
## So, mean, variance, and scale will be estimated here.
## Nugget is fixed and equals zero.
n <- 10
points <- 100
points <- 15
##
model <-"gencauchy"
x <- runif(points,0,3)
y <- runif(points,0,3) ## 100 random points in square [0, 3]^2
param <- c(0, 1, 0, 1, 1, 2)

estparam <- c(0, NA, 0, NA, NA, 2) ## NA means: "to be estimated"  
estparam <- param ## NA means: "to be estimated"
estparam[runif(length(estparam)) < 0.5] <- NA
print(estparam)
save(file="rs", .Random.seed)
d <- GaussRF(x=x, y=y, grid=FALSE, model=model, param=param, n=n)

if (FALSE)
  repeat{
  save(file="rs", .Random.seed)
  
  d <- GaussRF(x=x, y=y, grid=FALSE, model=model, param=param, n=n)
  str(fitvario(Print=3,
               x=cbind(x,y), data=d, model=model, param=estparam,
               lower=c(0.1, 0.1), upper=c(1.9, 5),
               lsq.methods=c("self", "plain", "sqrt.nr", "sd.inv"),
                                        #cross="cross.ign"
               )[c("variogram", "values")], vec.len=6, digits.d=5)
  str(try(Ofitvario(x=cbind(x,y), data=d, model=model, param=estparam,
                    lower=c(0.1, 0.1), upper=c(1.9, 5))[c("mle", "mle.value")]),
      vec.len=6, digits.d=5)
  print(estparam)
}


## The next two estimations give about the same result.
## For the first the sill is fixed to 1.5. For the second the sill
## is reached if the estimated variance is smaller than 1.5
estparam <-  c(0, NA, NA, NA, NA, NA) 
str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam, sill=1.5))
str(Ofitvario(x=cbind(x,y), data=d, model=model, param=estparam, sill=1.5))

estparam <-  c(0, NA, NaN, NA, NA, NA) 
parampositions(model=model, param=estparam)
f <- function(param) {
   param[5] <- max(0, 1.5 - param[1])
   return(param)
}
str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam, 
    sill=1, transform=f))
str(Ofitvario(x=cbind(x,y), data=d, model=model, param=estparam, 
    sill=1, transform=f))

## the next call gives a warning, since the user may programme
## strange things in this setup, and the program cannot check it.
estparam <- c(0, NA, NA, NA, NA, NaN) 
parampositions(model=model, param=estparam)
f <- function(param) {param[3] <- param[2]; param}
unix.time(str(fitvario(x=cbind(x,y), data=d, model=model,
      param=estparam, transform=f, standard.style=TRUE), vec.len=6))
unix.time(str(Ofitvario(x=cbind(x,y), data=d, model=model,
      param=estparam, transform=f, standard.style=TRUE), vec.len=6))

## much better programmed, but also much slower:
estmodel <- list(list(model="gencauchy", var=NA, scale=NA,
                      kappa=list(NA, function(m) m[[1]]$kappa[1])))
unix.time(str(fitvario(x=cbind(x,y), data=d, model=estmodel),
              vec.len=6))
unix.time(str(Ofitvario(x=cbind(x,y), data=d, model=estmodel),
              vec.len=6))

} # source exists