###############################################################################
#
# test_mvn.R 
#
# You are welcome to use this code for your own research but please 
# properly cite and attibute:
#
# Author: 
# Cameron Bracken, cameron.bracken@gmail.com
#
# Citation: 
# C. Bracken, K. Holman, B. Rajagopalan, H. Moradkhani, A Bayesian 
# hierarchical approach to multivariate nonstationary hydrologic frequency 
# analysis, Water Resources Research, 2017, Under Review. 
#
# This script performs tests for multivariate normality, the data must pass 
# in order to use the gaussian elliptical copula. 
#
###############################################################################

source('lib.R')

load_packages(c('ggplot2','data.table','reshape2','ismev','corpcor','evd','MVN'))

model_years = 1981:2015

snow_dt = fread('data/station_annmax_SWE_1979-2015.csv')
res = fread('data/annmax_WSE_Inflow_TPD_1963-2015.csv')

#snow_dt = fread('data/FD_station_annmax_SWE_1979-2016.csv')
#res = fread('data/FD_annmax_WSE_Inflow_1966-2016.csv')

snow_mat = data.matrix(dcast.data.table(snow_dt[WY %in% model_years],WY~id,value.var='mWESD')[,WY:=NULL])/1000 #km
snow = apply(snow_mat,1,mean,na.rm=TRUE)#snow_mat[,'USS0006K08S']

#snow_complete = snow_mat[,apply(snow_mat,2,function(x)!any(is.na(x)))]
#pca = prcomp(snow_complete)
#pca_covariates = scale(t(t(pca$rotation[,1])%*%t(snow_complete)))

flow = res[WY %in% model_years]$INmax/1000 # kcfs
elev = res[WY %in% model_years]$WSEmax/1000 #km

enso_winter = fread('data/mei_djf.csv')
pdo_winter = fread('data/pdo_djf.csv')
amo_winter = fread('data/amo_djf.csv')

indicies_winter = merge(merge(enso_winter,pdo_winter,by=c('year')),
    amo_winter,by=c('year'))[year %in% model_years]
#indicies = merge(enso,pdo,by=c('year'))[year %in% model_years]
covars_snow = cbind(1,scale(as.matrix(indicies_winter)))


enso_spring = fread('data/mei_mam.csv')
pdo_spring = fread('data/pdo_mam.csv')
amo_spring = fread('data/amo_mam.csv')

indicies_spring = merge(merge(enso_spring,pdo_spring,by=c('year')),
    amo_spring,by=c('year'))[year %in% model_years]
#indicies = merge(enso,pdo,by=c('year'))[year %in% model_years]
covars_flow = cbind(1,scale(as.matrix(indicies_spring)))

covars_elev = cbind(1,scale(model_years))

fit_snow = gev.fit(apply(snow_mat,1,mean,na.rm=TRUE),covars_snow,mul=2:5,show=F,method='BFGS')
fit_flow = gev.fit(flow,covars_flow,mul=2:5,show=F,method='BFGS')
fit_elev = gev.fit(elev,show=F,method='BFGS')


gev_normal = cbind(
    qnorm(with(fit_snow,pgev(snow,vals[,1],vals[1,2],vals[1,3]))),
    qnorm(with(fit_flow,pgev(flow,vals[,1],vals[1,2],vals[1,3]))),
    qnorm(with(fit_elev,pgev(elev,vals[,1],vals[1,2],vals[1,3]))))

print(roystonTest(gev_normal))
print(mardiaTest(gev_normal))
print(hzTest(gev_normal))

fit_snow2 = gev.fit(apply(snow_mat,1,mean,na.rm=TRUE),show=F,method='BFGS')
fit_flow2 = gev.fit(flow,show=F,method='BFGS')
fit_elev2 = gev.fit(elev,show=F,method='BFGS')

gev_normal2 = cbind(
    qnorm(with(fit_snow2,pgev(snow,vals[,1],vals[1,2],vals[1,3]))),
    qnorm(with(fit_flow2,pgev(flow,vals[,1],vals[1,2],vals[1,3]))),
    qnorm(with(fit_elev2,pgev(elev,vals[,1],vals[1,2],vals[1,3]))))

print(roystonTest(gev_normal2))
print(mardiaTest(gev_normal2))
print(hzTest(gev_normal2))