####################NCAP FUNCTION DEFINITIONS##############################
##########Nonlinear cannonical analysis of principle co-ordinates##########
###########################################################################
gradient.choice=function(type="vonB") { 
 gradientfunction=switch(type,
   vonB={function(b,x) return( pmax(0,(1-exp(-x%*%b))) )},
   hyperbolic={function(b,x) return( pmax(0,x%*%b/(1+x%*%b)) )},
   logistic={function(b,x) return( exp(x%*%b)/(1+exp(x%*%b)) )},
   stop("\n","Unrecognized gradient",
       "\n","Must be one of: vonB, hyperbolic, or logistic","\n") )
 return(gradientfunction)
}

model=function(formula.spec,fixed.intercept=T) {
  if(fixed.intercept)
    X=as.matrix(model.matrix(formula.spec)[,-1])
  else {
    cat("\n NOTE: Model includes intercept term in linear predictor.",
        "\n       This is NOT VALID for the vonB gradient","\n")
    X=as.matrix(model.matrix(formula.spec))
  }
  attr(X,"fixed.intercept")=fixed.intercept
  return(X)
}

NLCCor=function(b,x,Q,gradient,blow=NULL,bhigh=NULL,pwgt=0.001) {
 if(length(b)!=dim(x)[2])
   stop("\nERROR: Length of parameter vector does not equal model dimension")
 y=gradient(b,x)
 #y=y+0.001*(-1)^(1:length(y)) #Jiggle to avoid numerical problems
 penalty=0
 if(!is.null(blow) & !is.null(bhigh))
     penalty=pwgt*sum((pmax(0,b-bhigh)+pmin(0,b-blow))^2)
 fit1=lm(y~Q)
 return(-summary(fit1)$r.squared+penalty) }

NLCCorSeq=function(b0,x,Q,grad,...) {
  npar=length(b0)
  fitseq=matrix(0,nrow=ncol(Q),ncol=npar+2)
  parnames=paste("par",1:npar)
  colnames(fitseq)=c(parnames,"Max cor","Target cor")
  n=nrow(Q)
  best.m=0; best.rsq=0
  for(npco in 1:ncol(Q)) {
    target.rsq=1-exp(log(1-best.rsq)-(npco-best.m)*max(2,log(n))/n) 
    fitseq[npco,npar+2]=target.rsq
    MaxNLCor=nlm(NLCCor,b0,x=x,Q=Q[,1:npco],gradient=grad,iterlim=500,...)
    b=MaxNLCor$est
    fitseq[npco,1:npar]=b 
    fitseq[npco,npar+1]=-MaxNLCor$min
    if(fitseq[npco,npar+1]>target.rsq) {best.m=npco; best.rsq=fitseq[npco,npar+1]}
  }
  return(fitseq)
}

#Example of use
#NLCCorSeq(0.01,x=X,Q,gradient,blow=0.0,bhigh=1)

plot.NCAP=function(b,x,Q,gradient,resids=T,rescale=T,xpts=NULL,...) {
  y=gradient(b,x)
  fit=lm(y~Q)
  yfit=fit$fit; resid=fit$resid
  if(rescale) {
    gfit=lm(yfit~y)
    resid=gfit$resid/gfit$coef[2]
    yfit=y+resid #Reverse role of y and yfit
  }
  if (interactive() & ncol(x)==1) {
    if(is.null(xpts)) xpts=as.matrix(seq(min(x),max(x),length=100))
    plot(xpts,gradient(b,as.matrix(xpts)),type="l",ylab="Gradient",las=1,...)
    points(x,yfit) 
  }
  if (ncol(x)==2 & !attr(x,"fixed.intercept")) {
    x=as.matrix(x[,2])
    if(!is.null(xpts)) xpts=cbind(1,xpts)
      else xpts=cbind(1,as.matrix(seq(min(x),max(x),length=100)))   
    plot(xpts[,2],gradient(b,xpts),type="l",ylab="Gradient",las=1,...)
    points(x,yfit,ylab="Gradient",...)
  } 
  if(interactive() & resids) {
    cat("\nType  <Return>  for residual plot")
    readline() 
    plot(y,resid,xlab="Gradient",ylab="Residual",las=1)
    abline(h=0)
  }
}

BootNLCor=function(b0,x,D,npco,gradient,nsamps=100,coverage=0.95,...) {
  npar=length(b0)
  bootsumm=matrix(0,nrow=nsamps,ncol=npar+3)
  for(iter in 1:nsamps) {
    resampled=sample(1:nrow(x),replace=T)
    if(round(iter/1)==iter/1) 
      dput(list(iter=iter,resampled=resampled),file="BootProgress.txt")
    bootx=as.matrix(x[resampled,]); bootD=D[resampled,resampled]
    bootQ=pco(bootD,varplot=F)$vec[,1:npco]
    bootfit=nlm(NLCCor,b0,x=bootx,Q=bootQ,gradient=gradient,iterlim=500,...)
    LinCor=cancor(bootx,bootQ)$cor^2
    bootsumm[iter,]=c(bootfit$est,-bootfit$min,bootfit$code,LinCor[1]) }
  codecount=sum(bootsumm[,ncol(bootsumm)]>2)
  if(codecount>0) cat("\nWARNING:",codecount,"fits may not have converged \n")
  CIsamps=round(nsamps*(c((1-coverage)/2,1-(1-coverage)/2)))
  if(CIsamps[1]>0){ 
    cat("\n",100*coverage,"% bootstrap percentile CI",ifelse(npar<2,"is","are"),"\n")
    write(signif(apply(as.matrix(bootsumm[,1:npar]),2,sort)[CIsamps,],3),"",ncol=2) }
  attr(bootsumm,"m")=npco
  return(bootsumm)
}

PermNLCor=function(b0,x,D,npco,gradient,nsamps=100,...) {
  npar=length(b0)
  Q=pco(D,varplot=F)$vec[,1:npco]
  observedfit=nlm(NLCCor,b0,x=x,Q=Q,gradient=gradient,iterlim=500,...)
  observedF=(-observedfit$min-cancor(x,Q)$cor^2)/(1+observedfit$min)
  cat("\nObserved NCAP correlation and F-stat are",
      -observedfit$min,"and",observedF,"\n") 
  permsumm=matrix(0,nrow=nsamps,ncol=npar+3)
  for(iter in 1:nsamps) { 
    resampled=sample(1:nrow(x),replace=F)
    if(round(iter/1)==iter/1) 
      dput(list(iter=iter,resampled=resampled),file="PermProgress.txt")
    permx=as.matrix(x[resampled,])
    permfit=nlm(NLCCor,b0,x=permx,Q=Q,gradient=gradient,iterlim=500,...) 
    LinCor=cancor(permx,Q)$cor^2
    permsumm[iter,]=c(permfit$est,-permfit$min,permfit$code,LinCor[1]) }
  cor.pval=1-(rank(c(-observedfit$min,permsumm[,npar+1]))[1]-1)/(nsamps+1)
  F.rank=rank(c(observedF,(permsumm[,npar+1]-permsumm[,npar+3])/(1-permsumm[,npar+1])))
  F.pval=1-(F.rank[1]-1)/(nsamps+1)
  cat("\nNull model permutational p-value is",cor.pval,
      "\nLinear model hypothesis p-value is",F.pval,"\n")
  attr(permsumm,"m")=npco
  return(permsumm)
}