# To run procedure "monofit" in R, one needs libraries "mgcv" and "splines"

if(exists("is.R") && is.function(is.R) && is.R()){
library(mgcv)
library(splines)
}

monofit<-function(x,y,Z1=NULL,Z2=NULL,gridx=seq(min(x),max(x),len=101),
gridZ1=NULL,gridZ2=NULL,x0=NULL,Z0=NULL,maxknots=10,Nknots=21,degpol=3,
crit=0,levoff=0,incr=0,prin=1,alpha=0.05){

#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------

# x,y: data

# Z1: positive additional parameters included in the full fit and the monotone fit
# Z2: additional parameters included in the full fit

# gridx,gridZ1,gridZ2: grids related to x, Z1 and Z2 for the plots

# x0,Z0: points related to x and Z1 where the regression has to be calculated 

# Note: Z1, Z2, gridZ1, gridZ2 and Z0 should be matrices

# maxknots: maximum number of knots in the model
# Nknots  : number of candidate knots
# degpol  : degree of splines 
 
# crit=0: F-test approach (default)
# crit=1: aic criterion
# crit=2: bic criterion
# crit=3: gcv criterion

# levoff=0: without levelling-off (default)
# levoff=1: with levelling-off

# incr=0: monotone decreasing (default)
# incr=1: monotone increasing

# prin=0: without printing and ploting of the different steps
# prin=1: with printing and ploting of the different steps 
# prin=2: with printing and ploting of the different steps 
#         where the user has to click on the plots to pursue

# alpha: significance level used for the F-approach 

#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------

# Initialization of the knots

if((maxknots+1)>Nknots) Nknots<-maxknots+1
allknots<-seq(min(x),max(x),len=Nknots+2)[-c(1,Nknots+2)]

# Forward algorithm

if(prin==0) gridx<-NULL
if(prin>=1) cat("Begin forward procedure","\n")
finished<-0
knots<-NULL
opt0<-Inf

while(finished==0){

ssr0<-Inf
for(k in 1:length(allknots)){
ssr<-obj(c(knots,allknots[k]),x,y,Z1,incr,levoff,degpol)
if(ssr<ssr0){
ssr0<-ssr
k0<-k
}
}
if(prin>=1) cat("next candidate knot to be added:",round(allknots[k0],2),"\n")
nextknot<-allknots[k0]
allknots<-allknots[-k0]
if(crit==0) tmp<-testfit(x,y,Z1,Z2,gridx,gridZ1,gridZ2,x0,Z0,knots,nextknot,degpol,crit,levoff,incr,prin)
if(crit>=1) tmp<-testfit(x,y,Z1,Z2,gridx,gridZ1,gridZ2,x0,Z0,knots,NULL,degpol,crit,levoff,incr,prin)
if(crit==0 && (tmp$opt>=alpha || length(knots)==maxknots)){
finished<-1
finalknots<-knots
}
if(crit>=1 && tmp$opt<opt0){
opt0<-tmp$opt
finalknots<-knots
}
if(crit>=1 && length(knots)==maxknots) finished<-1
if(finished==0) knots<-c(knots,nextknot)

} # end while

if(length(finalknots)>0) finalknots<-sort(finalknots)
if(prin>=1 && length(finalknots)>0)
cat("knots selected after forward procedure:",round(finalknots,2),"\n")
if(prin>=1 && length(finalknots)==0) 
cat("no knots selected after forward procedure","\n")

# Backward algorithm

if(prin>=1 && length(finalknots)>0) cat("Begin backward procedure","\n")
finished<-0
knots<-finalknots
if(length(knots)==0) finished<-1

while(finished==0){

ssr0<-Inf
for(k in 1:length(knots)){
ssr<-obj(knots[-k],x,y,Z1,incr,levoff,degpol)
if(ssr<ssr0){
ssr0<-ssr
k0<-k
}
}
if(prin>=1) cat("next candidate knot to be removed:",round(knots[k0],2),"\n")
if(crit==0) tmp<-testfit(x,y,Z1,Z2,gridx,gridZ1,gridZ2,x0,Z0,knots[-k0],c(nextknot,knots[k0]),degpol,crit,levoff,incr,prin)
if(crit>=1) tmp<-testfit(x,y,Z1,Z2,gridx,gridZ1,gridZ2,x0,Z0,knots[-k0],NULL,degpol,crit,levoff,incr,prin)
if(crit==0 && tmp$opt>=alpha) finalknots<-knots[-k0]
if(crit==0 && tmp$opt<alpha) finished<-1
if(crit>=1 && tmp$opt<opt0){
opt0<-tmp$opt
finalknots<-knots[-k0]
}
if(finished==0){
nextknot<-c(nextknot,knots[k0])
knots<-knots[-k0]
}
if(length(knots)==0) finished<-1

} # end while

if(length(finalknots)>0) finalknots<-sort(finalknots)
if(prin>=1 && length(finalknots)>0)
cat("knots selected after backward procedure:",round(finalknots,2),"\n")
if(prin>=1 && length(finalknots)==0) 
cat("no knots selected after backward procedure","\n")

# Final fit

tmp<-testfit(x,y,Z1,Z2,gridx,gridZ1,gridZ2,x0,Z0,finalknots,NULL,degpol,crit,levoff,incr,prin)
yhat<-as.vector(tmp$yhat)
knots<-finalknots
if(!is.null(x0) && !is.null(Z0)) return(list(x0=x0,Z0=Z0,yhat=yhat,knots=knots))
if(!is.null(x0) && is.null(Z0)) return(list(x0=x0,yhat=yhat,knots=knots))
if(is.null(x0)) return(knots)
}

#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------

testfit<-function(x,y,Z1,Z2,gridx,gridZ1,gridZ2,x0,Z0,knots,nextknot,degpol,crit,levoff,incr,prin){

n<-length(x)
r<-0
s<-0
if(!is.null(Z1)) r<-dim(as.matrix(Z1))[2]
if(!is.null(Z2)) s<-dim(as.matrix(Z2))[2]

# constrained fit

tmp<-objmon(knots,x,y,Z1,gridx,gridZ1,x0,Z0,incr,levoff,degpol)
ssrc<-tmp$ssr
fitc<-tmp$fit
yhat<-tmp$yhat
q<-sum(tmp$betaint<=0.0001)
p<-length(knots)+degpol+1
# unconstrained fit

tmp<-objfull(c(knots,nextknot),x,y,Z1,Z2,gridx,gridZ1,gridZ2,degpol)
ssru<-tmp$ssr
fitu<-tmp$fit

# Goodness of constrained fit 

if(crit==0){
k<-length(nextknot)
opt<-1-pf((n-p-k-r-s)*(ssrc-ssru)/((q+k+s)*ssru),q+k+s,n-p-k-r-s)
}
if(crit==1) opt<-n*log(ssrc/n)+2*(p+r-q)
if(crit==2) opt<-n*log(ssrc/n)+log(n)*(p+r-q)
if(crit==3) opt<-log(ssrc/n)-2*log(1-(p+r-q)/n)

# plot of the fits

if(prin>=1 && !is.null(gridx)){
labtit<-paste(length(knots),"knots","  SSR=",round(ssrc,2))
if(crit==0 && length(nextknot)>0) labtit<-paste(labtit,"  p=",round(opt,3)) 
if(crit==1) labtit<-paste(labtit,"  AIC=",round(opt,2))
if(crit==2) labtit<-paste(labtit,"  BIC=",round(opt,2))
if(crit==3) labtit<-paste(labtit,"  GCV=",round(opt,2))
plot(x,y,type="n")
points(x,y)
i0<-1
i1<-1
while(i1<length(gridx)){
while(gridx[i1+1]>gridx[i1] && i1<length(gridx)) i1<-i1+1
if(length(nextknot)>0) lines(gridx[i0:i1],fitu[i0:i1],lty=2,lwd=2)
lines(gridx[i0:i1],fitc[i0:i1],lty=1,lwd=2)
i0<-i1+1
i1<-i1+1
}
if(length(knots)>0) for(i in 1:length(knots)) abline(v=knots[i],lty=2,lwd=2)
if(length(nextknot)>0) 
for(i in 1:length(nextknot)) abline(v=nextknot[i],lty=2,lwd=1)
title(labtit)
if(prin==2) locator(1)
}
return(list(opt=opt,yhat=yhat))
}

#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------

obj<-function(knots,x,y,Z1,incr,levoff,degpol){
B<-Bsplinebasis(knots,x,degpol,min(x),max(x))
Bint<-spline.int(B)
Bdim<-dim(Bint)[2]
if(incr==1) y<-(-y)
ymin<-min(min(y),0)
y<-y-ymin
if(levoff==0) X<-cbind(Bint,Z1)
if(levoff==1){
indlev<-(Bdim-degpol):(Bdim-1)
X<-cbind(Bint[,-indlev],Z1)
}
if(!exists("is.R")) res<-nnls.fit(X,y)$residuals
if(exists("is.R") && is.function(is.R) && is.R()){
d<-dim(X)[2]
M<-list(y=y,w=y*0+1,X=X,C=matrix(0,0,0),S=matrix(0,0,0),off=array(0,0),df=array(0,0),sp=array(0,0),p=rep(1,d),Ain=diag(1,d,d),bin=rep(0,d))
betaint<-pcls(M)
res<-y-X%*%betaint
}
return(sum(res^2))
}

objmon<-function(knots,x,y,Z1,gridx,gridZ1,x0,Z0,incr,levoff,degpol){

B<-Bsplinebasis(knots,x,degpol,min(x),max(x))
if(!is.null(gridx)) Bgrid<-Bsplinebasis(knots,gridx,degpol,min(x),max(x))
if(!is.null(x0)) Bx0<-Bsplinebasis(knots,x0,degpol,min(x),max(x))
Bint<-spline.int(B)
Bdim<-dim(Bint)[2]
if(incr==1) y<-(-y)
ymin<-min(min(y),0)
y<-y-ymin
if(levoff==0) X<-cbind(Bint,Z1)
if(levoff==1){
indlev<-(Bdim-degpol):(Bdim-1)
X<-cbind(Bint[,-indlev],Z1)
}
if(!exists("is.R")){
tmp<-nnls.fit(X,y)
betaint<-tmp$coef
res<-tmp$residuals
}
if(exists("is.R") && is.function(is.R) && is.R()){
d<-dim(X)[2]
M<-list(y=y,w=y*0+1,X=X,C=matrix(0,0,0),S=matrix(0,0,0),off=array(0,0),df=array(0,0),sp=array(0,0),p=rep(1,d),Ain=diag(1,d,d),bin=rep(0,d))
betaint<-pcls(M)
res<-y-X%*%betaint
}
ssr<-sum(res^2)
if(length(knots)>0 && levoff==1) betaint<-c(betaint[1:(Bdim-degpol-1)],rep(0,degpol),betaint[(Bdim-degpol):length(betaint)])
if(length(knots)==0 && levoff==1) betaint<-c(rep(0,degpol),betaint[(Bdim-degpol):length(betaint)])
beta<-coef.spline(betaint[1:Bdim])
if(length(betaint)>Bdim) beta<-c(beta,betaint[(Bdim+1):length(betaint)])

fit<-NULL
yhat<-NULL
if(!is.null(gridx)) fit<-cbind(Bgrid,gridZ1)%*%beta
if(!is.null(x0)) yhat<-cbind(Bx0,Z0)%*%beta
fit<-fit+ymin
yhat<-yhat+ymin
if(incr==1){
fit<--fit
yhat<--yhat
}

return(list(ssr=ssr,fit=fit,yhat=yhat,betaint=betaint))
}

#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------

objfull<-function(knots,x,y,Z1,Z2,gridx,gridZ1,gridZ2,degpol){
B<-Bsplinebasis(knots,x,degpol,min(x),max(x))
if(!is.null(gridx)) Bgrid<-Bsplinebasis(knots,gridx,degpol,min(x),max(x))
ls<-lsfit(cbind(B,Z1,Z2),y,intercept=F)
ssr<-sum(ls$residuals^2)
beta<-ls$coef
fit<-NULL
if(!is.null(gridx)) fit<-cbind(Bgrid,gridZ1,gridZ2)%*%beta
return(list(ssr=ssr,fit=fit))
}

#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------

Bsplinebasis<-function(knots,x,degpol,a,b){
bdryknots<-c(rep(a,degpol+1),rep(b,degpol+1))
B<-spline.des(sort(c(bdryknots,knots)),x,degpol+1)$design
return(B)
}

spline.int<-function(B){
Bint<-B
p<-dim(B)[2]
for(i in 2:p)
for(j in 1:(i-1))
Bint[,i]<-Bint[,i]+B[,j]
return(Bint)
}

coef.spline<-function(betaint){
p<-length(betaint)
beta<-rep(betaint[p],p)
for(i in 1:(p-1))
for(j in i:(p-1))
beta[i]<-beta[i]+betaint[j]
return(beta)
}