#+
# run_pca_on_arflib_866.r
# 	read in the entire sample library, run PCA on it, and save the eigenvalues and eigenvectors to text files
#
# INPUTS: edit this file to set these values
#	callib.dir	path name to directory containing the calibration library, a set of simulated ARFs
#	pattern.fil	pattern to match to find the files in callib.dir
#	output.dir	path name to directory to put the PCA results in
#	sqlamb.fil	output filename to hold the sqrt of the eigenvalues
#	pcomps.fil	output filename to hold the Principal Components
#	avgarf.fil	output filename to hold the average of the samples
#	numcomps	number of PCs to keep in the outputs (make sure this is .LE. size of callib.dir)
# 
# execute as source("run_pca_on_arflib_866.r")
# requires R v3+
#
# Vinay Kashyap (2014apr; based on Hyunsook Lee's /data/susfu/hyunsook/fumtu/BX/GIBBS/OLD/pca.r)
#-

require("FITSio")

callib.dir <- "store_arfs/"
pattern.fil <- "866_"
output.dir <- "./"
sqlamb.fil <- "sqlamb_866.txt"
pcomps.fil <- "pcomps_866.txt"
avgarf.fil <- "avgarf_866.txt"
fils.lst <- list.files(path=callib.dir,pattern=pattern.fil)
n.out <- length(fils.lst)
numcomps <- n.out

# read ARFs from calibration library
arfs <- numeric()
for (i in 1:n.out){
  tmp.fil <- paste(callib.dir,fils.lst[i],sep="")
  print(tmp.fil)
  tmp <- readFITS(file=tmp.fil,hdu=1)
  var <- tmp$col[[3]]
  if (i == 1) numcols <- length(var)
  arfs <- cbind(arfs,var)
}

# carry out PCA
ARF <- t(arfs)
pc.ARF <- prcomp(ARF,center=T,scale=F,retx=T)
Amean <- pc.ARF$center

#write outputs
write(as.numeric(pc.ARF$sdev[1:numcomps]),sqlamb.fil,ncol=1)
write(pc.ARF$rotation[,1:numcomps],pcomps.fil,ncol=numcols)
write(Amean,avgarf.fil,ncol=numcols)