##########################################################################
#
# amcmc:  A PACKAGE OF R FUNCTIONS FOR ADAPTIVE MCMC
#
# Copyright (c) 2007 by Jeffrey S. Rosenthal (probability.ca/jeff)
# 
# (Last modified March 2009.)
# 
# Instructions for use at:  probability.ca/amcmc
# 
# Licensed for general copying, distribution and modification according to
# the GNU General Public License (www.gnu.org/copyleft/gpl.html).
# 
##########################################################################


# DEFINE SOME PRELIMINARY FUNCTIONS.

indicator <- function(cond) {
  if (cond)
    return(1)
  else
    return(0)
}

argcount <- function(f) {
  thecount = 0
  for (iii in 1:length(formals(f))) {
    thecount = thecount + indicator( formals(f)[iii] == '' )
  }
  return(thecount)
}

firstcoordsq <- function(x) return(x^2)

# DEFINE THE MAIN FUNCTION.

amcmc <- function(densfn=dnorm, functional=firstcoordsq, numbatches=1000,
	   batchlength=10, logflag=FALSE, adapt=TRUE, burninfrac = 0.1,
	   initval = 1.0, verboselevel=0, cdensity=cfns, cfunctional=cfns,
	   cfns=FALSE, vectorlength=0, write=FALSE) {
  if ( argcount(functional) != argcount(densfn) )
    stop("Number of arguments must be the same for densfn and functional.")
  if (cdensity)
    K = -1
  else
    K = argcount(densfn)
  adaptval = indicator(adapt)
  logval = indicator(logflag)
  cdensityval = indicator(cdensity)
  cfunctionalval = indicator(cfunctional)
  writeval = indicator(write)
  if (vectorlength > 0) {
    listdensfn <- densfn
    listfunctional <- functional
  } else {
    listdensfn <- function(thelist) return( do.call(densfn,thelist) )
    listfunctional <- function(thelist) return( do.call(functional,thelist) )
  }
  .Call("amcmc", listdensfn, listfunctional, K, numbatches,
	batchlength, logval, adaptval, burninfrac, initval, verboselevel,
	cdensityval, cfunctionalval, vectorlength, writeval, parent.frame())
}

# DEFINE SOME EXTRA FUNCTIONS (NOT REALLY NEEDED).

firstcoord <- function(x) return(x)
identity <- function(x) return(x)
sq <- function(x) return(x^2)
lognormaldens <- function(x) return( log(dnorm(x)) )

lognormaldens2 <- function(x,y) return(
			log(dnorm(x)) + log(dnorm(y)) )
normaldens2 <- function(x,y) return( dnorm(x) * dnorm(y) )
firstcoordsq2 <- function(x,y) return(x^2)
ssq2 <- function(x,y) return(x^2+y^2)

listlognormaldens2 <- function(arglist) return(
			log(dnorm(arglist[1])+dnorm(arglist[2])) )
listfirstcoordsq2 <- function(arglist) return(arglist[1]^2)
listssq2 <- function(arglist) return(arglist[1]^2+arglist[2]^2)

lognormaldens5 <- function(x,y,z,v,w) return( log(
	dnorm(x))+log(dnorm(y))+log(dnorm(z))+log(dnorm(w))+log(dnorm(v)) )
normaldens5 <- function(x,y,z,v,w) return(
	dnorm(x)*dnorm(y)*dnorm(z)*dnorm(w)*dnorm(v) )
firstcoordsq5 <- function(x,y,z,v,w) return(x^2)
ssq5 <- function(x,y,z,v,w) return(x^2+y^2+z^2+v^2+w^2)

amcmc2 <- function(...) return( amcmc(normaldens2, firstcoordsq2, ...) )
amcmc5 <- function(...) return( amcmc(normaldens5, firstcoordsq5, ...) )

# FUNCTIONS SPECIFIC TO THE BASEBALL DATA OF MORRIS (1983).

baseballdata <- c(0.395, 0.375, 0.355, 0.334, 0.313, 0.313, 0.291, 0.269,
    0.247, 0.247, 0.224, 0.224, 0.224, 0.224, 0.224, 0.200, 0.175, 0.148);

baseballlogdens <- function(th1, th2, th3, th4, th5, th6, th7, th8, th9,
	th10, th11, th12, th13, th14, th15, th16, th17, th18, mu, A) {
    if (A<0)
      return(-Inf);
    s0 = 1; mu0=0; a1=a2=-1; b1=b2=2;
    V = 0.00434;  lV = log(V);
    theta <- c(th1, th2, th3, th4, th5, th6, th7, th8, th9, th10, th11,
				th12, th13, th14, th15, th16, th17, th18)
    result = -(mu-mu0)*(mu-mu0) / 2.0 / s0 / s0 -
                    b1/A - (a1+1) * log(A) -b2/V - (a2+1) * log(V);
    for (kk in 1:18) {
        result = result - 0.5 * log(A) - (theta[kk]-mu)^2 / 2.0 / A -
		0.5 * lV - (baseballdata[kk]-theta[kk])^2 / 2.0 / V;
    }
    return(result);
}

baseballfirstfunct <- function(th1, th2, th3, th4, th5, th6, th7, th8, th9,
	th10, th11, th12, th13, th14, th15, th16, th17, th18, mu, A) {
    return(th1);
}

baseballlogdensvect <- function(theargs) {
    theta <- theargs[1:18];
    mu <- theargs[19];
    A <- theargs[20];
    if (A<0)
      return(-Inf);
    s0 = 1; mu0=0; a1=a2=-1; b1=b2=2;
    V = 0.00434;  lV = log(V);
    result = -(mu-mu0)*(mu-mu0) / 2.0 / s0 / s0 -
                    b1/A - (a1+1) * log(A) -b2/V - (a2+1) * log(V);
    for (kk in 1:18) {
        result = result - 0.5 * log(A) - (theta[kk]-mu)^2 / 2.0 / A -
		0.5 * lV - (baseballdata[kk]-theta[kk])^2 / 2.0 / V;
    }
    return(result);
}

baseballfirstfunctvect <- function(theargs) {
    return(theargs[1]);
}

# ROUTINES TO PRECOMPILE THE USER'S OWN FUNCTIONS

currentcfns = NULL;

cfns <- function(stringname="dummyhold") {
    system("rm -f amcmc.so")
      system(paste("R CMD SHLIB amcmc.c ", stringname, ".c", sep=""))
    dyn.load("amcmc.so")
    currentcfns <<- stringname;
}

editcfns <- function(stringname=currentcfns) {
  if (is.null(stringname)) {
    cat('editcfns: must specify base filename\n');
  } else {
    thefilename = paste(stringname, ".c", sep="")
    thecommand = paste(' sh -c " if test -s ', thefilename,
		' ; then echo pre-existing ; else cp dummyhold.c ',
		thefilename, ' ; fi " ')
    system(thecommand)
    thecommand = paste(getOption('editor'), thefilename, sep=' ')
    system(thecommand)
    cfns(stringname)
  }
}

showcfns <- function() { currentcfns }

# CREATE A DUMMYHOLD FUNCTION

system("rm -f dummyhold.c")
system("echo '#include <math.h>' > dummyhold.c")
system("echo 'int mydim() { return(1); }' > dummyhold.c")
system("echo 'double mydensity() { return(1.0); }' >> dummyhold.c")
system("echo 'double myfunctional() { return(1.0); }' >> dummyhold.c")

# COMPILE THE C CODE IF NECESSARY ...

system(' sh -c " if test -s amcmc.so; then echo pre-compiled ; else R CMD SHLIB amcmc.c dummyhold.c ; fi " ')

# LOAD THE ASSOCIATED C SHARED OBJECT (.so) FILE.

dyn.load("amcmc.so")