support <-     ##- n=the number of balls, k=the number of boxes.
function(n, k) {  
   if(k == 1) return(as.matrix(n))
   if(n == 0) return(matrix(0, 1, k))
   X <- NULL
   for(i in 0:n) X <- rbind(X,cbind(i,Recall(n-i,k-1)))
   dimnames(X) <- NULL
   X
}
#-----------------------------------------------------------
dmultinom <- 
function(x, prob = stop("no arg")) {
   prob <- rep(prob, l=length(x))  
   prob <- prob / sum(prob)
      n <- sum(x) 
   if( any((prob<=0) & (x>0)) ) return(0) 
   logp <- ifelse( prob>0, log(prob), 0 )
   exp( lgamma(n+1) - sum(lgamma(x+1)) + sum(x*logp) )
}
#-----------------------------------------------------------
critical.value <- 
function(stat, prob, alpha = stop("no arg")) 
{
   TINY <- .Machine$single.eps     ##- be careful
   n    <- length(stat)
   idx  <- order(stat)
   stat <- stat[idx]
   prob <- prob[idx]
   cdf  <- cumsum(prob)
   cv   <- stat[n - sum(cdf > (1 - alpha)) + 1]
   alpha1 <- sum(prob[stat - cv >   TINY])   ## >
   alpha2 <- sum(prob[stat - cv > - TINY])   ## >=
   return(list(cv=cv, a1 = alpha1, a2 = alpha2))
}
#-----------
rpower <-   ##- randomized power of  size-a test
function(stat, prob, cv, alpha, alpha1, alpha2)
{
   TINY <- .Machine$single.eps     ##- be careful
   b1 <- sum(prob[stat - cv >   TINY])       ## > 
   b2 <- sum(prob[stat - cv > - TINY])       ## >= 
   return(b1 + (alpha - alpha1)/(alpha2 - alpha1) * (b2 - b1))
}
#-----------------------------------------------------------
PD1old <-
function(d,f,a,h) {   ##- a=tuning parameter, h=penalty
   HUGE <- .Machine$single.xmax
   d <- d/sum(d); f <- f/sum(f); 
   id1 <- (d > 0 )
   #-------------
   if (a < -1) {
      if ( any(!id1) ) return(HUGE)  ##- check empty cell
      tmp <- sum(d*(f/d)^(-a) - d) / (a*(a+1))
      return(tmp)
   }
   #-------------
   if (a == -1) {   ##- Kullback-Leibler
      if ( any(!id1) ) return(HUGE)  ##- check empty cell
      if ( any(f<=0) ) return(0)     ##- check empty cell of model
      return( sum( f*(log(f)-log(d)) ) )
   }
   #-------------
   if (a == 0)  {   ##- Likelihood Disparity 
      if(missing(h))  h <- 1
      d1 <- d[id1]; f1 <- f[id1]; 
      if ( any(f<=0) ) return(HUGE)  ##- check empty cell of model
      ## tmp <- sum( d1*(log(d1)-log(f1)) ) + h*(1-sum(f1))
      tmp <- sum( d1*(log(d1)-log(f1)) - d1 + f1 ) + h*(1-sum(f1))
      return( tmp )
   }
   #-------------
   if(missing(h))  h <- 1/(a+1)  
   d1 <- d[id1]; f1 <- f[id1]; 
   tmp <- sum(d1*(d1/f1)^a - (a+1)*d1 + a*f1) / (a*(a+1)) + h*(1-sum(f1))
   return(tmp)
}
#-----------------------------------------------------------
PD1 <-
function(d,f,a) {   ##- a=tuning parameter
   HUGE <- .Machine$single.xmax
   if ( length(d)==0 || length(f)==0 ) return(0)
   #-------------
   if (a < -1) return( sum(d*(f/d)^(-a)-(a+1)*d+a*f) / (a*(a+1)) )
   #-------------
   if (a == -1) {   ##- Kullback-Leibler
      if ( any(f<=0) ) return(0)     ##- check empty cell of model
      return( sum(f*(log(f)-log(d))+d-f) )
   }
   #-------------
   if (a == 0)  {   ##- Likelihood Disparity 
      if ( any(f<=0) ) return(HUGE)  ##- check empty cell of model
      return( sum(d*(log(d)-log(f))-d+f) )
   }
   #-------------
   if ( any(f<=0) ) return(HUGE)  ##- check empty cell of model
   return( sum(d*(d/f)^a-(a+1)*d+a*f)/(a*(a+1)) )
}
#-----------------------------------------------------------
PD2 <-
function(d,f, a, a.in, a.out, h) {   
   ##- a=tuning parameter, h=penalty
   ##- if either a or a.in is given, ordinary PD is performed
   ##- if both a.in and a.out is given, combined PD is performed
   ##- we can still use penalty with combined PD
   HUGE <- .Machine$single.xmax
   if(missing(a) && missing(a.in) )  stop("no arg in a or a.in") 
   if(!missing(a))     a.in <- a;
   if(missing(a.out))  a.out <- a.in;
   d <- d/sum(d);      f <- f/sum(f); 
   id0  <- d <= 0;     idin <- d < f;
   #-------------
   if (a.in <= -1) {
      if ( any(id0) ) return(HUGE)  ##- check empty cell
      din  <- d[idin];  dout <- d[!idin]
      fin  <- f[idin];  fout <- f[!idin]
      tmp1 <- PD1(din, fin, a.in) 
      tmp2 <- PD1(dout,fout,a.out) 
      return(tmp1+tmp2)
   }
   #-------------
   if(missing(h))  h <- 1/(a.in+1)  
   d1 <- d[!id0]; f1 <- f[!id0]; 
   idin <- d1 < f1
   din  <- d1[idin];  dout <- d1[!idin]
   fin  <- f1[idin];  fout <- f1[!idin]
   tmp1 <- PD1(din, fin, a.in) 
   tmp2 <- PD1(dout,fout,a.out) 
   return(tmp1+tmp2 + h*(1-sum(f1)) )
}
#-----------------------------------------------------------
PD <- PD2
#-----------------------------------------------------------
PBHM <-
function(d,f,a,h) {   ##- a=tuning parameter, h=penalty
   HUGE <- .Machine$single.xmax
   d <- d/sum(d); f <- f/sum(f);
   id1 <- (d > 0 )
   d1 <- d[id1]; f1 <- f[id1];
   #-------------
   if ( a==0 ) {
      if(missing(h))  h <- 0.5
      return( 0.5*sum((d1-f1)^2/f1) + h*(1-sum(f1)) )
   }
   #-------------
   if ( a==1 ) {    ##- a=1 and h=HUGE at empty cell
      if (missing(h)) {
         if( any(!id1) ) return(HUGE)  ##- check empty cell
         return( 0.5*sum((d-f)^2/f) )
      }
      return( 0.5*sum((d1-f1)^2/f1) + h*(1-sum(f1)) )
   }
   #-------------
   if (missing(h)) h <- 0.5*(a+1/(1-a))
   tmp <- 0.5*sum((d1-f1)^2 * (a/f1+(1-a)/(a*sqrt(d1)+(1-a)*sqrt(f1))^2))+
          h*(1-sum(f1))
   return(tmp)
   #-------------
}