#############################################################
#                                                           #
#	MLE.AIC function                                    #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

mle.aic <- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, var.full=0, alpha=2, contrasts = NULL) {

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$var.full <- mf$alpha <- mf$contrasts <- NULL
    mf$model <- mf$x <- mf$y <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

nrep <- 2^nvar-1

if (size<nvar+1) {stop("Number of observation must be at least equal to the number of predictors (including intercept) + 1")}
if (var.full<0) {
cat("mle.aic: the variance of the full model can not be negative, using default value \n")
var.full <- 0
}

  z <- .Fortran("mleaic",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(nrep),
	as.double(var.full),
	as.double(alpha),
	aic=mat.or.vec(nrep,nvar+1),
	param=mat.or.vec(nrep,nvar),
	var=double(nrep),
	resid=mat.or.vec(nrep,size),
	info=integer(1),
	PACKAGE="wle")
	
result$aic <- z$aic
result$coefficients <- z$param
result$scale <- sqrt(z$var)
result$residuals <- z$resid
result$call <- cl
result$info <- z$info
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

dn <- colnames(xdata)
dimnames(result$coefficients)_list(NULL,dn)
dimnames(result$aic)_list(NULL,c(dn,"aic"))

class(result) <- "mle.aic"

return(result)
}

summary.mle.aic <- function (object, num.max=20, ...) {

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'mle.aic\' object")
}

if (num.max<1) {
cat("summary.mle.aic: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

ans <- list()
aic <- z$aic
if (is.null(nmodel <- nrow(aic))) nmodel <- 1
num.max <- min(nmodel,num.max)
if (nmodel!=1) { 
    nvar <- ncol(aic)-1
    aic <- aic[order(aic[,(nvar+1)]),]
    aic <- aic[1:num.max,]
}

ans$aic <- aic
ans$num.max <- num.max
ans$call <- z$call

class(ans) <- "summary.mle.aic"
return(ans)
}

print.mle.aic <- function (x, digits = max(3, getOption("digits") - 3), ...) {
res_summary.mle.aic(object=x, num.max=nrow(x$aic), ...)
print.summary.mle.aic(res, digits=digits, ...)
}

print.summary.mle.aic <- function (x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")

    cat("\nAkaike Information Criterion (AIC):\n")
    if(x$num.max>1) {
    nvar <- ncol(x$aic)-1
    x$aic[,(nvar+1)] <- signif(x$aic[,(nvar+1)],digits)
    } else {
    nvar <- length(x$aic)-1
    x$aic[(nvar+1)] <- signif(x$aic[(nvar+1)],digits)
    }
    print(x$aic)
    cat("\n")

    cat("Printed the first ",x$num.max," best models \n") 
    invisible(x)
}






#############################################################
#                                                           #
#	MLE.CP function                                     #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

mle.cp <- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, var.full=0, contrasts=NULL)
{

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$var.full <- mf$contrasts <- NULL
    mf$model <- mf$x <- mf$y <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

nrep <- 2^nvar-1

if (size<nvar+1) {
stop("Number of observation must be at least equal to the number of predictors (including intercept) + 1")
}

if (var.full<0) {
cat("mle.cp: the variance of the full model can not be negative, using default value \n")
var.full <- 0
}

  z <- .Fortran("mlecp",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(nrep),
	as.double(var.full),
	cp=mat.or.vec(nrep,nvar+1),
	param=mat.or.vec(nrep,nvar),
	var=double(nrep),
	resid=mat.or.vec(nrep,size),
	info=integer(1),
	PACKAGE="wle")


result$cp <- z$cp
result$coefficients <- z$param
result$scale <- sqrt(z$var)
result$residuals <- z$resid
result$call <- cl
result$info <- z$info
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

dn <- colnames(xdata)
dimnames(result$coefficients)_list(NULL,dn)
dimnames(result$cp)_list(NULL,c(dn,"cp"))

class(result) <- "mle.cp" 

return(result)

}

summary.mle.cp <- function (object, num.max=20, ...) {

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'mle.cp\' object")
}

ans <- list()
cp <- z$cp

if (num.max<1) {
cat("summary.mle.cp: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

if(is.null(nmodel <- nrow(cp))) nmodel <- 1
num.max <- min(nmodel,num.max)
if (nmodel!=1) { 
    nvar <- ncol(cp)-1
    nparam <- apply(cp[,(1:nvar)],1,sum)
    cp <- cp[cp[,(nvar+1)]<=(nparam+0.00001),]
    if(!is.null(nrow(cp)) & nrow(cp)>1) {
	num.max <- min(nrow(cp),num.max)
    	cp <- cp[order(cp[,(nvar+1)]),]
    	cp <- cp[1:num.max,]
    } else num.max <- 1
}

ans$cp <- cp
ans$num.max <- num.max
ans$call <- z$call

class(ans) <- "summary.mle.cp"
return(ans)
}

print.mle.cp <- function (x, digits = max(3, getOption("digits") - 3), ...) {
res_summary.mle.cp(object=x, num.max=nrow(x$cp), ...)
print.summary.mle.cp(res, digits=digits, ...)
}

print.summary.mle.cp <- function (x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")

    cat("\nMallows Cp:\n")
    if(x$num.max>1) {
    nvar <- ncol(x$cp)-1
    x$cp[,(nvar+1)] <- signif(x$cp[,(nvar+1)],digits)
    } else {
    nvar <- length(x$cp)-1
    x$cp[(nvar+1)] <- signif(x$cp[(nvar+1)],digits)
    }
    print(x$cp)
    cat("\n")

    cat("Printed the first ",x$num.max," best models \n") 
    invisible(x)
}


#############################################################
#                                                           #
#	MLE.CV function                                     #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

mle.cv <- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, monte.carlo=500, split, contrasts=NULL)
{

if (missing(split)) {
split <- 0
}

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$monte.carlo <- mf$split <- mf$contrasts <- NULL
    mf$model <- mf$x <- mf$y <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

nrep <- 2^nvar-1

if(size<nvar+1){stop("Number of observation must be at least equal to the number of predictors (including intercept) + 1")}

if(split<nvar+2 | split>(size-2)){
split <- max(round(size^(3/4)),nvar+2)
cat("mle.cv: dimension of the split subsample set to default value = ",split,"\n")
}
maxcarlo <- sum(log(1:size))-(sum(log(1:split))+sum(log(1:(size-split))))
if(monte.carlo<1 | log(monte.carlo) > maxcarlo){
stop("MonteCarlo replication not in the range")
}

  z <- .Fortran("mlecv",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(nrep),
	as.integer(monte.carlo),
	as.integer(split),
	cv=mat.or.vec(nrep,nvar+1),
	info=integer(1))


result$cv <- z$cv
result$call <- cl
result$info <- z$info
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

dn <- colnames(xdata)
dimnames(result$cv) <- list(NULL,c(dn,"cv"))

class(result) <- "mle.cv" 

return(result)

}

summary.mle.cv <- function (object, num.max=20, ...) {

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'mle.cv\' object")
}

if (num.max<1) {
cat("summary.mle.cv: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

ans <- list()
cv <- z$cv
if(is.null(nmodel <- nrow(cv))) nmodel <- 1
num.max <- min(nmodel,num.max)
if (nmodel!=1) { 
nvar <- ncol(cv)-1
cv <- cv[order(cv[,(nvar+1)]),]
cv <- cv[1:num.max,]
}

ans$cv <- cv
ans$num.max <- num.max
ans$call <- z$call

class(ans) <- "summary.mle.cv"
return(ans)
}

print.mle.cv <- function (x, digits = max(3, getOption("digits") - 3), ...) {
res_summary.mle.cv(object=x, num.max=nrow(x$cv), ...)
print.summary.mle.cv(res, digits=digits, ...)
}

print.summary.mle.cv <- function (x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")

    cat("\nCross Validation selection criteria:\n")
    if(x$num.max>1) {
    nvar <- ncol(x$cv)-1
    x$cv[,(nvar+1)] <- signif(x$cv[,(nvar+1)],digits)
    } else {
    nvar <- length(x$cv)-1
    x$cv[(nvar+1)] <- signif(x$cv[(nvar+1)],digits)
    }
    print(x$cv)
    cat("\n")

    cat("Printed the first ",x$num.max," best models \n") 
    invisible(x)
}



#############################################################
#                                                           #
#	MLE.STEPWISE function                               #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

mle.stepwise <- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, type="Forward", f.in=4.0, f.out=4.0, contransts=NULL)
{

ntype <- switch(type,
	Forward = 1,
	Backward = 2,
	Stepwise = 3,
	stop("The type must be Forward, Backward or Stepwise"))

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$type <- mf$f.in <- mf$f.out <- NULL
    mf$max.iter <- mf$contrasts <- NULL
    mf$model <- mf$x <- mf$y <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

if (size<nvar+1) {
stop("Number of observation must be at least equal to the number of predictors (including intercept) + 1")
}

if (f.in<0 | f.out<0) {
stop("f.in and f.out can not be negative")
}

nrep_2^nvar-1

  z <- .Fortran("step",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(nrep),
	as.integer(ntype),
	as.double(f.in),
	as.double(f.out),
	step=mat.or.vec(nrep,nvar+1),
	info=integer(1),
	imodel=integer(1),
	PACKAGE = "wle")

result$step <- z$step[1:z$imodel,]
result$info <- z$info
result$call <- match.call()
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt
result$type <- type
result$f.in <- f.in
result$f.out <- f.out

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

dn <- colnames(xdata)

if (z$imodel==1) {
names(result$step) <- c(dn," ")
} else {
dimnames(result$step) <- list(NULL,c(dn," "))
}


class(result) <- "mle.stepwise"

return(result)

}


summary.mle.stepwise <- function (object, num.max=20, ...) {

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'mle.stepwise\' object")
}

if (num.max<1) {
cat("summary.mle.stepwise: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

ans <- list()
step <- z$step
if(is.null(nmodel <- nrow(step))) nmodel <- 1
num.max <- min(nmodel,num.max)
if (nmodel!=1) { 
    step <- step[(nmodel-num.max+1):nmodel,]
}

ans$step <- step
ans$num.max <- num.max
ans$type <- z$type
ans$f.in <- z$f.in
ans$f.out <- z$f.out
ans$call <- z$call

class(ans) <- "summary.mle.stepwise"
return(ans)
}

print.mle.stepwise <- function (x, digits = max(3, getOption("digits") - 3), ...) {
res_summary.mle.stepwise(object=x, num.max=nrow(x$step), ...)
print.summary.mle.stepwise(res, digits=digits, ...)
}

print.summary.mle.stepwise <- function (x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")

    cat("\n",x$type," selection procedure\n")
    if (x$type=="Forward" | x$type=="Stepwise") {
	cat("\nF.in: ",x$f.in)
    } 
    if (x$type=="Backward" | x$type=="Stepwise") {
	cat("\nF.out: ",x$f.out)
    }
    cat(" \n")
    cat("\nLast ",x$num.max," iterations:\n")

    if(x$num.max>1) {
    nvar <- ncol(x$step)-1
    x$step[,(nvar+1)] <- signif(x$step[,(nvar+1)],digits)
    } else {
    nvar <- length(x$step)-1
    x$step[(nvar+1)] <- signif(x$step[(nvar+1)],digits)
    }
    print(x$step)
    cat("\n")
    invisible(x)
}






#############################################################
#                                                           #
#	PLOT.MLE.CP function                                #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

plot.mle.cp <- function(object, base.line=0, num.max=20, plot.it=TRUE, log.scale=FALSE, xlab="Number of Predictors", ylab=NULL)
{

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'mle.cp\' object")
}

cp <- object$cp

if (num.max<1) {
cat("plot.mle.cp: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

if(is.null(nrow(cp)) | nrow(cp)==1) {
    num.model <- 1
} else {
    num.model <- nrow(cp) 
}

if(num.model<num.max)
{
cat("plot.mle.cp: The number of models is less than num.max \n")
num.max <- num.model
}

if(is.null(ncol(cp))) {
stop("No models to plot")
} else {
nvar <- ncol(cp)-1
}

good.model <- (apply(cp[,1:nvar],1,sum)+base.line>=cp[,nvar+1])
cp.good <- matrix(cp[good.model,],ncol=nvar+1)
cp.bad <- matrix(cp[!good.model,],ncol=nvar+1)
ordine.good <- order(cp.good[,nvar+1])
ordine.bad <- order(cp.bad[,nvar+1])
cp.good <- matrix(cp.good[ordine.good,],ncol=(nvar+1))
num.good <- dim(cp.good)[1]
cp.bad <- matrix(cp.bad[ordine.bad,],ncol=(nvar+1))
num.bad <- dim(cp.bad)[1]

label.good <- character()
for(i in 1:nvar){
label.good <- paste(label.good,cp.good[,i],sep="")
}
label.bad <- character()
for(i in 1:nvar){
label.bad <- paste(label.bad,cp.bad[,i],sep="")
}


xcoord.good <- apply(matrix(cp.good[,1:nvar],ncol=nvar),1,sum)[1:min(num.max,num.good)]
ycoord.good <- cp.good[,nvar+1][1:min(num.max,num.good)]

label.good <- label.good[1:min(num.max,num.good)]

xcoord.best <- xcoord.good[1]
ycoord.best <- ycoord.good[1]

label.best <- label.good[1]

if(length(xcoord.good)==1)
{
xcoord.good <- 0
ycoord.good <- 0
plot.good <- FALSE
}else
{
xcoord.good <- xcoord.good[-1]
ycoord.good <- ycoord.good[-1]
label.good <- label.good[-1]
plot.good <- TRUE
}

if(num.max>num.good)
{
xcoord.bad <- apply(matrix(cp.bad[,1:nvar],ncol=nvar),1,sum)[1:min(num.bad,num.max-num.good)]
ycoord.bad <- cp.bad[,nvar+1][1:min(num.bad,num.max-num.good)]
label.bad <- label.bad[1:min(num.bad,num.max-num.good)]
plot.bad <- TRUE
}else
{
xcoord.bad <- 0
ycoord.bad <- 0
plot.bad <- FALSE
}

xlim.min <- min(xcoord.good,xcoord.bad,xcoord.best)
xlim.max <- max(xcoord.good,xcoord.bad,xcoord.best)

yetichetta <- "Cp"

if(log.scale)
{
ycoord.good <- log10(ycoord.good+min(ycoord.good,ycoord.bad,ycoord.best)+1)
ycoord.bad <- log10(ycoord.bad+min(ycoord.good,ycoord.bad,ycoord.best)+1)
ycoord.best <- log10(ycoord.best+min(ycoord.good,ycoord.bad,ycoord.best)+1)
yetichetta <- "Cp log10 scale"
}

ylim.min <- min(ycoord.good,ycoord.bad,ycoord.best)
ylim.max <- max(ycoord.good,ycoord.bad,ycoord.best)

if (is.null(ylab)) {
ylab <- yetichetta
}

if(plot.it)
{
plot(xcoord.best,ycoord.best,xlim=c(xlim.min,xlim.max),ylim=c(ylim.min,ylim.max),xlab=xlab,ylab=ylab,type="n")
text(xcoord.best,ycoord.best,col=4,labels=label.best)

if(plot.good)
{
text(xcoord.good,ycoord.good,col=3,labels=label.good)
}

if(plot.bad)
{
text(xcoord.bad,ycoord.bad,col=2,labels=label.bad)
}


if(!log.scale)
{
abline(base.line,1,col=2)
abline(0,1)
}
else
{
vettx <- seq(xlim.min,xlim.max,0.5)
vetty <- log10(vettx+min(ycoord.good,ycoord.bad,ycoord.best)+1)
vetty.base.line <- log10(vettx+min(ycoord.good,ycoord.bad,ycoord.best)+1+base.line)
lines(vettx,vetty.base.line,col=2,type="l")
lines(vettx,vetty,type="l")
}

}

invisible(list(num.good=num.good,num.bad=num.bad,cp.good=cp.good, cp.bad=cp.bad))
}
#############################################################
#                                                           #
#	PLOT.WLE.CP function                                #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

plot.wle.cp <- function(object, base.line=0, num.max=20, plot.it=TRUE, log.scale=FALSE, xlab="Number of Predictors", ylab=NULL)
{

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'wle.cp\' object")
}

wcp <- object$wcp

if (num.max<1) {
cat("plot.wle.cp: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

if(is.null(nrow(wcp)) | nrow(wcp)==1) {
    num.model <- 1
} else {
    num.model <- nrow(wcp) 
}

if(num.model<num.max)
{
cat("plot.wle.cp: The number of models is less than num.max \n")
num.max <- num.model
}

if(is.null(ncol(wcp))) {
stop("No models to plot")
} else {
nvar <- ncol(wcp)-1
}

good.model <- (apply(wcp[,1:nvar],1,sum)+base.line>=wcp[,nvar+1])
wcp.good <- matrix(wcp[good.model,],ncol=nvar+1)
wcp.bad <- matrix(wcp[!good.model,],ncol=nvar+1)
ordine.good <- order(wcp.good[,nvar+1])
ordine.bad <- order(wcp.bad[,nvar+1])
wcp.good <- matrix(wcp.good[ordine.good,],ncol=(nvar+1))
num.good <- dim(wcp.good)[1]
wcp.bad <- matrix(wcp.bad[ordine.bad,],ncol=(nvar+1))
num.bad <- dim(wcp.bad)[1]

label.good <- character()
for(i in 1:nvar){
label.good <- paste(label.good,wcp.good[,i],sep="")
}
label.bad <- character()
for(i in 1:nvar){
label.bad <- paste(label.bad,wcp.bad[,i],sep="")
}


xcoord.good <- apply(matrix(wcp.good[,1:nvar],ncol=nvar),1,sum)[1:min(num.max,num.good)]
ycoord.good <- wcp.good[,nvar+1][1:min(num.max,num.good)]

label.good <- label.good[1:min(num.max,num.good)]

xcoord.best <- xcoord.good[1]
ycoord.best <- ycoord.good[1]

label.best <- label.good[1]

if(length(xcoord.good)==1)
{
xcoord.good <- 0
ycoord.good <- 0
plot.good <- FALSE
}else
{
xcoord.good <- xcoord.good[-1]
ycoord.good <- ycoord.good[-1]
label.good <- label.good[-1]
plot.good <- TRUE
}

if(num.max>num.good)
{
xcoord.bad <- apply(matrix(wcp.bad[,1:nvar],ncol=nvar),1,sum)[1:min(num.bad,num.max-num.good)]
ycoord.bad <- wcp.bad[,nvar+1][1:min(num.bad,num.max-num.good)]
label.bad <- label.bad[1:min(num.bad,num.max-num.good)]
plot.bad <- TRUE
}else
{
xcoord.bad <- 0
ycoord.bad <- 0
plot.bad <- FALSE
}

xlim.min <- min(xcoord.good,xcoord.bad,xcoord.best)
xlim.max <- max(xcoord.good,xcoord.bad,xcoord.best)

yetichetta <- "WCp"

if(log.scale)
{
ycoord.good <- log10(ycoord.good+min(ycoord.good,ycoord.bad,ycoord.best)+1)
ycoord.bad <- log10(ycoord.bad+min(ycoord.good,ycoord.bad,ycoord.best)+1)
ycoord.best <- log10(ycoord.best+min(ycoord.good,ycoord.bad,ycoord.best)+1)
yetichetta <- "WCp log10 scale"
}

ylim.min <- min(ycoord.good,ycoord.bad,ycoord.best)
ylim.max <- max(ycoord.good,ycoord.bad,ycoord.best)

if (is.null(ylab)) {
ylab <- yetichetta
}

if(plot.it)
{
plot(xcoord.best,ycoord.best,xlim=c(xlim.min,xlim.max),ylim=c(ylim.min,ylim.max),xlab=xlab,ylab=ylab,type="n")
text(xcoord.best,ycoord.best,col=4,labels=label.best)

if(plot.good)
{
text(xcoord.good,ycoord.good,col=3,labels=label.good)
}

if(plot.bad)
{
text(xcoord.bad,ycoord.bad,col=2,labels=label.bad)
}


if(!log.scale)
{
abline(base.line,1,col=2)
abline(0,1)
}
else
{
vettx <- seq(xlim.min,xlim.max,0.5)
vetty <- log10(vettx+min(ycoord.good,ycoord.bad,ycoord.best)+1)
vetty.base.line <- log10(vettx+min(ycoord.good,ycoord.bad,ycoord.best)+1+base.line)
lines(vettx,vetty.base.line,col=2,type="l")
lines(vettx,vetty,type="l")
}

}

invisible(list(num.good=num.good,num.bad=num.bad,wcp.good=wcp.good, wcp.bad=wcp.bad))
}
#############################################################
#                                                           #
#	PLOT.WLE.LM function                                #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

plot.wle.lm <- function(object, level.weight=0.5, ask = interactive() && .Device != "postscript") {

old.par <- par(no.readonly=TRUE)
on.exit(par(old.par))

if (!inherits(object, "wle.lm")) stop("Use only with 'wle.lm' objects")

if (ask) par(ask = TRUE)
 
if (!is.null(object$model)) {
ydata <- object$model[,1]
xdata <- object$model[,-1]
} else if (!is.null(object$x) & !is.null(object$y)) {
ydata <- object$y
xdata <- object$x
} else {
stop("Please, rerun wle.lm with model=TRUE (or x=TRUE and y=TRUE)")
}

if (level.weight<0 | level.weight >1) {
cat("plot.wle.lm: level.weight should be between zero and one, set to 0.5 \n")
level.weight <- 0.5
}

param <- as.matrix(object$coefficients)
res <- as.matrix(object$residuals)
y.fit <- as.matrix(object$fitted.values)
weight <- as.matrix(object$weights)
tot.weight <- object$tot.weights
tot.sol <- object$tot.sol

if (tot.sol>1) {
par(mfcol=c(tot.sol,tot.sol))

for (isol in 1:tot.sol) {
for (jsol in 1:tot.sol) {
y.fit.i <- y.fit[isol,]
res.i <- res[isol,]
weight.i <- weight[isol,]

y.fit.j <- y.fit[jsol,]
res.j <- res[jsol,]
weight.j <- weight[jsol,]

level.i <- weight.i>=level.weight
level.j <- weight.j>=level.weight

color <- rep(2,length(ydata))
color[level.i] <- 1
color[level.j] <- 3

color.res <- rep(2,length(ydata))
color.res[level.i & (res.i>res.j)] <- 1
color.res[level.j & (res.i<res.j)] <- 3

color.w <- rep(2,length(ydata))
color.w[level.i & (weight.i>weight.j)] <- 1
color.w[level.j & (weight.i<weight.j)] <- 3

if (isol==jsol) {
plot(weight.i,col=color,xlab="Observations",ylab="Weights",main=paste("Weights of the root: ",isol))
} else {
if (isol>jsol) {
plot(res.i,res.j,col=color.res,xlab=paste("Residuals of the root: ",isol),ylab=paste("Residuals of the root: ",jsol),main="Residuals")
abline(0,1)
} else {
plot(weight.i,weight.j,col=color.w,xlab=paste("Weights of the root: ",isol),ylab=paste("Weights of the root: ",jsol),main="Weights")
abline(0,1)
}
}
}
}


for (isol in 1:tot.sol) {
par(mfcol=c(2,2))
y.fit.temp <- y.fit[isol,]
res.temp <- res[isol,]
weight.temp <- weight[isol,]

level <- weight.temp>=level.weight
color <- rep(2,length(ydata))
color[level] <- 1

plot(y.fit.temp,res.temp,col=color,xlab="Fitted values",ylab="Residuals")

plot(y.fit.temp,res.temp*weight.temp,col=color,xlab="Fitted values",ylab="Weighted residuals")

qqnorm(res.temp,col=color)
qqline(res.temp)

qqnorm(res.temp*weight.temp,col=color)
qqline(res.temp*weight.temp)
}
} else {

par(mfcol=c(1,1))
level.i <- weight>=level.weight
color <- rep(2,length(ydata))
color[level.i] <- 3

plot(weight,col=color,xlab="Observations",ylab="Weights",main="Weights of the root")

par(mfcol=c(2,2))

level <- weight>=level.weight
color <- rep(2,length(ydata))
color[level] <- 1

plot(y.fit,res,col=color,xlab="Fitted values",ylab="Residuals")

plot(y.fit,res*weight,col=color,xlab="Fitted values",ylab="Weighted residuals")

qqnorm(res,col=color)
qqline(res)

qqnorm(res*weight,col=color)
qqline(res*weight)
}
}



#############################################################
#                                                           #
#	WLE.AIC function                                    #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.aic<- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, boot=30, group, var.full=0, num.sol=1, raf="HD", smooth=0.031, tol=10^(-6), equal=10^(-3), max.iter=500, min.weight=0.5, method="full", alpha=2, contrasts=NULL)
{

raf <- switch(raf,
	HD = 1,
	NED = 2,
	SCHI2 = 3,
	-1)

if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")

type <- switch(method,
	full = 0,
	reduced = 1,
	-1)

if (type==-1) stop("Please, choose the method: full=wieghts based on full model, reduced=weights based on the actual model")

if (missing(group)) {
group <- 0
}

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$boot <- mf$group <- mf$smooth <- NULL
    mf$tol <- mf$equal <- mf$num.sol <- NULL
    mf$min.weight <- mf$max.iter <- mf$raf <- NULL
    mf$var.full <- mf$alpha <- mf$contrasts <- NULL
    mf$model <- mf$x <- mf$y <- mf$method <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

nrep <- 2^nvar-1

if (size<nvar+1) {
stop("Number of observations must be at least equal to the number of predictors (including intercept) + 1")
}

if (!(group>1)) {
group <- max(round(size/4),nvar+2)
cat("wle.aic: dimension of the subsample set to default value = ",group,"\n")
}

maxboot <- sum(log(1:size))-(sum(log(1:group))+sum(log(1:(size-group))))

if (boot<1 | log(boot) > maxboot) {
stop("bootstrap replication not in the range")
}

if (!(num.sol>=1)) {
cat("wle.aic: number of solution to report set to 1 \n")
num.sol <- 1
}

if(max.iter<1) {
cat("wle.aic: max number of iteration set to 500 \n")
max.iter <- 500
}

if (smooth<10^(-5)) {
cat("wle.aic: the smooth parameter seems too small \n")
}

if (tol<0) {
cat("wle.aic: the accuracy can not be negative, using default value \n")
tol <- 10^(-6)
}

if (equal<0) {
cat("wle.aic: the equal parameter can not be negative, using default value \n")
equal <- 10^(-3)
}

if (var.full<0) {
cat("wle.aic: the variance of the full model can not be negative, using default value \n")
var.full <- 0
}

if (min.weight<0) {
cat("wle.aic: the minimum sum of the weights can not be negative, using default value \n")
min.weight <- 0.5
}

  z <- .Fortran("wleaic",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(boot),
	as.integer(group),
	as.integer(nrep),
	as.integer(raf),
	as.double(smooth),
	as.double(tol),
	as.double(equal),
	as.integer(max.iter),
	as.double(var.full),
	as.integer(num.sol),
	as.double(min.weight),
	as.integer(type),
	as.double(alpha),
	waic=mat.or.vec(nrep*num.sol,nvar+1),
	param=mat.or.vec(nrep*num.sol,nvar),
	var=double(nrep*num.sol),
	resid=mat.or.vec(nrep*num.sol,size),
	totweight=double(nrep*num.sol),
	weight=mat.or.vec(nrep*num.sol,size),
	same=integer(nrep*num.sol),
	info=integer(1),
	PACKAGE="wle")

delnull <- z$same==0

result$waic <- z$waic[!delnull,]
result$coefficients <- z$param[!delnull,]
result$scale <- sqrt(z$var[!delnull])
result$residuals <- z$resid[!delnull]
result$weights <- z$weight[!delnull,]
result$tot.weights <- z$totweight[!delnull]
result$freq <- z$same[!delnull]
result$call <- cl
result$info <- z$info
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

dn <- colnames(xdata)
dimnames(result$coefficients) <- list(NULL,dn)
dimnames(result$waic) <- list(NULL,c(dn,"waic"))

class(result) <- "wle.aic"

return(result)
}


summary.wle.aic <- function (object, num.max=20, ...) {

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'wle.aic\' object")
}

if (num.max<1) {
cat("summary.wle.aic: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

ans <- list()
waic <- z$waic
if (is.null(nmodel <- nrow(waic))) nmodel <- 1
num.max <- min(nmodel,num.max)
if (nmodel!=1) { 
    nvar <- ncol(waic)-1
    waic <- waic[order(waic[,(nvar+1)]),]
    waic <- waic[1:num.max,]
}

ans$waic <- waic
ans$num.max <- num.max
ans$call <- z$call

class(ans) <- "summary.wle.aic"
return(ans)
}

print.wle.aic <- function (x, digits = max(3, getOption("digits") - 3), ...) {
res_summary.wle.aic(object=x, num.max=nrow(x$waic), ...)
print.summary.wle.aic(res, digits=digits, ...)
}

print.summary.wle.aic <- function (x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")

    cat("\nWeighted Akaike Information Criterion (WAIC):\n")
    if(x$num.max>1) {
    nvar <- ncol(x$waic)-1
    x$waic[,(nvar+1)] <- signif(x$waic[,(nvar+1)],digits)
    } else {
    nvar <- length(x$waic)-1
    x$waic[(nvar+1)] <- signif(x$waic[(nvar+1)],digits)
    }
    print(x$waic)
    cat("\n")

    cat("Printed the first ",x$num.max," best models \n") 
    invisible(x)
}








#############################################################
#                                                           #
#	WLE.CP function                                     #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.cp <- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, boot=30, group, var.full=0, num.sol=1, raf="HD", smooth=0.031, tol=10^(-6), equal=10^(-3), max.iter=500, min.weight=0.5, method="full", alpha=2, contrasts=NULL)
{

raf <- switch(raf,
	HD = 1,
	NED = 2,
	SCHI2 = 3,
	-1)

if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")

type <- switch(method,
	full = 0,
	reduced = 1,
	-1)

if (type==-1) stop("Please, choose the method: full=wieghts based on full model, reduced=weights based on the actual model")

if (missing(group)) {
group <- 0
}

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$boot <- mf$group <- mf$smooth <- NULL
    mf$tol <- mf$equal <- mf$num.sol <- NULL
    mf$min.weight <- mf$max.iter <- mf$raf <- NULL
    mf$var.full <- mf$alpha <- mf$contrasts <- NULL
    mf$model <- mf$x <- mf$y <- mf$method <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

nrep <- 2^nvar-1

if (size<nvar+1) {
stop("Number of observations must be at least equal to the number of predictors (including intercept) + 2")
}

if (!(group>1)) {
group <- max(round(size/4),nvar+2)
cat("wle.cp: dimension of the subsample set to default value = ",group,"\n")
}

maxboot <- sum(log(1:size))-(sum(log(1:group))+sum(log(1:(size-group))))

if (boot<1 | log(boot) > maxboot) {
stop("Bootstrap replication not in the range")
}


if (!(num.sol>=1)) {
cat("wle.cp: number of solution to report set to 1 \n")
num.sol <- 1
}

if (max.iter<1) {
cat("wle.cp: max number of iteration set to 500 \n")
max.iter <- 500
}

if (smooth<10^(-5)) {
cat("wle.cp: the smooth parameter seems too small \n")
}

if (tol<0) {
cat("wle.cp: the accuracy can not be negative, using default value \n")
tol <- 10^(-6)
}
if (equal<0) {
cat("wle.cp: the equal parameter can not be negative, using default value \n")
equal <- 10^(-3)
}

if (var.full<0) {
cat("wle.cp: the variance of the full model can not be negative, using default value \n")
var.full <- 0
}

if (min.weight<0) {
cat("wle.cp: the minimum sum of the weights can not be negative, using default value \n")
min.weight <- 0.5
}

  z <- .Fortran("wlecp",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(boot),
	as.integer(group),
	as.integer(nrep),
	as.integer(raf),
	as.double(smooth),
	as.double(tol),
	as.double(equal),
	as.integer(max.iter),
	as.double(var.full),
	as.integer(num.sol),
	as.double(min.weight),
	as.integer(type),
	as.double(alpha),
	wcp=mat.or.vec(nrep*num.sol,nvar+1),
        param=mat.or.vec(nrep*num.sol,nvar),
	var=double(nrep*num.sol),
	resid=mat.or.vec(nrep*num.sol,size),
	totweight=double(nrep*num.sol),
	weight=mat.or.vec(nrep*num.sol,size),
	same=integer(nrep*num.sol),
	info=integer(1),
	PACKAGE = "wle")

delnull <- z$same==0

result$wcp <- z$wcp[!delnull,]
result$coefficients <- z$param[!delnull,]
result$scale <- sqrt(z$var[!delnull])
result$residuals <- z$resid[!delnull]
result$weights <- z$weight[!delnull,]
result$tot.weights <- z$totweight[!delnull]
result$freq <- z$same[!delnull]
result$call <- cl
result$info <- z$info
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

dn <- colnames(xdata)
dimnames(result$coefficients) <- list(NULL,dn)
dimnames(result$wcp) <- list(NULL,c(dn,"wcp"))

class(result) <- "wle.cp"

return(result)

}

summary.wle.cp <- function (object, num.max=20, ...) {

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'wle.cp\' object")
}

if (num.max<1) {
cat("summary.wle.cp: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

ans <- list()
wcp <- z$wcp
if(is.null(nmodel <- nrow(wcp))) nmodel <- 1
num.max <- min(nmodel,num.max)

if (nmodel!=1) { 
    nvar <- ncol(wcp)-1
    nparam <- apply(wcp[,(1:nvar)],1,sum)
    wcp <- wcp[wcp[,(nvar+1)]<=(nparam+0.00001),]
    if(!is.null(nrow(wcp)) & nrow(wcp)>1) {
	num.max <- min(nrow(wcp),num.max)
    	wcp <- wcp[order(wcp[,(nvar+1)]),]
    	wcp <- wcp[1:num.max,]
    } else num.max <- 1
}

ans$wcp <- wcp
ans$num.max <- num.max
ans$call <- z$call

class(ans) <- "summary.wle.cp"
return(ans)
}

print.wle.cp <- function (x, digits = max(3, getOption("digits") - 3), ...) {
res_summary.wle.cp(object=x, num.max=nrow(x$wcp), ...)
print.summary.wle.cp(res, digits=digits, ...)
}

print.summary.wle.cp <- function (x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")

    cat("\nWeighted Mallows Cp:\n")
    if(x$num.max>1) {
    nvar <- ncol(x$wcp)-1
    x$wcp[,(nvar+1)] <- signif(x$wcp[,(nvar+1)],digits)
    } else {
    nvar <- length(x$wcp)-1
    x$wcp[(nvar+1)] <- signif(x$wcp[(nvar+1)],digits)
    }
    print(x$wcp)
    cat("\n")

    cat("Printed the first ",x$num.max," best models \n") 
    invisible(x)
}



#############################################################
#                                                           #
#	WLE.CV function                                     #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.cv <- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, monte.carlo=500, split, boot=30, group, num.sol=1, raf="HD", smooth=0.031, tol=10^(-6), equal=10^(-3), max.iter=500, min.weight=0.5, contrasts=NULL)
{

raf <- switch(raf,
	HD = 1,
	NED = 2,
	SCHI2 = 3,
	-1)

if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")

if (missing(group)) {
group <- 0
}

if (missing(split)) {
split <- 0
}

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$monte.carlo <- mf$split <- NULL
    mf$boot <- mf$group <- mf$smooth <- NULL
    mf$tol <- mf$equal <- mf$num.sol <- NULL
    mf$min.weight <- mf$max.iter <- mf$raf <- NULL
    mf$contrasts <- NULL
    mf$model <- mf$x <- mf$y <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

nrep <- 2^nvar-1

if (size<nvar+1) {
stop("Number of observations must be at least equal to the number of predictors (including intercept) + 1")
}

if ((!(group>1))|group<(nvar+1)) {
group <- max(round(size/4),nvar+2)
cat("wle.cv: dimension of the subsample set to default value = ",group,"\n")
}

maxboot <- sum(log(1:size))-(sum(log(1:group))+sum(log(1:(size-group))))

if (boot<1 | log(boot) > maxboot) {
stop("Bootstrap replication not in the range")
}

if (split<nvar+2 | split>(size-2)) {
split <- max(round(size^(3/4)),nvar+2)
cat("wle.cv: dimension of the split subsample set to default value = ",split,"\n")
}

maxcarlo <- sum(log(1:size))-(sum(log(1:split))+sum(log(1:(size-split))))

if (monte.carlo<1 | log(monte.carlo) > maxcarlo) {
stop("MonteCarlo replication not in the range")
}

if (!(num.sol>=1)) {
cat("wle.cv:number of solution to report set to 1 \n")
num.sol <- 1
}

if (max.iter<1) {
cat("wle.cv: max number of iteration set to 500 \n")
max.iter <- 500
}

if (smooth<10^(-5)) {
cat("wle.cv: the smooth parameter seems too small \n")
}

if (tol<0) {
cat("wle.cv: the accuracy can not be negative, using default value \n")
tol <- 10^(-6)
}

if (equal<0) {
cat("wle.cv: the equal parameter can not be negative, using default value \n")
equal <- 10^(-3)
}

if (min.weight<0) {
cat("wle.cv: the minimum sum of the weights can not be negative, using default value \n")
min.weight <- 0.5
}

  z <- .Fortran("wlecv",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(boot),
	as.integer(group),
	as.integer(nrep),
	as.integer(monte.carlo),
	as.integer(split),
	as.integer(raf),
	as.double(smooth),
	as.double(tol),
	as.double(equal),
	as.integer(max.iter),
	as.integer(num.sol),
	as.double(min.weight),
	wcv=mat.or.vec(nrep,nvar+1),
	param=mat.or.vec(num.sol,nvar),
	var=double(num.sol),
	resid=mat.or.vec(num.sol,size),
	totweight=double(num.sol),
	weight=mat.or.vec(num.sol,size),
	same=integer(num.sol),
	index=integer(1),
	info=integer(1))

delnull <- z$same==0

result$wcv <- z$wcv[!delnull,]
result$coefficients <- z$param[!delnull,]
result$scale <- sqrt(z$var[!delnull])
result$residuals <- z$resid[!delnull]
result$weights <- z$weight[!delnull,]
result$tot.weights <- z$totweight[!delnull]
result$freq <- z$same[!delnull]
result$call <- cl
result$info <- z$info
result$index <- z$index
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

dn <- colnames(xdata)
if (is.null(nrow(result$coefficients))) {
names(result$coefficients) <- dn
} else {
dimnames(result$coefficients) <- list(NULL,dn)
}
dimnames(result$wcv) <- list(NULL,c(dn,"wcv"))

class(result) <- "wle.cv"

return(result)

}


summary.wle.cv <- function (object, num.max=20, ...) {

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'wle.cv\' object")
}

if (num.max<1) {
cat("summary.wle.cv: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

ans <- list()
wcv <- z$wcv
if(is.null(nmodel <- nrow(wcv))) nmodel <- 1
num.max <- min(nmodel,num.max)
if (nmodel!=1) { 
nvar <- ncol(wcv)-1
wcv <- wcv[order(wcv[,(nvar+1)]),]
wcv <- wcv[1:num.max,]
}

ans$wcv <- wcv
ans$num.max <- num.max
ans$call <- z$call

class(ans) <- "summary.wle.cv"
return(ans)
}

print.wle.cv <- function (x, digits = max(3, getOption("digits") - 3), ...) {
res_summary.wle.cv(object=x, num.max=nrow(x$wcv), ...)
print.summary.wle.cv(res, digits=digits, ...)
}

print.summary.wle.cv <- function (x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")

    cat("\nWeighted Cross Validation selection criteria:\n")
    if(x$num.max>1) {
    nvar <- ncol(x$wcv)-1
    x$wcv[,(nvar+1)] <- signif(x$wcv[,(nvar+1)],digits)
    } else {
    nvar <- length(x$wcv)-1
    x$wcv[(nvar+1)] <- signif(x$wcv[(nvar+1)],digits)
    }
    print(x$wcv)
    cat("\n")

    cat("Printed the first ",x$num.max," best models \n") 
    invisible(x)
}



#############################################################
#                                                           #
#	WLE.LM function                                     #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.lm <- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, boot=30, group, num.sol=1, raf="HD", smooth=0.031, tol=10^(-6), equal=10^(-3), max.iter=500, contrasts=NULL)
{

raf <- switch(raf,
	HD = 1,
	NED = 2,
	SCHI2 = 3,
	-1)

if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")

if (missing(group)) {
group <- 0
}

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$boot <- mf$group <- mf$smooth <- NULL
    mf$tol <- mf$equal <- mf$num.sol <- NULL
    mf$max.iter <- mf$raf <- mf$contrasts <- NULL
    mf$model <- mf$x <- mf$y <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

if (size<(nvar+1)) {
stop("Number of observations must be at least equal to the number of predictors (including intercept) + 1")
}

if ((!(group>1)) | group<(nvar+1)) {
group <- max(round(size/4),nvar+2)
cat("wle.lm: dimension of the subsample set to default value: ", group,"\n")
}

maxboot <- sum(log(1:size))-(sum(log(1:group))+sum(log(1:(size-group))))

if (boot<1 | log(boot) > maxboot) {
stop("Bootstrap replication not in the range")
}

if (!(num.sol>=1)) {
cat("wle.lm: number of solution to report set to 1 \n")
num.sol <- 1
}

if (max.iter<1) {
cat("wle.lm: max number of iteration set to 500 \n")
max.iter <- 500
}

if (smooth<10^(-5)) {
cat("wle.lm: the smooth parameter seems too small \n")
}

if (tol<0) {
cat("wle.lm: the accuracy can not be negative, using default value \n")
tol <- 10^(-6)
}

if (equal<0) {
cat("wle.lm: the equal parameter can not be negative, using default value \n")
equal <- 10^(-3)
}

  z <- .Fortran("wleregfix",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(nvar),
	as.integer(boot),
	as.integer(group),
	as.integer(num.sol),
	as.integer(raf),
	as.double(smooth),
	as.double(tol),
	as.double(equal),
	as.integer(max.iter),
	param=mat.or.vec(num.sol,nvar),
	var=double(num.sol),
	resid=mat.or.vec(num.sol,size),
	totweight=double(num.sol),
	weight=mat.or.vec(num.sol,size),
	same=integer(num.sol),
	nsol=integer(1),
	nconv=integer(1),
	PACKAGE = "wle")

if (z$nsol>0) {
z$var <- z$var[1:z$nsol]
z$totweight <- z$totweight[1:z$nsol]
z$same <- z$same[1:z$nsol]

if (num.sol==1)
{
z$param <- c(z$param)
z$resid <- c(z$resid)
z$weight <- c(z$weight)
}
else
{
z$param <- z$param[1:z$nsol,]
z$resid <- z$resid[1:z$nsol,]
z$weight <- z$weight[1:z$nsol,]
}

y.fit <- t(xdata%*%matrix(z$param,ncol=z$nsol,byrow=TRUE))

if(z$nsol==1)
{
devparam <- sqrt(z$var*diag(solve(t(xdata)%*%diag(z$weight)%*%xdata,tol=1e-100)))
y.fit <- as.vector(y.fit)
}
else
{
devparam <- sqrt(z$var[1]*diag(solve(t(xdata)%*%diag(z$weight[1,])%*%xdata,tol=1e-100)))
for(i in 2:z$nsol){
devparam <- rbind(devparam,sqrt(z$var[i]*diag(solve(t(xdata)%*%diag(z$weight[i,])%*%xdata,tol=1e-100))))
}
}

result$coefficients <- z$param
result$standard.error <- devparam
result$scale <- sqrt(z$var)
result$residuals <- z$resid
result$fitted.values <- y.fit
result$weights <- z$weight
result$tot.weights <- z$totweight
result$tot.sol <- z$nsol
result$not.conv <- z$nconv
result$freq <- z$same
result$call <- cl
result$info <- z$info
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

dn <- colnames(xdata)

if (is.null(nrow(result$coefficients))) {
names(result$coefficients) <- dn
} else {
dimnames(result$coefficients) <- list(NULL,dn)
}

class(result) <- "wle.lm"

return(result)
} else {
stop("No solutions are fuond, checks the parameters")
}

}

print.wle.lm <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("Coefficients:\n")
    print.default(format(coef(x), digits=digits),
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("Scale estimate: ",format(x$scale, digits=digits))
    cat("\n")
    cat("\nNumber of solutions ",x$tot.sol,"\n")
    cat("\n")
    invisible(x)
}

summary.wle.lm <- function(object, root="ALL", ...)
{
z <- .Alias(object)
if (is.null(z$terms)) stop("invalid \'wle.lm\' object")

tot.sol <- z$tot.sol

if (root!="ALL" & !is.numeric(root)) {
    stop("Please, choose one root, for print all root ALL")
} else if (root=="ALL") {
    root <- 1:tot.sol
} else if (tot.sol<root) {
    stop(paste("Root ",root," not found"))
}

ans <- list()
for (iroot in root) {
ans <- c(ans,list(summary.wle.lm.root(object=object, root=iroot)))
}
class(ans) <- "summary.wle.lm" 

return(ans)
}

print.summary.wle.lm <- function (x, digits = max(3, getOption("digits") - 3), signif.stars= getOption("show.signif.stars"),	...)
{
for (i in 1:length(x)) {
print.summary.wle.lm.root(x[[i]], digits = max(3, getOption("digits") - 3), signif.stars= getOption("show.signif.stars"))
}
}




summary.wle.lm.root <- function (object, root=1, ...) {
z <- .Alias(object)
if (is.null(z$terms)) stop("invalid \'wle.lm\' object")

tot.sol <- z$tot.sol
if (tot.sol<root) {
    stop(paste("Root ",root," not found"))
}
if (tot.sol!=1) {
n <- ncol(z$residuals)
p <- ncol(z$coefficients)
} else {
n <- length(z$residuals)
p <- length(z$coefficients)
}

rdf <- z$tot.weights[root]*n - p 
if (tot.sol>1) {   
    r <- z$residuals[root,]
    f <- z$fitted.values[root,]
    w <- z$weights[root,]
    est <- z$coefficients[root,]
    se <- z$standard.error[root,]
} else {
    r <- z$residuals
    f <- z$fitted.values
    w <- z$weights
    est <- z$coefficients
    se <- z$standard.error
}

    mss <- if (attr(z$terms, "intercept")) {
    		m <- sum(w * f /sum(w))
        	sum(w * (f - m)^2)
        	} else sum(w * f^2)
    rss <- sum(w * r^2)
  
    resvar <- rss/rdf
    tval <- est/se
    ans <- z[c("call", "terms")]
    ans$residuals <- sqrt(w)*r
    ans$coefficients <- cbind(est, se, tval, 2*(1 - pt(abs(tval), rdf)))
    dimnames(ans$coefficients)<-
	list(names(est),
	     c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
    ans$sigma <- sqrt(resvar)
    ans$df <- c(p, rdf, p)
    if (p != attr(z$terms, "intercept")) {
	df.int <- if (attr(z$terms, "intercept")) 1 else 0
	ans$r.squared <- mss/(mss + rss)
	ans$adj.r.squared <- 1 - (1 - ans$r.squared) *
	    ((n - df.int)/rdf)
	ans$fstatistic <- c(value = (mss/(p - df.int))/resvar,
			    numdf = p - df.int, dendf = rdf)
    }

    ans$root <- root
    return(ans)
}

print.summary.wle.lm.root <- function (x, digits = max(3, getOption("digits") - 3), signif.stars= getOption("show.signif.stars"),	...)
{
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")
    resid <- x$residuals
    df <- x$df
    rdf <- df[2]
 
   cat("Root ",x$root)
   cat("\n\nWeighted Residuals:\n", sep="")
    if (rdf > 5) {
	nam <- c("Min", "1Q", "Median", "3Q", "Max")
	rq <- if (length(dim(resid)) == 2)
	    structure(apply(t(resid), 1, quantile),
		      dimnames = list(nam, dimnames(resid)[[2]]))
	else  structure(quantile(resid), names = nam)
	print(rq, digits = digits, ...)
    }
    else if (rdf > 0) {
	print(resid, digits = digits, ...)
    } else { # rdf == 0 : perfect fit!
	cat("ALL", df[1], "residuals are 0: no residual degrees of freedom!\n")
    }
    if (nsingular <- df[3] - df[1])
	cat("\nCoefficients: (", nsingular,
	    " not defined because of singularities)\n", sep = "")
    else cat("\nCoefficients:\n")

    print.coefmat(x$coef, digits=digits, signif.stars=signif.stars, ...)
    ##
    cat("\nResidual standard error:",
	format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
    if (!is.null(x$fstatistic)) {
	cat("Multiple R-Squared:", formatC(x$r.squared, digits=digits))
	cat(",\tAdjusted R-squared:",formatC(x$adj.r.squared,d=digits),
	    "\nF-statistic:", formatC(x$fstatistic[1], digits=digits),
	    "on", x$fstatistic[2], "and",
	    x$fstatistic[3], "degrees of freedom,\tp-value:",
	    formatC(1 - pf(x$fstatistic[1], x$fstatistic[2],
			   x$fstatistic[3]), dig=digits),
	    "\n")
    }

    cat("\n")
    invisible(x)
}

fitted.wle.lm <- function(object, ...) object$fitted.values
coef.wle.lm <- function(object, ...) object$coefficients
weights.wle.lm <- .Alias(weights.default)
formula.wle.lm <- function(object, ...) formula(object$terms)

model.frame.wle.lm <-
    function(formula, data, na.action, ...) {
	if (is.null(formula$model)) {
	    fcall <- formula$call
	    fcall$method <- "model.frame"
	    fcall[[1]] <- as.name("wle.lm")
	    eval(fcall, sys.frame(sys.parent()))
	}
	else formula$model
    }

#############################################################
#                                                           #
#	WLE.NORMAL function                                 #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.normal <- function(x, boot=30, group, num.sol=1, raf="HD", smooth=0.003, tol=10^(-6), equal=10^(-3), max.iter=500)
{

raf <- switch(raf,
	HD = 1,
	NED = 2,
	SCHI2 = 3,
	-1)

if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")

if (missing(group)) {
group <- 0
}

x <- as.vector(x)
size <- length(x)
result <- list()

if (size<=2) {
stop("Number of observation must be at least equal to 2")
}

if (!(group>1)) {
group <- max(round(size/4),3)
cat("wle.normal: dimension of the subsample set to default value: ",group,"\n")
}

maxboot <- sum(log(1:size))-(sum(log(1:group))+sum(log(1:(size-group))))

if (boot<1 | log(boot) > maxboot) {
stop("Bootstrap replication not in the range")
}

if (!(num.sol>=1)) {
cat("wle.normal: number of solution to report set to 1 \n")
num.sol <- 1
}

if (max.iter<1) {
cat("wle.normal: max number of iteration set to 500 \n")
max.iter <- 500
}

if (smooth<10^(-5)) {
cat("wle.normal: the smooth parameter seems too small \n")
}

if (tol<0) {
cat("wle.normal: the accuracy can not be negative, using default value \n")
tol <- 10^(-6)
}

if (equal<0) {
cat("wle.normal: the equal parameter can not be negative, using default value \n")
equal <- 10^(-3)
}

  z <- .Fortran("wlenorm",
	as.double(x), 
	as.integer(size),
	as.integer(boot),
	as.integer(group),
	as.integer(num.sol),
	as.integer(raf),
	as.double(smooth),
	as.double(tol),
	as.double(equal),
	as.integer(max.iter),
	mean=double(num.sol),
	var=double(num.sol),
	totweight=double(num.sol),
	weight=mat.or.vec(num.sol,size),
	same=integer(num.sol),
	nsol=integer(1),
	nconv=integer(1),
	PACKAGE = "wle")

if (z$nsol>0) {
result$location <- z$mean[1:z$nsol]
result$scale <- sqrt(z$var[1:z$nsol])
result$tot.weights <- z$totweight[1:z$nsol]
result$weights <- z$weight[1:z$nsol,]
result$freq <- z$same[1:z$nsol]
result$tot.sol <- z$nsol
result$not.conv <- z$nconv
result$call <- match.call()

class(result) <- "wle.normal"

return(result)
} else{
stop("No solutions are fuond, checks the parameters")
}
}

print.wle.normal <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("Location:\n")
    print.default(format(x$location, digits=digits),
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("Scale:\n")
    print.default(format(x$scale, digits=digits),
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("\nNumber of solutions ",x$tot.sol,"\n")
    cat("\n")
    invisible(x)
}



#############################################################
#                                                           #
#	WLE.NORMAL.MULTI function                           #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.normal.multi <- function(x, boot=30, group, num.sol=1, raf="HD", smooth, tol=10^(-6), equal=10^(-3), max.iter=500)
{

raf <- switch(raf,
	HD = 1,
	NED = 2,
	SCHI2 = 3,
	-1)

if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")

if (missing(group)) {
group <- 0
}

if (is.null(size <- nrow(x)) | is.null(nvar <- ncol(x))) {
    if (is.vector(x)) {
	return(wle.normal(x=x, boot=boot, group=group, num.sol=num.sol, raf=raf, smooth=smooth, tol=tol, equal=equal, max.iter=max.iter)) 
    } else {
	stop("'x' must be a matrix or a vector")
    }
}

if (missing(smooth)) {
    smooth <- wle.smooth(dimension=nvar,costant=4,weight=0.5,interval=c(0.0001,20))$root
}

result <- list()

if (size<nvar*(nvar+1)) {
stop(paste("Number of observation must be at least equal to ",nvar*(nvar+1)))
}
if (!(group>nvar*(nvar+1))) {
group <- max(round(size/4),nvar*(nvar+1))
cat("wle.normal.multi: dimension of the subsample set to default value: ",group,"\n")
}

maxboot <- sum(log(1:size))-(sum(log(1:group))+sum(log(1:(size-group))))

if (boot<1 | log(boot) > maxboot) {
stop("Bootstrap replication not in the range")
}

if (!(num.sol>=1)) {
cat("wle.normal.multi: number of solution to report set to 1 \n")
num.sol <- 1
}

if (max.iter<1) {
cat("wle.normal.multi: max number of iteration set to 500 \n")
max.iter <- 500
}

if (smooth<10^(-5)) {
cat("wle.normal.multi: the smooth parameter seems too small \n")
}

if (tol<0) {
cat("wle.normal.multi: the accuracy can not be negative, using default value \n")
tol <- 10^(-6)
}

if (equal<0) {
cat("wle.normal.multi: the equal parameter can not be negative, using default value \n")
equal <- 10^(-3)
}

  z <- .Fortran("wlenormmulti",
	as.double(x), 
	as.integer(size),
	as.integer(nvar),
	as.integer(boot),
	as.integer(group),
	as.integer(num.sol),
	as.integer(raf),
	as.double(smooth),
	as.double(tol),
	as.double(equal),
	as.integer(max.iter),
	mean=mat.or.vec(num.sol,nvar),
	var=mat.or.vec(num.sol,nvar*nvar),
	totweight=double(num.sol),
	weight=mat.or.vec(num.sol,size),
	same=integer(num.sol),
	nsol=integer(1),
	nconv=integer(1))

if (z$nsol>0) {

dn <- colnames(x)

temp <- z$var[1:z$nsol,]
if (z$nsol>1) {
temp.a <- matrix(temp[1,],ncol=nvar)
dimnames(temp.a) <- list(dn,dn)
temp.b <- list(temp.a)
for (i in 2:z$nsol) {
temp.a <- matrix(temp[i,],ncol=nvar)
dimnames(temp.a) <- list(dn,dn)
temp.b <- c(temp.b,list(temp.a))
}
} else {
temp.a <- matrix(temp,ncol=nvar)
dimnames(temp.a) <- list(dn,dn)
temp.b <- list(temp.a)
}

result$location <- z$mean[1:z$nsol,]
result$variance <- temp.b
result$tot.weights <- z$totweight[1:z$nsol]
result$weights <- z$weight[1:z$nsol,]
result$freq <- z$same[1:z$nsol]
result$smooth <- smooth
result$tot.sol <- z$nsol
result$not.conv <- z$nconv
result$call <- match.call()


if (is.null(nrow(result$location))) {
names(result$location) <- dn
} else {
dimnames(result$location) <- list(NULL,dn)
}


class(result) <- "wle.normal.multi"

return(result)
} else {
stop("No solutions are fuond, checks the parameters")
}
}

print.wle.normal.multi <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("Location:\n")
    print.default(format(x$location, digits=digits),
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("\nVariance-Covariance matrix:\n")
    print.default(x$variance, digits=digits,
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("\nNumber of solutions ",x$tot.sol,"\n")
    cat("\n")
    invisible(x)
}







#############################################################
#                                                           #
#	WLE.ONESTEP function                                #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                             #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.onestep <- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, ini.param, ini.scale, raf="HD", smooth=0.0320018, num.step=1, contrasts=NULL)
{

raf <- switch(raf,
	HD = 1,
	NED = 2,
	SCHI2 = 3,
	-1)

if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$ini.param <- mf$ini.scale <- mf$smooth <- NULL
    mf$num.step <- mf$raf <- mf$contrasts <- NULL
    mf$model <- mf$x <- mf$y <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

if (size<(nvar+1)) {stop("Number of observation must be at least equal to the number of predictors (including intercept) + 1")}

if (!(ini.scale>=0)) {
stop("The initial scale error must be non negative")
}

if (!(num.step>=1)) {
cat("wle.onestep: number of steps can not be negative, set to 1 \n")
num.step <- 1
}

if (smooth<10^(-5)) {
cat("wle.onestep: the smooth parameter seems too small \n")
}

ini.var <- ini.scale^2

  z <- .Fortran("wleonestepfix",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(nvar),
	as.double(ini.param),
	as.double(ini.var),
	as.integer(raf),
	as.double(smooth),
	as.integer(num.step),
	param=double(nvar),
	var=double(1),
	resid=double(size),
	totweight=double(1),
	weight=double(size))

if(z$var>0) {

devparam <- sqrt(z$var*diag(solve(t(xdata)%*%diag(z$weight)%*%xdata)))

result$coefficients <- z$param
result$standard.error <- devparam
result$scale <- sqrt(z$var)
result$residuals <- z$resid
result$fitted.values <- as.vector(xdata%*%z$param)
result$weights <- z$weight
result$tot.weights <- z$totweight
result$call <- cl
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

if (is.null(names(ini.param))) {
dn <- colnames(xdata)
} else {
dn <- names(ini.param)
}

if (is.null(nrow(result$coefficients))) {
names(result$coefficients) <- dn
} else {
dimnames(result$coefficients) <- list(NULL,dn)
}

class(result) <- "wle.onestep"

return(result)
} else {
stop("The initial estimates do not seems very good: the total sum of the weights is less than number of independent variables")
}
}

print.wle.onestep <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("Coefficients:\n")
    print.default(format(coef(x), digits=digits),
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("Scale estimate: ",format(x$scale, digits=digits))
    cat("\n")

    invisible(x)
}
#############################################################
#                                                           #
#	WLE.SMOOTH function                                 #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: October, 10, 2000                             #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.smooth <- function(weight=0.31,costant=3,level=0.2,dimension=1,raf="HD",interval=c(0.00001,0.5),tol=10^-6,max.iter=1000)
{

raf <- switch(raf,
	HD = 1,
	NED = 2,
	SCHI2 = 3,
	-1)

if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")

delta <- function(smooth,costant,level,dimension){
level*(((smooth+1)/smooth)^(dimension/2)*exp(costant^2/(2*(dimension+1)))-1)}

if (raf==3) {w <- function(smooth,costant,level,dimension,weight){(1-((delta(smooth,costant,level,dimension)**2)/((delta(smooth,costant,level,dimension)**2) + 2)))-weight}
} else {
if (raf==2) {
adelta <- function(d) {2-(2+d)*exp(-d)} 
} else {
adelta <- function(d) {2*(sqrt(d+1)-1)}
}
w <- function(smooth,costant,level,dimension,weight){
(adelta(delta(smooth,costant,level,dimension))+1)/(delta(smooth,costant,level,dimension)+1)-weight
}
}

result <- uniroot(w,interval=interval,costant=costant,level=level,dimension=dimension,weight=weight,maxiter=max.iter,tol=tol)

result$call <- match.call()

class(result) <- "wle.smooth"

return(result)
}

print.wle.smooth <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("\nBandwidth: ",format(x$root, digits=digits))
    cat("\n")
    invisible(x)
}



#############################################################
#                                                           #
#	WLE.STEPWISE function                               #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: December, 19, 2000                            #
#	Version: 0.3                                        #
#                                                           #
#	Copyright (C) 2000 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.stepwise <- function(formula, data=list(), model=TRUE, x=FALSE, y=FALSE, boot=30, group, num.sol=1, raf="HD", smooth=0.031, tol=10^(-6), equal=10^(-3), max.iter=500, min.weight=0.5, type="Forward", f.in=4.0, f.out=4.0, method="WLE", contrasts=NULL)
{

raf <- switch(raf,
	HD = 1,
	NED = 2,
	SCHI2 = 3,
	-1)

if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")

ntype <- switch(type,
	Forward = 1,
	Backward = 2,
	Stepwise = 3,
	-1)

if (ntype==-1) stop("The type must be Forward, Backward or Stepwise")

nmethod <- switch(method,
		WLE = 0,
	        WLS = 1,
		-1)

if (nmethod==-1) stop("The method must be WLE, or WLS the default value is WLE")

if (missing(group)) {
group <- 0
}

    ret.x <- x
    ret.y <- y
    result <- list()	
    mt <- terms(formula, data = data)
    mf <- cl <- match.call()
    mf$boot <- mf$group <- mf$smooth <- NULL
    mf$tol <- mf$equal <- mf$num.sol <- NULL
    mf$max.iter <- mf$raf <- mf$contrasts <- NULL
    mf$min.weight <- NULL
    mf$type <- mf$f.in <- mf$f.out <- NULL
    mf$model <- mf$x <- mf$y <- mf$method <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    xvars <- as.character(attr(mt, "variables"))[-1]
    inter <- attr(mt, "intercept")
    if((yvar <- attr(mt, "response")) > 0) xvars <- xvars[-yvar]
    xlev <-
	if(length(xvars) > 0) {
	    xlev <- lapply(mf[xvars], levels)
	    xlev[!sapply(xlev, is.null)]
	}
    ydata <- model.response(mf, "numeric")
    if (is.empty.model(mt)) 
	stop("The model is empty")
    else 
	xdata <- model.matrix(mt, mf, contrasts)

if (is.null(size <- nrow(xdata)) | is.null(nvar <- ncol(xdata))) stop("'x' must be a matrix")
if (length(ydata)!=size) stop("'y' and 'x' are not compatible")

if (size<nvar+1) {
stop("Number of observations must be at least equal to the number of predictors (including intercept) + 1")
}

if (f.in<0 | f.out<0) {
stop("f.in and f.out can not be negative")
}

if (!(group>1)) {
group <- max(round(size/4),nvar+2)
cat("wle.stepwise: dimension of the subsample set to default value: ",group,"\n")
}

maxboot <- sum(log(1:size))-(sum(log(1:group))+sum(log(1:(size-group))))

if (boot<1 | log(boot) > maxboot) {
stop("Bootstrap replication not in the range")
}

if (!(num.sol>=1)) {
cat("wle.stepwise: number of solution to report set to 1 \n")
num.sol <- 1
}

if (max.iter<1) {
cat("wle.stepwise: max number of iteration set to 500 \n")
max.iter <- 500
}

if (smooth<10^(-5)) {
cat("wle.stepwise: the smooth parameter seems too small \n")
}

if (tol<0) {
cat("wle.stepwise: the accuracy can not be negative, using default value \n")
tol <- 10^(-6)
}

if (equal<0) {
cat("wle.stepwise: the equal parameter can not be negative, using default value \n")
equal <- 10^(-3)
}

if (min.weight<0) {
cat("wle.stepwise: the minimum sum of the weights can not be negative, using default value \n")
min.weight <- 0.5
}

nrep <- 2^nvar-1

  z <- .Fortran("wstep",
	as.double(ydata),
	as.matrix(xdata),
	as.integer(0), 
	as.integer(size),
	as.integer(nvar),
	as.integer(boot),
	as.integer(group),
	as.integer(nrep),
	as.integer(raf),
	as.double(smooth),
	as.integer(ntype),
	as.double(tol),
	as.double(equal),
	as.integer(max.iter),
	as.integer(num.sol),
	as.double(min.weight),
	as.double(f.in),
	as.double(f.out),
	as.integer(nmethod),
	wstep=mat.or.vec(nrep,nvar+1),
	param=mat.or.vec(num.sol,nvar),
	var=double(num.sol),
	resid=mat.or.vec(num.sol,size),
	totweight=double(num.sol),
	weight=mat.or.vec(num.sol,size),
	same=integer(num.sol),
	indice=integer(1),
	info=integer(1),
	imodel=integer(1),
	nsol=integer(1))

result$wstep <- z$wstep[1:z$imodel,]
result$coefficients <- z$param[1:z$nsol,]
result$scale <- sqrt(z$var[1:z$nsol])
result$residuals <- z$resid[1:z$nsol,]
result$tot.weights <- z$totweight[1:z$nsol]
result$weights <- z$weight[1:z$nsol,]
result$freq <- z$same[1:z$nsol]
result$index <- z$indice
result$info <- z$info
result$call <- cl
result$contrasts <- attr(xdata, "contrasts")
result$xlevels <- xlev
result$terms <- mt
result$type <- type
result$method <- method
result$f.in <- f.in
result$f.out <- f.out

if (model)
    result$model <- mf
if (ret.x)
    result$x <- xdata
if (ret.y)
    result$y <- ydata

dn <- colnames(xdata)

if (is.null(nrow(result$coefficients))) {
names(result$coefficients) <- dn
} else {
dimnames(result$coefficients) <- list(NULL,dn)
}

if (z$imodel<=1) {
names(result$wstep) <- c(dn," ")
} else {
dimnames(result$wstep) <- list(NULL,c(dn," "))
}

class(result) <- "wle.stepwise"

return(result)

}

summary.wle.stepwise <- function (object, num.max=20, ...) {

z <- .Alias(object)
if (is.null(z$terms)) {
    stop("invalid \'wle.stepwise\' object")
}

if (num.max<1) {
cat("summary.wle.stepwise: num.max can not less than 1, num.max=1 \n")
num.max <- 1
}

ans <- list()
wstep <- z$wstep
if(is.null(nmodel <- nrow(wstep))) nmodel <- 1
num.max <- min(nmodel,num.max)
if (nmodel!=1) { 
    wstep <- wstep[(nmodel-num.max+1):nmodel,]
}

ans$wstep <- wstep
ans$num.max <- num.max
ans$type <- z$type
ans$f.in <- z$f.in
ans$f.out <- z$f.out
ans$call <- z$call

class(ans) <- "summary.wle.stepwise"
return(ans)
}

print.wle.stepwise <- function (x, digits = max(3, getOption("digits") - 3), ...) {
res_summary.wle.stepwise(object=x, num.max=nrow(x$wstep), ...)
print.summary.wle.stepwise(res, digits=digits, ...)
}

print.summary.wle.stepwise <- function (x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")

    cat("\n",x$type," selection procedure\n")
    if (x$type=="Forward" | x$type=="Stepwise") {
	cat("\nF.in: ",x$f.in)
    } 
    if (x$type=="Backward" | x$type=="Stepwise") {
	cat("\nF.out: ",x$f.out)
    }
    cat(" \n")
    cat("\nLast ",x$num.max," iterations:\n")

    if(x$num.max>1) {
    nvar <- ncol(x$wstep)-1
    x$wstep[,(nvar+1)] <- signif(x$wstep[,(nvar+1)],digits)
    } else {
    nvar <- length(x$wstep)-1
    x$wstep[(nvar+1)] <- signif(x$wstep[(nvar+1)],digits)
    }
    print(x$wstep)
    cat("\n")
    invisible(x)
}







.First.lib <- function(lib, pkg) {
  library.dynam("wle", pkg, lib)
}
