#
# fonction pour faire xvalid sur du GLSM
#

`xvalid.binom` <-
function (geodata, coords = geodata$coords, data = geodata$data, 
    model, units.m = "default",mcmc.input, variog.obj = NULL, 
    locations.xvalid = "all", data.xvalid = NULL, trend.data, trend.loc, messages, nwin="full", ask=T,borders=NULL, logitp.data,trend0=T,resid=F,  ...) 
{
    call.fc <- match.call()
    if (missing(messages)) 
        messages.screen <- as.logical(ifelse(is.null(getOption("geoR.messages")), 
            TRUE, getOption("geoR.messages")))
    else messages.screen <- messages
    if (any(class(model) == "eyefit")) {
        model <- model[[1]]
        model$method <- "eyefit"
    }
    if (missing(geodata)) 
        geodata <- list(coords = coords, data = data, units.m =units.m)
    if (any(units.m == "default")) {
        if (!is.null(geodata$units.m)) 
            units.m <- geodata$units.m
        else units.m <- rep(1, n)
    }
    n <- nrow(coords)
    data <- as.vector(data)
    if (length(data) != n) 
        stop("coords and data have incompatible sizes")
    if(trend0){
      xmat <- trend.spatial(trend = trend.data, geodata = geodata)
      if (nrow(xmat) != n) 
        stop("coords and trend have incompatible sizes")
    }
    else xmat=NULL
    if (missing(mcmc.input)) 
        stop("xvalid.binom : argument mcmc.input must be given")
    mcmc.input <- .mcmc.check.aux(mcmc.input, fct = "binom.krige")
#    if (is.null(model$aniso.pars)) 
#        model$aniso.pars <- c(0, 1)
    if (is.null(model$kappa)) 
        model$kappa <- 0.5
    if (all(locations.xvalid == "all") | is.vector(locations.xvalid)) {
        autocross <- TRUE
        if (all(locations.xvalid == "all")) 
            locations.xvalid <- 1:n
        else if (any(locations.xvalid > n) | mode(locations.xvalid) != 
            "numeric") 
            stop("\nxvalid: vector indicating locations to be validated is not a numeric vector and/or has element(s) with value greater than the number of data loccations")
        crossvalid <- TRUE
    }
    else {
        autocross <- FALSE
        if (is.matrix(locations.xvalid) | is.data.frame(locations.xvalid)) 
            if (dim(locations.xvalid)[2] <= 2) {
                if (dim(locations.xvalid)[2] == 1) {
                  locations.xvalid <- is.vector(locations.xvalid)
                  crossvalid <- TRUE
                  if (any(locations.xvalid) > n | length(locations.xvalid) > 
                    n) 
                    stop("incorrect value to the argument locations.xvalid.\nThis must be a numeric vector with numbers indicating the locations to be cross-validated")
                }
                else {
                  if (messages.screen) 
                    cat("xvalid: cross-validation to be performed on locations provided by the user\n")
                  if (is.null(data.xvalid)) 
                    stop("the argument \"data.xvalid\" must be provided in order to perform validation on a set of locations different from the original data")
                  crossvalid <- FALSE
                }
            }
            else if (!is.vector(locations.xvalid) | mode(locations.xvalid) != 
                "numeric") 
                stop("\nargument locations.xvalid must be either:\n a numeric vector with numbers indicating the locations to be cross-validated\n a matrix with coordinates for the locations to be cross-validated.")
            else if (any(locations.xvalid) > n | length(locations.xvalid) > 
                n) 
                stop("incorrect value to the argument locations.xvalid.\nThis must be a numeric vector with numbers indicating the locations to be cross-validated")
    }
    if (crossvalid == FALSE) 
        n.pt.xv <- dim(locations.xvalid)[[1]]
    else n.pt.xv <- length(locations.xvalid)
    if (messages.screen) {
        cat(paste("xvalid: number of data locations       =", 
            n))
        cat("\n")
        cat(paste("xvalid: number of validation locations =", 
            n.pt.xv))
        cat("\n")
        if (crossvalid) 
            cat("xvalid: performing cross-validation")
        else cat("xvalid: performing validation at the locations provided")
        cat("\n")

    }
    if (crossvalid) {

#
# pb d'allocation je mets le calcul qui était sous forme apply
#        res <- as.data.frame(t(apply(matrix(locations.xvalid),1, cv.f,data.mcmc=Tres.mcmc)))
#  dans une boucle
#
#
      if(!ask) Sscale=mcmc.input$S.scale
      repeat{
      if(ask){
	print("valeur pour S.scale (0 pour continuer)")
	Sscale=scan(n=1)
      }
	if(Sscale<=0) break
#	S.scale=0.15
	mcmc.input$S.scale=Sscale
	beta.prior <- "flat"
	Tres.mcmc=binom.krige.simul1(data, coords, units.m,mcmc.input, model$cov.model,model$kappa,model$beta,model$cov.pars,model$nugget,model$aniso.pars,trend.data,beta.prior,messages.screen,nwin=nwin)$Sdata
        if(!ask) Sscale=0
    }
	res=matrix(0,n.pt.xv,5)
	for(ndata in 1:n.pt.xv){
           if (messages.screen) 
                cat(paste(ndata, ", ", sep = ""))
            coords.out <- coords[ndata, , drop = FALSE]
            data.out <- data[ndata]/units.m[ndata]
            if(trend0) xmat.out <- xmat[ndata, , drop = FALSE]
	    if(nwin!="full"){
		cv.coords=coords[-ndata,]
		dm0 <- loccoords(coords = cv.coords, locations = coords.out)
		cv.coords <- cv.coords[order(dm0)[1:nwin], ]
		cv.data.mcmc = Tres.mcmc[-ndata,][order(dm0)[1:nwin],]
		if(trend0) cv.xmat <- xmat[-ndata, , drop = FALSE][order(dm0)[1:nwin], , drop = FALSE]
              }
	    else{
	         cv.coords <- coords[-ndata, ]
        	 cv.data.mcmc <- Tres.mcmc[-ndata,]
           	 if(trend0) cv.xmat <- xmat[-ndata, , drop = FALSE]
               }
           if(trend0){
	    Tkgc=krige.control(type.krige="ok",trend.d = ~cv.xmat + 0, trend.l = ~xmat.out + 0, cov.model = model$cov.model, beta=model$beta,cov.pars = model$cov.pars, nugget = model$nugget, kappa = model$kappa, lambda = model$lambda, aniso.pars = model$aniso.pars)
          }
           else{
	    Tkgc=krige.control(type.krige="ok", cov.model = model$cov.model, cov.pars = model$cov.pars, nugget = model$nugget, kappa = model$kappa, lambda = model$lambda, aniso.pars = model$aniso.pars)
           }
            kr <- glm.krige.aux1(data=cv.data.mcmc, coords = cv.coords, 
		loc = coords.out, borders=borders, krige = Tkgc, 
                output = output.control(n.pred=0,mess = FALSE))
	    res[ndata,1]=data.out
	    res[ndata,2]=mean(kr$pred)
	    res[ndata,3]=kr$krige.var[1] + var(as.double(kr$pred))
	    res[ndata,4]=mean(kr$npred)
	    res[ndata,5]=mean(kr$nvar)+var(as.double(kr$npred))          	    	
         }
	res=as.data.frame(res)
    }
    else {
	Tkgc=krige.glm.control(trend.d = ~xmat , trend.l = trend.loc,
		cov.model = model$cov.model, cov.pars = model$cov.pars, 
                nugget = model$nugget, kappa = model$kappa, 
                lambda = model$lambda, aniso.pars = model$aniso)
        res <- binom.krige1(coords = coords, data = data, loc = locations.xvalid, 
            krige = Tkgc,output = output.glm.control(mess = FALSE),nwin=nwin)[1:2]
        res <- data.frame(data.xvalid, res$pred, res$krige.var, res$npred, res$nvar)
    }
    if (messages.screen) 
        cat("\nxvalid: end of cross-validation\n")
    names(res) <- c("pdata", "predicted", "krige.var","npred","nvar")
    res$nerr <- res$pdata - res$npred
    res$std.nerr <- res$nerr/sqrt(res$nvar)
    res$data=logitp.data
    res$error <- res$data - res$pred
    res$std.error <- res$err/sqrt(res$krige.var)    
    res$prob <- pnorm(res$data, mean = res$pred, sd = sqrt(res$krige.var))
    attr(res, "row.names") <- NULL
    if (autocross) {
        if (!is.null(call.fc$geodata)) 
            attr(res, "geodata.xvalid") <- call.fc$geodata
        else attr(res, "locations.xvalid") <- call.fc$locations.xvalid
    }
    else{ if (!is.null(locations.xvalid)) 
        attr(res, "locations.xvalid") <- call.fc$locations.xvalid
	}
    attr(res, "class") <- "xvalid"

    if(resid){
	resn=list(err=res$error,nerr=res$nerr,nerrs=res$std.nerr)
    }
 	else{
    resn=list(Bias=mean(res$error),BiasR=mean(res$error/abs(res$data)),
       MSE=mean(res$error^2),
      MSEStd=mean(res$std.err^2),BiasT=mean(res$nerr),
      BiasTR=mean(res$nerr/abs(res$pdata)),
      MSET=mean(res$nerr^2),MSETStd=mean(res$std.nerr^2))
	}
	
    resn
}