####### Make a plot of the best MR designs for static and dynamic 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 10 Jan 2014
#Runs with new scripts for compiling results

setwd('mydir/Mod_TMP/results')
dat =read.csv('results_base.csv', header = T)

names(dat)

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 = 1.8
statcol = 'red'
dyncol = 'blue'
cexlab = 1.9
axwd = 3
lnwd=3
cexax = 1.5
cexleg=1.5
flabx = -5
flaby = 110
flabcex = cexax

#arrows
iarvals = 7
arvals2 = 75
arrlen = 0.15
arrwd = 4
arcol = 'grey20'

##### Figure
dev.new(width = 10, 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[,3:4]),-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

plot(datn$constarg, datn$statsave_bhigh_long, xlim = xlims, ylim = alims, pch = 17, cex = cexpts, col = statcol,
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)
axis(side = 2, labels = ytmplabs, at = ytmpticks, lwd= axwd, cex.axis = cexax)


lines(datn$constarg, datn$statsave_bhigh_long, col = statcol, lwd = lnwd)

points(datn$constarg, datn$dynsave_bhigh_long, pch = 22, cex = cexpts, col = dyncol, bg ='transparent')
lines(datn$constarg, datn$dynsave_bhigh_long, col = dyncol, lwd = lnwd)

arrows(datn$constarg[iarvals], datn$statsave_bhigh_long[iarvals], datn$constarg[iarvals], datn$dynsave_bhigh_long[iarvals]-1.5,
length = arrlen, lwd = arrwd, col = arcol)

arrows(datn$constarg[iarvals], datn$statsave_bhigh_long[iarvals], 57, datn$statsave_bhigh_long[iarvals],
length = arrlen, lwd = arrwd, col = arcol)

legend(x=27, y = 115, legend = c('Realistic static scenario', 'Best-case scenario','Value of dynamic model'), lty = 1, col = c(statcol, dyncol, arcol), pch = c(17,22,NA), bty ='n', cex= cexleg, lwd = lnwd) 

text(-2, 110, 'A', cex = cexlab)

# text(10,60,'Larger reserves', cex = cexlab, srt = 320, pos = 4, col = arcol)
# arrows(11,50,41,22, length = arrlen, lwd = arrwd, col = arcol, lty =1)


# arrows(15, 99.2, 50, 99.2,
# length = arrlen, lwd = lnwd, col = arcol)

#### TMP short-term
### Parameters
datn =dat[1:nrow(dat),]*100 #new database
ylims = c(0, round(max(datn[,6:7]),-1))
ytmpticks = seq(ylims[1], ylims[2], by = 20)
ytmplabs = ytmpticks


##Figure

plot(datn$constarg, datn$statsave_bhigh_short, xlim = xlims, ylim = alims, pch = 17, cex = cexpts, col = statcol,
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)
axis(side = 2, labels = ytmplabs, at = ytmpticks, lwd= axwd, cex.axis = cexax)


lines(datn$constarg, datn$statsave_bhigh_short, col = statcol, lwd = lnwd)

points(datn$constarg, datn$dynsave_bhigh_short, pch = 22, cex = cexpts, col = dyncol, bg ='transparent')
lines(datn$constarg, datn$dynsave_bhigh_short, col = dyncol, lwd = lnwd)

text(-2, 110, 'B', cex = cexlab)

# text(20,60,'Larger reserves', cex = cexlab, srt = 320, pos = 4, col = arcol)
# arrows(21,51,50,25, length = arrlen, lwd = arrwd, col = arcol, lty =1)

#### TMP optimal fishing , long-term
### Parameters
ylims = c(0, 100)
ytmpticks = seq(ylims[1], ylims[2], by = 20)
ytmplabs = ytmpticks


##Figure

plot(datn$constarg, datn$statsave_bopt_long, xlim = xlims, ylim = alims, pch = 17, cex = cexpts, col = statcol,
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)
axis(side = 2, labels = ytmplabs, at = ytmpticks, lwd= axwd, cex.axis = cexax)


lines(datn$constarg, datn$statsave_bopt_long, col = statcol, lwd = lnwd)

points(datn$constarg, datn$dynsave_bopt_long, pch = 22, cex = cexpts, col = dyncol, bg ='transparent')
lines(datn$constarg, datn$dynsave_bopt_long, col = dyncol, lwd = lnwd)

text(-2, 110, 'C', cex = cexlab)

# text(20,70,'Larger reserves', cex = cexlab, srt = 320, pos = 4, col = arcol)
# arrows(21,61,50,34, length = arrlen, lwd = arrwd, col = arcol, lty =1)

#### TMP optimal fishing , short-term
### Parameters
ylims = c(0, 100)
ytmpticks = seq(ylims[1], ylims[2], by = 20)
ytmplabs = ytmpticks


##Figure

plot(datn$constarg, datn$statsave_bopt_short, xlim = xlims, ylim = alims, pch = 17, cex = cexpts, col = statcol,
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)
axis(side = 2, labels = ytmplabs, at = ytmpticks, lwd= axwd, cex.axis = cexax)


lines(datn$constarg, datn$statsave_bopt_short, col = statcol, lwd = lnwd)

points(datn$constarg, datn$dynsave_bopt_short, pch = 22, cex = cexpts, col = dyncol, bg ='transparent')
lines(datn$constarg, datn$dynsave_bopt_short, col = dyncol, lwd = lnwd)

text(-2, 110, 'D', cex = cexlab)

# text(20,70,'Larger reserves', cex = cexlab, srt = 320, pos = 4, col = arcol)
# arrows(21,61,50,34, length = arrlen, lwd = arrwd, col = arcol, lty =1)