---
# System prepended metadata

title: Liora's Project

---

# Liora's Project
## Notes from meeting (8/21)
* we changed the arrow on the interactions diagram (now it goes from medusa to small fish --> this affects the derivative function of the small fish)
* we finished the interactions diagram
* we began making tables on here

## Next Meeting Agenda
* final checkover of functions
* finalize title
* go over subheadings for intro
* begin coding

* ****
* ~~**Define new schedule**~~

## Voluntary tasks for next meeting
* Send functions (Liora).
* Update diagram (Liora).
* Review and comment functions (Carlos).
* **Regeneron**
    * Write Regeneron recommendation (Carlos). **(POSTPONED)**

## State variables table

| Variable | Symbol | Description |
| -------- | -------- | -------- |
| Polyp     | $$P$$     |Polyp Population |
|Ephyra| $$E$$| Ephyra Population
|Medusa|$$M$$|Medua Population
|Plankton|$$PL$$|Plankton Population
|Small Fish|$$F_S$$|Small Fish Population
|Larvae|$$F_L$$|Larvae Population
|Egg|$$F_E$$|Fish Egg Population
|Adult Fish|$$F_A$$|Adult Fish Population

## Parameter table
| Parameter | Symbol | Description |
| -------- | -------- | -------- |
| Predation rate     | $$\alpha_{b,a}$$     | Predation rate of species a on species b     |
| Transition Rate    | $$\eta_{j,i}$$     | Rate at which i becomes j     |
| Death Rate     | $$\delta$$     | $$F_a (\delta_{F_a}+\lambda)$$ with $\lambda=$ Fishing rate     |
| Efficiency Rate| $$\epsilon_i$$ | Effect of species j consuming species i. Which means, how many of species j is produced per unit of consumed species i |



### Regeneron important dates (and docs)

* Project is due November 9th.
    * Our own date is October 16th (Everything but conclusions must be ready).
        * Methods (Carlos and Liora) **Sept 4th**:
            * Diagram (Liora) **Ago 28th**.
            * Functions () (solved to equilibrium?) **Sept 4th** (Carlos and Liora).
            * Code (Liora and Carlos) **Sept 4th**.
            * Debugging (Liora and Carlos) **Sept 11th**
        * Introduction (Liora and Carlos) **SET DATE (Liora)**.
        * Results (Liora) **Sept 25th [tops]**.
        * Discussion (Liora and Carlos) **TBD**.
        * Conclusions (Liora and Carlos) **TBD**.
* Ms. York (will supervise the process).



$$P$$
$$M_i \in \{p, e, m\} => M_p, M_e, M_m$$
$$F_j \in \{h, l, s, a\} => F_h, F_l, F_s, F_a$$
$$ \alpha_{M_p F_s} $$

# Chart

![LiorasDiagram](https://i.imgur.com/4NkoGx6.jpg)


## Key for Chart
| Number | Symbol | Description |
| -------- | -------- | -------- |
| 1 + 5 + 6 + 17 + 18 + 19   |    i = small fish population, j = adult fish population, l = medusa population, $$g \in \{ polyp, ephyra, medusa ,plankton\}$$| $$\frac{dF_j}{dt}=\eta_{ji} F_i (1-\frac{F_j+\alpha_{F_jJ_l}J_l}{K_j})-\delta_{F_j}F_j-F_j\sum_{g}\alpha_{jg}J_g$$   |
| 1 + 5 + 6 + 17 + 18 + 19   |    i = small fish population, j = adult fish population, l = medusa population, $$g \in \{ polyp, ephyra, medusa ,plankton\}$$| $$\frac{dJ_j}{dt}=\eta_{ji} J_i (1-\frac{J_j+\alpha_{J_jF_l}F_l}{K_j})-\delta_{J_j}J_j-J_j\sum_{g}\alpha_{jg}F_g$$   |
|  | $$g \in \{ polyp, ephyra, medusa ,plankton \}$$ c = saturation constant  $$\eta_{ji_{max}}$$ = the max rate at which stage i transitions to stage j |$$\eta_{ji}= \eta_{ji_{max}} e^{ \left(  \frac{(\epsilon_iF_i\sum_{g}  \alpha_{gi} J_g) -c}{\epsilon_iF_i\sum_{g}  \alpha_{gi} J_g +c } +1\right)}$$ |
| 7 + 8 + 17   |    i = ephyra population, j = medusa population, l = adult fish population|$$\frac{dF_j}{dt}=\eta_{ji} F_i (1-\frac{F_j+\alpha_{F_jJ_l}J_l}{K_j})-\delta_{F_j}F_j$$   |
| 1 + 2 + 12 + 13   |   i = larvae population, j = small fish population, l = adult fish population,   $$ g \in\{m,n \}$$ |$$\frac{dF_j}{dt}=\eta_{ji} F_i(1-\frac{F_j}{Kj})-\eta_{lj}F_j-F_j\sum_{g}\alpha_{jg}J_g$$      |
|3 + 4 + 20 |i = adult fish population, j = egg population, h = larvae population, b = medusa population | $$\frac{dF_j}{dt}=\eta_{ji} F_i(1-\frac{F_j}{Kj})-\eta_{hj}F_j-\alpha_{jb}F_jJ_b$$ |
|2 + 3 + 21 |i = egg population, j = larvae population, s = small fish population, b = medusa population|$$\frac{dF_j}{dt}=\eta_{ji} F_i(1-\frac{F_j}{Kj})-\eta_{sj}F_j-\alpha_{jb}F_jJ_b$$ |
|9 + 10 + 11 + 12 + 14 + 19 |i = medusa population, j = polyp population, e = ephyra population population $$ g \in\{m,p \}$$|$$\frac{dF_j}{dt}=\sum_{g}\eta_{ji}F_i(1-\frac{F_j}{Kj})-\eta_{ej}F_j$$|
|8 + 9 + 13 + 15 + 18 |i = polyp population, j = ephyra population| $$\frac{dF_j}{dt}=\sum_{g}\eta_{ji}F_i(1-\frac{F_j}{Kj})-\eta_{bj}F_j$$ |


## General Equations:
- $$\frac{dF_j}{dt}=(\sum_{i}\eta_{ji}F_i)(1-\frac{F_j +\sum_l \beta_{jl} J_l}{Kj})-\eta_{bj}F_j-F_j\sum_{g}\sigma_{jg}J_g$$
- $$\frac{dJ_j}{dt}=(\sum_{i}\gamma_{ji}J_i)(1-\frac{J_j+\sum_l \beta_{jl} F_l}{Bj})-\gamma_{bj}J_j-J_j\sum_{g}\alpha_{jg}F_g$$

# Template model

```r
# logistic ####
# predict jellyfish population
library(deSolve)

Parms=list(a=15,e=1.5,K=30000,r2=1.5)
State=c(N=10000, R=30000)
Times=seq(from=0,to=5,by=0.1)

logistic=function(Time,State,Parms){# a,e,K,r2,N,R
  # resource ####
  # small fish
  Parms_t=unlist(Parms)
  a=Parms_t[1]
  e=Parms_t[2]
  K=Parms_t[3]
  r2=Parms_t[4]
  R=State[2]
  N=State[1]
  R=r2*R*(1-R/K)-a*N*R # r2 is growth of the resource, R is the resource population
  # K is carrying capacity, a is the consumption rate of the consumer, N is consumer population
  
  
  # consumer ####
  # jellyfish
  N=a*e*R*N*(1-N/R) # a is the consumption rate of the consumer, e is the efficiency rate
  # N is the concumer population, R is the resource population
  return(list(c(N=N,R=R)))
}
sim=ode(y=State,times=Times,func=logistic,parms=Parms) 
plot(x=sim[,1],y=sim[,3],col="blue",type="l")
plot(x=sim[,1],y=sim[,2],col="red",type="l")
```



# New code version


```
# logistic ####
# predict jellyfish population
library(deSolve)


State=c(initialpopsvalues)
Times=seq(from=0,to=10,by=0.1)

parameters=list(
  k = ,# Carrying capacity -> vector
  eta , # rate at which fish stage i becomes j  -> Matrix
  gamma_ = ,# rate at which jellyfish stage i becomes j  -> matrix
  alpha = , # Predation fish over jellyfish -> matrix
  sigma = , # Predation jellyfish over fish -> matrix
  constan = , # Saturation constant -> scalar
  beta_ = , # Competition coefficient between fish and jellyfish -> matrix
  efficiency =  # efficiency (rate of biomass conversion) -> scalar
    
)

logistic=function(Time,State,Parms){# a,e,K,r2,N,R
  # resource ####
  # small fish
  Parms_t=unlist(Parms)
  a=Parms_t[1] # call each variable a word
  e=Parms_t[2]
  K=Parms_t[3]
  r2=Parms_t[4]
  R=State[2]
  N=State[1]
  
  dFdt= (eta%*%F_i) * (1 - (F_i+ beta_%*%J_i)/(K_j)) - # new individuals per stage in Fish population
    eta%*%F_i - # Individuals leaving stage i
      F_i*(sigma_%*%J_i) # Number of fish that are consumed by any jelly fish stage
  
  
  # dJdt =  
  
  outcome = c(dFdt, dJdt) 
  names(dFdt)<-c('names for each stage')
  
  return(list(c(N=N,R=R)))
}
sim=ode(y=State,times=Times,func=logistic,parms=parameters)

```

**we are assuming that efficiency is the same
```
library(deSolve)
library(deldir)
library(tidyverse)
library(tensor)
library(deldir)
library(MASS)
library(abind)
library(future)
library(purrr)
library(furrr)
library(parallel)
library(tictoc)
library(future.apply)
trabajadores=detectCores()-1

source("Liora_Code2.R")
source("parallelparameters_Liora.R")


###############Parameter tibble #############

# parameter_tibble = # This line inputs the parameters for estimating other parameters for Liora model
#   expand_grid(
#     efficiencyf=0.9*c(0.1, 5.5, 10),
#     efficiencyj=0.05*c(0.1, 5.5, 10),
#     competitionfactorf=c(0.1, 5.5, 10),
#     competitionfactorj=c(0.1, 5.5, 10),
#     constant=1*c(0.1, 5.5, 10),
#     k=1000*c(0.1, 5.5, 10),
#     b=1000*c(0.1, 5.5, 10),
#     predation_rate_factorf=c(0.1, 5.5, 10),
#     predation_rate_factorj=c(0.1, 5.5, 10)
#   )

parameter_tibble = # This line inputs the parameters for estimating other parameters for Liora model
  expand_grid(
    efficiencyf=0.9,
    efficiencyj=0.005,
    competitionfactorf=seq(0, 1, 0.25),
    competitionfactorj=seq(0.001, 0.08, 0.02),
    constant=0.1,
    k=1000,
    b=1000,
    predation_rate_factorf=seq(from=0.4, to=0.9, by=0.1),
    predation_rate_factorj=seq(from=0.05, to=0.1, by=0.02)
  )



#############################################


#######################################################################################################

# efficiencyf_=efficiencyf,
# efficiencyj_=efficiencyj_,
# competitionfactorf_=competitionfactorf_,
# competitionfactorj_=competitionfactorj,
# constant_=constant,
# k_=k,
# b_=b,
# predation_rate_factorf_=predation_rate_factorf,
# predation_rate_factorj_=predation_rate_factorj,
# transitionf=Fish, # matrix
# transitionj=Jellyfish,
# predation_fish_=predation_fish,
# predation_jellyfish_=predation_jellyfish,
# competition_fish_to_medusa_=competition_fish_to_medusa,
# competition_medusa_to_fish_=competition_medusa_to_fish

#######################################################################################################

sim = 
  function(efficiencyf_=efficiencyf,
           efficiencyj_=efficiencyj,
           competitionfactorf_=competitionfactorf,
           competitionfactorj_=competitionfactorj,
           constant_=constant,
           k_=k,
           b_=b,
           predation_rate_factorf_=predation_rate_factorf,
           predation_rate_factorj_=predation_rate_factorj,
           state, t){
    
    
    settings=parallel_parameters(efficiencyf_=efficiencyf_,
                                 efficiencyj_=efficiencyj_,
                                 competitionfactorf_=competitionfactorf_,
                                 competitionfactorj_=competitionfactorj_,
                                 constant_=constant_,
                                 k_=k_,
                                 b_=b_,
                                 predation_rate_factorf_=predation_rate_factorf_,
                                 predation_rate_factorj_=predation_rate_factorj_
    )
    Params = settings$parameters
    State=settings$State
    out = ode(y = State, times = times, func = logistic, parms=Params)
    
    res_df = 
      as.data.frame(out) |>
      as_tibble()
    
    #names(res_df) = names_df$name
    
    return(list(outcome = res_df, parameters = Params))# as_tibble(
    
  }

times <- seq(0, 0.5, by =0.1)



results = 
  pmap(
    .l = parameter_tibble,
    .f = sim,
    t = timeseq
  )


source("sortingresults.R")

df=as_tibble(bind_rows(lapply(results, results_sorting)))

dtf=df |>
  mutate(totalfish=rowSums(dplyr::select(df, contains('fish'))))|>
  mutate(totalmedusa=rowSums(dplyr::select(df, c("polyp", "ephyra", "medusa"))))



plot_total_fish = 
  dtf |>
  #filter(totalfish>10)|>
  mutate(competitionfactor_ratio=as.factor(competitionfactorf/competitionfactorj))|>
  ggplot(aes(predation_rate_factorj,totalfish, group = predation_rate_factorf, color = as.factor(predation_rate_factorf))) +
  geom_line() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))  +
  facet_grid(competitionfactorf ~ competitionfactorj, labeller = label_both)

plot_total_fish |>
  show()


```