## See a short tutorial at: 
## "http://membres-timc.imag.fr/Olivier.Francois/tutorial_spfa.html"

spfa <- function(data,coord,K=3,scale=0.01, noise = 0.001) {
 
  
  library(fields);
  ####### INPUT ################
  # data matrix of size N individuals x M loci
  # coord, matrix of spatial coordinates of size Nx2
  # requires longitude,latitude distinct for all individuals
  # N, the number of individuals
  # M, the number of loci
  # K the number of factors
  # scale, the parameter to tune
  
  ####### OUTPUT ##############
  # U of size NxK
  # V of size MxK
  # sigma, the covariance matrix used.
  
  # dimensions
  data = as.matrix(data);
  dd = dim(data);
  N = dd[1];
  M = dd[2];

  #impute missing data (naive)
  for (j in 1:M) data[is.na(data[,j]) ,j] = sample( data[!is.na(data[,j]) ,j], sum(is.na(data[,j])), rep = T)
  
  #print(class(data))

  #print(data)

  #centering
  data = data - matrix( rep( apply(data,MARGIN=2,mean), N), nrow = N, byrow = T  );

  geo = as.matrix(coord + cbind(rnorm(N,sd = noise),rnorm(N,sd = noise)) ); 
  dg = dim(geo);
  D = dg[2];
  
  # calculate m
  # for euclidian distance
  sigma = as.matrix(rdist.earth(geo,miles=FALSE))
  m = scale*mean(sigma);
  
  # calculate sigma
  sigma = exp(-sigma/m)
  
  # cholesky of the inverse of sigma
  Ch = chol(solve(sigma));
  
  # first svd
  res = svd(Ch%*%data,K,K);
  
  # second svd
  res2 = svd(solve(Ch) %*% res$u %*% diag(res$d[1:K]) %*% t(res$v),K,K)
  
  U = res2$u %*% diag(res2$d[1:K]);
  V = res2$v;
  l =list("U"=U,"sigma"=sigma,"V"=V,"mean"=m/scale,"diag"=res2$d[1:K],"diag2"=res$d[1:K])
  
  #The new factors are in U: plot(object$U) 
  return(l)
}