#############################################################
#                                                           #
#	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) {
stop("Number of observations must be at least equal to the number of predictors (including intercept)")
}

if (group<nvar) {
group <- max(round(size/4),nvar)
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.BINOMIAL function                               #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: February, 13, 2001                            #
#	Version: 0.1                                        #
#                                                           #
#	Copyright (C) 2001 Claudio Agostinelli              #
#                                                           #
#############################################################

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

result <- list()

if (raf!="HD" & raf!="NED" & raf!="SCHI2") 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)
nsize <- length(x)
result <- list()

if (nsize<1) {
stop("Number of observation must be at least equal to 1")
}

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

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

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

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

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

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

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

tot.sol <- 0
not.conv <- 0
iboot <- 0

while (tot.sol < num.sol & iboot < boot) {
   iboot <- iboot + 1
   x.boot <- x[round(runif(group,0.501,nsize+0.499))]
   p <- sum(x.boot)/(size*group)

   ff <- rep(0,nsize)
   x.diff <- tol + 1
   iter <- 0
   while (x.diff > tol & iter < max.iter) {
   iter <- iter + 1
   p.old <- p 
       tff <- table(x)/nsize
       nff <- as.numeric(names(tff))
       for (i in 1:nsize) {
           ff[i] <- tff[nff==x[i]] 
       }
       mm <- dbinom(x,size=size,prob=p)
       dd <- ff/mm - 1
       
       ww <- switch(raf,
                 HD =  2*(sqrt(dd + 1) - 1) ,
	         NED =  2 - (2 + dd)*exp(-dd) ,
	         SCHI2 =  1-(dd^2/(dd^2 +2)) )       

       if (raf=="HD" | raf=="NED") {
            ww <- (ww + 1)/(dd + 1)
       }

       ww[ww > 1] <- 1
       ww[ww < 0] <- 0

       p <- ww%*%x/(sum(ww)*size)

       x.diff <- abs(p - p.old)
   }
#### end of while (x.diff > tol & iter < max.iter)

   if (iter < max.iter) {

   if (tot.sol==0) {
      p.store <- p
      w.store <- ww
      tot.sol <- 1
   } else {
      if (min(abs(p.store-p))>equal) {
          p.store <- c(p.store,p)
          w.store <- rbind(w.store,ww)
          tot.sol <- tot.sol + 1
      }
   }

   } else not.conv <- not.conv + 1
   

}
##### end of while (tot.sol < num.sol & iboot < boot)

if (tot.sol) {
result$p <- p.store
result$tot.weights <- sum(ww)/nsize
result$weights <- w.store
result$tot.sol <- tot.sol
result$not.conv <- not.conv
result$call <- match.call()

class(result) <- "wle.binomial"

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

print.wle.binomial <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("p:\n")
    print.default(format(x$p, digits=digits),
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("\nNumber of solutions ",x$tot.sol,"\n")
    cat("\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) {
stop("Number of observations must be at least equal to the number of predictors (including intercept)")
}

if (group<nvar) {
group <- max(round(size/4),nvar)
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) {
stop("Number of observations must be at least equal to the number of predictors (including intercept)")
}

if (group<nvar) {
group <- max(round(size/4),nvar)
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.GAMMA function                                  #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: February, 16, 2001                            #
#	Version: 0.1                                        #
#                                                           #
#	Copyright (C) 2001 Claudio Agostinelli              #
#                                                           #
#############################################################

wle.gamma <- function(x, boot=30, group, num.sol=1, raf="HD", smooth=0.008, tol=10^(-6), equal=10^(-3), max.iter=500, shape.int=c(0.01,100), shape.tol=10, use.smooth=TRUE, tol.int) {

wsolve <- function (o, media, medialog) {
   medialog + log(o/media) - digamma(o)
}

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<2) {
group <- max(round(size/4),2)
cat("wle.gamma: 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.gamma: number of solution to report set to 1 \n")
num.sol <- 1
}

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

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

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

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

if (!is.logical(use.smooth)) {
cat("wle.gamma: the use.smooth must be a logical value, using default value \n")
use.smooth <- TRUE
}

if (length(shape.int)!=2) stop("shape.int must be a vector of length 2 \n")

shape.int <- rev(sort(shape.int))
shape.first <- shape.int

if (shape.int[2] <= 0) {
stop("the elements of shape.int must be positive \n")
}

if (shape.int[1] <= 0) {
cat("wle.gamma: the elements of shape.int must be positive, using default value \n")
shape.int[1] <- tol
}

if (shape.tol <=0) {
cat("wle.gamma: shape.tol must be positive, using default value \n")
shape.tol <- 10
} 

if (missing(tol.int)) {
   tol.int <- tol*10^(-4)
} else {
   if (tol.int <=0) {
     cat("wle.gamma: tol.int must be positive, using default value \n")
     tol.int <- tol*10^(-4) 
   } 
}

tot.sol <- 0
not.conv <- 0
iboot <- 0

xlog <- log(x)

while (tot.sol < num.sol & iboot < boot) {
   iboot <- iboot + 1
   x.boot <- x[round(runif(group,0.501,size+0.499))]
   xlog.boot <- log(x.boot)

   media <- sum(x.boot)/group
   medialog<- sum(xlog.boot)/group

   o <- uniroot(wsolve,interval=shape.first, media=media, medialog=medialog)$root

   if (o < tol) o <- 2*tol

   l <- media/o

xdiff <- tol + 1
iter <- 0
while (xdiff > tol & iter < max.iter) {

iter <- iter + 1
shape <- o
lambda <- 1/l
temp <- shape/lambda^2
dsup <- max(x)+ 3*smooth*temp

   z <- .Fortran("wlegamma",
	as.double(x), 
	as.integer(size),
	as.integer(raf),
	as.double(smooth*temp),
        as.integer(1*use.smooth),
        as.double(dsup),
	as.double(tol),
        as.double(tol.int),
	as.double(lambda),
	as.double(shape),
	weights=double(size),
	density=double(size),
	model=double(size))

ww <- z$weights
wsum <- sum(ww)
wmedia <- ww%*%x/wsum
wmedialog <- ww%*%xlog/wsum

shape.int <- c(max(tol,(o-shape.tol)),(o+shape.tol))

o <- uniroot(wsolve, interval=shape.int, media=wmedia, medialog=wmedialog)$root

if (o < tol) o <- 2*tol

l <- wmedia/o

xdiff <- max(abs(c(o-shape,l-1/lambda)))
}

   if (iter < max.iter) {

   if (tot.sol==0) {
      o.store <- o
      l.store <- l
      w.store <- ww
      tot.sol <- 1
   } else {
      if (min(abs(o.store-o))>equal & min(abs(l.store-l))>equal) {
          o.store <- c(o.store,o)
          l.store <- c(l.store,l)
          w.store <- rbind(w.store,ww)
          tot.sol <- tot.sol + 1
      }
   }

   } else not.conv <- not.conv + 1
   

}
##### end of while (tot.sol < num.sol & iboot < boot)



if (tot.sol) {
   result$scale <- c(l.store)
   result$shape <- o.store
   if (tot.sol>1) {
       tot.w <- apply(w.store,1,sum)/size
   } else tot.w <- sum(w.store)/size
  
   result$tot.weights <- tot.w
   result$weights <- w.store
   result$tot.sol <- tot.sol
   result$not.conv <- not.conv
   result$call <- match.call()

   class(result) <- "wle.gamma"

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


print.wle.gamma <- function(x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("Scale:\n")
    print.default(format(x$scale, digits=digits),
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("Shape:\n")
    print.default(format(x$shape, digits=digits),
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("\nNumber of solutions ",x$tot.sol,"\n")
    cat("\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) {
stop("Number of observations must be at least equal to the number of predictors (including intercept)")
}

if (group<nvar) {
group <- max(round(size/4),nvar)
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<2) {
group <- max(round(size/4),2)
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])
     if (z$nsol>1) {
          result$residuals <- matrix(rep(x,z$nsol),nrow=z$nsol,byrow=TRUE) - matrix(rep(z$mean[1:z$nsol],size),nrow=z$nsol,byrow=FALSE)
     } else {
          result$residuals <- x - result$location
     }
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)/2+nvar)) {
stop(paste("Number of observation must be at least equal to ",nvar*nvar))
}
if (group<(nvar*(nvar+1)/2+nvar)) {
group <- max(round(size/4),(nvar*(nvar+1)/2+nvar))
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.POISSON function                                #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: February, 16, 2001                            #
#	Version: 0.1                                        #
#                                                           #
#	Copyright (C) 2001 Claudio Agostinelli              #
#                                                           #
#############################################################

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

result <- list()

if (raf!="HD" & raf!="NED" & raf!="SCHI2") 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<1) {
stop("Number of observation must be at least equal to 1")
}

if (group<1) {
group <- max(round(size/4),1)
cat("wle.poisson: 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.poisson: number of solution to report set to 1 \n")
num.sol <- 1
}

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

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

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

tot.sol <- 0
not.conv <- 0
iboot <- 0

while (tot.sol < num.sol & iboot < boot) {
   iboot <- iboot + 1
   x.boot <- x[round(runif(group,0.501,size+0.499))]
   p <- sum(x.boot)/group

   ff <- rep(0,size)
   x.diff <- tol + 1
   iter <- 0
   while (x.diff > tol & iter < max.iter) {
   iter <- iter + 1
   p.old <- p 
       tff <- table(x)/size
       nff <- as.numeric(names(tff))
       for (i in 1:size) {
           ff[i] <- tff[nff==x[i]] 
       }
       mm <- dpois(x,lambda=p)
       dd <- ff/mm - 1
       
       ww <- switch(raf,
                 HD =  2*(sqrt(dd + 1) - 1) ,
	         NED =  2 - (2 + dd)*exp(-dd) ,
	         SCHI2 =  1-(dd^2/(dd^2 +2)) )       

       if (raf=="HD" | raf=="NED") {
            ww <- (ww + 1)/(dd + 1)
       }

       ww[ww > 1] <- 1
       ww[ww < 0] <- 0

       p <- ww%*%x/sum(ww)

       x.diff <- abs(p - p.old)
   }
#### end of while (x.diff > tol & iter < max.iter)

   if (iter < max.iter) {

   if (tot.sol==0) {
      p.store <- p
      w.store <- ww
      tot.sol <- 1
   } else {
      if (min(abs(p.store-p))>equal) {
          p.store <- c(p.store,p)
          w.store <- rbind(w.store,ww)
          tot.sol <- tot.sol + 1
      }
   }

   } else not.conv <- not.conv + 1
   

}
##### end of while (tot.sol < num.sol & iboot < boot)

if (tot.sol) {
result$lambda <- c(p.store)
result$tot.weights <- sum(ww)/size
result$weights <- w.store
result$tot.sol <- tot.sol
result$not.conv <- not.conv
result$call <- match.call()

class(result) <- "wle.poisson"

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

print.wle.poisson <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("lambda:\n")
    print.default(format(x$lambda, digits=digits),
		  print.gap = 2, quote = FALSE)
    cat("\n")
    cat("\nNumber of solutions ",x$tot.sol,"\n")
    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) {
stop("Number of observations must be at least equal to the number of predictors (including intercept)")
}

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)
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)
}







#############################################################
#                                                           #
#	WLE.T.TEST function                                 #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: Febraury 9, 2001                              #
#	Version: 0.1                                        #
#                                                           #
#	Based on t.test function in                         #
#       ctest package version 1.2.0                         #
#                                                           #
#############################################################

wle.t.test <- function(x, y=NULL, alternative = c("two.sided", "less", "greater"), mu=0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, boot=30, group, num.sol=1, raf="HD", smooth=0.003, tol=10^(-6), equal=10^(-3), max.iter=500) {

    alternative <- match.arg(alternative)
    wy <- NULL
    x.out <- x
    y.out <- y

    if(!missing(mu) && (length(mu) != 1 || is.na(mu)))
        stop("mu must be a single number")
    if(!missing(conf.level) &&
       (length(conf.level) != 1 || !is.finite(conf.level) ||
        conf.level < 0 || conf.level > 1))
        stop("conf.level must be a single number between 0 and 1")
    if( !is.null(y) ) {
	dname <- paste(deparse(substitute(x)),"and",
		       deparse(substitute(y)))
	if(paired)
	    xok <- yok <- complete.cases(x,y)
	else {
	    yok <- !is.na(y)
	    xok <- !is.na(x)
	}
	y <- y[yok]
    }
    else {
	dname <- deparse(substitute(x))
	if( paired ) stop("y is missing for paired test")
	xok <- !is.na(x)
	yok <- NULL
    }
    x <- x[xok]
    if( paired ) {
	x <- x-y
	y <- NULL
    }

    nx <- length(x)
    if(nx < 2) stop("not enough x observations")
    x.est <- wle.normal(x=x,boot=boot, group=group, num.sol=num.sol, raf=raf, smooth=smooth, tol=tol, equal=equal, max.iter=max.iter)
    x.tot.sol <- x.est$tot.sol
    y.tot.sol <- 1
    y.root <- 1
     
    if( !is.null(y) ) {
    ny <- length(y)
    if(ny < 2) stop("not enough y observations")
    y.est <- wle.normal(x=y,boot=boot, group=group, num.sol=num.sol, raf=raf, smooth=smooth, tol=tol, equal=equal, max.iter=max.iter)
    y.tot.sol <- y.est$tot.sol
    }

    wtstat <- matrix(0,ncol=y.tot.sol,nrow=x.tot.sol)
    wdf <- matrix(0,ncol=y.tot.sol,nrow=x.tot.sol)
    num.col <- 1
    if (!is.null(y) & paired==FALSE) num.col <- 2
    westimate <- matrix(0,nrow=x.tot.sol*y.tot.sol,ncol=num.col) 

    for (x.root in 1:x.tot.sol) {

    if (x.tot.sol==1) {
    mx <- x.est$location
    vx <- x.est$scale^2
    wx <- x.est$weights
    if (is.nan(mx) | is.nan(vx)) stop("no solutions are found for location and scale of 'x'")
    } else {
    mx <- x.est$location[x.root]
    vx <- (x.est$scale^2)[x.root]
    wx <- x.est$weights[x.root,]    
    }

    sumwx <- sum(wx)

    estimate <- mx
    if(is.null(y)) {
	df <- sumwx-1
	stderr <- sqrt(vx/sumwx)
	tstat <- (mx-mu)/stderr
	method <- ifelse(paired,"Paired wt-test for normal distributed data","One Sample wt-test for normal distributed data")
	names(estimate) <- ifelse(paired,"mean of the differences","mean of x")
    wtstat[x.root,y.root] <- tstat
    wdf[x.root,y.root] <- df
    westimate[x.root*y.root,] <- estimate

    } else {

        for (y.root in 1:y.tot.sol) {

        if (y.tot.sol==1) {
            my <- y.est$location
            vy <- y.est$scale^2
            wy <- y.est$weights
            if (is.nan(my) | is.nan(vy)) stop("no solutions are found for location and scale of 'y'")
        } else {
            my <- y.est$location[y.root]
            vy <- (y.est$scale^2)[y.root]
            wy <- y.est$weights[y.root,]    
        }

        sumwy <- sum(wy)

	method <- paste(if(!var.equal) "Welch", "Two Sample wt-test for normal distributed data")
	estimate <- c(mx,my)
	names(estimate) <- c("mean of x","mean of y")
	if(var.equal) {
	    df <- sumwx+sumwy-2
	    v <- ((sumwx-1)*vx + (sumwy-1)*vy)/df
	    stderr <- sqrt(v*(1/sumwx+1/sumwy))
	} else {
	    stderrx <- sqrt(vx/sumwx)
	    stderry <- sqrt(vy/sumwy)
	    stderr <- sqrt(stderrx^2 + stderry^2)
	    df <- stderr^4/(stderrx^4/(sumwx-1) + stderry^4/(sumwy-1))
	}
        tstat <- (mx - my - mu)/stderr
        wtstat[x.root,y.root] <- tstat
        wdf[x.root,y.root] <- df
        westimate[x.root*y.root,] <- estimate
    }
    }
    }

    result <- list()

    for (x.root in 1:x.tot.sol) {

    result.x <- list()
    for (y.root in 1:y.tot.sol) {

    tstat <- wtstat[x.root,y.root]
    df <- wdf[x.root,y.root]
    estimate <- westimate[x.root*y.root,]

    if (alternative == "less") {
	pval <- pt(tstat, df)
	cint <- c(NA, tstat + qt(conf.level, df) )
    }
    else if (alternative == "greater") {
	pval <- pt(tstat, df, lower = FALSE)
	cint <- c(tstat - qt(conf.level, df), NA)
    }
    else {
	pval <- 2 * pt(-abs(tstat), df)
	alpha <- 1 - conf.level
        cint <- qt(1 - alpha/2, df)
	cint <- tstat + c(-cint, cint)
    }
    cint <- mu + cint * stderr
    names(tstat) <- "wt"
    names(df) <- "df"
    names(mu) <- if(paired || !is.null(y)) "difference in means" else "mean"
    attr(cint,"conf.level") <- conf.level
    if (x.tot.sol>1) {
       wx <- x.est$weights[x.root,]
    } else {
       wx <- x.est$weights
    }


    if (!is.null(y)) { 
       if (y.tot.sol>1) {
           wy <- y.est$weights[y.root,]
       } else {
           wy <- y.est$weights
       }
    }

    rval <- list(statistic = tstat, parameter = df, p.value = pval,
	       conf.int=cint, estimate=estimate, null.value = mu,
	       alternative=alternative,
	       method=method, data.name=dname, 
               x.weights=wx, y.weights=wy, 
               x.root=x.root, y.root=y.root)

    class(rval) <- "htest"

    result.x <-  c(result.x,list(rval))
    }
    result$test <- c(result$test,list(result.x))
    result.x <- list()
    }

    result$x.tot.sol <- x.tot.sol
    result$y.tot.sol <- y.tot.sol
    result$call <- match.call()   
    result$paired <- paired
    result$x <- x.out
    result$y <- y.out

    class(result) <- "wle.t.test"

    return(result)
}


print.wle.t.test <- function(x, x.root="ALL", y.root="ALL", digits = 4, quote = TRUE, prefix = "", ...) {

x.tot.sol <- x$x.tot.sol
y.tot.sol <- x$y.tot.sol

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

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

cat("\nCall:\n")
cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="") 
cat("\n\nWeighted t test:\n", sep="")

for (xx.root in x.root) {
for (yy.root in y.root) {
    cat("\n'x' Root ",xx.root)
    if (!is.null(x$y) & x$paired==FALSE) cat (" 'y' Root ",yy.root)

    print.htest(x$test[[xx.root]][[yy.root]], digits=digits, quote=quot, prefix=prefix, ...)
}
}
 
    cat("\nNumber of 'x' solutions ",x.tot.sol,"\n")
    if (!is.null(x$y) & x$paired==FALSE) cat("\nNumber of 'y' solutions ",y.tot.sol,"\n")
    cat("\n")
    invisible(x)

}


#############################################################
#                                                           #
#	WLE.VAR.TEST function                               #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@stat.unipd.it                       #
#	Date: Febraury 21, 2001                             #
#	Version: 0.1                                        #
#                                                           #
#	Based on var.test function in                       #
#       ctest package version 1.2.0                         #
#                                                           #
#############################################################

wle.var.test <- function(x, y, ratio = 1, alternative = c("two.sided", "less", "greater"), conf.level = 0.95, x.root=1, y.root=1) {

    if (!((length(ratio) == 1) && is.finite(ratio) && (ratio > 0)))
        stop("ratio must be a single positive number")

    alternative <- match.arg(alternative)

    if (!((length(conf.level) == 1) && is.finite(conf.level) &&
          (conf.level > 0) && (conf.level < 1)))
        stop("conf.level must be a single number between 0 and 1")

    DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))

    if (inherits(x, "wle.lm") && inherits(y, "wle.lm")) {

        x.tot.sol <- x$tot.sol
        if (x.tot.sol<x.root) {
            stop(paste("'x' Root ",root," not found"))
        }
        if (x.tot.sol!=1) {
            x.res <- x$residuals[x.root,]
            x.c <- x$coefficients[x.root,]
            x.w <- x$weights[x.root,]
        } else {
            x.res <- x$residuals
            x.c <- x$coefficients
            x.w <- x$weights
        }

            x.n <- length(x.res)
            x.p <- length(x.c)
        
        DF.x <- x$tot.weights[x.root]*x.n - x.p 

        y.tot.sol <- y$tot.sol
        if (y.tot.sol<y.root) {
            stop(paste("'y' Root ",root," not found"))
        }

        if (y.tot.sol!=1) {
            y.res <- y$residuals[y.root,]
            y.c <- y$coefficients[y.root,]
            y.w <- y$weights[y.root,]
        } else {
            y.res <- y$residuals
            y.c <- y$coefficients
            y.w <- y$weights
        }

            y.n <- length(y.res)
            y.p <- length(y.c)

        DF.y <- y$tot.weights[y.root]*y.n - y.p 

    } else {

    if (inherits(x, "wle.normal") && inherits(y, "wle.normal")) {


       x.tot.sol <- x$tot.sol
        if (x.tot.sol<x.root) {
            stop(paste("'x' Root ",root," not found"))
        }
        if (x.tot.sol!=1) {
            x.res <- x$residuals[x.root,]
            x.w <- x$weights[x.root,]
        } else {
            x.res <- x$residuals
            x.w <- x$weights
        }
        
        DF.x <- x$tot.weights[x.root] - 1 

        y.tot.sol <- y$tot.sol
        if (y.tot.sol<y.root) {
            stop(paste("'y' Root ",root," not found"))
        }

        if (y.tot.sol!=1) {
            y.res <- c(y$residuals[y.root,])
            y.w <- y$weights[y.root,]
        } else {
            y.res <- y$residuals
            y.w <- y$weights
        }

        DF.y <- y$tot.weights[y.root] - 1
    
      } else {
          stop("'x' and 'y' must be of class wle.lm or wle.normal")
      }
      }
      
    V.x <- c(x.w%*%(x.res^2) / DF.x)
    V.y <- c(y.w%*%(y.res^2) / DF.y)

    ESTIMATE <- V.x / V.y
    STATISTIC <- ESTIMATE / ratio
    PARAMETER <- c(DF.x, DF.y)

    PVAL <- pf(STATISTIC, DF.x, DF.y)
    if (alternative == "two.sided") {
        PVAL <- 2 * min(PVAL, 1 - PVAL)
        BETA <- (1 - conf.level) / 2
        CINT <- c(ESTIMATE / qf(1 - BETA, DF.x, DF.y),
                  ESTIMATE / qf(BETA, DF.x, DF.y))
    }
    else if (alternative == "greater") {
        PVAL <- 1 - PVAL
        CINT <- c(ESTIMATE / qf(conf.level, DF.x, DF.y), Inf)
    }
    else
        CINT <- c(0, ESTIMATE / qf(1 - conf.level, DF.x, DF.y))
    names(STATISTIC) <- "WF"
    names(PARAMETER) <- c("num df", "denom df")
    names(ESTIMATE) <- names(ratio) <- "ratio of variances"
    attr(CINT, "conf.level") <- conf.level
    RVAL <- list(statistic = STATISTIC,
                 parameter = PARAMETER,
                 p.value = PVAL,
                 conf.int = CINT,
                 estimate = ESTIMATE,
                 null.value = ratio,
                 alternative = alternative,
                 method = "WF test to compare two variances",
                 data.name = DNAME)
    attr(RVAL, "class") <- "htest"
    return(RVAL)
}


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