#Integration over [0,1]^p

#Quasi-Zufallszahlen (Halton)

f.halton <- function(n,p)
{
  ## Purpose: Generating the Halton sequence with base  p 
  ## -------------------------------------------------------------------------
  ## Arguments: n = length of the sequence, p = base (must be prime)
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 29 Oct 2001, 10:53
nn <- 1:n
k <- ceiling(log(n)/log(p))
a <- matrix(0,nrow=n,ncol=k)
for (i in (1:k)){
  a[,i] <- nn%%p
  nn <- nn%/%p
}
a%*%(p^(-(1:k)))
}

f.halton(20,2)
plot(c(0,f.halton(50,2)),type="b")
f.halton(20,3)
f.halton(20,5)
plot(f.halton(100,7),type="b")

#comparing grid points, quasi- und pseudo-random points in 10 dimensions
xq <- matrix(0,nrow=1024,ncol=10) #quasi-random points
xq[,1] <- ((1:1024)-0.5)/1024
pp <- c(2,3,5,7,11,13,17,19,23)
for (i in (2:10)) xq[,i] <- f.halton(1024,pp[i-1])
round(xq[1:20,],3)

xp <- matrix(runif(10240),nrow=1024,ncol=10) # pseudo-random points


n <- 1:1024 # generating a regular grid 
xg <- matrix(0,nrow=1024,ncol=10)
for (i in (1:10)){
  xg[,i] <- n%%2
  n <- n%/%2
}
xg <- 0.25+0.5*xg
xg[1:20,]

# Visualisation of some coordinate pairs
pairs(xq[,1:5],cex=0.5)
pairs(xq[,6:10],cex=0.5)
pairs(xp[,1:5],cex=0.5)
pairs(xp[,6:10],cex=0.5)



# Projection on two random orthogonal directions
par(mfcol=c(3,2))
u1 <- rnorm(10)
u2 <- rnorm(10)
u1 <- u1/sqrt(sum(u1^2))
u2 <- u2-sum(u1*u2)*u1
u2 <- u2/sqrt(sum(u2^2))
plot(xg%*%u1,xg%*%u2,cex=0.5)
plot(xp%*%u1,xp%*%u2,cex=0.5)
plot(xq%*%u1,xq%*%u2,cex=0.5)



# Integration of the function exp(-sum_i c*|x_i-0.5|)
f.int <- function(x,c)
{
  ## Purpose: Approximation of the integral of exp(-sum_i c*|x_i-0.5|)
  ##          over x aus [0,1]^10
  ## -------------------------------------------------------------------------
  ## Arguments: x=integration points, c=Parameter of the integrand
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  5 Nov 2001, 10:10
  h <- apply(abs(x-0.5),1,sum)
  wert <- (2*(1-exp(-0.5*c))/c)^10
  mean(exp(-c*h))/wert -1 #relative error
}

f.int(xp,0.25)
f.int(xp,0.5)
f.int(xp,1)
f.int(xp,2)

f.int(xq,0.25)
f.int(xq,0.5)
f.int(xq,1)
f.int(xq,2)

f.int(xg,0.25)
f.int(xg,0.5)
f.int(xg,1)
f.int(xg,2)

#How well do the points fill space ?

#distances to nearest neighbor and to boundary.
par(mfrow=c(2,2))
dnn <- as.matrix(dist(xp)) #all pairwise distances
dnn <- dnn + diag(sqrt(10),1024)
dnn <- apply(dnn,1,min)
mean(dnn)
hist(dnn,breaks=seq(0.10,0.90,0.05)) 
dbd <- abs(apply(xp,1,range)-matrix(c(0,1),nrow=2,ncol=1024))
dbd <- apply(dbd,2,min) #distances to the boundary
mean(dbd)
hist(dbd)
dnn <- as.matrix(dist(xq)) #all pairwise distances
dnn <- dnn + diag(sqrt(10),1024)
dnn <- apply(dnn,1,min)
mean(dnn)
hist(dnn,breaks=seq(0.10,0.90,0.05))
dbd <- abs(apply(xq,1,range)-matrix(c(0,1),nrow=2,ncol=1024))
dbd <- apply(dbd,2,min) #distances to the boundary
mean(dbd)
hist(dbd)




# distance of an additional random point to the nearest member of the sequence
f.abst <- function(x,nrep=100)
{
  ## Purpose: Distance of a random point in [0,1]^p to the nearest member of x
  ## -------------------------------------------------------------------------
  ## Arguments: x=sequence in [0,1]^p, nrep=number of replicates.
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  5 Nov 2001, 09:43
  n <- nrow(x)
  p <- ncol(x)
  dnn <- rep(0,nrep)
  for (i in (1:nrep)) {
    d <- apply((x - matrix(runif(p),nrow=n,ncol=p,byrow=T))^2,1,sum)
    dnn[i] <- sqrt(min(d))
  }
  dnn
}

par(mfrow=c(2,2))
abst <- f.abst(xp)
mean(abst)
hist(abst,breaks=seq(0.20,0.80,0.05))
abst <- f.abst(xq)
mean(abst)
hist(abst,breaks=seq(0.20,0.80,0.05))
abst <- f.abst(xg)
mean(abst)
hist(abst,breaks=seq(0.20,0.80,0.05))