Main Model

---

Check for the existance of needed datafiles.

If the files are not present, then run their corresponding scripts and generate them.

The corresponding R Markdown scripts associated with these files can be found on RPub:

• Solar Data Cleaning

• Thermocline Model

• Mysocline Depth Model

setwd("/Users/Nick/mysisModeling") #Locate ourselves in the computer. 

#Check for solar data
if(!file.exists("data/light_day_moon_hour.csv")){ source("dataGen/solarData.r") } 

#Check for thermocline data
if(!file.exists("data/Depth_Thermocline_Hour.csv")){ source("dataGen/thermoclineLevels.r") } 

#Mysocline data
if(!file.exists("data/mysocline_hour.csv")){ source("dataGen/depthModel.r") }

#Mysocline data
if(!file.exists("data/FoodAvail_Hour.csv")){ source("dataGen/foodAvailability.r") }

#Read in the data now. 
depthData = read.csv("data/mysocline_hour.csv")$x
foodCurve = read.csv("data/FoodAvail_hour.csv")
foodAvail = foodCurve$foodAvail
foodVar   = foodCurve$variability

Class and method declarations:

I have chosen to use the S4 object oriented programing style. The details are minor but more info can be found here.

Set up the mysis class:

This basically is a digital mysis that simplistically holds a couple of peices of information and upon which the model will act.

setClass("mysis",
         representation(
           energy    = "numeric",    # The energy reserves are a numeric value
           migrating = "logical",    # If mysis is migrating is a logical value
           alive     = "logical",    # Alive still?
           depth     = "numeric"),   
         prototype(
           energy    = 0,            # The default instantiated mysis starts with zero energy...
           migrating = FALSE,        # and not migrating...
           alive     = TRUE,         # alive
           depth     = 100),         # and at the bottom.  
         )

#Set the show method, basically how we want the program to display the info about the mysis on calling. 
setMethod("show", "mysis", 
          function(object){
            if (object@alive){
              print(object@energy)
              print(object@migrating)
              print(object@depth)
            } else {
              print("dead")
            }
          })

Real quick we need to make a model for decisions to migrate based upon condition:

The form of this model is a logistic curve: \[f(x) = \frac{1}{1 + e^{-k(x - m)}}\] Where \(k\) is an indicator of steepness, and \(m\) is the midepoint of the curve.

migrationProb = function(condition){ #I am choosing to leave the random roll outside of script for continuity, could be changed
  m = 120 #value of curve midpoint
  k = 0.03 #steepness of curve
  x = condition 
  
  dist = 1/(1 + exp(-k* (x - m)))   
  return(1 - dist)
}

The decision tree model:

This is a stochastic model, meaning it utilizes randomness to simulate the real world.

The function of the decision tree is such:

• Draw random number to decide if migrating:
• If the number drawn is above the ratio needed to migrate (food ratio) then the mysis migrates.
  • After migration is decided, another draw is done to see if the mysis is killed by predation.

rewardUnits = .68

setGeneric( "nextTime", function(object, ...){standardGeneric("nextTime")})
setMethod("nextTime","mysis", 
          function(object, foodAvail, foodVar, time){       #Takes in the mysis object and the ratio of food quality at a given time. 
            
            
            #Run the draws: 
            migrationDraw   = runif(1) #random number between 0 and 1, this will be used for the migration decision
            migrationDraw_2 = runif(1) #draw for second migration decision, this time compared to logistic curve
            predationDraw   = runif(1) #" " to see if killed by predation
            
            #Lets set some constants real quick: 
            migrationCost = 21       #How many energy units the mysis use up migrating
            migrationRisk = 0  #We don't really need to simulate predation.
            stayRisk      = 0 
#             migrationRisk = 0.0001   #Chance of being eaten if migrating
#             stayRisk      = 0.00001  #Chance of being eaten if they stay on the bottom
            
            #Energy rewards:
            migrationReward = rnorm(1, (rewardUnits*(1 + foodAvail)), foodVar)
            stayReward      = migrationReward * .2
            
            threshold = 0.2  #Let's make sure that the benthic reward can never drop below a threshold.
            if(stayReward < threshold){
              stayReward = threshold
            }
            
#             migrationReward = rnorm(1,rewardUnits, foodVar)
#             stayReward      = rnorm(1,(rewardUnits * ( 1 - foodAvail)), foodVar) 
            #this makes the rewards responsive to conditions
            
            sunset  = 18 #hardcode sunset, fill this with real data later. 
            sunrise = 7
            conditionCurve = migrationProb(object@energy)
            
            
            if (object@energy <= 0){ #did the mysis starve?
              object@alive = FALSE
            }     
            
            #if the mysis is in good condition use standard migration chance, if bad use 1 - chance
            if ((time %% 24 == sunset + 1)&&(object@alive)){
              
              if (migrationDraw < foodAvail) {
                object@migrating = migrationDraw_2 > conditionCurve 
                object@energy = object@energy - migrationCost 
              } else {
                 object@migrating = FALSE
              }
              
            } else if (time %% 24 == sunrise){ #bring them back down at sunrise. 
              object@migrating = FALSE
            }    
            

            #Here is the decision tree:
            if (object@alive){ #If the mysis is alive let's run the decision tree
  
              if (object@migrating){ 
                object@depth = depthData[time] #Grab the mysocline limit at this hour 
                
                if (predationDraw > migrationRisk){ #The mysis evades predation
                  object@energy = object@energy + migrationReward 
                
                } else { #Mysis is eaten
                  object@alive = FALSE
                }
                
              } else { #The mysis didn't migrate
                
                object@depth = 100
                
                if (predationDraw > stayRisk){ #The mysis evades predation
                  object@energy = object@energy +  stayReward
                  
                } else { #Mysis is eaten
                  object@alive = FALSE
                }
              }
            }
            object} #Return the mysis
          )

Testing it:

Now that the model is all set up can run it over a small subset of the eventual numbers (~10,000 mysids over 8,700 hours) and see what results we get.

mysids = NULL
numOfMysids = 50

for (i in 1:numOfMysids){
  initialenergy = rnorm(1,150,25) #draw initial energy from normal dist centered at 30 with stdev of 5. 
  mysids = c(mysids, new("mysis", energy = initialenergy) )
}

… and now we can run the model on those mysids:

migrations = NULL
conditions = NULL
hours = 1:(24*365) #150 days
for (mysid in mysids){ #loop through the mysis
  migration = NULL 
  condition = NULL
  for (i in hours){

    mysid    = nextTime(mysid, foodAvail[i], foodVar[i] , i)
    migration = c(migration, mysid@depth)
    condition = c(condition, mysid@energy)
  }
  migrations = cbind(migrations, migration) # add to the object of migrations. 
  conditions = cbind(conditions, condition)
} 

What do we have?

Let’s plot this small run to see what kind of behavior we are getting.

Points to address:

This obviously isn’t right.

## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.

Let’s take a look at the energy reward distributions, monte carlo style:

We will set the food quality at 0.7

# #rewardUnits = 1
# foodAvail   = 0.7
# foodVar     = 0.3
# migrateVals = NULL
# stayVals = NULL
# for (i in 1:5000){
#   
# #   rnorm(1,rewardUnits, rewardUnits*stdMult)
# #             stayReward      = rnorm(1,(rewardUnits * ( 1 - foodAvail)), rewardUnits*stdMult) 
# migrateVals = c(migrateVals, rnorm(1,rewardUnits, foodVar))
# stayVals = c(stayVals, rnorm(1,(rewardUnits * ( 1 - foodAvail)), foodVar )  )
# }
# 
# par(mfrow=c(2,1))
# hist(migrateVals, main = "Migrate", xlim=c(-3,6))
# hist(stayVals, main = "Stay", xlim = c(-3,6))?

It would appear that we need to be a litle more ‘mean’ with our taking away of energy.