library('rstan')
library('fields')
library('parallel')

stan_model_file = 'daily_precip_occurrence.stan'
iterations = 1500
warmup = 500
chains = 1
cores = 1#min(detectCores(), 4)

# first day
data_url = 'http://cires1.colorado.edu/~aslater/CVEN_6833/colo_pcp_daily_1997_01_11.dat'
# second day
# data_url = 'http://cires1.colorado.edu/~aslater/CVEN_6833/colo_pcp_daily_1997_01_12.dat'
# third day
# data_url = 'http://cires1.colorado.edu/~aslater/CVEN_6833/colo_pcp_daily_1997_01_13.dat'

# read into data table 
df = na.omit(read.table(data_url,na.strings='-999.999'))[1:100,]
names(df) <- c('lat','lon','elev','precip')


# first column is intercept
# scale the other covariates to reduce roundoff error
covars = cbind(1,scale(df[,c('lon','lat','elev')]))
y = (df$precip > 0) + 0

S = length(y)
P = ncol(covars) 

# distance matrix in kilometers
dist = rdist.earth(df[,c('lon','lat')],miles=FALSE)

stan_data = list(S=S,P=P,y=y,covars=covars)

# initial guesses
initf = function()list(
    sill = 10,
    range = max(dist),
    nugget = 10,
    betas = coefficients(glm(y~covars+0, family=binomial(logit)))
    )

# note this model uses gaussian process on logit residuals, which are 
# likely non-Gaussian.  consider transforming the residuals??
#
# still, the model converges nicely.  

stanfit = stan(stan_model_file, data=stan_data, chains=chains, 
    iter=iterations, warmup=warmup, cores=cores, init=initf,
    refresh=1)