######################################################################
### Uses a for loop to calculate result for each sample in xs

mydnorm1 <- function(xs,
                     mu = matrix(0,nrow(xs),1),
                     sigma = diag(1,nrow(xs))) {

  xs <- as.matrix(xs)                   # to handle args like c(0,0,0)
  nDimensions <- nrow(xs)                         # or dim(xs)[1]
  nSamples <- ncol(xs)

  if (length(mu) == 1)
    detsigma <- sigma
  else
    detsigma <- det(sigma)
  normConstant <- 1/((2*pi)^(nDimensions/2) * sqrt(detsigma))
  invSigma <- solve(sigma)

  results <- c()                        # or NULL
  for (coli in 1:nSamples) {
    x <- xs[,coli,drop=FALSE]           # see ?drop
    results <- c(results,
                 normConstant * exp(-1/2 * t(x-mu) %*% invSigma %*% (x-mu)))
  }
  results
}

######################################################################
## For loop replaced with matrix operations.

mydnorm2 <- function(xs,
                     mu = matrix(0,nrow(xs),1),
                     sigma = diag(1,nrow(xs))) {

  xs <- as.matrix(xs)                   # to handle args like c(0,0,0)
  nDimensions <- nrow(xs)                         # or dim(xs)[1]
  nSamples <- ncol(xs)

  if (length(mu) == 1)
    detsigma <- sigma
  else
    detsigma <- det(sigma)
  normConstant <- 1/((2*pi)^(nDimensions/2) * sqrt(detsigma))
  invSigma <- solve(sigma)
  xsZeroMean <- xs - matrix(rep(mu,nSamples),nDimensions,nSamples)

  normConstant * exp(-1/2 * colSums( t(t(xsZeroMean) %*% invSigma) * xsZeroMean ))
}

######################################################################
mydnorm <- function(xs,
                    mu = matrix(0,nrow(xs),1),
                    sigma = diag(1,nrow(xs))) {
  if (missing(xs)) {
    cat(" mydnorm(xs, mu, sigma)\n",
        "    xs: nDimensions x nSamples\n",
        "    mu: nDimensions x 1\n",
        " sigma: nDimensions x nDimensions\n")
    cat(" Example: 
    > samples <- matrix(seq(0,1,length=15),3,5)
    > samples
               [,1]      [,2]      [,3]      [,4]      [,5]
    [1,] 0.00000000 0.2142857 0.4285714 0.6428571 0.8571429
    [2,] 0.07142857 0.2857143 0.5000000 0.7142857 0.9285714
    [3,] 0.14285714 0.3571429 0.5714286 0.7857143 1.0000000
    > mu <- matrix(c(0.2,0.3,0.1))
    > mu
         [,1]
    [1,]  0.2
    [2,]  0.3
    [3,]  0.1
    > sigma <- matrix(c(1,0.2,0,0.2,1,0,0,0,0.5),3,3)
    > sigma
         [,1] [,2] [,3]
    [1,]  1.0  0.2  0.0
    [2,]  0.2  1.0  0.0
    [3,]  0.0  0.0  0.5
    > mydnorm(samples, mu, sigma)
    [1] 0.04224973 0.04116455 0.03389237 0.02358084 0.01386423\n")

    return()
  }
  xs <- as.matrix(xs)                   # to handle args like c(0,0,0)
  nDimensions <- nrow(xs)                         # or dim(xs)[1]
  nSamples <- ncol(xs)

  if (length(mu) == 1)
    detsigma <- sigma
  else
    detsigma <- det(sigma)
  normConstant <- 1/((2*pi)^(nDimensions/2) * sqrt(detsigma))
  invSigma <- solve(sigma)
  xsZeroMean <- xs - matrix(rep(mu,nSamples),nDimensions,nSamples)

  normConstant * exp(-1/2 * colSums( t(t(xsZeroMean) %*% invSigma) * xsZeroMean ))

  
}

######################################################################

nSamples <- 10000
nDimensions <- 10
samples <- 0.01 * matrix(1:(nSamples*nDimensions), nDimensions, nSamples)

print("With for loop:")
print(system.time(a <- mydnorm1(samples)))
print(a[1:10])

print("Without for loop:")
print(system.time(b <- mydnorm(samples)))
print(b[1:10])