library(gdata)


readO<-function() {
obs<-NULL
obs<-cbind(obs,read.table("../data/csiro/rain.cy.murrdarl.5pc.area.dat",skip=1)$V2)
obs<-cbind(obs,read.table("../data/csiro/rain.cy.nsw.5pc.area.dat",skip=1)$V2)
obs<-cbind(obs,read.table("../data/csiro/rain.cy.nwquad.5pc.area.dat",skip=1)$V2)
obs<-cbind(obs,read.table("../data/csiro/rain.cy.qld.5pc.area.dat",skip=1)$V2)
obs<-cbind(obs,read.table("../data/csiro/rain.cy.swaust.5pc.area.dat",skip=1)$V2)
obs<-cbind(obs,read.table("../data/csiro/rain.cy.swquad.5pc.area.dat",skip=1)$V2)
obs<-cbind(obs,read.table("../data/csiro/rain.cy.victas.5pc.area.dat",skip=1)$V2)
obs
}
readE<-function() {
exp<-NULL
exp<-c(exp,list(read.xls("../data/csiro/Percent_area_below_5thpercentile_Rainfall_MDB.xls",sheet=1)[1:141,]))
exp<-c(exp,list(read.xls("../data/csiro/Percent_area_below_5thpercentile_Rainfall_NSW.xls",sheet=1)[1:141,]))
exp<-c(exp,list(read.xls("../data/csiro/Percent_area_below_5thpercentile_Rainfall_NWAust.xls",sheet=1)[1:141,]))
exp<-c(exp,list(read.xls("../data/csiro/Percent_area_below_5thpercentile_Rainfall_Qld.xls",sheet=1)[1:141,]))
exp<-c(exp,list(read.xls("../data/csiro/Percent_area_below_5thpercentile_Rainfall_SW-WA.xls",sheet=1)[1:141,]))
exp<-c(exp,list(read.xls("../data/csiro/Percent_area_below_5thpercentile_Rainfall_SWAust.xls",sheet=1)[1:141,]))
exp<-c(exp,list(read.xls("../data/csiro/Percent_area_below_5thpercentile_Rainfall_VicTas.xls",sheet=1)[1:141,]))
exp
}
readOt<-function() {
	obs<-NULL
	obs<-cbind(obs,read.table("../data/csiro/tmean.cy.murrdarl.95pc.area.dat",skip=1)$V2)
	obs<-cbind(obs,read.table("../data/csiro/tmean.cy.nsw.95pc.area.dat",skip=1)$V2)
	obs<-cbind(obs,read.table("../data/csiro/tmean.cy.nwquad.95pc.area.dat",skip=1)$V2)
	obs<-cbind(obs,read.table("../data/csiro/tmean.cy.qld.95pc.area.dat",skip=1)$V2)
	obs<-cbind(obs,read.table("../data/csiro/tmean.cy.swaust.95pc.area.dat",skip=1)$V2)
	obs<-cbind(obs,read.table("../data/csiro/tmean.cy.swquad.95pc.area.dat",skip=1)$V2)
	obs<-cbind(obs,read.table("../data/csiro/tmean.cy.victas.95pc.area.dat",skip=1)$V2)
	obs
}
readEt<-function() {
	exp<-NULL
	exp<-c(exp,list(read.xls("../data/csiro/Percent_area_above_95thpercentile_Temperature_MDB.xls",sheet=1)[1:141,]))
	exp<-c(exp,list(read.xls("../data/csiro/Percent_area_above_95thpercentile_Temperature_NSW.xls",sheet=1)[1:141,]))
	exp<-c(exp,list(read.xls("../data/csiro/Percent_area_above_95thpercentile_Temperature_NWAust.xls",sheet=1)[1:141,]))
	exp<-c(exp,list(read.xls("../data/csiro/Percent_area_above_95thpercentile_Temperature_Qld.xls",sheet=1)[1:141,]))
	exp<-c(exp,list(read.xls("../data/csiro/Percent_area_above_95thpercentile_Temperature_SW-WA.xls",sheet=1)[1:141,]))
	exp<-c(exp,list(read.xls("../data/csiro/Percent_area_above_95thpercentile_Temperature_SWAust.xls",sheet=1)[1:141,]))
	exp<-c(exp,list(read.xls("../data/csiro/Percent_area_above_95thpercentile_Temperature_VicTas.xls",sheet=1)[1:141,]))
	exp
}


#O<-readO()
#E<-readE()
#Ot<-readOt()
#Et<-readEt()

regions<-c("MDB","NSW","NWAust","Qld","SW-WA","SWAust","VicTas")
models<-c("cgcm-t47", "cgcm-t63", "csiro-mk3.0", "csiro-mk3.5", "giss-aom", "giss-er",   "iap",  "inmcm",   "ipsl", "miroc-h", "miroc-m",   "mri", "ncar-ccsm")


# Fit trend line

linear<- function(obs,exp) {
t<-1:108
#Get correlations
l<-lm(obs ~ exp)
s<-summary(l)
r<-c(s$r.squared, s$coefficients[2,4])
}

# Return periods test

returns<-function(m1,m2) {
r1<-diff(which(m1>0))-1
r2<-diff(which(m2>0))-1
t.test(r1,r2)$p.value
}


# Slopes test

test_slopes<-function(x,y,t) {
fit1 <- lm(x~t)
s1 <- summary(fit1)$coefficients
fit2 <- lm(y~t)
s2 <- summary(fit2)$coefficients
db <- (s2[2,1]-s1[2,1])
sd <- sqrt(s2[2,2]^2+s1[2,2]^2)
df <- (fit1$df.residual+fit2$df.residual)
td <- db/sd
c(s1[2,1],s2[2,1],2*pt(-abs(td), df)) 
}

#NS-Efficiency

efficiency<-function(obs,exp) {
1-sum((obs-exp)^2)/sum((obs-mean(obs))^2)	
}

# Table of stats for a region

do_region<-function(region=1) {
	results<-NULL
	for (j in 2:14) {
		r<-NULL
		obs<-O[,region]
		exp<-as.vector(E[[region]][j])[1:108,1]
		r<-c(linear(obs,exp),efficiency(obs,exp),returns(obs,exp))
		results<-rbind(results,r)
	}
	rownames(results)<-models
	colnames(results)<-c("r2","corr-sig","ns-eff","returnp-p")
	results
}

# Table of stats for all regions

do_all<-function() {
	result<-NULL
	for(i in 1:7) {
		d<-do_region(i) 
		result<-c(result,list(as.data.frame(d)))
	}
	result
}

#Distribution of trends for a region

do_trendt<-function(region=1) {
	results<-NULL
	for (j in 2:14) {
		r<-NULL
		obs<-O[,region]
		exp<-as.vector(E[[region]][j])[1:108,1]
		r<-test_slopes(obs,exp,1:108)
		results<-rbind(results,r)
	}
	rownames(results)<-models
	colnames(results)<-c("trend-o","trend-e","slopes-p")
	t<-t.test(results[,2],mu=results[1,1])
	c(results[1,1],mean(results[,2]),sd(results[,2]),t$statistic,t$p.value)
}

#Table with all trends

do_all_trends<-function() {
	result<-NULL
	for(i in 1:7) {
		d<-do_trendt(i) 
		result<-rbind(result,d)
	}
	rownames(result)<-regions
	colnames(result)<-c("trend-o","mean-trend-e","sd-trend-e","t-test","p")
	result
}

#Figure 1

do_dist_plot<-function() {
	d<-read.xls("../data/csiro/Percent_area_below_5thpercentile_Rainfall_SW-WA.xls",sheet=1)[1:141,]
	pr<-as.matrix(d[1:107,2:14])
	dim(pr)<-NULL
	fr<-as.matrix(d[108:141,2:14])
	dim(fr)<-NULL
	p<-hist(fr,plot=F)
	plot(p$breaks[1:(length(p$breaks)-1)],p$density,col="red",xlab="Percentage area low rainfall",type="l",ylab="Proportion")
	f<-hist(pr,plot=F,breaks=p$breaks)
    lines(p$breaks[1:(length(p$breaks)-1)],f$density,col="blue")
}

go<-function() {
	print(do_all_trends())
	print(do_all())
	do_dist_plot()
}

#Distribution of trends for a region

freq<-function(obs) {
	length(obs[obs>0])/length(obs)
}


do_rankt<-function(Ot,Et,region=1) {
	results<-NULL
	for (j in 2:14) {
		r<-NULL
		obs<-Ot[,region]
		exp<-as.vector(Et[[region]][j])[1:98,1]
		r1<-freq(obs[1:49])-freq(exp[1:49])
		r2<-freq(obs[50:98])-freq(exp[50:98])
		results<-rbind(results,c(r1,r2))
	}
	rownames(results)<-models
	colnames(results)<-c("freq1","freq2")
	results
}

do_rank_all<-function(O,E,region=1) {
	result<-NULL
	for(i in 1:7) {
		r<-do_rankt(O,E,region=i)
		l<-lm(r[,2]~r[,1])
		s<-summary(l)$coefficients[2,4]
		result<-c(result,s)
	}
	result
}

rank_table<-function()
{
 table<-NULL
 t1<-do_rank_all(Ot,Et)
 t2<-do_rank_all(O,E)
 table<-cbind(table,t1)
 table<-cbind(table,t2)
 rownames(table)<-regions
 colnames(table)<-c("95%Temp","5%Rain")
 table
}