####### Make a plot of the best MR designs for static naive vs informed models of Tun Mustapha Park Fishery model
#
#Supplementary file for:
# Christopher J. Brown, Crow White, Maria Beger, Hedley S. Grantham, Benjamin S. Halpern, Carissa J. Klein, Peter J. Mumby, Vivitskaia J.D. Tulloch, Mary Ruckelshaus and Hugh P. Possingham. Fisheries and biodiversity benefits of using static versus dynamic models for designing marine reserve networks. Ecosphere
#
#CJ Brown 11 Jan 2014

setwd('mydir/Mod_TMP/results')

dat =read.csv('results_base.csv', header = T)
nvals = 10
biovals = seq(5, 95, length.out = nvals)

profvals.stat = 99 + -1*biovals

profvals.dyn = profvals.stat  + 1*biovals - (0.01*biovals^2)

#addvals = seq(0,40, length.out = nvals/2)
#profvals.dyn = profvals.stat + c(addvals, rev(addvals))

#plot parameters
alims = c(0,100)

cexpts = 2
statcol = 'grey70' #colour for short-term informed
dyncol = 'black' #actual colour for naive
cexlab = 1.9
axwd = 3
lnwd=3
cexax = 1.5
cexleg=1.4
flabx = -2
flaby = 108
flabcex = 1.8

#arrows
iarvals = 5
arvals2 = 67
arrlen = 0.15
arrwd = 4
arcol = 'grey20'

##### Figure
dev.new(width = 11, height = 10)
par(mar = c(5,6,4,2), xpd=T, las = 1, mfrow = c(2,2))


##### TMP figure

### Parameters
datn =dat[1:nrow(dat),]*100 #new database
ylims = c(0, round(max(datn[,5:6]),-1))
xlims = c(0, round(max(datn[,1]),-1))
xtmpticks = seq(xlims[1], xlims[2], by = 10)
xtmplabs = xtmpticks

ytmpticks = seq(ylims[1], ylims[2], by = 20)
ytmplabs = ytmpticks


##Figure A

plot(datn$constarg, datn$naivesave_bhigh_long, xlim = xlims, ylim = ylims, pch = 18, cex = cexpts, col = dyncol,
xaxt='n', yaxt ='n', bty ='n', xlab = 'Habitat reservation target (%)', ylab ='Fishery profits (%)', cex.lab = cexlab, xaxs='i', yaxs='i')

axis(side = 1, labels = xtmplabs, at = xtmpticks, lwd= axwd, cex.axis = cexax, pos = 0)
axis(side = 2, labels = ytmplabs, at = ytmpticks, lwd= axwd, cex.axis = cexax)

lines(datn$constarg, datn$naivesave_bhigh_long, col = dyncol, lwd = lnwd)

points(datn$constarg, datn$statsave_bhigh_long, pch = 17, cex = cexpts, col = statcol)
lines(datn$constarg, datn$statsave_bhigh_long, col = statcol, lwd = lnwd)

arrows(datn$constarg, datn$statsave_bhigh_long, datn$constarg, datn$naivesave_bhigh_lon, length = 0, lwd = 1.3, col = 'grey50')


text(flabx, flaby, 'A: Malaysia long-term', cex = flabcex, pos = 4)

legend(x=32, y = 105, legend = c('Naive', 'Realistic'), lty = 1, col = c(dyncol, statcol), pch = c(18,17), bty ='n', cex= cexleg, lwd = lnwd) 

## Reserve size arrow
text(15,55,'Greater area in reserves', cex = cexlab, srt = 325, pos = 4, col = arcol)
arrows(15, 47, 39, 25, length = arrlen, lwd = arrwd, col = arcol, lty =1)

##Figure B

plot(datn$constarg, datn$naivesave_bhigh_short, xlim = xlims, ylim = ylims, pch = 18, cex = cexpts, col = dyncol,
xaxt='n', yaxt ='n', bty ='n', xlab = 'Habitat reservation target (%)', ylab ='Fishery profits (%)', cex.lab = cexlab, xaxs='i', yaxs='i')

axis(side = 1, labels = xtmplabs, at = xtmpticks, lwd= axwd, cex.axis = cexax, pos = 0)
axis(side = 2, labels = ytmplabs, at = ytmpticks, lwd= axwd, cex.axis = cexax)

lines(datn$constarg, datn$naivesave_bhigh_short, col = dyncol, lwd = lnwd)

points(datn$constarg, datn$statsave_bhigh_short, pch = 17, cex = cexpts, col = statcol)
lines(datn$constarg, datn$statsave_bhigh_short, col = statcol, lwd = lnwd)

arrows(datn$constarg, datn$statsave_bhigh_short, datn$constarg, datn$naivesave_bhigh_short, length = 0, lwd = 1.3, col = 'grey50')

text(flabx, flaby, 'B: Malaysia short-term', cex = flabcex, pos = 4)

###############################################################
################# 					SCB MODEL 				#############
###############################################################

##### Figure

##### SCB figure
setwd('/Users/uqcbrow9/Documents/CT project/static dynamic comparison/Jmodels/SCB results')
scbdat = read.csv('table_save_equil_v7.csv', header=F)
scbdat[,2:7] = scbdat[,2:7]*100
# alpha, 1
# naive static profit then biomass, 2 3 
# informed static profit then biomass, 4 5  
# dynamic profit then biomass 6 7

### Parameters
datn =dat[1:10,]*100 #new database
ylims = c(0, 100)
xlims = c(50, 100)
xtmpticks = seq(xlims[1], xlims[2], by = 10)
xtmplabs = xtmpticks

ytmpticks = seq(ylims[1], ylims[2], by = 20)
ytmplabs = ytmpticks


##Figure

plot(scbdat[,3], scbdat[,2], xlim = xlims, ylim = alims, pch = 18, cex = cexpts, col = dyncol,
xaxt='n', yaxt ='n', bty ='n', xlab = 'Biomass (%)', ylab ='Fishery profits (%)', cex.lab = cexlab, xaxs='i', yaxs='i')

# text(5,-38,'Biomass (%)', cex = cexlab, pos= 4)

axis(side = 1, labels = xtmplabs, at = xtmpticks, lwd= axwd, cex.axis = cexax, pos = 0)
axis(side = 2, labels = ytmplabs, at = ytmpticks, lwd= axwd, cex.axis = cexax)

## Add connecting lines
arrows(scbdat[,3], scbdat[,2], scbdat[,5], scbdat[,4], length = 0, lwd = 1.3, col = 'grey50')

points(scbdat[,5], scbdat[,4], pch = 17, cex = cexpts, col = statcol)

lines(scbdat[,3], scbdat[,2], col = dyncol, lwd = lnwd)

lines(scbdat[,5], scbdat[,4], col = statcol, lwd = lnwd)

text(40, flaby, 'C: California long-term', cex = flabcex, pos = 4)

######### NPV results

scbdat = read.csv('table_save_npv_v7.1.csv', header=F)
scbdat[,2:7] = scbdat[,2:7]*100
# alpha, 1
# naive static profit then biomass, 2 3 
# informed static profit then biomass, 4 5  
# dynamic profit then biomass 6 7

### Parameters
datn =dat[1:10,]*100 #new database
ylims = c(0, 100)
xlims = c(50, 100)
xtmpticks = c(seq(xlims[1], xlims[2], by = 10), 100)
xtmplabs = xtmpticks

ytmpticks = seq(ylims[1], ylims[2], by = 20)
ytmplabs = ytmpticks


##Figure

plot(scbdat[,3], scbdat[,2], xlim = xlims, ylim = alims, pch = 18, cex = cexpts, col = dyncol,
xaxt='n', yaxt ='n', bty ='n', xlab = 'Biomass (%)', ylab ='Fishery profits (%)', cex.lab = cexlab, xaxs='i', yaxs='i')

# text(5,-38,'Biomass (%)', cex = cexlab, pos= 4)

axis(side = 1, labels = xtmplabs, at = xtmpticks, lwd= axwd, cex.axis = cexax, pos = 0)
axis(side = 2, labels = ytmplabs, at = ytmpticks, lwd= axwd, cex.axis = cexax)

## Add connecting lines
arrows(scbdat[,3], scbdat[,2], scbdat[,5], scbdat[,4], length = 0, lwd = 1.3, col = 'grey50')

points(scbdat[,5], scbdat[,4], pch = 17, cex = cexpts, col = statcol)

lines(scbdat[,3], scbdat[,2], col = dyncol, lwd = lnwd)

lines(scbdat[,5], scbdat[,4], col = statcol, lwd = lnwd)

# legend(x=80, y = 105, legend = c('Naive', 'Informed'), lty = 1, col = c(dyncol, statcol), pch = c(18,17), bty ='n', cex= cexleg, lwd = lnwd) 


text(40, flaby, 'D: California short-term', cex = flabcex, pos = 4)