#### >>>>> #### (Fourth Check, confirming uniformity on a box)

## Find a box around rep(1,7) contained in the region. If it is not
   too small, then:

### Divide this box into 2^7 = 128 pieces by conditions that each
    coord is <1 or > 1, using "sign" function

### Now for each of a very large number of MCMC generated vectors
####   (1) check that all coord's within .074 of 1, discard if not.
####   (2) for those within range, create vector of signs,then takes
###        (s+1)/2 in each coord and inner product with 2^(0:6) and
###         add 1. This gives numerical category 1:128 which should 
###         all be equally likely.

#### This calculation will do two things:
####      (i)  check uniformity within the box
####      (ii) if we compare the total number in the box compared to
####           the total number N simulated, we get an estimate of 
####           the volume of the region, whose reciprocal is the
####           normalizing constant for the density !

Counts = numeric(128)
a0 = 2^(0:6)
## assume the box half-lengths (not nec equal) are given in vector "box"

## We start with the last MCMC iterate X1:  then

for(i in 1:1e5) {
        X1= Xtrans(X1,carr,lowbd)$Xout
        if(max(abs(X1-1)/box)<1) {
           m = sum((sign(X1-1)+1)*a0)/2 + 1
           Counts[m] = Counts[m]+1 } } 

### more details to follow

===========================================
xbd = c(carr%*%rep(1,7))-lowbd
boxmat = array(0, c(1000,8))
for (i in 1:1000) {                 #### choose direction at random
    box = runif(7)            #### find largest box in region centered about
    box = box/sum(box)           #### rep(1,7) with these relative dimensions
    t = min(xbd/abs(carr %*% box))
    boxmat[i,1:7] = t*box
    boxmat[i,8] = 128*prod(t*box) }
> summary(boxmat[,8])       ### incredibly large variation in vlumes of boxes !!
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
1.084e-11 3.482e-08 3.723e-07 2.251e-05 3.690e-06 3.516e-03 
> (1:1000)[boxmat[,8]>3e-3]
[1] 438 983
> round(boxmat[438,],6)
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
[1,] 0.323946 0.090960 0.183911 0.589822 0.218544 0.142199 0.276536 0.003516
              #### half-widths in coordinate directions of largest
            ### box found centered at rep(1,7) lying within region

X1 = rep(1,7)
for(i in 1:4e4) X1 = Xtrans(X1,carr,lowbd)$Xout    #### 40,000 pre-iterates

for(i in 1:1e6) {
        X1= Xtrans(X1,carr,lowbd)$Xout
        if(max(abs(X1-1)/boxmat[438,1:7])<1) {
           m = sum((sign(X1-1)+1)*a0)/2 + 1
           Counts[m] = Counts[m]+1 } } 

sum(Counts)
[1] 366   

## Here simulated N=1e6 numbers of which 366 fell into this box of vol= .003516.
## Thus we estimate the total volume is : .003516*1e6/366 = 9.606557
## To check equal distribution in cells:

  [1] 13  0  8  0  4  3  3  0  1  0  1  0  0  0  1  0 10  1  0  1  2  0  1  0  0  0
 [27]  1  0  0  0  4  1  3  0  1  1  0  1  3  3  3  2  3 11  0  0  0  4  2  1  0  0
 [53]  0  0 11  1  2  2 11  1  2  0  2  5  9  2  3  2  6  4  3  7  1  0  8  0  6  1
 [79]  6  4  9  0  4  3  3  2  5  7  2  1  4  2  6  1  4  2  2  6  3  4  3  0  5  0
[105]  7  3  3  3  8  6  4  3  6  4  2  2  2  1  0  1  2  2  8  5  5  4  8  7

### Could be, but it is hard to check ... might need larger sample sizes,
###   but here are some indications that convergence is not so good !

> SevnArr = array(Counts,c(rep(2,7)))

> apply(SevnArr,1,sum)
[1] 239 127
> apply(SevnArr,4,sum)
[1] 183 183
> apply(SevnArr,5,sum)
[1] 191 175
> apply(SevnArr,2,sum)
[1] 166 200
> apply(SevnArr,3,sum)
[1] 191 175
> apply(SevnArr,6,sum)
[1] 172 194
> apply(SevnArr,7,sum)
[1] 130 236
#### So there is bad imbalance in a couple of coordinate directions
###   within the box. But recall:

> round(boxmat[438,1:7],4)
[1] 0.3239 0.0910 0.1839 0.5898 0.2185 0.1422 0.2765
### So the dimensions of the box in those two directions are not tiny.

###====================================================================


### One could do this same kind of exercise with other feasible 
###    points X1 as centers of boxes in the region  !


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

BigXout = array(7e5, c(1e5,7))
for(i in 1:1e5) {
        X1= Xtrans(X1,carr,lowbd)$Xout
        BigXout[i,] = X1 }
> round(apply(BigXout,2,range),4)
        [,1]    [,2]    [,3]    [,4]    [,5]    [,6]   [,7]
[1,] -1.0070 -0.6680 -0.2874 -0.3866 -0.6064 -1.2731 0.2661
[2,]  1.7775  2.4612  2.5809  3.6276  2.7269  1.8072 3.5983

### other mathematical methods for finding regions over which 
    the generated points will be uniform:

Radon transform (heavily used in 2-dimensional and I think
    3-dimensional image reconstruction ("tomography")

> det(BigXout[1:7,]-array(1,c(7,7)))
[1] 4.476586e-13
> det(BigXout[101:107,]-array(1,c(7,7)))
[1] 8.776188e-09
> det(BigXout[100*(1:7),]-array(1,c(7,7)))
[1] 0.2387285
### these are 7! * the volumes of the regions bounded by 
###   rep(1,7) and the seven indicated MCMC generated vectors

> dmat = cbind(diag(7),rep(-1,7))
>  det(dmat %*% BigXout[100*(1:8),])
[1] -0.09907137
### 
maxvol=0
for(i in 1:1e4) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > maxvol) maxvol=tmp }
> maxvol
[1] 1.947853

for(i in 8e4+(1:1e4)) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > maxvol) maxvol=tmp }

for(i in 6e4+(1:1e4)) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > maxvol) maxvol=tmp }

for(i in 5e4+(1:1e4)) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > maxvol) maxvol=tmp }

for(i in 1e4+(1:1e4)) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > maxvol) maxvol=tmp }
                                  ### now improved to 2.1455
for(i in 2e4+(1:1e4)) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > maxvol) maxvol=tmp }

for(i in 3e4+(1:1e4)) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > maxvol) maxvol=tmp }

for(i in 7e4+(1:1e4)) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > maxvol) maxvol=tmp }

for(i in 9e4-800+(1:1e4)) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > maxvol) maxvol=tmp }

### So the maximum volume seen so far for the convex body spanned
###    by 8 outputted MCMC vectors (at lag 1 or 100)  is  
###    2.1455/gamma(8) =  0.000397
### This is not really any more promising than the way we did it before
###    using rectangular boxes.

##-------------------- Wrapping up. ------------------------------

### Here is another calculation/estimate of the overall volume of the 
###   region in which we are MCMC simulating uniformly distributed
###   random points.

### First create a test for falling in the largest-volume convex hull

for(i in 1e4+(1:1e4)) {
    tmp =det(dmat %*% BigXout[i+100*(1:8),])
    if(tmp > 2.145) cat(i, "\n") }
14490

HullPts = BigXout[14490+100*(1:8),]  ### this is the set of 8 vectors
### whose convex hull falls in the desired region and has volume
###  det(dmat %*% HullPts)=2.1455  divided by 7 factorial (ie /5040)

## Need to code the criterion that a vector falls in this convex hull.
## Start by finding a set of vectors w_j respectively orthogonal to
## the span of all HullPts[i,]-HullPts[8,], over i in (1:7)\{j}
Wmat = array(0, c(7,7))
HullAux = dmat %*% HullPts

for(j in 1:7) {
         Wmat[j,] = HullAux[j,] - t(HullAux[-j,]) %*% solve(HullAux[-j,]%*%
                      t(HullAux[-j,])) %*% HullAux[-j,]%*% HullAux[j,]
         Wmat[j,] = Wmat[j,]/sqrt(sum(Wmat[j,]^2)) }
> pupr = diag(Wmat %*% t(HullAux))
> pupr
[1] 0.9208133 0.5156522 0.4399978 0.7149128 0.9588133 0.4358079 1.2058547
> Wmat = (1/pupr)*Wmat

### net result is that HullAux %*% t(Wmat) is essentially Identity
### and also: inclusion in the convex-hull is equivalent to 
### Wmat %*% (X-HullPts[8,]) having nonneg components summing to <=1
### This is our test for inclusion.

  HullTst = function(X, v0=HullPts[8,], w=Wmat)  {
           rv= c(w %*% (X-v0))
           (min(rv) > -1e-10) & (sum(rv) < 1+1e-10) }

HullTst(rep(1,7))  ## F
HullTst(HullPts[3,])  ## T    OK

> { u0 = runif(8); HullTst(t(HullPts) %*% (u0/sum(u0)))}
[1] TRUE
> { u0 = runif(8); HullTst(t(HullPts) %*% (u0/sum(u0)))}
[1] TRUE
> { u0 = runif(8); HullTst(t(HullPts) %*% (u0/sum(u0)))}
[1] TRUE
> { u0 = runif(8); HullTst(t(HullPts) %*% (u0/sum(u0)))}
[1] TRUE
> { u0 = runif(8); HullTst(t(HullPts) %*% (u0/sum(u0)))}
[1] TRUE
> { u0 = runif(8); HullTst(t(HullPts) %*% (1.1*u0/sum(u0)))}
[1] FALSE
> { u0 = runif(8); HullTst(t(HullPts) %*% (1.1*u0/sum(u0)))}
[1] FALSE
> { u0 = runif(8); HullTst(t(HullPts) %*% (0.9*u0/sum(u0)))}
[1] FALSE                  ########## So HullTst seems to work.


Hcount = 0
for(i in 1:1e6)  { 
      X1 = Xtrans(X1,carr,lowbd)$Xout
      Hcount = Hcount+ HullTst(X1) }

> Hcount
[1] 130      ### Far too small to be a good test

IN CONCLUSION: IF WE REALLY WANTED TO CHECK UNIFORMITY, MAYBE THE
BEST WAY WOULD BE TO FIND A SET OF DISJOINT RECTANGULAR BOXES ALL
LYING INSDE THE DESIRED REGION AND IF THE TOTAL VOLUME WAS NOT TOO
SMALL THEN WE COULD DISCARD MCMC GENERATED X1 VECTORS WHICH DO NOT
FALL IN ANY OF THEM, AND THEN CHECK THAT THE RELATIVE FREQUENCIES 
OF FALLING IN INDIVIDUAL ONES OF THEM ARE COMPATIBLE WITH THE RELATIVE
VOLUMES. RECTANGULAR SUBSETS ARE MUCH EASIER TO HANDLE IN THIS WAY 
THAN CONVEX HULLS F GENERAL SETS OF 8 VECTORS.