#####################################################################################
###  PROGRAM FOR GRAPHING PRED-PREY SIMULATION RESULTS
###         by: Andre Buchheister
###             Created: 2-10-13
###             Updated: 8-29-13
#####################################################################################

# DESCRIPTION
#   This file is for graphing results from the predator prey simulation, PredPreySim_ce_9-4-13.r

# NOTES
#  3/19/13  - Modifying this graphing code to include the STEP scenarios
#  8/29/13  - Updating code to generate nicer figures for the manuscript.
#  9/30/13  - Adding plots for % increase in recruitment.
# 10/5/13   - Made the SYR plots into a 3 panel plot.  For individual plots, (e.g. for seminar graph, i 
#             may need to go to version 9-30-13 of this file.
# 11/23/13  - Slight modifications to make plots for defense powerpoint presentation
#  7/27/14  - When running 1000 simulations, i had to break the files up into 3 different ones.  
#           - So, here i will need to merge those files together to then graph things.


library(reshape2)    #for melt()
library(gdata)        #for drop.levels()
library(Hmisc)        #for minor.tick()
library(plyr)         #for ddply()
library(abind)        #for abind()



### LOADING RESULTS AND DATA FILES
  
    #Importing LW and age data for flounder from ChesMMAP
      LWdata=read.csv("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\Data\\Flounder_LW_Age_6-11-12.csv")
        
  #LOAD DATA FOR SIMULATION MODEL - don't run this if plotting sensitivity results    
    #importing model results
    #OLD: load(file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_9-30-13.RData")
    #Load the 3 files that have the different sets of simulations
    load(file="C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_600_7-25-14.RData")
      SIM.MOD.1 = SIM.MOD
    load(file="C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_200a_7-25-14.RData")
      SIM.MOD.2 = SIM.MOD
    load(file="C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_200b_7-25-14.RData")
      SIM.MOD.3 = SIM.MOD

    summary(SIM.MOD.1); summary(SIM.MOD.2); summary(SIM.MOD.3)


summary(SIM.MOD); ls()
    
    #Check dimensions of files to bind them by simulation run.
    dim(SIM.MOD.1$S.all)
    dim(SIM.MOD.1$S.diff.mean)
    dim(SIM.MOD.1$S.diff.p.mean)
    dim(SIM.MOD.1$Y.all)
    dim(SIM.MOD.1$Y.diff.mean)
    dim(SIM.MOD.1$Y.diff.p.mean)
    dim(SIM.MOD.1$R.all)
    dim(SIM.MOD.1$R.diff.mean)
    dim(SIM.MOD.1$R.diff.p.mean)
    dim(SIM.MOD.1$w.all)
    dim(SIM.MOD.1$w.diff.mean)
    dim(SIM.MOD.1$w.diff.p.mean)
    dim(SIM.MOD.1$time.target.df)
    dim(SIM.MOD.1$pulses.all)
    dim(SIM.MOD.1$function.call$scen.names)


    #Bind, and rename the various arrays to facilitate coding.
      S.all = abind(SIM.MOD.1$S.all, SIM.MOD.2$S.all, SIM.MOD.3$S.all, along=2)
      S.diff.mean = abind(SIM.MOD.1$S.diff.mean, SIM.MOD.2$S.diff.mean, SIM.MOD.3$S.diff.mean, along=1)
      S.diff.p.mean = abind(SIM.MOD.1$S.diff.p.mean, SIM.MOD.2$S.diff.p.mean, SIM.MOD.3$S.diff.p.mean, along=1)
      Y.all = abind(SIM.MOD.1$Y.all, SIM.MOD.2$Y.all, SIM.MOD.3$Y.all, along=2)
      Y.diff.mean = abind(SIM.MOD.1$Y.diff.mean, SIM.MOD.2$Y.diff.mean, SIM.MOD.3$Y.diff.mean, along=1)
      Y.diff.p.mean = abind(SIM.MOD.1$Y.diff.p.mean, SIM.MOD.2$Y.diff.p.mean, SIM.MOD.3$Y.diff.p.mean, along=1)
      R.all = abind(SIM.MOD.1$R.all, SIM.MOD.2$R.all, SIM.MOD.3$R.all, along=2)
      R.diff.mean = abind(SIM.MOD.1$R.diff.mean, SIM.MOD.2$R.diff.mean, SIM.MOD.3$R.diff.mean, along=1)
      R.diff.p.mean = abind(SIM.MOD.1$R.diff.p.mean, SIM.MOD.2$R.diff.p.mean, SIM.MOD.3$R.diff.p.mean, along=1)
      w.all = abind(SIM.MOD.1$w.all, SIM.MOD.2$w.all, SIM.MOD.3$w.all, along=2)
        dimnames(w.all)[[2]] = as.character(c(1:1000)) #rename sims from 1:1000 to avoid problems with plotting
      w.diff.mean = abind(SIM.MOD.1$w.diff.mean, SIM.MOD.2$w.diff.mean, SIM.MOD.3$w.diff.mean, along=1)
        dimnames(w.diff.mean)[[1]] = paste("sim",c(1:1000),sep="") #rename sims from 1:1000 to avoid problems with plotting
      w.diff.p.mean = abind(SIM.MOD.1$w.diff.p.mean, SIM.MOD.2$w.diff.p.mean, SIM.MOD.3$w.diff.p.mean, along=1)
        dimnames(w.diff.p.mean)[[1]] = paste("sim",c(1:1000),sep="") #rename sims from 1:1000 to avoid problems with plotting
      time.target.df = rbind(SIM.MOD.1$time.target.df, SIM.MOD.2$time.target.df, SIM.MOD.3$time.target.df)

      pulses.all = abind(SIM.MOD.1$pulses.all, SIM.MOD.2$pulses.all, SIM.MOD.3$pulses.all, along=1)
      
      scen.names = SIM.MOD.1$function.call$scen.names #these are same for all files
  
        #Check Dimensions of the binded arrays
        dim(S.all)
        dim(S.diff.mean)
        dim(S.diff.p.mean)
        dim(Y.all)
        dim(Y.diff.mean)
        dim(Y.diff.p.mean)
        dim(R.all)
        dim(R.diff.mean)
        dim(R.diff.p.mean)
        dim(w.all)
        dim(w.diff.mean)
        dim(w.diff.p.mean)
        dim(time.target.df)
        dim(pulses.all)
        dim(scen.names)
    
      

#Renaming arrays to facilitate coding
      if(FALSE) { #Use this if all the simulations are in one file
        S.all = SIM.MOD$S.all
        S.diff.mean = SIM.MOD$S.diff.mean
        S.diff.p.mean = SIM.MOD$S.diff.p.mean
        Y.all = SIM.MOD$Y.all
        Y.diff.mean = SIM.MOD$Y.diff.mean
        Y.diff.p.mean = SIM.MOD$Y.diff.p.mean
        R.all = SIM.MOD$R.all
        R.diff.mean = SIM.MOD$R.diff.mean
        R.diff.p.mean = SIM.MOD$R.diff.p.mean
        w.all = SIM.MOD$w.all
        w.diff.mean = SIM.MOD$w.diff.mean
        w.diff.p.mean = SIM.MOD$w.diff.p.mean
        time.target.df = SIM.MOD$time.target.df
        pulses.all = SIM.MOD$pulses.all
        scen.names = SIM.MOD$function.call$scen.names
      }
 
  #LOAD DATA FOR MONTE CARLO SENSITIVITY RESULTS - only run this when plotting sensitivity results    
    if(FALSE) {#RUN THIS ONLY IF INTERESTED IN PLOTTING MC SENSITIVITY RESULTS
    #Import SENSITIVITY results (Monte Carlo Method) - S.sens[sims,scenarios,sens.run], Y.sens, w.sens[sim,age,scen,sens], t.sens[sim,scen,PvNP,sens]
      #Sensitivity results run in 3 parts (with different seeds) due to memory issues.  
      #Part 1 (100 sims)
      load("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_MC_100a_seed500_10-5-13.RData")
      ls()
      R.sens.1=R.sens; S.sens.1=S.sens; t.sens.1=t.sens; w.sens.1=w.sens; Y.sens.1=Y.sens
        
      #Part 2 (200 sims)
      load("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_MC_200a_seed200_7-30-14.RData")
      R.sens.2=R.sens; S.sens.2=S.sens; t.sens.2=t.sens; w.sens.2=w.sens; Y.sens.2=Y.sens
      
      #Part 3 (700 sims)
      load("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_MC_700a_seed301_7-30-14.RData")
      R.sens.3=R.sens; S.sens.3=S.sens; t.sens.3=t.sens; w.sens.3=w.sens; Y.sens.3=Y.sens

      #Merge the 3 parts together
      dim(R.sens); dim(S.sens); dim(t.sens); dim(w.sens); dim(Y.sens)
      R.sens = abind(R.sens.1, R.sens.2, R.sens.3, along=3); dim(R.sens)
      S.sens = abind(S.sens.1, S.sens.2, S.sens.3, along=3); dim(S.sens)
      t.sens = abind(t.sens.1, t.sens.2, t.sens.3, along=4); dim(t.sens)
      w.sens = abind(w.sens.1, w.sens.2, w.sens.3, along=4); dim(w.sens)
      Y.sens = abind(Y.sens.1, Y.sens.2, Y.sens.3, along=3); dim(Y.sens)
      
      #Reformat data to be the same as SIM.MOD outputs (to utilize plotting functions for normal simulation analysis)
        S.diff.p.mean = t(S.sens[1,,])  #Pick the first of the 2 simulation runs, and treat the sensitivity runs as the rows (to match SIM.MOD$S.diff.p.mean)
        Y.diff.p.mean = t(Y.sens[1,,])
        R.diff.p.mean = t(R.sens[1,,])
        w.diff.p.mean = aperm(w.sens[1,,,],c(3,1,2))
        
        #Create Junk placeholders
        S.diff.mean = S.diff.p.mean;  Y.diff.mean = Y.diff.p.mean; R.diff.mean = R.diff.p.mean; w.diff.mean = w.diff.p.mean;

        #importing function.call information from regular simulation results
        load(file="C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_200a_7-25-14.RData")
        summary(SIM.MOD); ls()
        scen.names = SIM.MOD$function.call$scen.names
     }
  
      
  #FORMAT DATA              
    #Melt arrays into data frames to facilitate plotting
      S.diff.p.mean.df = melt(S.diff.p.mean, varnames=c("Sim","Scenario"), value.name="S.diff.p.mean")
      Y.diff.p.mean.df = melt(Y.diff.p.mean, varnames=c("Sim","Scenario"), value.name="Y.diff.p.mean")
      R.diff.p.mean.df = melt(R.diff.p.mean, varnames=c("Sim","Scenario"), value.name="R.diff.p.mean")
      w.diff.p.mean.df = melt(w.diff.p.mean, varnames=c("Sim","Age","Scenario"), value.name="w.diff.p.mean")  #melt array to facilitate plotting by group

      S.diff.mean.df = melt(S.diff.mean, varnames=c("Sim","Scenario"), value.name="S.diff.mean")
      Y.diff.mean.df = melt(Y.diff.mean, varnames=c("Sim","Scenario"), value.name="Y.diff.mean")
      R.diff.mean.df = melt(R.diff.mean, varnames=c("Sim","Scenario"), value.name="R.diff.mean")
      w.diff.mean.df = melt(w.diff.mean, varnames=c("Sim","Age","Scenario"), value.name="w.diff.mean")  #melt array to facilitate plotting by group

    #Merge data.frames together
      SY.diff.df = merge(S.diff.mean.df, S.diff.p.mean.df)
      SY.diff.df = merge(SY.diff.df, Y.diff.mean.df)
      SY.diff.df = merge(SY.diff.df, Y.diff.p.mean.df);    
      SY.diff.df = merge(SY.diff.df, R.diff.mean.df)
      SY.diff.df = merge(SY.diff.df, R.diff.p.mean.df);    head(SY.diff.df)
      
      w.diff.df = merge(w.diff.mean.df, w.diff.p.mean.df); head(w.diff.df)
      
    #Bring in ancillary descriptive scenario data from scen.names
      SY.diff.df  = merge(SY.diff.df, scen.names, by.x="Scenario", by.y="scen.name"); head(SY.diff.df)
      w.diff.df   = merge(w.diff.df,  scen.names, by.x="Scenario", by.y="scen.name"); head(w.diff.df)

    #Re-sort dataframes to facilitate plotting of factor levels later
      SY.diff.df$F.scen = factor(SY.diff.df$F.scen, levels=c("highF","decrF","dec2F","lowF","mor.F"))  #Set F.scen level order
      SY.diff.df = SY.diff.df[order(SY.diff.df$F.scen,SY.diff.df$I.scen,SY.diff.df$P.scen),]
          SY.diff.df$Scenario = factor(SY.diff.df$Scenario, levels=SY.diff.df[!duplicated(SY.diff.df$Scenario),"Scenario"]) #re-specify factor levels so that they are in a good order
          
    gc()

### SIMULATION BOUNDS/CONDITIONS
    #Numbers of predators, ages, prey, seasons, years, and regions in simulation
    predators = SIM.MOD$function.call$predators
    ages  = SIM.MOD$function.call$ages          #represents ages 0 to 7+ to be consistent with stock assessment
    preys = SIM.MOD$function.call$preys         #prey 1 = Crustaceans (mysid), prey 2 = Small fishes (anchovy), 3 = larger fishes (spot)
    years = SIM.MOD$function.call$years
    years.2 = years                             #This is for a coding nuance... necessary when estimating the target SSB
    seasons = SIM.MOD$function.call$seasons
    regions = SIM.MOD$function.call$regions     #Region 1 is the Chesapeake Bay
    #sims = SIM.MOD$function.call$sims
    sims = SIM.MOD.1$function.call$sims + SIM.MOD.2$function.call$sims + SIM.MOD.3$function.call$sims
    n.pulse = SIM.MOD$function.call$n.pulse

### PLOTTING PARAMETERS
    #plot.mar = c(5.1,5.1,4.1,2.1)
    plot.mar = c(10.1,5.1,4.1,2.1)    #Use this when using default X-axis labels that are long
    plot.mar.w= c(5.1,5.1,4.1,2.1)    #For weight at age plots
    notch.gr=T             #Notch for boxplots
    prey.col = c("gray", "red", "blue")
    cexlab  = 1.2    #cexlab  = 1.6
    cexaxis = 1.1    #cexaxis = 1.4 
    cexmain = 1.2    #cexmain = 1.8

### DELETE FILES TO CONSERVE MEMORY
    SIM.MOD.1 = NULL; SIM.MOD.2 = NULL; SIM.MOD.3 = NULL; SIM.MOD=NULL;   
    memory.size(); gc(); memory.size()

### SAVE OR LOAD SESSION (to avoid running code above)
  library(session) #to save session
  #save.session(file="C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\PredPreySim_graphing_session_8-27-14.RSession")
  restore.session(file="C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\PredPreySim_graphing_session_8-27-14.RSession")
 

####################################################################
### RUNNING PLOTTING FUNCTIONS AND STORING GRAPHS (Need to run functions below first)
####################################################################
#  pdf("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SimFigs_7-27-14.pdf", height=7,width=7) 
  #pdf("C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SensitivityFigs_9-11-13.pdf") 
    #Note: when exporting graphs to PDF, the "windows()" statements in the plotting functions have to be commented out
    
  #WEIGHT AT AGE LINEPLOTS
      #For Pulse and STEP scenarios
      w.age.lineplot()
      par(mfrow=c(1,1))  #Reset Par to be one graph per page.

  #WEIGHT AT AGE BOXPLOTS
      #For Pulse scenarios
      w.age.boxplot(I.scens = "mig", P.scens.type = "Pulse", F.scens = "decrF")    #NOTE: F scenarios don't affect the W@age plots.
      w.age.boxplot(I.scens = "mix", P.scens.type = "Pulse", F.scens = "decrF")    #NOTE: F scenarios don't affect the W@age plots.
  
      #For Step scenarios
      #w.age.boxplot(I.scens = "mig", P.scens.type = "Step", F.scens = "decrF")    #NOTE: F scenarios don't affect the W@age plots.

  #TIME TO TARGET PLOTS
      #DecrF, mig scenarios, both types of prey scenarios:
      t.plots(F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse")  
      t.prob( F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse")
      #t.plots.sens(F.scens = "decrF", I.scens = "mig", P.scens.type = "Both")  #Plot of MC sensitivity results  
      par(mfrow=c(1,1))  #Reset Par to be one graph per page.

  #SY PLOTS
      #All Scenarios
      SY.plots()

   dev.off()


### S, Y, AND R PLOTS IN 3 PANELS
#  pdf("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SimFigs_SYR panels_7-27-14.pdf", height=7,width=5) 

  SY.plots(Resp.Type="PulseYears")
  SY.plots(Resp.Type="AllYears")
  
  dev.off()
  
  
  
##################################################################################
### STORE PLOTS FOR POWERPOINT PRESENTATION (BIGGER FONTS, ETC).  11/23/13
##################################################################################
  pdf("C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\PPT_SimFigs_11-23-13.pdf", height=7,width=7) 
    #Note: when exporting graphs to PDF, the "windows()" statements in the plotting functions have to be commented out
    
  #WEIGHT AT AGE LINEPLOTS
      #For Pulse and STEP scenarios
      w.age.lineplot()
      par(mfrow=c(1,1))  #Reset Par to be one graph per page.

  #WEIGHT AT AGE BOXPLOTS
      #For Pulse scenarios
      w.age.boxplot(I.scens = "mig", P.scens.type = "Pulse", F.scens = "decrF")    #NOTE: F scenarios don't affect the W@age plots.
      #w.age.boxplot(I.scens = "mix", P.scens.type = "Pulse", F.scens = "decrF")    #NOTE: F scenarios don't affect the W@age plots.
  
      #For Step scenarios
      #w.age.boxplot(I.scens = "mig", P.scens.type = "Step", F.scens = "decrF")    #NOTE: F scenarios don't affect the W@age plots.

  #TIME TO TARGET PLOTS
      #DecrF, mig scenarios, both types of prey scenarios:
      #t.plots(F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse")  
      t.prob( F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse", LABEL=F)
      par(mfrow=c(1,1))  #Reset Par to be one graph per page.

  #SY PLOTS
      #Subset of Scenarios for presentation
      SY.plots.ppt()

   dev.off()





####################################################################
### SENSITIVITY PLOTS (Need to run functions below first)
####################################################################
  #pdf("C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SimFigs_9-26-13.pdf", height=7,width=7) 
  #pdf("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SensitivityFigs_MC_8-4-14.pdf", height=7,width=7) 
    #Note: when exporting graphs to PDF, the "windows()" statements in the plotting functions have to be commented out
    
  #WEIGHT AT AGE BOXPLOTS
      #For Pulse scenarios
      w.age.boxplot(I.scens = "mig", P.scens.type = "Pulse", F.scens = "decrF")    #NOTE: F scenarios don't affect the W@age plots.
      w.age.boxplot(I.scens = "mix", P.scens.type = "Pulse", F.scens = "decrF")    #NOTE: F scenarios don't affect the W@age plots.
  
  #TIME TO TARGET PLOTS
      #DecrF, mig scenarios, both types of prey scenarios:
      t.plots.sens(F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse")  #Plot of MC sensitivity results  
      #t.prob( F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse")
      par(mfrow=c(1,1))  #Reset Par to be one graph per page.

  #SY PLOTS
      #All Scenarios
      SY.plots()

   dev.off()


### S, Y, AND R PLOTS IN 3 PANELS
  #pdf("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SensitivityFigs_SYR panels_MC_8-4-14.pdf", height=7,width=5) 

  SY.plots(Resp.Type="PulseYears")
  SY.plots(Resp.Type="AllYears")
  
  dev.off()

### TIME TO TARGET PLOTS
  #pdf("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SensitivityFigs_time to target_MC_8-4-14.pdf", height=7,width=7) 
      #DecrF, mig scenarios, both types of prey scenarios:
      t.plots.sens(F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse")  #Plot of MC sensitivity results  
      #t.prob( F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse")
  dev.off()

   
##########################################   
### FUNCTION FOR SSB, Y, and R PLOTS
##########################################
  head(SY.diff.df)
  
  ## Function for plotting %Change in SSB and Yield.  Specify entries (can be multiples) for:
  #    I.scen ("mig" or "mix"), P.scens.type ("Pulse","Step", or "Both"), F.scens ("highF", "lowF", "decrF")  
  
SY.plots = function(I.scens = c("mig","mix"), P.scens.type = "Pulse", F.scens = c("highF","lowF","decrF"), 
                    Resp.Type="PulseYears", Title="All Scen. Comp:") {
      #I.scens = c("mig"); P.scens.type = "Pulse"; F.scens = c("highF","decrF","lowF")
      
      if(P.scens.type == "Pulse") P.scens = c("Pulse1","Pulse2","Pulse3")
      if(P.scens.type == "Step")  P.scens = c("Step.1","Step.2","Step.3")
      if(P.scens.type == "Both")  P.scens = c("Pulse1","Pulse2","Pulse3", "Step.1","Step.2","Step.3")
      if(P.scens.type == "All")   P.scens = unique(scen.names$P.scen)

  #GENERAL PLOTS - NEW ORDER
    #Create subset of data based on specifications
    SY.sub.0 = SY.diff.df[SY.diff.df$I.scen %in% I.scens & SY.diff.df$P.scen %in% P.scens & SY.diff.df$F.scen %in% F.scens,];   #(Just change SY.sub.0 above, and change the Y variable and title)
        SY.sub.0 = SY.sub.0[order(SY.sub.0$I.scen,SY.sub.0$P.scen,SY.sub.0$F.scen),]  #Use this for I>P>F ordering      
        SY.sub.0$Scenario = factor(SY.sub.0$Scenario, levels=SY.sub.0[!duplicated(SY.sub.0$Scenario),"Scenario"])


    #MANUSCRIPT FIGURES (9-25-13)
      if(Resp.Type=="PulseYears"){
        SY.var = c("S.diff.p.mean","Y.diff.p.mean","R.diff.p.mean")   #Variables to be plotted
        SY.ylab= c(expression(paste(Delta,"S (%, for pulse years)")),expression(paste(Delta,"Y (%, for pulse years)")),expression(paste(Delta,"R (%, for pulse years)")) ) #Y labels to be used for each plot
        }
      if(Resp.Type=="AllYears"){
        SY.var = c("S.diff.mean","Y.diff.mean","R.diff.mean")   #Variables to be plotted
        SY.ylab= c(expression(paste(Delta,"S (%, for all years)")), expression(paste(Delta,"Y (%, for all years)")), expression(paste(Delta,"R (%, for all years)"))  )   #Y labels to be used for each plot
         }         
                 
        par(mfrow=c(3,1))
        #par(mar=c(9.1,5.1,2,2))
        par(mar=c(1,5,1,1))
        par(oma=c(8,1,1,1))
        
      for(i in 1:3) {
        SY.plot.form = as.formula(paste(SY.var[i],"~Scenario"))
        Panel = c("A.", "B.", "C.")
        #windows()
        boxplot(SY.plot.form, data=SY.sub.0, main="", 
            ylab=SY.ylab[i] , las=2, col=rep(prey.col,each=length(F.scens)), notch=notch.gr, # "Mean % diff. in YIELD (from No-pulse)"
            xlab="", cex.lab=cexlab, 
            xaxt="n"  )  #, cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis
            abline(v=c(3.5,6.5,9.5,12.5,15.5),lty=3,lwd=0.5)
            abline(v=9.5)
            minor.tick(ny=5,nx=0)
            axis(1,at=c(1:18),labels=F, las=1)    #c("High", "Decr", "Low")  cex.axis=cexaxis, 
            mtext(Panel[i],line=-2,at=0.5)
            
            if(i==1){  #Insert Legend in top figure
            legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
                  col=prey.col,pch=15,bg="white",pt.cex=1.5)  # , cex=cexlab
                     }
            if(i==3){   #Insert axes labels for bottom figure only.
            axis(1,at=c(1:18),labels=rep(c("H", "D", "L"),6), las=1)    #c("High", "Decr", "Low")  cex.axis=cexaxis, 
            axis(1,at=c(2,5,8,11,14,17),tick=F,line=2, labels=rep(c("Prey 1","Prey 2","Prey 3"),2))
            axis(1,at=c(5,14),tick=F,line=4, labels=c("Mig","Mix"))
            axis(1,at=c(9.5),tick=F,line=6, labels="Scenario", cex.axis=cexlab)
            
            axis(1,at=0,labels=expression(italic("Fishing                    ")), tick=F, line=0, padj=0)  #Labels for the type of scenario    
            axis(1,at=0,labels=expression(italic("Pulse                      ")), tick=F, line=2, padj=0)     
            axis(1,at=0,labels=expression(italic("Migration                   ")), tick=F, line=4, padj=0)     
                  }
            } #END PLOT LOOP

  #Summarize medians by scenario
  SY.summary = ddply(SY.sub.0,.(Scenario,P.scen,I.scen,F.scen),summarise, 
    S.diff.p.med=median(S.diff.p.mean,na.rm=T), 
    S.diff.med=median(S.diff.mean,na.rm=T), 
    Y.diff.p.med=median(Y.diff.p.mean,na.rm=T), 
    Y.diff.med=median(Y.diff.mean,na.rm=T), 
    R.diff.p.med=median(R.diff.p.mean,na.rm=T), 
    R.diff.med=median(R.diff.mean,na.rm=T)) 
  #SY.summary = SY.summary[order(SY.summary$F.scen,SY.summary$P.scen,SY.summary$I.scen),]  #Order as needed: e.g. Order so that the mig/mix are next to one another
   SY.summary

                 }  #End SY.plots()

SY.plots()


### SYR PLOTS FOR PRESENTATION (11-23-13)
SY.plots.ppt = function(I.scens = c("mig"), P.scens.type = "Pulse", F.scens = c("highF","lowF"), 
                    Resp.Type="PulseYears", Title="All Scen. Comp:") {
      #I.scens = c("mig"); P.scens.type = "Pulse"; F.scens = c("highF","decrF","lowF")
      
      if(P.scens.type == "Pulse") P.scens = c("Pulse1","Pulse2","Pulse3")
      if(P.scens.type == "Step")  P.scens = c("Step.1","Step.2","Step.3")
      if(P.scens.type == "Both")  P.scens = c("Pulse1","Pulse2","Pulse3", "Step.1","Step.2","Step.3")
      if(P.scens.type == "All")   P.scens = unique(scen.names$P.scen)

  #GENERAL PLOTS - NEW ORDER
    #Create subset of data based on specifications
    SY.sub.0 = SY.diff.df[SY.diff.df$I.scen %in% I.scens & SY.diff.df$P.scen %in% P.scens & SY.diff.df$F.scen %in% F.scens,];   #(Just change SY.sub.0 above, and change the Y variable and title)
        #SY.sub.0 = SY.sub.0[order(SY.sub.0$I.scen,SY.sub.0$P.scen,SY.sub.0$F.scen),]  #Use this for I>P>F ordering      
        SY.sub.0 = SY.sub.0[order(SY.sub.0$I.scen,SY.sub.0$F.scen,SY.sub.0$P.scen),]  #Use this for I>F>P ordering      
        SY.sub.0$Scenario = factor(SY.sub.0$Scenario, levels=SY.sub.0[!duplicated(SY.sub.0$Scenario),"Scenario"])


    #DEFENSE PRESENTATION PPT FIGURES (11-23-13)
      if(Resp.Type=="PulseYears"){
        SY.var = c("S.diff.p.mean","Y.diff.p.mean","R.diff.p.mean")   #Variables to be plotted
        SY.ylab= c(expression(paste(Delta,"S (%, for pulse years)")),expression(paste(Delta,"Y (%, for pulse years)")),expression(paste(Delta,"R (%, for pulse years)")) ) #Y labels to be used for each plot
        }
      if(Resp.Type=="AllYears"){
        SY.var = c("S.diff.mean","Y.diff.mean","R.diff.mean")   #Variables to be plotted
        SY.ylab= c(expression(paste(Delta,"S (%, for all years)")), expression(paste(Delta,"Y (%, for all years)")), expression(paste(Delta,"R (%, for all years)"))  )   #Y labels to be used for each plot
         }         
                 
        #par(mfrow=c(1,1))
        #par(mar=c(4,5,1,1))
        #par(oma=c(1,1,1,1))
        par(mar=plot.mar.w)
        
      for(i in 1:3) {
        SY.plot.form = as.formula(paste(SY.var[i],"~Scenario"))
        Panel = c("A.", "B.", "C.")
        #windows()
        boxplot(SY.plot.form, data=SY.sub.0, main="", 
            ylab=SY.ylab[i] , las=2, col=rep(prey.col,length(F.scens)), notch=notch.gr, #  
            xlab="Scenario", cex.lab=cexlab, cex.axis=cexaxis,
            xaxt="n"  )  #, cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis
            abline(v=c(3.5),lty=3,lwd=0.5)
            #abline(v=9.5)
            minor.tick(ny=2,nx=0)
            axis(1,at=c(2,5),labels=c("High F", "Low F"), las=1, cex.axis=cexaxis )    #c("High", "Decr", "Low")  cex.axis=cexaxis, 
            #mtext(Panel[i],line=-2,at=0.5)
            
            #if(i==1){  #Insert Legend in top figure
            legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
                  col=prey.col,pch=15,bg="white",pt.cex=1.5, cex=cexaxis)  # , cex=cexlab
                    # }
            } #END PLOT LOOP

  #Summarize medians by scenario
  SY.summary = ddply(SY.sub.0,.(Scenario,P.scen,I.scen,F.scen),summarise, 
    S.diff.p.med=median(S.diff.p.mean,na.rm=T), 
    S.diff.med=median(S.diff.mean,na.rm=T), 
    Y.diff.p.med=median(Y.diff.p.mean,na.rm=T), 
    Y.diff.med=median(Y.diff.mean,na.rm=T), 
    R.diff.p.med=median(R.diff.p.mean,na.rm=T), 
    R.diff.med=median(R.diff.mean,na.rm=T)) 
  #SY.summary = SY.summary[order(SY.summary$F.scen,SY.summary$P.scen,SY.summary$I.scen),]  #Order as needed: e.g. Order so that the mig/mix are next to one another
   SY.summary

                 }  #End SY.plots.ppt()

SY.plots.ppt()



##################################
## WEIGHT-AT-AGE BOXPLOT FUNCTION
##################################
  ## Function for plotting %Change in weight at age.  Specify single entries for:
  #    I.scen ("mig" or "mix"), P.scens.type ("Pulse" or "Step"), F.scens ("highF", "lowF", "decrF")  
  #    Note that the F scenario does not affect the plots.
  
  
  w.age.boxplot = function(I.scens = "mig", P.scens.type = "Pulse", F.scens = "highF") {

      if(P.scens.type == "Pulse") P.scens = c("Pulse1","Pulse2","Pulse3")
      if(P.scens.type == "PulsB") P.scens = c("PulsB1","PulsB2","PulsB3")
      if(P.scens.type == "Step")  P.scens = c("Step.1","Step.2","Step.3")
      
    #Create subset of data
    w.sub.0 = w.diff.df[w.diff.df$I.scen %in% I.scens & w.diff.df$P.scen %in% P.scens & w.diff.df$F.scen %in% F.scens,];
        w.sub.0$Scenario = factor(w.sub.0$Scenario, levels=w.sub.0[!duplicated(w.sub.0$Scenario),"Scenario"])
        head(w.sub.0)
        
    #Summarize medians by scenario
    w.summary = ddply(w.sub.0,.(Scenario,Age,P.scen,I.scen,F.scen),summarise, 
      w.diff.p.med=median(w.diff.p.mean,na.rm=T), 
      w.diff.med=median(w.diff.mean,na.rm=T)) 
    #w.summary = w.summary[order(w.summary$Age,w.summary$P.scen,w.summary$F.scen,w.summary$I.scen),]  #Order as needed: e.g. Order so that the mig/mix are next to one another
     w.summary

  
      #windows()
      par(mar=plot.mar.w)
  
    #Weight at age BOXPLOT - ALL years - F scenarios have no impact on weight at age
      boxplot(w.diff.mean~Scenario*Age, data=w.sub.0, main=paste("Weight-at-age change (all years, ",I.scens,", ",F.scens,")",sep=""), cex.main=cexmain,
            ylab=expression(paste(Delta,"w (%, for all years)")), col=prey.col, las=2, 
            xlab="Age", xaxt="n", cex.lab=cexlab, cex.axis=cexaxis, notch=notch.gr )  #, cex.main=cexmain, cex.axis=cexaxis
              for(i in 0:7) {                        #make "bars" to help identify the ages.
                lines(c(1:3)+i*3,rep(0,3), lwd=3)  }
            minor.tick(ny=2,nx=0)
            axis(1,at=seq(2,24, by=3), labels=c(0:(ages-1)), las=1, cex.axis=cexaxis)   #Make new x-axis
            legend("topright",c("Prey 1","Prey 2","Prey 3"), title=paste(substr(P.scens[1],1,5),"s in:",sep=""),
                  col=prey.col,pch=15,bg="white",pt.cex=1.5, cex=cexaxis)   # , cex=cexlab
  
      #windows()
      par(mar=plot.mar.w)

    #Weight at age BOXPLOT - PULSE years - F scenarios have no impact on weight at age
      boxplot(w.diff.p.mean~Scenario*Age, data=w.sub.0, main=paste("Weight-at-age change (pulse years, ",I.scens,", ",F.scens,")",sep=""), cex.main=cexmain,
            ylab=expression(paste(Delta,"w (%, for pulse years)")), col=prey.col, las=2, 
            xlab="Age", xaxt="n", cex.lab=cexlab, cex.axis=cexaxis, notch=notch.gr )  #, cex.main=cexmain, cex.axis=cexaxis, 
              for(i in 0:7) {                        #make "bars" to help identify the ages.
                lines(c(1:3)+i*3,rep(-0.5,3), lwd=3)  }
            minor.tick(ny=5,nx=0)
            axis(1,at=seq(2,24, by=3), labels=c(0:(ages-1)), las=1, cex.axis=cexaxis)   #Make new x-axis
            legend("topright",c("Prey 1","Prey 2","Prey 3"), title=paste(substr(P.scens[1],1,5),"s in:",sep=""),
                  col=prey.col,pch=15,pt.cex=1.5,bg="white", cex=cexaxis)  #, cex=cexlab
           
           }   #end w.age.boxplot

   w.age.boxplot(I.scens = "mig", P.scens.type = "Pulse", F.scens = "highF")


##################################
## WEIGHT-AT-AGE LINE PLOT FUNCTION
##################################


#w.all[scen,sim,yr,seas,age, NP v Pulse]
  #w.all[1,85,40:50,,,2];  head(w.all)  dim(w.all)  dimnames(w.all)

  #Melt weight at age data into data frame (exclude NoPulse runs, because they are identical)
      w.all.df = melt(w.all[,,26:60,,,2], varnames=c("Scenario", "Sim","Year","Season","Age"), value.name="Weight")
      head(w.all.df); tail(w.all.df) 
      w.all.df[1:100,]

    #Identify fish by their cohorts (ie year in which they were born)
      w.all.df$Age = as.numeric(substr(w.all.df$Age,4,4))
      w.all.df$Age.2 = w.all.df$Age+0.25*as.numeric(substr(w.all.df$Season,2,2))
      w.all.df$Cohort= w.all.df$Year - w.all.df$Age
      w.all.df$P.scen = substr(w.all.df$Scenario,1,6) 
           
    #Subset data for migration and highF scenarios
      w.all.df = w.all.df[substr(w.all.df$Scenario,8,13)=="mig-hi",]
      summary(w.all.df)
       
  #Calculate MEAN cohort trajectories across simulations and cohorts (for years 26+)
    #Melt weight at age data into data frame (exclude NoPulse runs, because they are identical)
        w.all.mean = apply(w.all,c(1,3,4,5,6),mean)  # w.all.mean[scen,yr,seas,age, NP v Pulse]
        
        w.all.df.mean = melt(w.all.mean[,,,,], varnames=c("Scenario","Year","Season","Age","NPvP"), value.name="Weight")
        head(w.all.df.mean)
  
      #Identify fish by their cohorts (ie year in which they were born)
        w.all.df.mean$Age = as.numeric(substr(w.all.df.mean$Age,4,4))
        w.all.df.mean$Age.2 = w.all.df.mean$Age+0.25*as.numeric(substr(w.all.df.mean$Season,2,2))
        w.all.df.mean$Cohort= w.all.df.mean$Year - w.all.df.mean$Age
        w.all.df.mean$P.scen = substr(w.all.df.mean$Scenario,1,6) 
             
      #Subset data for migration and highF scenarios and cohorts that were born after year 25
        w.all.df.mean.sub = w.all.df.mean[substr(w.all.df.mean$Scenario,8,13)=="mig-hi" & w.all.df.mean$Cohort>25,]
        summary(w.all.df.mean.sub)
      
      #Calculate means
        w.mean = aggregate(Weight~Scenario + NPvP + Age.2 + P.scen, data=w.all.df.mean.sub, mean)
        head(w.mean)  
        
      xlim.1 = c(0,8)  #c(0,8)
      ylim.1 = c(0,2.5)  #c(0,2.5)
      plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="gray", xlim=xlim.1, ylim=ylim.1, xlab="Age", ylab="Weight (kg)", 
          main="Mean weight at age across pulse simulations"); minor.tick(ny=2, nx=2, tick.ratio=0.5)  
          lines(w.mean[w.mean$NPvP=="Pulse" & w.mean$P.scen=="Pulse1",c("Age.2","Weight")],col="black", lwd=2)
          lines(w.mean[w.mean$NPvP=="Pulse" & w.mean$P.scen=="Pulse2",c("Age.2","Weight")],col="red", lwd=2)
          lines(w.mean[w.mean$NPvP=="Pulse" & w.mean$P.scen=="Pulse3",c("Age.2","Weight")],col="blue", lwd=2)
          lines(w.mean[w.mean$NPvP=="NoPulse" & w.mean$P.scen=="Pulse1",c("Age.2","Weight")],col="black", lty=2,lwd=1)
                legend("topleft",c("Pulse 1", "Pulse 2", "Pulse 3", "No Pulse"), col=c("black","red","blue"), lty=c(1,1,1,2), lwd=c(2,2,2,1))
       
  #Trajectory of a given cohort
    xlim.1 = c(0,8)  #c(0,8)
    ylim.1 = c(0,2.5)  #c(0,2.5)
    sim.pl = 500  #Originally 85
    Title.P = "Weight at age for a pulse simulation"
    Title.S = "Weight at age for a step simulation"
    years.pl = 60  #Max year of "experimental period"

  #OLD: Data for fastest growing cohorts:  a plot of mean No Pulse runs with the single cohorts from each pulse scenario that had the largest growth.
    if(FALSE){
     head(w.all.df)
     temp = aggregate(Weight~Scenario+Age.2+P.scen,data=w.all.df,max)  #Calculate max weights at each age by scenario
     temp.2 = temp[(temp$P.scen=="Pulse1" & temp$Age.2==0.75) | (temp$P.scen=="Pulse2" & temp$Age.2==1.75) | (temp$P.scen=="Pulse3" & temp$Age.2==3.75),] #subset for the max weights that occur at the age of peak pulse effects      
     w.all.df.max = w.all.df[w.all.df$Weight==temp.2[,"Weight"],]   #Identify the Sim and Cohort exhibiting the maximum weight at age.
     w.all.df.max
    }

  #Identify fastest growing cohort from all simulations
    #Calculate the residual sum of squares of each sim relative to the NoPulse base model    
    w.base = w.mean[w.mean$NPvP=="NoPulse" & w.mean$P.scen=="Pulse1",c("Age.2","Weight")] #Base model (NoPulse) weight at age
    head(w.all.df)
    w.all.df$Weight_base = w.base$Weight[match(w.all.df$Age.2,w.base$Age.2)] #"vlookup" to get base model weights at each age.
    w.all.df$Resid = w.all.df$Weight - w.all.df$Weight_base  #Calculate residual
    w.all.df$SquareResid = (w.all.df$Resid)**2
    gc()
    temp = aggregate(Weight~Scenario+P.scen+Sim+Cohort,data=w.all.df,sum)  #Calculate residual sum of squares
    head(temp)
    temp.2 = aggregate(Weight~Scenario+P.scen,data=temp,max)  #Identify max RSS
    head(temp.2)
    w.all.df.max = temp[temp$Weight %in% temp.2[,"Weight"],]   #Identify the Sim and Cohort exhibiting the maximum weight at age.
    w.all.df.max = w.all.df.max[order(w.all.df.max$P.scen),]  #Order by P.scen
    #ENTER THE IDENTIFIED COHORTS BELOW IN FUNCTION

w.age.lineplot = function(){ 
  
  #UPDATE COHORTS BELOW IF NEEDED!!!

  #FASTEST GROWING SINGLE COHORTS IN PULSE SCENARIOS:
     par(mar=plot.mar.w)
     plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="gray", xlim=xlim.1, ylim=ylim.1, xlab="Age", ylab="Weight (kg)", cex=1.5,
          cex.lab=cexlab, cex.axis=cexaxis, main="Fastest growing pulse cohorts", xaxt="n", cex.main=cexmain)
          axis(1,at=c(0:8), cex.axis=cexaxis) 
          minor.tick(ny=2, nx=0, tick.ratio=0.5)  
          lines(w.mean[w.mean$NPvP=="NoPulse" & w.mean$P.scen=="Pulse1",c("Age.2","Weight")],col="black", lty=2,lwd=2)
          lines(w.all.df[w.all.df$Sim==734 & w.all.df$Cohort==46 & w.all.df$P.scen=="Pulse1",c("Age.2","Weight")],col="black", lwd=2)
          lines(w.all.df[w.all.df$Sim==734 & w.all.df$Cohort==45 & w.all.df$P.scen=="Pulse2",c("Age.2","Weight")],col="red", lwd=2)
          lines(w.all.df[w.all.df$Sim==875 & w.all.df$Cohort==40 & w.all.df$P.scen=="Pulse3",c("Age.2","Weight")],col="blue", lwd=2)
                legend("topleft",c("Base", "Pulse 1", "Pulse 2", "Pulse 3"), col=c("black","black","red","blue"), lty=c(2,1,1,1), lwd=c(1,2,2,2), cex=cexaxis)
          #lines(x.ages,w, lwd=3)   #This is the "mean from the wild... need to run PredPreySim() with default parameters to get this to work.

  #PULSE GRAPHS - multiple lines tracking all cohorts        
    if(FALSE){
    par(mfrow=c(2,2)) #Plots of each prey pulse separate, and then all combined
      plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="gray", xlim=xlim.1, ylim=ylim.1, xlab="Age", ylab="Weight (kg)", main=Title.P); minor.tick(ny=2, nx=2, tick.ratio=0.5) 
        for(y in 25:years.pl){
          lines(w.all.df[w.all.df$Cohort==y & w.all.df$Sim==sim.pl & w.all.df$P.scen=="Pulse1",c("Age.2","Weight")],col="black")
          }
      plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="gray", xlim=xlim.1, ylim=ylim.1, xlab="Age", ylab="Weight (kg)", main=Title.P); minor.tick(ny=2, nx=2, tick.ratio=0.5)     #
        for(y in 25:years.pl){
          lines(w.all.df[w.all.df$Cohort==y & w.all.df$Sim==sim.pl & w.all.df$P.scen=="Pulse2",c("Age.2","Weight")],col="red")
          }
      plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="gray", xlim=xlim.1, ylim=ylim.1, xlab="Age", ylab="Weight (kg)", main=Title.P); minor.tick(ny=2, nx=2, tick.ratio=0.5)     #
        for(y in 25:years.pl){
          lines(w.all.df[w.all.df$Cohort==y & w.all.df$Sim==sim.pl & w.all.df$P.scen=="Pulse3",c("Age.2","Weight")],col="blue")
          }
      plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="gray", xlim=xlim.1, ylim=ylim.1, xlab="Age", ylab="Weight (kg)", main=Title.P); minor.tick(ny=2, nx=2, tick.ratio=0.5)  
        for(y in 25:years.pl){
          lines(w.all.df[w.all.df$Cohort==y & w.all.df$Sim==sim.pl & w.all.df$P.scen=="Pulse1",c("Age.2","Weight")],col="black")
          }
        for(y in 25:years.pl){
          lines(w.all.df[w.all.df$Cohort==y & w.all.df$Sim==sim.pl & w.all.df$P.scen=="Pulse2",c("Age.2","Weight")],col="red")
          }
        for(y in 25:years.pl){
          lines(w.all.df[w.all.df$Cohort==y & w.all.df$Sim==sim.pl & w.all.df$P.scen=="Pulse3",c("Age.2","Weight")],col="blue")
          }
          }

     } #End w.age.lineplot()

     
w.age.lineplot()  




##################################
## TIME TO TARGET SSB PLOTS
##################################

  #FUNCTION TO CALCULATE DIFF IN TIME TO TARGET (relative to no-pulse)
    #summary(time.target.df)
    #head(time.target.df)

  t.plots = function(F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse") {     
    #F.scens = c("decrF","exp.F","ass.F","knifF","mor.F")
    #I.scens = "mig"
    #P.scens.type = "Both"  
      if(P.scens.type == "Pulse") P.scens = c("Pulse1","Pulse2","Pulse3")
      if(P.scens.type == "Step")  P.scens = c("Step.1","Step.2","Step.3")
      if(P.scens.type == "Both") P.scens = c("Pulse1","Pulse2","Pulse3", "Step.1","Step.2","Step.3")
    
    t.diff = dcast(time.target.df, Scenario + Sim + P.scen+I.scen+F.scen~ PvNP, value.var="Time.target", mean)
    t.diff$t.diff = t.diff$NoPulse - t.diff$Pulse   #Calculate diff in time to target SSB
    t.diff$t.diff.perc = round(t.diff$t.diff/t.diff$NoPulse*100,2)   #Calculate diff in time to target SSB
    t.diff.sub = na.omit(subset(t.diff,F.scen %in% F.scens & I.scen%in% I.scens & P.scen %in% P.scens))  #Subset the dataset based on conditions
    #t.diff.sub = t.diff.sub[order(t.diff.sub$F.scen, t.diff.sub$P.scen),]   #Order the data set in the way you want the scenarios ordered for plotting    
    t.diff.sub = t.diff.sub[order(t.diff.sub$P.scen, t.diff.sub$F.scen),]   #New order (9-25-13)    
    t.diff.sub$Scenario = factor(t.diff.sub$Scenario, levels=t.diff.sub[!duplicated(t.diff.sub$Scenario),"Scenario"])   #Re-set factor levels
        levels(t.diff.sub$Scenario)
        head(t.diff.sub)
        dim(t.diff.sub)        
      #Calculate specific values that are plotted
      aggregate(data=t.diff.sub, cbind(t.diff.perc,t.diff)~Scenario, median)  #Calculate the median values shown in the figure
      aggregate(data=t.diff.sub, cbind(t.diff.perc,t.diff)~Scenario, max)     #Calculate the median values shown in the figure
      aggregate(data=t.diff.sub, cbind(t.diff.perc,t.diff)~Scenario, quantile)     #Calculate the median values shown in the figure
  
    #PLOT OF DIFFERENCES (PULSE - NOPULSE)   
      par(mar=c(7.1,5.1,4,2))
        boxplot(t.diff.perc~Scenario, data=t.diff.sub, main="Reduction in time to reach target SSB", 
            ylab=expression(paste(Delta,"t (%)")), las=2, col=rep(prey.col,each=length(F.scens)), notch=notch.gr, # "Mean % diff. in YIELD (from No-pulse)"
            xlab="", cex.lab=cexlab, 
            xaxt="n"  )  #, cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis
            abline(v=c(3.5,6.5),lty=3,lwd=0.5)
            minor.tick(ny=5,nx=0)
            axis(1,at=c(1:18),labels=rep(c("D", "D2", "M"),6), las=1)    #c("High", "Decr", "Low")  cex.axis=cexaxis, 
            axis(1,at=c(2,5,8),tick=F,line=2, labels=rep(c("Prey 1","Prey 2","Prey 3"),1))
            axis(1,at=c(5),tick=F,line=4, labels=c("Scenario"), cex.axis=cexlab)
            
            axis(1,at=0.5,labels=expression(italic("Fishing                         ")), tick=F, line=0, padj=0)  #Labels for the type of scenario    
            axis(1,at=0.5,labels=expression(italic("Pulse                             ")), tick=F, line=2, padj=0)     
            #axis(1,at=0,labels=expression(italic("Migration                   ")), tick=F, line=4, padj=0)     
            legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
                  col=prey.col,pch=15,bg="white",pt.cex=1.5)  # , cex=cexlab

      if (FALSE){
      par(mar=plot.mar)
      #windows(); par(mar=plot.mar)
      boxplot(t.diff~Scenario, data=t.diff.sub, main="Reduction in years to reach target SSB", 
          ylab="Reduction in years to reach target SSB (years) (from No-pulse)", las=2, col=prey.col, notch=T)
          minor.tick(ny=5, nx=0)

      #windows()
      par(mar=plot.mar)
      boxplot(t.diff.perc~Scenario, data=t.diff.sub, main="Reduction in time to reach target SSB", 
          ylab="Reduction in time to reach target SSB (%)", las=2, col=prey.col, notch=T)  #,
          #xaxt="n")
          minor.tick(ny=5, nx=0)
          #axis(1,at=c(1:length(P.scens)),labels=P.scens,cex.axis=cexaxis, las=2)
      }
       
    #PLOT OF TIME TO TARGET BY SCENARIO
      #Organize data 
      t.diff.sub.base = subset(t.diff.sub, Scenario == "Pulse1-mig-decrF")  #Specify one scenario to extract the "NoPulse" baseline data
      t.diff.sub.base$Pulse = t.diff.sub.base$NoPulse                       #Copy the nopulse values into the pulse column for plotting later
      t.diff.sub.base$t.diff = 0                                            #Set t.diff to 0 to avoid confusion
      t.diff.sub.base[,c("Scenario","P.scen")] = "NoPulse"                  #Rename the scenario
      t.diff.sub.2 = rbind(t.diff.sub.base,t.diff.sub)                      #combine the baseline data with the rest.
      t.diff.sub.2$Scenario = factor(t.diff.sub.2$Scenario, levels=t.diff.sub.2[!duplicated(t.diff.sub.2$Scenario),"Scenario"])   #Re-set factor levels
      levels(t.diff.sub.2$Scenario)
      head(t.diff.sub.2)

      if(FALSE){
      #Plot time to target 
        #windows()
        par(mar=plot.mar)  # c(7,5,5,2)
        #Box plot for multiple F.scenarios
        boxplot(Pulse~Scenario, data=t.diff.sub.2, main=paste("Pulse Effects on Management Goals (",I.scens,")",sep=""), 
                col=c("white",rep(c("gray","red","blue"),length(F.scens))),
                ylab="Time to reach target SSB (yr)", xlab="", las=2, cex.main=cexmain, cex.lab=cexlab, notch=notch.gr )   
                #xlab="", xaxt="n"  )
                #axis(1,at=c(1:(1+length(P.scens))),labels=c("No Pulse", P.scens),cex.axis=cexaxis, las=2)
                minor.tick(ny=5, nx=0, tick.ratio=0.5)                                                             
                abline(v=c(1.5,4.5), lty=2, col="red")
        #OLD: For one F.scenario only:
          boxplot(Pulse~P.scen, data=t.diff.sub.2, main=paste("Pulse Effects on Management Goals (",F.scens,", ",I.scens,")",sep=""), 
                  col=c("white",rep(c("gray","red","blue"),length(P.scens)/3)),
                  ylab="Time to reach target SSB (yr)", xlab="", las=2, cex.main=cexmain, cex.lab=cexlab, notch=notch.gr ,#)   
                  xlab="", xaxt="n"  )
                  axis(1,at=c(1:(1+length(P.scens))),labels=c("No Pulse", P.scens),cex.axis=cexaxis, las=2)
                  minor.tick(ny=5, nx=0, tick.ratio=0.5)                                                             
                  abline(v=c(1.5,4.5), lty=2, col="red")
                }
         } #END t.plots
   
   
   t.plots()


#################################################
## SENSITIVITY PLOT: TIME TO TARGET SSB PLOTS
#################################################

###SENSITIVITY PLOT: FUNCTION TO CALCULATE DIFF IN TIME TO TARGET (relative to no-pulse)
  t.plots.sens = function(F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse") {     
    #F.scens = c("decrF","dec2F","mor.F")
    #I.scens = "mig"
    #P.scens.type = "Both"  
      if(P.scens.type == "Pulse") P.scens = c("Pulse1","Pulse2","Pulse3")
      if(P.scens.type == "Step")  P.scens = c("Step.1","Step.2","Step.3")
      if(P.scens.type == "Both") P.scens = c("Pulse1","Pulse2","Pulse3", "Step.1","Step.2","Step.3")

      #Calculations for Sensitivity analysis plots (slightly different from t.plots())
            t.sens.df = melt(t.sens[1,,,], varnames=c("Scenario","NPvP","sensitivity.run"), value.name="t")
            t.diff = dcast(t.sens.df, Scenario + sensitivity.run~ NPvP, value.var="t", mean)
            t.diff = merge(scen.names,t.diff, by.x="scen.name", by.y="Scenario", all.y=T); colnames(t.diff)[1] = "Scenario"
            t.diff$t.diff = t.diff$NoPulse - t.diff$Pulse   #Calculate diff in time to target SSB
            t.diff$t.diff.perc = round(t.diff$t.diff/t.diff$NoPulse*100,2)   #Calculate diff in time to target SSB
            t.diff.sub = na.omit(subset(t.diff,F.scen %in% F.scens & I.scen%in% I.scens & P.scen %in% P.scens))  #Subset the dataset based on conditions
            t.diff.sub = t.diff.sub[order(t.diff.sub$P.scen, t.diff.sub$F.scen),]   #New order (9-25-13)    
            t.diff.sub$Scenario = factor(t.diff.sub$Scenario, levels=t.diff.sub[!duplicated(t.diff.sub$Scenario),"Scenario"])   #Re-set factor levels
              levels(t.diff.sub$Scenario)
              #t.diff.sub[,-c(2,6)]
            
    
    #PLOT OF DIFFERENCES (PULSE - NOPULSE)   
      par(mar=plot.mar)
      #windows(); par(mar=plot.mar)
      boxplot(t.diff~Scenario, data=t.diff.sub, main="SENSITIVITY: Reduction in years to reach target SSB", 
          ylab="Reduction in years to reach target SSB (years) (from No-pulse)", las=2, col=prey.col, notch=T)
          minor.tick(ny=5, nx=0)

      #windows()
      #par(mar=plot.mar.w)
      par(mar=c(7.1,5.1,4,2))
      boxplot(t.diff.perc~Scenario, data=t.diff.sub, main="SENSITIVITY: Reduction in time to reach target SSB", 
          ylab=expression(paste(Delta,"t (%)")), las=2, col=rep(prey.col,each=length(F.scens)), notch=T,
          xaxt="n", cex.lab=cexlab)
            abline(v=c(3.5,6.5),lty=3,lwd=0.5)
            minor.tick(ny=5,nx=0)
            axis(1,at=c(1:18),labels=rep(c("D", "D2", "M"),6), las=1)    #c("High", "Decr", "Low")  cex.axis=cexaxis, 
            axis(1,at=c(2,5,8),tick=F,line=2, labels=rep(c("Prey 1","Prey 2","Prey 3"),1))
            axis(1,at=c(5),tick=F,line=4, labels=c("Scenario"), cex.axis=cexlab)
            
            axis(1,at=0.5,labels=expression(italic("Fishing                         ")), tick=F, line=0, padj=0)  #Labels for the type of scenario    
            axis(1,at=0.5,labels=expression(italic("Pulse                             ")), tick=F, line=2, padj=0)     
            #axis(1,at=0,labels=expression(italic("Migration                   ")), tick=F, line=4, padj=0)     
            legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
                  col=prey.col,pch=15,bg="white",pt.cex=1.5)  # , cex=cexlab
      aggregate(data=t.diff.sub, t.diff.perc~Scenario, median)  #Calculate the median values shown in the figure

       
    #PLOT OF TIME TO TARGET BY SCENARIO
      if(FALSE){
      #Organize data 
      t.diff.sub.base = subset(t.diff.sub, Scenario == "Pulse1-mig-decrF")  #Specify one scenario to extract the "NoPulse" baseline data
      t.diff.sub.base$Pulse = t.diff.sub.base$NoPulse                       #Copy the nopulse values into the pulse column for plotting later
      t.diff.sub.base$t.diff = 0                                            #Set t.diff to 0 to avoid confusion
      t.diff.sub.base[,c("Scenario","P.scen")] = "NoPulse"                  #Rename the scenario
      t.diff.sub.2 = rbind(t.diff.sub.base,t.diff.sub)                      #combine the baseline data with the rest.
      

      #Plot time to target 
        #windows()
        par(mar=plot.mar.w)  # c(7,5,5,2)
        boxplot(Pulse~P.scen, data=t.diff.sub.2, main=paste("SENSITIVITY: Pulse Effects on Management Goals (",F.scens,", ",I.scens,")",sep=""), 
                col=c("white",rep(c("gray","red","blue"),length(P.scens)/3)),
                ylab="Time to reach target SSB (yr)", xlab="", las=2, cex.main=cexmain, cex.lab=cexlab, notch=notch.gr ,#)   
                xlab="", xaxt="n"  )
                axis(1,at=c(1:(1+length(P.scens))),labels=c("No Pulse", P.scens),cex.axis=cexaxis, las=2)
                minor.tick(ny=5, nx=0, tick.ratio=0.5)                                                             
                abline(v=c(1.5,4.5), lty=2, col="red")
            }
       
       } #END t.plots.sens
   
   
   t.plots.sens()





######################################
#SSB TRAJECTORIES THROUGH TIME (9-9-13) - Figuring out why sim#12,13 "NoPulse" are anomalous for time to target - yielding a large effect of "Pulses" on time to target SSB.  
######################################
  #See manuscript outline> results> time to target, for notes.
  #S.all[scen,sim,year,NPvP]
  #dimnames(S.all)
  #head(t.diff.sub)
  
  
  scen.pl = "Pulse3-mig-decrF"
  j = 24  #sim=24 was the biggest outlier, particularly for Pulse3-mig-decrF
  #Plot the SSB trajectories, and highlight some specific cases
  for(j in c(24))  {   # Big outliers 24, 40 ),  Normal: 61, 55, 54, 28, 23, 15
  #windows()
  plot(1:years, S.all[scen.pl,j,,"NoPulse"]/1E6, col="red", type="l", xlab="Year", ylab="SSB (kmt)", ylim=c(0,120),
        main = paste("SSB trajectories -",scen.pl,"-",j))
    abline(v=25, lty=3)
    abline(v=50, lty=3)
    abline(h=91.4, lty=2) #This is the S.target with Fmsy=0.31
    minor.tick(nx=4)
      lines(S.all[scen.pl,j,,"Pulse"]/1E6, col="black", type="l")
  par(new=T)
     }
  par(new=F)
    


  #PLOT OF ONE OF THE LARGEST OUTLIERS
  scen.pl = "Pulse3-mig-decrF"
  j = 24  #sim=24 was the biggest outlier, particularly for Pulse3-mig-decrF

  #pdf("C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SimFigs_SSB trajectory_outlier_10-8-13.pdf", height=7,width=7) 
  
  plot(1:years, S.all[scen.pl,j,,"NoPulse"]/1E6, col="red", type="l", xlab="Year", ylab="S (kmt)",
        #main = paste("SSB trajectories -",scen.pl,"-",j),
        #ylim=c(0,120), 
        cex.lab=cexlab) 
    abline(v=25, lty=3)
    abline(v=50, lty=3)
    abline(h=91.4, lty=2) #This is the S.target with Fmsy=0.31
    minor.tick(nx=4)
      lines(S.all[scen.pl,j,,"Pulse"]/1E6, col="black", type="l")
   dev.off()
   

    #OLD
    for (i in 1:14) {
      lines(S.all[1,i,,"NoPulse"]/1E6, col="black", type="l")
      }
      lines(S.all[1,12,,"NoPulse"]/1E6, col="red", type="l", lwd=3)   #Took 33 yrs to reach target (sim=12)
      lines(S.all[1,13,,"NoPulse"]/1E6, col="blue", type="l", lwd=3)  #Took 34 yrs to reach target (sim=13)
      lines(S.all[1,9,,"NoPulse"]/1E6, col="green", type="l", lwd=3)  #Never reached the target (sim=9)



    #JUNK - looking at outlier time to target...
      #Particularly problematic  simulation runs: Sim=12, 13, 84.
    junk = subset(t.diff.sub, t.diff>10)
    junk

    
    t.nopulse = subset(t.diff, Scenario=="Pulse1-mig-decrF")
    t.nopulse$Pulse = t.nopulse$NoPulse
    t.nopulse
    

######################################
# PROBABILITY OF REACHING S.TARGET (9-9-13)   
######################################

#Function for plotting the probability of reaching S.target by year.
    
t.prob = function(F.scens = c("decrF","dec2F","mor.F"), I.scens = "mig", P.scens.type = "Pulse", LABEL=T) { 
      if(P.scens.type == "Pulse") P.scens = c("Pulse1","Pulse2","Pulse3")
      if(P.scens.type == "Step")  P.scens = c("Step.1","Step.2","Step.3")
      if(P.scens.type == "Both") P.scens = c("Pulse1","Pulse2","Pulse3", "Step.1","Step.2","Step.3")

    #Troubleshooting: running code outside of function
      #F.scens = c("decrF","dec2F","mor.F"); I.scens = "mig"; P.scens.type = "Pulse"; LABEL=T
      #P.scens = c("Pulse1","Pulse2","Pulse3")
    
    #Re-arrange data and subset for scenarios of interest
    t.diff = dcast(time.target.df, Scenario + Sim + P.scen+I.scen+F.scen~ PvNP, value.var="Time.target", mean)
    t.diff$t.diff = t.diff$NoPulse - t.diff$Pulse   #Calculate diff in time to target SSB
    t.diff$t.diff.perc = round(t.diff$t.diff/t.diff$NoPulse*100,2)   #Calculate diff in time to target SSB
    t.diff.sub = na.omit(subset(t.diff,F.scen %in% F.scens & I.scen%in% I.scens & P.scen %in% P.scens))  #Subset the dataset based on conditions
    t.diff.sub$Scenario = factor(t.diff.sub$Scenario)   #Re-set factor levels
        head(t.diff.sub)
        dim(t.diff.sub)        

    head(time.target.df)

    #Make table to store plotted values
    CumProbTable = data.frame(Years = seq(0,50, by=1), 
                              M.P1=NA, M.P2=NA, M.P3=NA, M.NP=NA,
                              D2.P1=NA, D2.P2=NA, D2.P3=NA, D2.NP=NA,
                              D.P1=NA, D.P2=NA, D.P3=NA, D.NP=NA                             
                              )
    
  par(mar=plot.mar.w)
  #par(mfrow=c(3,1))  
  for (i in F.scens) {   #Run a loop to plot all of the probability plots for each F.scenario
    F2plot = i  #Pick the F.scenario that you want to plot below
    Tx = time.target.df[time.target.df$P.scen=="Pulse1" & time.target.df$I.scen=="mig" & 
                        time.target.df$F.scen==F2plot & time.target.df$PvNP=="NoPulse", "Time.target"]
        Tx.breaks = seq(0,50, by=1)
        Tx.cut = cut(Tx, Tx.breaks, right=FALSE)
        Tx.freq = table(Tx.cut)
        Tx.cumfreq = c(0,cumsum(Tx.freq))

    T1 = time.target.df[time.target.df$P.scen=="Pulse1" & time.target.df$I.scen=="mig" & 
                        time.target.df$F.scen==F2plot & time.target.df$PvNP=="Pulse", "Time.target"]
        T1.breaks = seq(0,50, by=1)
        T1.cut = cut(T1, T1.breaks, right=FALSE)
        T1.freq = table(T1.cut)
        T1.cumfreq = c(0,cumsum(T1.freq))
    
    T2 = time.target.df[time.target.df$P.scen=="Pulse2" & time.target.df$I.scen=="mig" & 
                        time.target.df$F.scen==F2plot & time.target.df$PvNP=="Pulse", "Time.target"]
        T2.breaks = seq(0,50, by=1)
        T2.cut = cut(T2, T2.breaks, right=FALSE)
        T2.freq = table(T2.cut)
        T2.cumfreq = c(0,cumsum(T2.freq))

    T3 = time.target.df[time.target.df$P.scen=="Pulse3" & time.target.df$I.scen=="mig" & 
                        time.target.df$F.scen==F2plot & time.target.df$PvNP=="Pulse", "Time.target"]
      T3.breaks = seq(0,50, by=1)
      T3.cut = cut(T3, T3.breaks, right=FALSE)
      T3.freq = table(T3.cut)
      T3.cumfreq = c(0,cumsum(T3.freq))

    #Fill in table of values
      if(i=="decrF") CumProbTable[,c("D.P1","D.P2","D.P3","D.NP")] = cbind(T1.cumfreq/sims,T2.cumfreq/sims,T3.cumfreq/sims,Tx.cumfreq/sims)
      if(i=="dec2F") CumProbTable[,c("D2.P1","D2.P2","D2.P3","D2.NP")] = cbind(T1.cumfreq/sims,T2.cumfreq/sims,T3.cumfreq/sims,Tx.cumfreq/sims)
      if(i=="mor.F") CumProbTable[,c("M.P1","M.P2","M.P3","M.NP")] = cbind(T1.cumfreq/sims,T2.cumfreq/sims,T3.cumfreq/sims,Tx.cumfreq/sims)
    
    
  plot(Tx.breaks, Tx.cumfreq/sims, lty=2, type ="l", lwd=ifelse(i=="decrF",3,1), xlim=c(0,40), ylim=c(0,1),
    #main=ifelse(i=="decrF",expression("Time to"~italic("S"["target"])),""), 
    xlab=ifelse(i=="decrF","Years",""),  
    ylab=ifelse(i=="decrF",expression("Cum. probability of acheiving S "["MSY"]),""),
    cex.lab=cexlab, cex.axis=cexaxis, xaxt="n", yaxt="n", las=2 )
    abline(h=0.5,lty=3)      
    if(i=="decrF") {  #Draw axes only once
      axis(1,at=seq(0,40,by=1),tck=-0.01, labels=F);   axis(1,at=seq(0,40,by=5), cex.axis=cexaxis)
      axis(2,at=seq(0,1,by=0.1),tck=-0.01, labels=F); axis(2,at=seq(0,1,by=0.2),las=2, cex.axis=cexaxis)  }
  
    lines(T1.breaks, T1.cumfreq/sims, col="black", lwd=ifelse(i=="decrF",3,1))
    lines(T2.breaks, T2.cumfreq/sims, col="red",  lwd=ifelse(i=="decrF",3,1))
    lines(T3.breaks, T3.cumfreq/sims, col="blue", lwd=ifelse(i=="decrF",3,1))
    if(i=="decrF") points(c(0,12),c(0,0), cex=1.5, pch=25, bg="black")
    if(i=="dec2F") points(c(0,6),c(0,0), cex=1.5, pch=25, bg="black")
    if(i=="mor.F") points(c(0,0),c(0,0), cex=1.5, pch=25, bg="black") 
    if(i=="decrF") {legend("bottomright",c("Base", "Pulse 1","Pulse 2","Pulse 3"), title="Scenario",
            col=c("black","black","red","blue"), lty=c(2,1,1,1), lwd=3, cex=cexaxis)  }
    par(new=T)
       
       } #END LOOP
    par(new=F)
  
    #Add F scenario labels
    if(LABEL==T){
    text(0,0.01,labels="M",srt=0, adj=0, pos=3, cex=cexlab)
    text(6,0.01,labels="D2",srt=0, adj=0, pos=3, cex=cexlab)
    text(12,0.01,labels="D",srt=0, adj=0, pos=3, cex=cexlab)
    #text(0,0.03,labels="Mo-F",srt=70, adj=0, cex=cexlab)
    #text(6,0.03,labels="D2-F",srt=70, adj=0, cex=cexlab)
    #text(11,0.03,labels="De-F",srt=70, adj=0, cex=cexlab)
      }

  #Display table with the cumulative probability values
    print(CumProbTable)
  
  
  
      } #END t.prob()

t.prob()

  #pdf("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SimFigs_CumProbPlot_7-27-14.pdf", height=7,width=7)    
   t.prob()
   dev.off()
 




########################################################   
### BASIC GROWTH PLOTS (PredPreySim() Output) 
########################################################
  #NOTE: to run the code below, first run the PredPreySim() function (using default parameters)
  # to generate the necessary tables and objects.
  
  #Run Function to get outputs
  PredPreySim(seed=350); N.prey.pl.base = N.prey
  PredPreySim(P.scenario="Pulse1", seed=350); N.prey.pl.pulse = N.prey

 #pdf("C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SimFigs_Growth&Cons_10-7-13.pdf", height=8,width=4) 
   
  ### GRAPHING GROWTH OUTPUT
    #windows() #open new graph window for plotting    width=7,height=12
    #Temporary plotting to investigate Growth trajectories and Cmax
      par(oma=c(4,0,0,0))
      par(mfrow=c(3,1))
      par(mar=c(0.5,5.1,2.1,2.1))
      x.ages = seq(0.25,ages,by=0.25)

      #Relative diet by age
      plot(x.ages,rep(-50,ages*seasons), ylim=c(0,100),ylab="Diet (%W)",xlab="Age", main="", cex.lab=cexlab)    #Create empty plot
        for (y in (ages+1):years){
          lines(x.ages,100*FxnR[y,,1,,1]/apply(FxnR[y,,1,,],c(1,2),sum),col="black")
          lines(x.ages,100*FxnR[y,,1,,2]/apply(FxnR[y,,1,,],c(1,2),sum),col="red")    
          lines(x.ages,100*FxnR[y,,1,,3]/apply(FxnR[y,,1,,],c(1,2),sum),col="blue")    
           }
        lines(x.ages, diet.1.exp*100, col="gray40", lwd=4)
        lines(x.ages, diet.2.exp*100, col="red4", lwd=4)               
        lines(x.ages, diet.3.exp*100, col="midnightblue", lwd=4)               
        minor.tick(nx=2,ny=2)
        legend("right",c("Prey 1","Prey 2","Prey 3"), col=c("black","red","blue"), lty=1, lwd=2, bg="white")
      
      #Consumption relative to Cmax by age
      plot(x.ages,Cmax, type="l",lwd=3, col="green", main="",xlab="Age", cex.lab=cexlab, 
      ylab=expression("Consumption (kg  fish"^" -1"~"season"^" -1"~")") )   
      for (y in (ages+1):years){
        lines(x.ages,apply(FxnR[y,,1,,],c(1,2),sum))  }    #calculates sum of prey consumed by season and age (for each year y, and region 1)
        minor.tick(nx=2,ny=2)
        legend("topleft",c(expression("C"^"max"),"C"), col=c("green","black"), lty=1, lwd=c(3,1))
      
      #par(mar=c(4.5,5.1,2.1,2.1) )
      #Weight at age plot
      plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="gray", main="",ylab="Weight (kg)",xlab="Age",xlim=c(0,8),ylim=c(0,2.5), cex.lab=cexlab)    
      for (y in (ages+1):years){
        lines(x.ages,w.pred[y,,1,], col="blue")  }
        lines(x.ages,w, col="red", lwd=2)
        minor.tick(nx=2,ny=2)
        legend("topleft",c("Mean from wild","Simulated"), col=c("red","blue"), lty=1, lwd=c(2,1))
        
        axis(1,at=c(4),tick=F,line=2, labels="Age", cex.axis=cexlab) #Add label to X-axis of bottom plot
      
      dev.off()
      
  ### PREY BIOMASS TIME-SERIES OUTPUT

    #pdf("C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SimFigs_PreyBiomass_10-5-13.pdf", height=3.5,width=4) 

      #Time series of prey biomass                                N.prey.1.mean/1E6/1E3
      par(mfrow=c(1,1))
      par(mar=c(4.5,4.5,1,1) )
      plot(seq(1,by=0.25,length=years*seasons),log(t(N.prey.pl.pulse[,,1,1]/1E6/1E3)), type="l", lwd=1, col="red",     #Converting N.prey.1.mean (in #prey) to kmt, assuming 1g/fish
        xlim=c(20,60), cex.lab=cexlab,
        main="",xlab="Year",ylab="log(Prey Biomass)")
        lines(seq(1,by=0.25,length=years*seasons),log(t(N.prey.pl.base[,,1,1]/1E6/1E3)), type="l", lwd=1, col="black")
        abline(v=c(25,50), lty=3)
        minor.tick(nx=5, ny=0)
          
      dev.off()
      
      
      PredPreySim(P.scenario = "Pulse1")








































### ORIGINAL (Pre 3/19). MAKING PLOTS 

##SSB PLOTS
  #Create subsets for different plots
  S.sub.1 = S.diff.p.mean.df[S.diff.p.mean.df$I.scen!="mix",];
      S.sub.1$Scenario = factor(S.sub.1$Scenario, levels=S.sub.1[!duplicated(S.sub.1$Scenario),"Scenario"])
  S.sub.2 = S.diff.p.mean.df[S.diff.p.mean.df$F.scen!="decrF",];
      S.sub.2$Scenario = factor(S.sub.2$Scenario, levels=S.sub.2[!duplicated(S.sub.2$Scenario),"Scenario"])
  S.sub.3 = S.diff.p.mean.df[S.diff.p.mean.df$F.scen!="decrF" & S.diff.p.mean.df$I.scen!="mix",];
      S.sub.3$Scenario = factor(S.sub.3$Scenario, levels=S.sub.3[!duplicated(S.sub.3$Scenario),"Scenario"])
  S.sub.4 = S.diff.p.mean.df[S.diff.p.mean.df$F.scen=="highF",];
      S.sub.4$Scenario = factor(S.sub.4$Scenario, levels=S.sub.4[!duplicated(S.sub.4$Scenario),"Scenario"])
      levels(S.sub.4$Scenario)
      
  #Plot of all scenarios
    windows()
    par(mar=c(10,6,5,2))
    boxplot(S.diff.p.mean, main="SSB change due to prey pulses - Pulse years",
            cex.lab=1.4,xlab="",ylab="mean SSB %diff", notch=notch.gr, las=2)

  #SSB vs. F.scen | Prey  (Migration only)
    par(mar=plot.mar)
    boxplot(S.diff.p.mean~Scenario, data=S.sub.1, main="SSB change due to prey pulses",
        xlab="Scenario",ylab="Mean % Difference (from No-pulse)", notch=notch.gr, col=prey.col,
        xaxt="n", cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis )
        abline(v=c(3.5,6.5), lty=2, lwd=0.5)
        axis(1,at=c(2,5,8),labels=c("High F", "Low F", "Decr. F"), cex.axis=cexaxis)
        legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
              col=prey.col,pch=15,bg="white", cex=cexlab)
  #SSB vs. F.scen | Prey  (Migration only, exclude decreasing F)
    par(mar=plot.mar)
    boxplot(S.diff.p.mean~Scenario, data=S.sub.3, main="SSB change due to prey pulses",
        xlab="Scenario",ylab="Mean % Difference (from No-pulse)", notch=notch.gr, col=prey.col,
        xaxt="n", cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis )
        abline(v=c(3.5), lty=2, lwd=0.5)
        axis(1,at=c(2,5),labels=c("High F", "Low F"), cex.axis=cexaxis)
        legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
              col=prey.col,pch=15,bg="white", cex=cexlab)

  #SSB vs. I.scen | Prey  (Migration only)
    par(mar=plot.mar)
    boxplot(S.diff.p.mean~Scenario, data=S.sub.2, main="SSB change due to prey pulses",
        xlab="",ylab="Mean % Difference (from No-pulse)", notch=notch.gr, col=prey.col, las=2,
        xaxt="n",cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis )         #
        abline(v=c(3.5,6.5,9.5), lty=2, lwd=0.5)
        abline(v=c(6.5), lty=1, lwd=2)
        axis(1,at=c(2,5,8,11),labels=rep(c("Mig","Mix"),2), las=1, cex.axis=cexaxis)
        axis(1,at=c(3.5,9.5), labels=c("High F", "Low F"), tick=F, line=2, cex.axis=cexaxis)
        legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
              col=prey.col,pch=15,bg="white", cex=cexlab)
  #SSB vs. I.scen | Prey  (Migration only; High F only)
    par(mar=plot.mar)
    boxplot(S.diff.p.mean~Scenario, data=S.sub.4, main="SSB change due to prey pulses",
        xlab="Scenario",ylab="Mean % Difference (from No-pulse)", notch=notch.gr, col=prey.col, las=2,
        xaxt="n",cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis )         #
        abline(v=c(3.5), lty=2, lwd=0.5)
        #abline(v=c(6.5), lty=1, lwd=2)
        axis(1,at=c(2,5),labels=rep(c("Mig","Mix"),1), las=1, cex.axis=cexaxis)
        #axis(1,at=c(3.5,9.5), labels=c("High F", "Low F"), tick=F, line=2, cex.axis=cexaxis)
        legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
              col=prey.col,pch=15, bg="white", cex=cexlab)

##YIELD PLOTS
  #Create subsets for different plots
  Y.sub.1 = Y.diff.p.mean.df[Y.diff.p.mean.df$I.scen!="mix",];
      Y.sub.1$Scenario = factor(Y.sub.1$Scenario, levels=Y.sub.1[!duplicated(Y.sub.1$Scenario),"Scenario"])
  Y.sub.2 = Y.diff.p.mean.df[Y.diff.p.mean.df$F.scen!="decrF",];
      Y.sub.2$Scenario = factor(Y.sub.2$Scenario, levels=Y.sub.2[!duplicated(Y.sub.2$Scenario),"Scenario"])

  #Plot of all scenarios
    windows()
    par(mar=c(10,6,5,2))
    boxplot(Y.diff.p.mean, main="YIELD change due to prey pulses - Pulse years",
            cex.lab=1.4,xlab="",ylab="mean YIELD %diff", notch=notch.gr, las=2)

  #YIELD vs. F.scen | Prey  (Migration only)
    par(mar=plot.mar)
    boxplot(Y.diff.p.mean~Scenario, data=Y.sub.1, main="YIELD change due to prey pulses",
        xlab="Scenario",ylab="Mean % Difference (from No-pulse)", notch=notch.gr, col=prey.col,
        xaxt="n", cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis )
        abline(v=c(3.5,6.5), lty=2, lwd=0.5)
        axis(1,at=c(2,5,8),labels=c("High F", "Low F", "Decr. F"), cex.axis=cexaxis)
        legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
              col=prey.col,pch=15,bg="white", cex=cexlab)

  #YIELD vs. I.scen | Prey  (Migration only)
    par(mar=plot.mar)
    boxplot(Y.diff.p.mean~Scenario, data=Y.sub.2, main="YIELD change due to prey pulses",
        xlab="",ylab="Mean % Difference (from No-pulse)", notch=notch.gr, col=prey.col, las=2,
        xaxt="n",cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis )         #
        abline(v=c(3.5,6.5,9.5), lty=2, lwd=0.5)
        abline(v=c(6.5), lty=1, lwd=2)
        axis(1,at=c(2,5,8,11),labels=rep(c("Mig","Mix"),2), las=1, cex.axis=cexaxis)
        axis(1,at=c(3.5,9.5), labels=c("High F", "Low F"), tick=F, line=2, cex.axis=cexaxis)
        legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
              col=prey.col,pch=15,bg="white", cex=cexlab)


##WEIGHT-AT-AGE PLOTS

  #Create subsets for different plots
  w.sub.1 = w.diff.p.mean.df[w.diff.p.mean.df$I.scen=="mig" & w.diff.p.mean.df$F.scen=="highF",];
      w.sub.1$Scenario = factor(w.sub.1$Scenario, levels=w.sub.1[!duplicated(w.sub.1$Scenario),"Scenario"])
      levels(w.sub.1$Scenario)
      summary(w.sub.1)

  #Weight at age BOXPLOT (migration only and highF only)  - F scenarios have no impact on weight at age
    par(mar=plot.mar)
    boxplot(w.diff.p.mean~Scenario*Age, data=w.sub.1, main="Weight-at-age change due to prey pulses",
          xlab="Age", ylab="Mean % Difference (from No-pulse)", notch=notch.gr, col=prey.col,
          xaxt="n", cex.lab=cexlab, cex.main=cexmain, cex.axis=cexaxis)
            for(i in 0:7) {                        #make "bars" to help identify the ages.
              lines(c(1:3)+i*3,rep(0,3), lwd=5)  }
          minor.tick(ny=2,nx=0)
          axis(1,at=seq(2,24, by=3), labels=c(0:(ages-1)), las=1, cex.axis=cexaxis)
          legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
                col=prey.col,pch=15,bg="white", cex=cexlab)

  #ORIGINAL CODE: Change in WEIGHT-AT-AGE
      for (i in unique(paste(scen.names$I.scen,scen.names$F.scen[1],sep="-")))  {
        w.diff.p.mean.df.2 = subset(w.diff.p.mean.df, Scenario %in% paste(unique(scen.names$P.scen),i,sep="-"))  #subset only for the migration
        w.diff.p.mean.df.2$Scenario = factor(w.diff.p.mean.df.2$Scenario)                                         #re-"factorize" Scenario to remove any unused factor levels
        windows()
        boxplot(w.diff.p.mean~Scenario*Age, data=w.diff.p.mean.df.2, main=paste(i,"scenario"), col=c("gray","red","blue"),
              ylab="Weight-at-age increase (%) during prey pulse", xlab="Age", xaxt="n", notch=notch.gr )
              axis(1,at=c(1:(preys*ages)),labels=rep(0:(ages-1),each=preys))
              legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Pulses in:",
                  col=c("gray","red","blue"),pch=15,bg="white")
                   }  #end i loop

  #Weight at age plot         #w.all[scen,sim,years,seas,ages,PulseNoPulse]
    #Make vector of fractional ages
    x.ages = seq(0.25,ages,by=0.25)

    #Create an array with only the weight data for the years in which a pulse occurred
    w.pulses = w.all[,,c(1:n.pulse),,,2]   # make a smaller array with only n.pulse # of years
      w.pulses[,,,,]  =0                   #turn all entries into zero to help avoid errors
      for (s in 1:sims){                   #fill in the array with data from PULSE years only.
          w.pulses[,s,,,] = w.all[,s,pulses.all[s,],,,2]
               }
    #Make Plot
    windows()
    plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="darkgray", main=paste("Mean weights at age",scen.names$I.scen[1],sep=" - "),
      ylab="Weight (kg)",xlab="Age",xlim=c(0,10),ylim=c(0,3))
      lines(x.ages, apply(w.all[1,,,,,1],c(3,4),mean), col="black", lwd=3, lty=2)  #Scen=1 (NO PULSE); averaging across years and simulations. NOTE: The NO-pulse values are identical across Scenarios as long as the I.scen is the same.
      lines(x.ages, apply(w.pulses[1,,,,],c(3,4),mean), col="black", lwd=6)  #Scen=1 (PULSE YEARS); averaging across years and simulations
      lines(x.ages, apply(w.pulses[2,,,,],c(3,4),mean), col="red", lwd=3)  #Scen=2 (PULSE YEARS); averaging across years and simulations
      lines(x.ages, apply(w.pulses[3,,,,],c(3,4),mean), col="blue", lwd=5)  #Scen=3 (PULSE YEARS); averaging across years and simulations
        legend("topleft",c("No Pulse","Prey-1 Pulse","Prey-2 Pulse","Prey-3 Pulse"), col=c("black","black","red","blue"),
              lty=c(2,1,1,1), lwd=c(2,6,3,5))
          #JUNK: calculate %diff in w@age due to pulse, to compare to boxplots:
          #(apply(w.pulses[1,,,,],c(3,4),mean) - apply(w.all[1,,,,,1],c(3,4),mean)) /  apply(w.all[1,,,,,1],c(3,4),mean)

## TIME TO TARGET SSB
  #PLOTS OF DECREASING F SUBSCENARIOS:
    #Organize data
      time.target.df$Scen = paste(time.target.df$Scenario,time.target.df$PvNP,sep="-")     #Create new variable to use to subset the data
      time.target.sub = time.target.df[-which(time.target.df[,"PvNP"]=="NoPulse" & time.target.df[,"P.scen"] !="Pulse1"),]
         #time.target.sub$Scen = factor(time.target.sub$Scen, levels =   #Specify order of factor levels:
          c("Pulse1-mig-decrF.24-NoPulse","Pulse1-mig-decrF.22-NoPulse", "Pulse1-mig-decrF.21-NoPulse","Pulse1-mig-decrF-NoPulse","Pulse1-mig-decrF.19-NoPulse","Pulse1-mig-decrF.18-NoPulse","Pulse1-mig-decrF.16-NoPulse",
            "Pulse1-mig-decrF.24-Pulse","Pulse1-mig-decrF.22-Pulse", "Pulse1-mig-decrF.21-Pulse","Pulse1-mig-decrF-Pulse","Pulse1-mig-decrF.19-Pulse","Pulse1-mig-decrF.18-Pulse","Pulse1-mig-decrF.16-Pulse",
            "Pulse2-mig-decrF.24-Pulse","Pulse2-mig-decrF.22-Pulse", "Pulse2-mig-decrF.21-Pulse","Pulse2-mig-decrF-Pulse","Pulse2-mig-decrF.19-Pulse","Pulse2-mig-decrF.18-Pulse","Pulse2-mig-decrF.16-Pulse",
            "Pulse3-mig-decrF.24-Pulse","Pulse3-mig-decrF.22-Pulse", "Pulse3-mig-decrF.21-Pulse","Pulse3-mig-decrF-Pulse","Pulse3-mig-decrF.19-Pulse","Pulse3-mig-decrF.18-Pulse","Pulse3-mig-decrF.16-Pulse")
            )
         time.target.sub$Scen = factor(time.target.sub$Scen, levels =   #Specify order of factor levels:
          c("Pulse1-mig-decrF.24-NoPulse","Pulse1-mig-decrF.24-Pulse","Pulse2-mig-decrF.24-Pulse","Pulse3-mig-decrF.24-Pulse",
            "Pulse1-mig-decrF.22-NoPulse", "Pulse1-mig-decrF.22-Pulse","Pulse2-mig-decrF.22-Pulse", "Pulse3-mig-decrF.22-Pulse",
            "Pulse1-mig-decrF.21-NoPulse","Pulse1-mig-decrF.21-Pulse","Pulse2-mig-decrF.21-Pulse","Pulse3-mig-decrF.21-Pulse",
            "Pulse1-mig-decrF-NoPulse","Pulse1-mig-decrF-Pulse","Pulse2-mig-decrF-Pulse","Pulse3-mig-decrF-Pulse",
            "Pulse1-mig-decrF.19-NoPulse","Pulse1-mig-decrF.19-Pulse","Pulse2-mig-decrF.19-Pulse","Pulse3-mig-decrF.19-Pulse",
            "Pulse1-mig-decrF.18-NoPulse","Pulse1-mig-decrF.18-Pulse","Pulse2-mig-decrF.18-Pulse","Pulse3-mig-decrF.18-Pulse",
            "Pulse1-mig-decrF.16-NoPulse", "Pulse1-mig-decrF.16-Pulse","Pulse2-mig-decrF.16-Pulse", "Pulse3-mig-decrF.16-Pulse",
            "Pulse1-mig-decrF.14-NoPulse", "Pulse1-mig-decrF.14-Pulse","Pulse2-mig-decrF.14-Pulse", "Pulse3-mig-decrF.14-Pulse",
            "Pulse1-mig-decrF.10-NoPulse", "Pulse1-mig-decrF.10-Pulse","Pulse2-mig-decrF.10-Pulse", "Pulse3-mig-decrF.10-Pulse")
            )
    #Plot all F subscenarios
      windows()
      par(mar=c(14,5,5,2))
      boxplot(Time.target~Scen, data=time.target.sub, main="Pulse Effects on Management Goals", #col=c("gray","red","blue"),
              ylab="Time to reach target SSB (yr)", xlab="", las=2, notch=notch.gr, cex.lab=1.2, ylim=c(10,42) )   # xaxt="n", 
              #axis(1,at=c(1:4),labels=c("None", "Prey 1", "Prey 2", "Prey 3"))
              minor.tick(ny=5, nx=0, tick.ratio=0.5)                                                             
              abline(v=c(4.5,8.5,12.5,16.5,20.5,24.5,28.5,32.5), lty=2, col="red")

    #Plot single, select subscenario
    F.plot = "decrF.14"       #SPECIFY SCENARIO OF INTEREST
    F.list = c(paste("Pulse1-mig-",F.plot,"-NoPulse",sep=""), paste("Pulse1-mig-",F.plot,"-Pulse",sep=""),   #Build pertinent scenario subsets
               paste("Pulse2-mig-",F.plot,"-Pulse",sep=""), paste("Pulse3-mig-",F.plot,"-Pulse",sep=""))
      time.target.sub.1 = subset(time.target.df, Scen %in% F.list)
        time.target.sub.1$Scen = factor(time.target.sub.1$Scen )      #Re-factorize variable to get rid of unwanted factor levels.
      
      windows()
      par(mar=plot.mar)
      boxplot(Time.target~Scen, data=time.target.sub.1, main="Pulse Effects on Management Goals", col=c("white", prey.col),
              ylab="Time to reach target SSB (yr)", xlab="Pulse Scenario", xaxt="n", las=2, notch=notch.gr, cex.lab=cexlab, cex.main=cexmain,
              ylim=c(14,39) )
              axis(1,at=c(1:4),labels=c("None", "Prey 1", "Prey 2", "Prey 3"),cex.axis=cexaxis)
              minor.tick(ny=5, nx=0, tick.ratio=0.5)