DemGenTrout 1.3

;; DemGenTrout 1.3 (2012)
;;
;; Developed by
;; Beatrice M. Frank
;; Earth and Life Institute
;; Universite catholique de Louvain
;; Beatrice.Frank@gmail.com
;; 
;; Last updated: September 2012

breed[ trout ]
breed[ dead ]
breed[ gone ]
breed[ stocked ]
breed[ hybrid ]

;; trout state variables
trout-own[ sex genotype natal-stream week-of-birth week-of-birth2 ;; fixed
  current-stream age stage status body-length body-weight condition-factor num-offspring
  spawned? moved-to-spawn? returned? ]

;; gone trout state variables
gone-own[ sex genotype natal-stream week-of-birth week-of-birth2 ;; fixed
  current-stream age stage status body-length body-weight condition-factor num-offspring
  spawned? moved-to-spawn? returned? ]

;; dead trout state variables
dead-own[ sex genotype natal-stream week-of-birth week-of-birth2 ;; fixed
  current-stream age stage status body-length body-weight condition-factor num-offspring
  spawned? moved-to-spawn? returned? ]

;; stocked state variables
stocked-own[ sex genotype natal-stream week-of-birth week-of-birth2 ;; fixed
  current-stream age stage status body-length body-weight condition-factor num-offspring
  spawned? moved-to-spawn? returned? newly-stocked? ]

;; hybrid state variables
hybrid-own[ sex genotype natal-stream week-of-birth week-of-birth2 ;; fixed
  current-stream age stage status body-length body-weight condition-factor num-offspring
  spawned? moved-to-spawn? returned? ]
  
;; stream state variables
patches-own [ stream flow temperature ] 

globals[ week week-of-year month-of-year year
 ;; trout initialisation
 prop-C prop-L init-N
 prop-age0-C prop-age1-C prop-age2-C prop-age3-C
 prop-age0-L prop-age1-L prop-age2-L
 meanl-age0-C meanl-age0-L meanl-age1-C meanl-age1-L meanl-age2-C meanl-age2-L meanl-age3-C meanl-age3-L
 sdl-age0-C sdl-age0-L sdl-age1-C sdl-age1-L sdl-age2-C sdl-age2-L sdl-age3-C sdl-age3-L 
 trdown-age0-C trdown-age0-L trdown-age1-C trdown-age1-L trdown-age2-C trdown-age2-L trdown-age3-C trdown-age3-L
 trup-age0-C trup-age0-L trup-age1-C trup-age1-L trup-age2-C trup-age2-L trup-age3-C trup-age3-L
 meanl-birth-C sdl-birth-C meanl-birth-L sdl-birth-L trdown-C trup-C trdown-L trup-L 
 mean-bcf-C sd-bcf-C mean-bcf-L sd-bcf-L
 propC-inL 
 ;; stream input data
 streamC-discharge streamL-discharge streamC-temperature streamL-temperature 
 mean-dischargeC sd-dischargeC mean-dischargeL sd-dischargeL
 mean-temperatureC sd-temperatureC mean-temperatureL sd-temperatureL
 max-dischargeC max-dischargeL

 ;; survival parameters 
 streamC-age0-survival streamC-age1-survival streamC-age2-survival streamC-age3-survival
 streamL-age0-survival streamL-age1-survival streamL-age2-survival streamL-age3-survival 
 streamC-age0-survival-list streamC-age1-survival-list streamC-age2-survival-list streamC-age3-survival-list
 streamL-age0-survival-list streamL-age1-survival-list streamL-age2-survival-list streamL-age3-survival-list  
 C-age0-survival C-age1-survival C-age2-survival C-age3-survival
 L-age0-survival L-age1-survival L-age2-survival L-age3-survival 
 predation-factor
 ;; growth parameters
 streamC-parK streamL-parK
 streamC-max-length streamC-parA streamC-parB streamL-max-length streamL-parA streamL-parB   
 ;; spawning parameters
 spawn-start spawn-end moved-prop spawn-mean-length spawn-mean-cond
 spawn-sd-length spawn-sd-cond spawn-min-age
 spawn-mean-temperature spawn-sd-temperature spawn-mean-flow spawn-sd-flow
 streamC-capacity offprod offprod-min-length offprod-max-length offprod-min
 offspringC-varA offspringC-varB offspringL-varA offspringL-varB
 offprodC-max offprodL-max 
 length-heritability 
 ;; hatching parameters
 hatch-start hatch-end
 ;; downstream movement parameters
 move-start move-end move-min-age move-max-age
 move-mean-length move-sd-length 
 move-age1-prob move-age1-varB move-age1-varA
 move-age2-prob move-age2-varB move-age2-varA 
 move-mean-temperature move-sd-temperature move-mean-flow move-sd-flow
 ;; leaving stream L forever parameters
 leaving-propC leaving-propL
 ;; stocking
; stocking-start stocking-end stocking-coefficient survival-factor
 num-stocked hybrid-spawning-prob ; trout-spawning-prob stocked-spawning-prob 

 out-year out-month out-week
]  

to setup
  clear-all      
  no-display     
  ;; controlling the random numbers
  random-seed 1223251200  
  
  ;; first week of the year = 1st to 7th October included
  set week 1
  set week-of-year 1
  set month-of-year 1
  set year 1
  
  import-world "DemGenTrout-world.csv"
  ask trout [set genotype lput "000/000" genotype]
  
  reset-ticks
  restore-defaults
end 

to restore-defaults
  ;; restore the default parameters linked to stocking
  set stocking-start 1
  set stocking-end 10
  set stocking-coefficient 0.50     
  set survival-factor 10  
  set trout-spawning-prob 1.00
  set stocked-spawning-prob 0.10  
  ;; compute spawning probability for hybrids
  set hybrid-spawning-prob (trout-spawning-prob * (1 - ((1 - stocked-spawning-prob) / 2)))    
end 

to set-parameters  
  ;; compute spawning probability for hybrids
  set hybrid-spawning-prob (trout-spawning-prob * (1 - ((1 - stocked-spawning-prob) / 2)))   
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;; SCHEDULING ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to go
  no-display 

  tick 
  ;; end at year 35
  if ticks = 1820 [ user-message "Input data are only available for 35 years."  stop ]
  ;; a year = 52 weeks  
  set week-of-year (week mod 52) + 1
  set month-of-year (int(week-of-year / 4.5) mod 52) + 1 
  set year (int(week / 52)) + 1 ; mod 52
  
  ;; each week

;-- 1. Update stream hydrological conditions ----------------------------------------------------------
  
  ask patches [ update-stream-conditions ] 
   
;-- New. Introduce stocked trout in the system --------------------------------------------------------

  ;; introduction of stocked trout
  if week-of-year = 27 and year >= stocking-start and year <= stocking-end [ intro-stocked-trout ]   
  
let selected-breed (list trout stocked hybrid)
foreach selected-breed 
[    
  
;-- 2. Kill trout in each stream ----------------------------------------------------------------------

  ask ? [ trout-die ]         
    
;-- 3. Update trout length, weight and condition factor -----------------------------------------------
  
  ask ? [ trout-grow ]      

;-- 4. Reproduce trout in each stream ----------------------------------------------------------------- 

  ;; spawning time window: duration = 11 weeks
  if week-of-year >= spawn-start and week-of-year <= spawn-end
  [   
    ;--- 4.1. Identify candidate spawners ---;
    ask ? [ trout-become-candidate-spawners ] 
    ;; select some of the candidate spawners born in stream L in order to make them reproduce in C
    let candidate-spawners-streamL ? with [ status = "candidate-spawner" and current-stream = "L" 
      and natal-stream = "L" and not moved-to-spawn? ]
    while [ (count ? with [status = "selected-candidate-spawner"]) < (moved-prop * count candidate-spawners-streamL) ]
    [ ask one-of candidate-spawners-streamL [ set status "selected-candidate-spawner" ] ]      
    
    ;--- 4.2. Move candidate spawners upstream ---;   
    ask ? [ trout-move-to-spawn ]         
    
    ;; adapt offspring production if abundance in C is strictly higher than streamC-capacity
    ifelse count ? with [ current-stream = "C" ] > streamC-capacity 
    [ set offprodC-max offprod / 10 ]
    [ set offprodC-max offprod ]        

    ;--- 4.3. Produce offspring in each stream ---;  
    trout-spawn-in-streamC
    trout-spawn-in-streamL
  ]
   
;-- 5. Increment trout age and update stage -----------------------------------------------------------

  if week-of-year = (hatch-start - 1) ;; one week before hatching
  [ 
    ask ? [ trout-update-age ]
    ask ? with [age > 6] [die!] ;; trout die at age 7    
  ]       

;-- 6. Reveal offspring -------------------------------------------------------------------------------   
  
  ;; hatching time window: duration = 11 weeks
  if week-of-year >= hatch-start and week-of-year <= hatch-end
  [ (foreach sort ? with [age = -1]
      [ ask ?
        [
          if (week-of-year - week-of-birth) = (hatch-start - spawn-start) ;; delay = 10 weeks
          [ set hidden? false set age 0 ]
        ]])
  ]

;-- 7. Move trout of stream C downstream ----------------------------------------------------------
    
  ;; migration time window: duration = 25 weeks
  if week-of-year >= move-start and week-of-year <= move-end
  [             
    ;--- 7.1. Identify candidate migrants among juveniles ---; 
    ask ? [ trout-become-candidate-migrants ]          
     
    ;--- 7.2. Update movement probabilities ---;     
    update-move-age1-prob   
    update-move-age2-prob      
            
    ;--- 7.3. Move candidate migrants ---;
    ask ? [ trout-move-downstream ]                  
  ]        

;-- 8. Move upswimming spawners back to stream L ----------------------------------------------------- 

  if week-of-year = (spawn-end + 1) ;; the week after the spawning process
  [ 
    ask ? [ trout-move-back ]
  ]               
  
;-- 9. Remove young trout of stream L from the system ------------------------------------------------    
  
  ;; at the end of the year
  if week-of-year = 52
  [ 
    ;; some juveniles born in C and living in stream L do not settle and go elsewhere 
    let bornC-in-streamL ? with [current-stream = "L" and natal-stream = "C" and status = "migrant"]    
    let bornC-leaving (leaving-propC * count bornC-in-streamL)
    ask n-of bornC-leaving bornC-in-streamL [ leave-streamL ]  
    ;; some juveniles born in L and living in stream L do not settle and go elsewhere  
    let bornL-in-streamL ? with [current-stream = "L" and natal-stream = "L" and (age = 1 or age = 2)]
    let bornL-leaving (leaving-propL * count bornL-in-streamL)
    ask n-of bornL-leaving bornL-in-streamL [ leave-streamL ]   
  ]                
  
; Trout reset --------------------------------------------------------------------------------------
  
  ;; at the end of the year
  if week-of-year = 52
  [      
    ;; reset trout status and boolean state variables    
    ask ? [ set status "non-spawner" set moved-to-spawn? false set spawned? false set returned? false ] 
    ;; kill all hidden turtles  
    ask gone [die]     
    ask dead [die]          
  ]
]

; Increment the week -------------------------------------------------------------------------------
   
  set week week + 1
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;; STREAM PROCEDURE ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; update water temperature and flow of each stream

to update-stream-conditions
  ifelse ticks < (6 * 52)
  ;; first 6 years (2004 to 2009)
  [
    ifelse stream = "C"
    [
      set flow item week streamC-discharge
      set temperature item week streamC-temperature 
    ]
    [
      set flow item week streamL-discharge
      set temperature item week streamL-temperature
    ]  
  ]
  ;; next years (2010 to 2038)
  [
    ifelse stream = "C"
    [
      set flow min (list max-dischargeC (item week streamC-discharge + random-lognormal mean-dischargeC sd-dischargeC))
      set temperature (item week streamC-temperature + random-normal mean-temperatureC sd-temperatureC) 
    ]
    [
      set flow min (list max-dischargeL (item week streamL-discharge + random-lognormal mean-dischargeL sd-dischargeL))
      set temperature (item week streamL-temperature + random-normal mean-temperatureL sd-temperatureL)
    ]  
  ]
end 

;;;;;;;;;;;;;;;;;;;;;;;;;; TROUT PROCEDURES ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; trout mortality

to trout-die
  ;; in stream C
  if (breed = trout or breed = stocked) and stage = "fry" and current-stream = "C" [if random-float 1 > ((item week-of-year streamC-age0-survival) ^ (1 / 52)) [die!]]  
  if breed = hybrid and stage = "fry" and current-stream = "C" [if random-float 1 > (((item week-of-year streamC-age0-survival) * (1 - ((1 - (item week-of-year streamC-age0-survival)) / 2))) ^ (1 / 52)) [die!]]
  if (breed = trout or breed = stocked) and stage = "juvenile" and age = 1 and current-stream = "C" [if random-float 1 > ((item week-of-year streamC-age1-survival) ^ (1 / 52)) [die!]]
  if breed = hybrid and stage = "juvenile" and age = 1 and current-stream = "C" [if random-float 1 > (((item week-of-year streamC-age1-survival) * (1 - ((1 - (item week-of-year streamC-age1-survival)) / 2))) ^ (1 / 52)) [die!]]  
  if (breed = trout or breed = stocked) and stage = "juvenile" and age = 2 and current-stream = "C" [if random-float 1 > ((item week-of-year streamC-age2-survival) ^ (1 / 52)) [die!]]
  if breed = hybrid and stage = "juvenile" and age = 2 and current-stream = "C" [if random-float 1 > (((item week-of-year streamC-age2-survival) * (1 - ((1 - (item week-of-year streamC-age2-survival)) / 2))) ^ (1 / 52)) [die!]]    
  if (breed = trout or breed = stocked) and stage = "adult" and age < 6 and current-stream = "C" [if random-float 1 > ((item week-of-year streamC-age3-survival) ^ (1 / 52)) [die!]] 
  if breed = hybrid and stage = "adult" and age < 6 and current-stream = "C" [if random-float 1 > (((item week-of-year streamC-age3-survival) * (1 - ((1 - (item week-of-year streamC-age3-survival)) / 2))) ^ (1 / 52)) [die!]]    
  if week-of-year >= spawn-start and week-of-year <= spawn-end and (breed = trout or breed = stocked) and moved-to-spawn? = true and current-stream = "C" 
   [if random-float 1 > (((item week-of-year streamC-age3-survival) / predation-factor) ^ (1 / (spawn-end - spawn-start))) [die!]] 
  if week-of-year >= spawn-start and week-of-year <= spawn-end and breed = hybrid and moved-to-spawn? = true and current-stream = "C" 
   [if random-float 1 > ((((item week-of-year streamC-age3-survival) * (1 - ((1 - (item week-of-year streamC-age3-survival)) / 2))) / predation-factor) ^ (1 / (spawn-end - spawn-start))) [die!]]       
  ;; in stream L
  if (breed = trout or breed = stocked) and stage = "fry" and current-stream = "L" [if random-float 1 > ((item week-of-year streamL-age0-survival) ^ (1 / 52)) [die!]]
  if breed = hybrid and stage = "fry" and current-stream = "L" [if random-float 1 > (((item week-of-year streamL-age0-survival) * (1 - ((1 - (item week-of-year streamL-age0-survival)) / 2))) ^ (1 / 52)) [die!]]  
  if (breed = trout or breed = stocked) and stage = "juvenile" and age = 1 and current-stream = "L" [if random-float 1 > ((item week-of-year streamL-age1-survival) ^ (1 / 52)) [die!]]
  if breed = hybrid and stage = "juvenile" and age = 1 and current-stream = "L" [if random-float 1 > (((item week-of-year streamL-age1-survival) * (1 - ((1 - (item week-of-year streamL-age1-survival)) / 2))) ^ (1 / 52)) [die!]]  

  if breed = stocked and newly-stocked? = true and stage = "juvenile" and age = 2 and current-stream = "L" [if random-float 1 > (((item week-of-year streamL-age2-survival) / survival-factor) ^ (1 / 52)) [die!]]  
  if (breed = trout or breed = stocked) and stage = "juvenile" and age = 2 and current-stream = "L" [if random-float 1 > ((item week-of-year streamL-age2-survival) ^ (1 / 52)) [die!]]
  if breed = hybrid and stage = "juvenile" and age = 2 and current-stream = "L" [if random-float 1 > (((item week-of-year streamL-age2-survival) * (1 - ((1 - (item week-of-year streamL-age2-survival)) / 2))) ^ (1 / 52)) [die!]]
  
  if (breed = trout or breed = stocked) and stage = "adult" and age < 6 and current-stream = "L" [if random-float 1 > ((item week-of-year streamL-age3-survival) ^ (1 / 52)) [die!]]  
  if breed = hybrid and stage = "adult" and age < 6 and current-stream = "L" [if random-float 1 > (((item week-of-year streamL-age3-survival) * (1 - ((1 - (item week-of-year streamL-age3-survival)) / 2))) ^ (1 / 52)) [die!]]    
end 

;; individual growth

to trout-grow
  ifelse current-stream = "C"
  [
    ;; von Bertalanffy growth equation for trout in stream C
    set body-length (body-length + (streamC-parK * (streamC-max-length - body-length)))
    let healthy-weight (streamC-parA * ((body-length) ^ streamC-parB)) ;; weight-length regression 
    set body-weight (condition-factor * healthy-weight)      
  ]
  [
    ;; von Bertalanffy growth equation for trout in stream L
    set body-length (body-length + (streamL-parK * (streamL-max-length - body-length)))  
    let healthy-weight (streamL-parA * ((body-length) ^ streamL-parB)) ;; weight-length regression
    set body-weight (condition-factor * healthy-weight)      
  ]
  set condition-factor max (list 0.5 (condition-factor + (- 0.008 + random-float 0.016)))
  ;; test for negative length and negative weight
  if body-length < 0 [ user-message "Negative fish length" ]
  if body-weight < 0 [ user-message "Negative fish weight" ]    
end 

to die!
  set breed dead
  hide-turtle
end 

;; identification of candidate spawners

to trout-become-candidate-spawners
  if age >= spawn-min-age and status != "spawner" and not moved-to-spawn? 
  [ if random-normal spawn-mean-length spawn-sd-length < body-length and 
    random-normal spawn-mean-cond spawn-sd-cond < condition-factor
    [ set status "candidate-spawner" ] ]  
end 

;; movement to stream C for reproduction

to trout-move-to-spawn
 let temperature-streamL [temperature] of one-of patches with [ stream = "L" ]      
 let flow-streamL [flow] of one-of patches with [ stream = "L" ] 
 ;; spawners born in C or in L move to C
 if (status = "candidate-spawner" and not moved-to-spawn? and current-stream = "L" and natal-stream = "C") 
   or (status = "selected-candidate-spawner" and not moved-to-spawn?)     
 [ if random-lognormal spawn-mean-flow spawn-sd-flow < flow-streamL   
   [ set current-stream "C" set xcor 1 set moved-to-spawn? true set status "candidate-spawner" ] ]  
end 

;; update trout age and stage

to trout-update-age
  set age age + 1
  trout-update-stage 
end   

;; update trout stage: fry, juvenile, adult 

to trout-update-stage
  if age < 1 [ set stage "fry" stop ]
  if age >= 1 and age <= 2 [ set stage "juvenile"stop ]
  if age > 2 [ set stage "adult" stop ] 
end 

;; reproduction

to trout-spawn-in-streamC
  let trout-spawners trout with [ status = "candidate-spawner" or status = "spawner" and current-stream = "C" ] 
  let stocked-spawners-temp stocked with [ status = "candidate-spawner" or status = "spawner" and current-stream = "C" ]  
  let hybrid-spawners-temp hybrid with [ status = "candidate-spawner" or status = "spawner" and current-stream = "C" ]
  let num-stocked-spawners (floor ((count stocked-spawners-temp) * stocked-spawning-prob))
  let num-hybrid-spawners (floor ((count hybrid-spawners-temp) * hybrid-spawning-prob))  
  let stocked-spawners n-of num-stocked-spawners stocked-spawners-temp
  let hybrid-spawners n-of num-hybrid-spawners hybrid-spawners-temp  
  let spawners (turtle-set trout-spawners stocked-spawners hybrid-spawners) 
  ifelse count spawners < 2
  [ stop ]
  [  
    ;; length of parents when born (VBLF backwards)
    let born-length-parents VBLF-reverse spawners    
    ;; mean and variance of length of parents when born
    if length born-length-parents < 2 [ stop ]
    let mean-parents mean born-length-parents
    let var-parents variance born-length-parents
    let num-crossover ((count trout-spawners + count stocked-spawners + count hybrid-spawners) / 2)  
    ;; polygamous mating with satellite males (not in good condition)    
    repeat num-crossover
    [
      let num-males 1 + random 4 ; number of males drawn randomly from an uniform distribution from 1 to 4
      let all-parent-M spawners with [ sex = "M" ]
      if count all-parent-M < num-males [ stop ]         
      let no-parent-M max-n-of num-males spawners with [ sex = "M" ] [condition-factor]
      let sel-parent-M all-parent-M with [ not member? self no-parent-M and not spawned? ]  
      if count spawners with [ sex = "F" and not spawned? ] <= 1 [ stop ] 
      if count sel-parent-M < num-males [ stop ]  
      ;; select the parents: first, the spawners that have moved to spawn, then the other spawners  
      let parent-M n-of num-males sel-parent-M
      let parent-F max-one-of spawners with [ sex = "F" and not spawned? ] [condition-factor]     
      if (count spawners with [ sex = "M" and not spawned? and moved-to-spawn? ] - num-males) >= num-males  
      [ 
        if any? spawners with [ sex = "M" and not spawned? and moved-to-spawn? ] 
        [
          set all-parent-M spawners with [ sex = "M" and moved-to-spawn? ]
          set no-parent-M max-n-of num-males spawners with [ sex = "M" and moved-to-spawn? ] [condition-factor]
          set sel-parent-M all-parent-M with [ not member? self no-parent-M and not spawned? ]           
          set parent-M n-of num-males sel-parent-M
        ]
      ]
      if any? spawners with [ sex = "F" and not spawned? and moved-to-spawn?]
      [ set parent-F max-one-of spawners with [ sex = "F" and not spawned? and moved-to-spawn?] [condition-factor] ]
      ;; length deviances of parents M and F from mean-parents
      let dummy VBLF-reverse parent-M
      let dev-M 0
      let temp (random-normal meanl-birth-C sdl-birth-C)
      if temp < trdown-C [set temp trdown-C]   
      if temp > trup-C [set temp trup-C]                
      if length dummy = 0 [ set dev-M (mean-parents - temp) ]      
      if length dummy = 1 [ set dev-M (mean-parents - (item 0 dummy)) ]        
      if length dummy > 1 [ set dev-M (mean-parents - mean dummy) ]
      let dev-F (mean-parents - (VBLF-reverse parent-F))       
      ask parent-F
      [
        set spawned? true
        set status "spawner"
        set xcor 1
        let num-fry update-num-offspringC [body-length] of parent-F   
        hatch-trout num-fry
        [
          set hidden? true ;; hide trout until hatching
          set age -1 set stage "fry"  
          ;; genotype of parent-M is a mixture of the several genotypes of males
          let genotype-mix mixture ([genotype] of parent-M)
          ;; define genotypes for eggs
          let egg-genotype cross-genotype (genotype-mix) ([genotype] of parent-F)  
          set genotype egg-genotype   
          set breed item 0 assign-breed egg-genotype                                           
          set xcor 1
          set current-stream "C"
          set natal-stream "C"
          set num-offspring 0
          set week-of-birth week-of-year      
          set week-of-birth2 (week + 1)          
          ifelse random-float 1 > 0.5 [ set sex "M" ] [ set sex "F" ]
          set status "non-spawner" set spawned? false
          set moved-to-spawn? false set returned? false          
          ;; set body-length as function of parents' length and environment
          let genetic-dev (sqrt length-heritability) * (dev-M + dev-F)    
          let envi-var ((1 - length-heritability) ^ 2) * var-parents
          let envi-dev random-normal 0 (sqrt envi-var) 
          set body-length mean-parents + genetic-dev + envi-dev
          let temp2 (random-normal meanl-birth-C sdl-birth-C)  
          if temp2 < trdown-C [set temp2 trdown-C]   
          if temp2 > trup-C [set temp2 trup-C]                
          if body-length < trdown-C or body-length > trup-C [ set body-length temp2 ]           
          ;; set body-weight and condition-factor
          set condition-factor (random-normal mean-bcf-C sd-bcf-C)         
          set body-weight (condition-factor * ((body-length) ^ 3) / 100000)       
        ] 
        set num-offspring num-fry
      ]
      ask parent-M
      [
        ;set spawned? true ;; shut because males can reproduce more than once
        set status "spawner" 
        set xcor 1        
        set num-offspring [num-offspring] of parent-F               
      ] 
    ]
  ]               
end  

to trout-spawn-in-streamL
  let trout-spawners trout with [ status = "candidate-spawner" or status = "spawner" and current-stream = "L" ]
  let stocked-spawners-temp stocked with [ status = "candidate-spawner" or status = "spawner" and current-stream = "L" ]
  let hybrid-spawners-temp hybrid with [ status = "candidate-spawner" or status = "spawner" and current-stream = "L" ]
  let num-stocked-spawners (floor ((count stocked-spawners-temp) * stocked-spawning-prob))
  let num-hybrid-spawners (floor ((count hybrid-spawners-temp) * hybrid-spawning-prob))  
  let stocked-spawners n-of num-stocked-spawners stocked-spawners-temp
  let hybrid-spawners n-of num-hybrid-spawners hybrid-spawners-temp    
  let spawners (turtle-set trout-spawners stocked-spawners hybrid-spawners) 
  ifelse count spawners < 2
  [ stop ]
  [  
    let born-length-parents VBLF-reverse spawners 
    ;; mean and variance of length of parents when born
    if length born-length-parents < 2 [ stop ]  
    let mean-parents mean born-length-parents
    let var-parents variance born-length-parents  
    let num-crossover ((count trout-spawners + count stocked-spawners + count hybrid-spawners) / 2)  
    ;; polygamous mating with satellite males (not in good condition)    
    repeat num-crossover  
    [
      let num-males 1 + random 4 ; number of males drawn randomly from an uniform distribution from 1 to 4
      let all-parent-M spawners with [ sex = "M" ]
      if count all-parent-M < num-males [ stop ]         
      let no-parent-M max-n-of num-males spawners with [ sex = "M" ] [condition-factor]
      let sel-parent-M all-parent-M with [ not member? self no-parent-M and not spawned? ]  
      if count spawners with [ sex = "F" and not spawned? ] <= 1 [ stop ] 
      if count sel-parent-M < num-males [ stop ]  
      ;; select the parents: first, the spawners that have moved to spawn, then the other spawners  
      let parent-M n-of num-males sel-parent-M
      let parent-F max-one-of spawners with [ sex = "F" and not spawned? ] [condition-factor]     
      if (count spawners with [ sex = "M" and not spawned? and moved-to-spawn? ] - num-males) >= num-males  
      [ 
        if any? spawners with [ sex = "M" and not spawned? and moved-to-spawn? ] 
        [
          set all-parent-M spawners with [ sex = "M" and moved-to-spawn? ]
          set no-parent-M max-n-of num-males spawners with [ sex = "M" and moved-to-spawn? ] [condition-factor]
          set sel-parent-M all-parent-M with [ not member? self no-parent-M and not spawned? ]           
          set parent-M n-of num-males sel-parent-M
        ]
      ]
      if any? spawners with [ sex = "F" and not spawned? and moved-to-spawn?]
      [ set parent-F max-one-of spawners with [ sex = "F" and not spawned? and moved-to-spawn?] [condition-factor] ]
      ;; length deviances of parents M and F from mean-parents
      let dummy VBLF-reverse parent-M
      let dev-M 0
      let temp (random-normal meanl-birth-L sdl-birth-L)
      if temp < trdown-L [set temp trdown-L]   
      if temp > trup-L [set temp trup-L]                
      if length dummy = 0 [ set dev-M (mean-parents - temp) ]            
      if length dummy = 1 [ set dev-M (mean-parents - (item 0 dummy)) ]        
      if length dummy > 1 [ set dev-M (mean-parents - mean dummy) ]    
      let dev-F (mean-parents - (VBLF-reverse parent-F))      
      ask parent-F
      [
        set spawned? true
        set status "spawner"
        set xcor 0
        let num-fry update-num-offspringL [body-length] of parent-F 
        hatch-trout num-fry
        [
          set hidden? true ;; hide trout until hatching
          set age -1 set stage "fry"
          ;; genotype of parent-M is a mixture of the several genotypes of males
          let genotype-mix mixture ([genotype] of parent-M)
          ;; define genotypes for eggs
          let egg-genotype cross-genotype (genotype-mix) ([genotype] of parent-F)  
          set genotype egg-genotype   
          set breed item 0 assign-breed egg-genotype                                 
          set xcor 0
          set current-stream "L"
          set natal-stream "L"  
          set num-offspring 0   
          set week-of-birth week-of-year      
          set week-of-birth2 (week + 1)              
          ifelse random-float 1 > 0.5 [ set sex "M" ] [ set sex "F" ]    
          set status "non-spawner" set spawned? false
          set moved-to-spawn? false set returned? false                   
          ;; set body-length as function of parents' length and environment
          let genetic-dev (sqrt length-heritability) * (dev-M + dev-F)     
          let envi-var ((1 - length-heritability) ^ 2) * var-parents
          let envi-dev random-normal 0 (sqrt envi-var) 
          set body-length mean-parents + genetic-dev + envi-dev
          let temp2 (random-normal meanl-birth-L sdl-birth-L)  
          if temp2 < trdown-L [set temp2 trdown-L]   
          if temp2 > trup-L [set temp2 trup-L]                         
          if body-length < trdown-L or body-length > trup-L [ set body-length temp2 ]       
          ;; set body-weight and condition-factor
          set condition-factor (random-normal mean-bcf-L sd-bcf-L)
          set body-weight (condition-factor * ((body-length) ^ 3) / 100000)              
        ]
        set num-offspring num-fry
      ]
      ask parent-M
      [
        ;set spawned? true ;; shut because males can reproduce more than once
        set status "spawner"
        set xcor 0      
        set num-offspring [num-offspring] of parent-F
      ]       
    ]  
  ]      
end   

;; stream C

to-report update-num-offspringC [ female-length ] ;; linear function
  ;; compute the intermediate variables for the linear offspring production function for trout of stream C
  set offspringC-varA (offprodC-max - offprod-min) / (offprod-max-length - offprod-min-length)
  set offspringC-varB offprod-min - (offspringC-varA * offprod-min-length) 
  report ((offspringC-varA * female-length) + offspringC-varB)
end 

;; stream L

to-report update-num-offspringL [ female-length ] ;; linear function
  ;; compute the intermediate variables for the linear offspring production function  for trout of stream C
  set offspringL-varA (offprodL-max - offprod-min) / (offprod-max-length - offprod-min-length)
  set offspringL-varB offprod-min - (offspringL-varA * offprod-min-length)    
  report ((offspringL-varA * female-length) + offspringL-varB)
end 

to-report cross-genotype [ genotype1 genotype2 ]
  let new-genotype []
  let bit 0
  repeat length genotype1
  [
    let choice random 4
    if choice = 0
    [ set new-genotype (lput (word(substring (item bit genotype1) 0 4)(substring (item bit genotype2) 0 3)) new-genotype)]
    if choice = 1
    [ set new-genotype (lput (word(substring (item bit genotype1) 0 3) (substring (item bit genotype2) 3 7)) new-genotype)]
    if choice = 2
    [ set new-genotype (lput (word(substring (item bit genotype2) 0 3) (substring (item bit genotype1) 3 7)) new-genotype)]
    if choice = 3
    [ set new-genotype (lput (word(substring (item bit genotype1) 4 7) (substring (item bit genotype2) 3 7)) new-genotype)]      
    set bit bit + 1
  ]
  report new-genotype
end 

to-report mixture [ genotype-temp ]
  let genotype-mix []
  let bit 0
  repeat length one-of genotype-temp
  [
    let choice random (length genotype-temp) ; choose randomly one genotype among n males
    set genotype-mix lput item bit (item choice genotype-temp) genotype-mix
    set bit bit + 1
  ] 
  report genotype-mix
end 

to-report VBLF-reverse [ trout-indiv ]  
  ifelse is-agentset? trout-indiv
  [
    let length-when-born []
    (foreach sort trout-indiv
      [
        let time 0
        let length-diff 0
        let length-when-born-temp 0
        ask ?
        [
          set time (week - week-of-birth2)
          ifelse natal-stream = "C" 
          [ 
            set length-diff (streamC-max-length - body-length)
            set length-when-born-temp (streamC-max-length - (length-diff / ((1 - streamC-parK) ^ (time))))            
          ]
          [ 
            set length-diff (streamL-max-length - body-length)
            set length-when-born-temp (streamL-max-length - (length-diff / ((1 - streamL-parK) ^ (time))))            
          ]
        ]     
        set length-when-born lput length-when-born-temp length-when-born                
      ])
    report length-when-born 
  ]  
  [
    let time (week - ([week-of-birth2] of trout-indiv))      
    ifelse [natal-stream] of trout-indiv = "C"
    [ 
      let length-when-born (streamC-max-length - ((streamC-max-length - [body-length] of trout-indiv) / ((1 - streamC-parK) ^ (time)))) 
      report length-when-born
    ]
    [ 
      let length-when-born (streamL-max-length - ((streamL-max-length - [body-length] of trout-indiv) / ((1 - streamL-parK) ^ (time)))) 
      report length-when-born  
    ]   
  ]   
end 

to trout-become-candidate-migrants
  if (age = move-min-age or age = move-max-age) and current-stream = "C" and natal-stream = "C" 
  [ if random-normal move-mean-length move-sd-length < body-length
    [ set status "candidate-migrant" ] ]
end 

;; movement from C to L

to trout-move-downstream
  let temperature-streamC [temperature] of one-of patches with [ stream = "C" ]  
  let flow-streamC [flow] of one-of patches with [ stream = "C" ]    
  if status = "candidate-migrant"
  [ if random-normal move-mean-temperature move-sd-temperature < temperature-streamC 
    and random-lognormal move-mean-flow move-sd-flow < flow-streamC
    [ ifelse age = 1
      ;; juveniles age 1
      [ if random-float 1 < move-age1-prob 
        [ set xcor 0 set current-stream "L" set status "migrant" ] ]
      ;; juveniles age 2
      [ if random-float 1 < move-age2-prob
        [ set xcor 0 set current-stream "L" set status "migrant" ] ]
    ] 
  ]   
end     

to update-move-age1-prob ;; logistic function
  let nb1 count trout with [current-stream = "C" and body-length > 70]
  let Z1 (exp((move-age1-varA * nb1) + move-age1-varB))
  set move-age1-prob (Z1 / (1 + Z1))  
end 

to update-move-age2-prob ;; logistic function
  let nb2 count trout with [current-stream = "C" and body-length > 70]
  let Z1 (exp((move-age2-varA * nb2) + move-age2-varB))
  set move-age2-prob (Z1 / (1 + Z1))  
end 

to trout-move-back
  ;; candidate spawners return to previous stream
  if (status = "candidate-spawner" and current-stream = "C" and moved-to-spawn?)
  ;; effective spawners return to previous stream
    or (status = "spawner" and current-stream = "C" and moved-to-spawn?)
  [ set current-stream "L" set xcor 0 set returned? true ]  
end 

to leave-streamL
  set breed gone
  hide-turtle
end 

to intro-stocked-trout
  if year = 1
  [
    let num-trout count trout with [current-stream = "L"]
    set num-stocked round (stocking-coefficient * num-trout)
    ;print num-stocked
  ]
  create-stocked num-stocked
  [ 
    set breed stocked
    set newly-stocked? true
    setxy random-xcor -15
    set shape "box" set size 0.8
    set current-stream "L"
    set natal-stream "X"
    ifelse random-float 1 > 0.5
      [ set sex "M" ]
      [ set sex "F" ]
    ;; attribute trout genotypes at 12 microsatellites markers (loci)     
    let markers []
    file-open "input-genotypes-stocking.txt"
    set markers (sentence file-read (list ["999/999"]))  
    file-close
    let freq []
    file-open "input-genofreq-stocking.txt"
    set freq (sentence file-read (list [1]))
    file-close    
    set genotype []          
    set genotype assign-genotype (markers) (freq)
    set num-offspring 0     
    set status "non-spawner"
    set spawned? false set moved-to-spawn? false set returned? false      
    set age 2 
    trout-update-stage ;; a trout procedure setting stage from age  
    ;; attribute body-length, body-weight and condition-factor
    set body-length random-normal 248.47 26.71 set condition-factor 1.1 
    set body-weight (condition-factor * ((body-length) ^ 3) / 100000)     
    ;; test for negative length and negative weight
    if body-length < 0 [ user-message "Negative fish length" ]
    if body-weight < 0 [ user-message "Negative fish weight" ]               
  ] ;; end of create-stocked in stream L  
end 

to-report assign-genotype [markers freq]
  let locus []
  let bit 0
  repeat length markers
  [
    set locus assign-locus (bit) (markers) (freq)     
    set genotype lput (first locus) genotype  
    set bit bit + 1
  ] 
  report genotype
end 

to-report assign-locus [bit markers freq]
    let locus []
    let q 0
    let ran random-float 1.0
    (foreach (item bit markers) (item bit freq)
    [
      set q q + ?2
      if ran <= q [set locus lput ?1 locus] 
      if empty? locus and ?1 = last (item bit markers)
      [ set locus assign-locus (bit) (markers) (freq) ]
    ])    
    report locus
end 

;; breed attribution to offspring in function of 13th loci in genotype

to-report assign-breed [ genotype-temp ]
  ;; wild trout
  if (last genotype-temp) = "000/000" [ report (list trout)  ]
  ;; stocked trout 
  if (last genotype-temp) = "999/999" [ report (list stocked)]
  ;; hybrids
  if (last genotype-temp) = "000/999" [ report (list hybrid) ]
  if (last genotype-temp) = "999/000" [ report (list hybrid) ]
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to-report random-lognormal [mu sigma]
  let S sqrt(ln(((sigma / mu) ^ 2) + 1))
  let M ln (mu) - ((S ^ 2) / 2)
  report exp(random-normal M S)
end