# mvnLL.R : Defines the function mvLL
#           source("mvnLL.R"), probably in a program defining another function.
#
mvnLL <- function(Sigma,Ymat,check=T,centered=F)
#        -2 LL for a multivariate normal: Function of the covariance matrix 
    {
    p <- dim(Ymat)[2]
    if(check)
        {
        p1 <- dim(Sigma)[1] ; p2 <- dim(Sigma)[2]
        if(p1 != p2)
            {cat("     Sigma = \n")
            print(Sigma)
            stop("Sigma must be square!")
            }
        if(p1 != p) stop("Num cols of Y must = num cols of Sigma.")
        } # End checking of input
    deter <- det(Sigma) 
#    print(Sigma) ; print(deter) # For debugging
    # If determinant is neg or very near zero, return +Infinity
    # Use the default tolerence value for solve
    if(deter < 1e-7) mvnLL <- Inf else
       {
       siginv <- solve(Sigma)
       nn <- dim(Ymat)[1]
       const <- nn*p*log(2*3.1415926)
       sumq <- 0 # Initializing sum of quadratic forms.
       for(i in 1:nn) sumq <- sumq + 
            as.numeric( Ymat[i,]%*%siginv%*%cbind(Ymat[i,]) )
       if(! centered)
           {
           ybar <- cbind(apply(Ymat,2,mean)) # Y-bar is a column vector
           sumq <- sumq - nn*t(ybar)%*%siginv%*%ybar
           }
       mvnLL <- as.numeric( nn*log(deter) + const + sumq ) # Make it a number
       } # End if the determinant is okay
#    cat("    mvnLL = \n") ; print(mvnLL)
    mvnLL
    } # End of function mvnLL