Adding plotting code to sense_beta

[kaydin]

Sense_beta branch is synced with master (March 10 2021). Add plotting code.

[kaydin]

Also include other bits (i.e. from Andy and Kerim)

[sgaichas]

@kaydin @aw2113 looking at my old plotting code, it is highly dependent on how I processed the sense outputs, which used 5+ other scripts I had deployed in various folder, etc. I had a combination of .sh and .R scripts.

I'm guessing you have similar processing scripts somewhere since ecosense.R outputs a single perturbed system, so you have to run it N times, figure out which systems passed the initial burnin (I have some diagnostic plots for this too), bin results into quantiles by perturbed species or scenario (which I may have embedded in C code before because otherwise the output was overwhelming), etc.

If you want to post your scripts somewhere I could try turning them into a wrapper function, then the plotting code would work from that.

not a rush, just thought I would post this note now since I just looked. thanks!

[kaydin]

`

Non-stanzas model test

Ebase <- "models/sensetests/EBS_base_Aug2016_AGG.csv"
Ediet <- "models/sensetests/EBS_diet_Aug2016_AGG.csv"
Eped <- "models/sensetests/EBS_ped_Aug2016_AGG.csv"
unbal <- read.rpath.params(Ebase, Ediet, Eped)

Setup Base Ecopath and Base Rsim scenario

bal <- rpath(unbal) # balanced
all_years <- 1991:2090
scene <- rsim.scenario(bal, unbal, years = all_years) # Ecosim params

UNCHANGED FROM Andy's test loop EXCEPT set seed.

Should come up with identical output.

scene$params$BURN_YEARS <- 50
NUM_RUNS <- 1000
parlist<-as.list(rep(NA,NUM_RUNS))
kept<-rep(NA,NUM_RUNS)

set.seed(12345) # for testing

for (i in 1:NUM_RUNS){
EBSsense <- scene

parlist[[i]]<- scene$params # Base ecosim params

parlist[[i]]<- rsim.sense.orig(scene, bal, unbal) # Replace the base params with Ecosense params

New rsim.sense should give same accept/rejects as Andy's rsim.sense.orig given the same seed

parlist[[i]]<- rsim.sense(scene, unbal, Vvary=c(-4.5,4.5)) # Replace the base params with Ecosense params
EBSsense$start_state$BB <- parlist[[i]]$B_BaseRef
parlist[[i]]$BURN_YEARS <- 50 # Set Burn Years to 50
EBSsense$params <- parlist[[i]]
EBStest <- rsim.run(EBSsense, method="AB", years=all_years)
failList <- which(is.na(EBStest$end_state$BB))
{if (length(failList)>0)
{cat(i,": fail in year ",EBStest$crash_year,": ",failList,"\n"); kept[i]<-F; flush.console()}
else
{cat(i,": success!\n"); kept[i]<-T; flush.console()}}
parlist[[i]]$BURN_YEARS <- 1
}

`

[sgaichas]

for completeness, sense scripts used with AK-ecosim inC are posted at https://github.com/sgaichas/model-scriptshare/tree/main/ecosense and other repositories linked in the readme there