# Author: Martin Tingley, Joe Werner (R adaption)
###############################################################################

setwd(BARCAST.wd);

Prior.Pars <- function(BARCAST.datdir, BARCAST.infile){
load(BARCAST.infile);
		
PRIORS<-list()
#Prior for alpha:
PRIORS$alpha<-c(0.5,.01); #[mu, sigma^2 of normal (trunc'ed normal, actually, cf. Theta_Updater)
		
#Prior for mu
#Normal (mean, variance). Take the mean to be mean(All_Inst)
#Once more, use a largish variance:
PRIORS$mu<-c( 0, .1)

#Prior for sigma2
#Inverse gamma (alpha, beta)
PRIORS$sigma2<-c(0.5,0.5);
#This is equivalent to 1 observations with an average square deviation of 2
# Very general, broad prior. 
# note, caliing the pars a and b, and the pars of the scaled inverse chi^2
# nu and s^2, we have nu=2a and s^2=b/a. These give, in the bayesian sense,
# the number of prior samples and the mean squared deviation.
# We have about 100 instrumental time series -> 100 samples -> a=50
# Those are distributed (max likelihood estimate for s^2=n/Sum(1/x_i)=0.84
# -> b = a*s^2 = 35
#PRIORS$sigma2<-c(5.0,42);

#Prior for phi. 
#RECALL we are using a log transformation in the Metropolis sampler for
#this step, so the prior is log normal. A largely uninformative prior would specify the
#range to be somewhere between 50 and 3000 km, so the log of the 
#inverse range should be between -8 and -3.5
# The log of the inverse range of a sensible temperature field should be
# in the order of -7 (= 1000 km) which is about 2 grid cells.
#PRIORS$phi<-c(-7, 1.2);
PRIORS$phi<-c(-7, 2);
#Mean then variance paramters; prior on log(phi) is normal with these
#parameters
CovMat.model <- 1

#Prior for lambda
#Normal (mu, sigma^2)
PRIORS$lambda <- c(0, .5^2)

#Prior for tau2_I
#Inverse gamma (alpha, beta)
#PRIORS$tau2.I<-c(0.5,0.5);
#This is equivalent to 1 observations (nu=2*a) with an average square
#deviation of 1 (s^2=b/a)
PRIORS$tau2.I<-c(10,0.5);

# it actually makes more sense to do everything for each proxy by itself,
# so I should flatten the loop and also set the proxy type here.

ParNum <- length(PRIORS)
for( proxIdx in grep("Prox", names(BARCAST.INPUT)) ){
  thisproxpar <- proxIdx - min(grep("Prox", names(BARCAST.INPUT)) ) + 1
  PRIORS[[ParNum + thisproxpar]] <- list();
  PRIORS[[ParNum + thisproxpar]]$type <- 0;
  PRIORS[[ParNum + thisproxpar]]$tau2.P.1<-c(0.5,.5);
  PRIORS[[ParNum + thisproxpar]]$Beta.0.1 <-c(0 , .2^2);
  PRIORS[[ParNum + thisproxpar]]$Beta.1.1 <- c( (1-PRIORS[[ParNum + thisproxpar]]$tau2.P.1[2]/(PRIORS[[ParNum + thisproxpar]]$tau2.P.1[1]+1)*(1-mean(PRIORS$alpha)^2)/(PRIORS$sigma2[2]/(PRIORS$sigma2[1]-1)))^(1/2), .2^2);
  PRIORS[[ParNum + thisproxpar]]$Beta.2.1 <- c(0, .1^2);
  names(PRIORS)[ParNum + thisproxpar] <- paste("PROXY",thisproxpar,sep=".")
}

PRIORS[[ParNum + 1]]$type <- 20;

#Prior for the T(0) vector: Normal 
#with mean given by mean of all inst and standard deviation 2
#times the std of all Inst.
#So the prior parameters (mean, var) are: 
PRIORS$T.0<-PRIORS$mu
#Note that this assumes a constant mean and diagonal cov mat. 
		
		
# SAVE the prior structure:
save(PRIORS,file = paste(BARCAST.datdir, "/PRIORSvNewMeth1.R", sep="") );
#save PRIORS_vNewMeth1 PRIORS

## ALSO set the MCMC jumping parameters
#The phi step.
# The jumping distribution of log(phi) is normal with mean zero; this sole paramter is
# the VARIANCE. We expect the posterior variance to be far smaller than the
# prior variance, so the jumping variance is set low. this can be adjusted
# as needed. 
#ALSO Specify the number of iterations of the Metropolis step for each iteration
#of the Gibbs: the Metroplis sampler needs to (come close to) converging each
#time. Until the algorithm settles down, this will take a little while.
#NOTE that this step is not time consuming. 
#First par is variance, second is number.
MHpars<-list()
MHpars$log.phi<-c(.1^2, 500);

save(MHpars, file=paste( BARCAST.datdir, "/MHparsvNewMeth1.R", sep="") )
}