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

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

data_url = 'http://cires1.colorado.edu/~aslater/CVEN_6833/colo_monthly_precip_01.dat'
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 = as.vector(scale(df$precip,center=F))

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(lm(y~covars+0))
    )

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