library(lmomRFA)
library(parallel)
library(data.table)

source("findLCV_h_05.R")
outfile = "output.csv"

cl = makeCluster(detectCores(),outfile='')

### READ IN ALL GRIDS FOR TVA LMOMENTS
mean=na.omit(fread("MLC_48hr_AtSiteMean_P2_08282014.csv"))
lcv=na.omit(fread("MLC_48hr_LCv_09032014.csv"))$value
lsk=na.omit(fread("MLC_48hr_LSkew_09032014.csv"))$value
latlon=with(mean,cbind(lon,lat))
mean = mean$value

lcv2 = lcv*mean
lkurt = 0.17
lmom = data.table(cbind(mean,lcv2,lsk))

# comment this out to do the whole dataset
#s = sample(1:nrow(lmom),1000)
#lmom = lmom[s]

lkurt = parApply(cl,lmom,1,findLK05_optim)
lmom[,lkurt:=lkurt]
kappar = t(parApply(cl,lmom,1,pelkap))

diff05 = abs(kappar[,4]-0.05)
aa = which(diff05<0.005)

quakap2=function(x){
	qqq=c(0.90,0.99,0.999,0.9999,0.99999)
	quakap(qqq,x)
}

quants = t(parApply(cl,kappar,1,quakap2))

final = data.table(cbind(latlon,kappar,quants))
colnames(final) = c("lon","lat","xi","alp","kap","hon","r10","r100","r1000","r10000","r100000")
write.csv(final,file=outfile,row.names=F,quote=F)