svdinv=function(a,tol=1e-6){
 # psuedo inverse of  a matrix.
 sva=svd(a)
 maxsv=max(sva$d)
 dd=ifelse(sva$d<tol*maxsv,0,1/sva$d)
 t(sva$u %*% (t(sva$v)*dd))
}

# generate multivariate Normals
rmultnorm=function(n,mu,sigma){
   # returns n rows of multivariate normal with specified mu and sigma
   # ie: returns a n x length(mu) matrix
  p=length(mu)
  z=matrix(rnorm(n * p),nrow=n)
  ch=chol(sigma,pivot=T)
  piv=attr(ch,"pivot")
  zz=(z%*%ch)
  zzz=0*zz
  zzz[,piv]=zz
  zzz + matrix(mu,nrow=n,ncol=p,byrow=T)
}

# generate a vectorized distance matrix
dist=function(x){
  # input is a matrix (typically a data matrix). The distances
  #   computed will be among the rows of x.
  # returns a vector of Eucledian distances among the rows of x (rather
  #   than a symmetric matrix).  For i < j, the distance between 
  #   row i and row j is element nrow(x)*(i-1) - i*(i-1)/2 + j-i of the result
  x=as.matrix(x)
  dists=numeric(nrow(x)*(nrow(x)-1)/2)
  for(i in 1:(nrow(x)-1)) {
    for(j in (i+1):nrow(x)) {
      dists[nrow(x)*(i-1) - i*(i-1)/2 + j-i]=sqrt(t(x[i,]-x[j,])%*%(x[i,]-x[j,]))
    }
  }  
  return(dists)
}

covpow=function(locs,pow=2,scale=5){
  # browser()
  d1=dist(locs)
  n=dim(locs)[1]
  mat=matrix(0,n,n)
  mat[lower.tri(mat)]=d1
  mat=mat+t(mat)
  cc=exp(-(mat/scale)^pow)
  return(cc)
}