######Nonlinear CAP FUNCTION DEFNS for multiple environmental variables####
#NCAP program for multiple environmental variables
###########################################################################
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, 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 CAUTION: Returning design matrix for 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)
}

LinCCor=function(X,pcoD,m,stat="Rsquare",verbose=T) {
  if(toupper(substr(stat,1,3))=="RDA") {
    if(dim(X)[2]>1) {
      dim(X)[2]=npar
      linear=function(c,X) return(1+X[,-1]%*%c)
      LinCor=-nlm(NLCCor,rep(1,npar-1),X=X,pcoD=pcoD,gradient=linear,m=m,stat="RDA")$min }
    else {
      m.max=sum(cumsum(abs(pcoD$val))<=sum(pcoD$val))
      wgts=pcoD$val[1:m]/sum(pcoD$val[1:m.max])
      LinCor=0
      for(i in 1:m) LinCor=LinCor+wgts[i]*summary(lm(pcoD$vec[,i]~X))$r.sq }
  }
  else 
    LinCor=cancor(X,pcoD$vec[,1:m])$cor[1]^2
  if(verbose) cat("\nLinear",stat,"value",LinCor,"\n")
  return(LinCor)
}

NLCCor=function(b,X,pcoD,gradient,m,stat="Rsquare",
                blow=NULL,bhigh=NULL,pwgt=0.001) {
 y=gradient(b,X)
 penalty=0
 if(toupper(substr(stat,1,3))=="RDA") {
     m.max=sum(cumsum(abs(pcoD$val))<=sum(pcoD$val))
     wgts=pcoD$val[1:m]/sum(pcoD$val[1:m.max]) }
   else wgts=1
 if(!is.null(blow) & !is.null(bhigh))
     penalty=pwgt*sum((pmax(0,b-bhigh)+pmin(0,b-blow))^2)
 Q=pcoD$vec[,1:m]
 S2yQ=t(y-mean(y))%*%t(t(Q)-apply(as.matrix(Q),2,mean))
 S2yy=max(sum((y-mean(y))^2),1e-12)
 return(-sum(wgts*S2yQ^2/S2yy)+penalty) }

NLCCorSeq=function(b0,X,pcoD,grad,m=NULL,stat="Rsquare",plots=T,...) {
  npar=length(b0)
 if(toupper(substr(stat,1,3))=="RDA") {
    cat("\nUsing RDA trace statistic\n"); Rsq=F 
    }
  else {
    cat("\nUsing total Rsquare statistic\n"); Rsq=T
    }
  m.max=sum(cumsum(abs(pcoD$val))<=sum(pcoD$val))
  if(is.null(m)) m=m.max
  if(m>m.max) {
    m=m.max
    cat("\nNumber of pco's reduced to",m,"to avoid negative variation \n")
  }
  fitseq=matrix(0,nrow=m,ncol=npar+1+(Rsq==T))
  if(Rsq)   
    attr(fitseq,"dimnames")=
      list(NULL,c(paste(rep("b",npar),1:npar,sep=""),"TotalRsq","Target"))
  else 
    attr(fitseq,"dimnames")=
      list(NULL,c(paste(rep("b",npar),1:npar,sep=""),"RDATrace"))
  n=nrow(pcoD$vec)
  best.m=0; best.rsq=0
  for(npco in 1:m) { 
    if(Rsq==T) {
      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,pcoD=pcoD,gradient=grad,m=npco,stat=stat,
                 iterlim=500,...)
    b=MaxNLCor$est
    fitseq[npco,1:npar]=b 
    fitseq[npco,npar+1]=-MaxNLCor$min
    if(Rsq==T)
      if(fitseq[npco,npar+1]>target.rsq) 
        {best.m=npco; best.rsq=fitseq[npco,npar+1]}
  }
  if(plots==T) {
    for(i in 1:npar)
      plot(fitseq[1:m,i],xlim=c(0,m),las=1,cex.axis=0.8,
            xlab="Number of pco's axes used",ylab=expression(hat(b)))
    if(Rsq==T) {
      plot(fitseq[1:m,npar+1],xlim=c(0,m),ylim=c(0,max(fitseq[,npar+c(1,2)])),
           las=1,cex.axis=0.8,xlab="Number of pco axes used",ylab="Total R squared")
      points(fitseq[1:m,npar+2],type="l") }
    else
      plot(fitseq[1:m,npar+1],xlim=c(0,m),ylim=c(0,max(fitseq[,npar+1])),
           las=1,cex.axis=0.8,xlab="Number of pco axes used",ylab="RDA trace statistic")
  }
  if(Rsq==T) cat("\nBest m between 1 and",m,"is indicated to be",
                 max((1:m)[fitseq[1:m,npar+1]>fitseq[1:m,npar+2]],na.rm=T),"\n")
  return(fitseq)
}

plot.NCAP=function(b,x,pcoD,gradient,m=5,stat="Rsquare",resids=T,rescale=T,
                   xpts=NULL,...) {
  Rsq=T; if(toupper(substr(stat,1,3))=="RDA") Rsq=F 
  y=gradient(b,x)
  fit=lm(y~pcoD$vec[,1:m])
  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,...)
    if(Rsq) 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,...)
    if(Rsq) points(x,yfit,...)
  } 
  if(interactive() & resids * Rsq) {
    cat("\nType  <Return>  for residual plot")
    readline() 
    plot(y,resid,xlab="Gradient",ylab="Residual",las=1)
    abline(h=0)
  }
}

BootNLCor=function(b0,X,D,gradient,m,stat="Rsquare",nsamps=100,coverage=0.95,...) {
  if(toupper(substr(stat,1,3))=="RDA")    
    cat("\nBootstrapping using RDA statistic\n") 
  else
    cat("\nBootstrapping using Rsquare statistic\n")
  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]
    bootpco=pco(bootD,varplot=F)
    bootfit=nlm(NLCCor,b0,X=bootX,pcoD=bootpco,gradient=gradient,m=m,stat=stat,
                iterlim=500,...)
    LinCor=LinCCor(bootX,pcoD,m,stat=stat,verbose=F)
    bootsumm[iter,]=c(bootfit$est,-bootfit$min,bootfit$code,LinCor) }
  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) }
  else 
    cat("\nInsufficient number of bootstraps to compute",100*coverage," % CI \n")
  attr(bootsumm,"m")=m
  return(bootsumm)
}

PermNLCor=function(b0,X,D,gradient,m,stat="Rsq",nsamps=100,...) {
  if(toupper(substr(stat,1,3))=="RDA")    
    cat("\nPermuting using RDA statistic\n") 
  else
    cat("\nPermuting using Rsquare statistic\n")
  npar=length(b0)
  pcoD=pco(D,varplot=F)
  observedfit=nlm(NLCCor,b0,X=X,pcoD=pcoD,gradient=gradient,m=m,stat=stat,
                  iterlim=500,...)
  observedF=(-observedfit$min-LinCCor(X,pcoD,m,stat=stat,verbose=F))/(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,pcoD=pcoD,gradient=gradient,m=m,stat=stat,
                iterlim=500,...) 
    LinCor=LinCCor(permX,pcoD,m,stat=stat,verbose=F)
    permsumm[iter,]=c(permfit$est,-permfit$min,permfit$code,LinCor) }
  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")=m
  return(permsumm)
}

###############################################################################
#Function defintions for obtaining principal co-ordinates
###############################################################################

#########Calculates dissimilarity matrix from abundance data##########
distance=function(N,measure="BC",trans="none"){
 n=dim(N)[1]; p=dim(N)[2]
 N=switch(trans,none=N,sqrt=sqrt(N),fourthroot=N^0.25,pa=(N>0),
          rowpropns=N/apply(N,1,sum),
  stop("\n","Unrecognized transformation",
       "\n","Must be one of: none, sqrt, fourthroot, pa, rowpropns","\n"))
 if(n<2) stop("Must be at least 2 sites")
 D=matrix(0,nrow=n,ncol=n)
 for(i in 1:(n-1)) {
  for(j in (i+1):n) {
   D[i,j]=switch(measure,
   BC=(sum(abs(N[i,]-N[j,]))/sum(N[i,]+N[j,])),
   sqrtBC=(sum(abs(N[i,]-N[j,]))/sum(N[i,]+N[j,]))^0.5,
   Can=sum(ifelse(N[i,]+N[j,]>0,abs(N[i,]-N[j,])/(N[i,]+N[j,]),0)),
   sqrtCan=sum(ifelse(N[i,]+N[j,]>0,abs(N[i,]-N[j,])/(N[i,]+N[j,]),0))^0.5,
   HornM=1-2*sum(N[i,]*N[j,])/( ( sum(N[i,]^2)/(sum(N[i,]))^2 + sum(N[j,]^2)/
               (sum(N[j,]))^2)*sum(N[i,])*sum(N[j,]) ),
   Eucl=sqrt(sum((N[i,]-N[j,])^2)),
   stop("\n","Unrecognized distance measure","\n",
      "Must be one of ","BC, sqrtBC, Can, sqrtCan, HornM, or Eucl","\n") )}}
 D=D+t(D)
 return(D)
}


################Returns PCO's and eigenvalues#############################
pco=function(D,varplot=T) {
 A=-0.5*D^2 
 G=centre.matrix(A) 
 eigenG=eigen(G,symmetric=T)
 if(varplot==T) {
 evals=eigenG$val
 propnvar=cumsum(abs(evals))/sum(evals)
 propnvar=propnvar[propnvar<=1]
 m=length(propnvar)
 cat("Total variation is explained by the first",m,"eigenvalues\n")
 plot(0:m,c(0,propnvar),type="l",las=1,
     ylim=c(0,1),xlab="# pco's",ylab="Propn of variation explained")
 abline(h=1,v=seq(0,m,5),lty=3) }
 return(eigenG)
}

###########Biplots of PCO's and PCO's vs environmental variable###########
plot.pco=function(pcoD,npco=3,x=NULL,ask=FALSE,...) {
  Q=pcoD$vec
  if(interactive() & is.null(x)) {
    for(i in 1:(npco-1)) {
      for(j in (i+1):npco) {
        if((i+j)>3 & ask) {
          cat("\nType  <Return>  for plot of pco",i,"vs pco",j)
          readline()
        } 
        plot(Q[,i],Q[,j],xlab=paste("Q[,",i,"]"),ylab=paste("Q[,",j,"]"),...) 
      }
    }
  } 
  if (interactive() & !is.null(x)) { 
    for(j in 1:npco) {
      if(ask) {
        cat("\nType  <Return>  for plot of enviromental variable vs pco",j)
        readline()
      }
      plot(x,Q[,j],xlab=paste(deparse(substitute(x))),
           ylab=paste("Q[,",j,"]"),...) 
    }
  }
  invisible()
}

#################Centres a matrix#########################################
centre.matrix=function(A) {n=dim(A)[1]
 return(A-apply(A,1,mean)%o%rep(1,n)-rep(1,n)%o%apply(A,2,mean)+mean(A)) }

##########Plots of species abundance data vs environmental variable#######
lattice.plot=function(x,N,label=NULL,mfrow=NULL,...) {
  old.par=par(no.readonly=TRUE) # all par settings which could be changed.
  on.exit(par(old.par))
  ratio=7/8
  n=nrow(N); nvars=ncol(N)
  if(is.null(mfrow)) 
     mfrow=c(ceiling(sqrt(nvars)/ratio),ceiling(sqrt(nvars)*ratio))
  par(mfrow=mfrow,mar=c(0,0,0,0),oma=c(0,0,2,0))
  for(i in 1:nvars) plot(x,N[,i],xaxt="n",yaxt="n",pch=20)
  if(is.null(label)) label=paste(deparse(substitute(x)))
  mtext(label,outer=T,line=0.5)
  invisible()
}