Social Network of Sensors

extensions [profiler]
breed [ sensors sensor ]
breed [ events event ]
globals [
  current-showing ; current sensor layer showed
  nmessages ; number of exchanged messages
  dist ; The distance of spreading from the center of the environment
  excor ; Save x cordinate of the sensors that get the event
  eycor ; Save y cordinate of the sensors that get the event
  n-degree
  total ; for test purposes
  CDT ;; Cumulative Dynamic Threshold
  strong-estimate ; The estimated number of strong ties
  weak-estimate ; The estimated number of weak ties
  reported ; flag to stop the simulation if the event has been reported to the sink
  training ; flag to let the simulation run if we are under the training time
]
patches-own [
  ;; keep track of how many times the patch has been covered
  ;; that is, how many times fell in sensor radius (but might not be visited!)
  overlapping
]
;; these variables are available for sensors and events
turtles-own [
  ;; virtual coordinates
  home-x home-y current-x current-y
]
;; sensor specific variables
sensors-own [
  jump-size ;; length of the jump
  waiting-time ;; time to wait before the next jump
  mobile ;; True if the sensor is mobile
  sink ;; True if the sensor is the sink
  found ;; True if has knowledge of the event
  sensor-radius ;; radius of the sensor
  displacement ;; sum of displacements from the home position
  Pnew ;; probability of a new jump. We need this for plotting
  visited-locations ;; list of visited locations with the frequency
  S ;; keep track of the number of disctinct visited locations
  sum-of-frequencies ;; keep track of the sum of the visited location (number of jumps)
  allmeet
  ;; Time of Last Encounter list
  ;; in Gradient (FRESH) routing is the time of last encounter of the sink
  ;; in DoNothing is the time of last encounter of other sensors
  TLE
  ;; message copies to spray in [Binary] Spray & Wait routing
  nMC
  ;; This list should be used to store probabilities or related quantities
  ;; Delivery Predictability metric in PROPHET routing
  ;; Encounter frequency in YAP and DoNothing
  P
  ;; Cumulative Score of Device d;
  ;; Barry Lavelle, Daragh Byrne, Cathal Gurrin, Alan F. Smeaton, Gareth J.F. Jones
  ;; "Bluetooth Familiarity: Methods of Calculation, Applications and Limitations."
  CS_d
  ;; total probability delivery of a sensor in YAP
  Pdelivery
  ;; OPTIMIZATION: keep track of the sum of encounter frequencies (for YAP)
  F
  ;; Total intervals where device d is present.
  I_d
  ;; Dunbar's Number
  D
  Ti-list; The history of encounter of sensor i
  CSTi   ; The list of Strong Ties
  CWTi   ; The list of Weak Ties
  strong-ties  ;Strong Ties Estimate
  weak-ties    ;Weak Ties Estimate

]

to setup
  clear-all
  set-default-shape events "x"
  set-default-shape sensors "triangle"
  set current-showing -1 ;; show-next display sensor 0
  set reported false ;; event has not been reported yet

  ;============================= Deploy Sensors =============================
  ;; deploy sensor first so their agent number start from 0 to avoid problems with sensor-index
  deploy-sensors mobile-sensor-distribution n-of-mobile-sensors
  ask sensors [
    set mobile true
    set S 1 ;; the first location is the sensor's "home"
  ]
  deploy-sensors static-sensor-distribution n-of-static-sensors
  ask sensors with [ mobile = 0 ] [ set mobile false ]  ;; we need to do this since NetLogo defaults variables to 0
  ;; WATCH-OUT!!! If you deploy sensors in a lattice then the actual number of them might have been changed by deploy-sensors procedure
  ;; so we set n-of-mobile-sensors and n-of-static-sensors to the actual value
  set n-of-mobile-sensors count sensors with [mobile]
  set n-of-static-sensors count sensors with [not mobile]
  ;; WATCH-OUT!!! Sinks MUST be deployed before events since the who variable is used as index in P
  deploy sink-location "sink" 1
  deploy event-location "event" 1

  ;============================= Intialize Sensors' State =============================
  ask sensors [
    set size 2 ;; size 2 is cached by netlogo, so it will run faster!
    set found false
    set home-x current-x
    set home-y current-y
    setxy current-x current-y ;; send sensors to "home"

    ifelse (mobile = true) [
      set shape "person"
      set color green
      set sensor-radius mobile-radius
      set sum-of-frequencies 0
      set displacement 0
      ;; [initial-patch-index frequence visit-time] set frequence = 0 because it is updated to 1 at first move call.
      set visited-locations n-values 1 [(list (p-index pxcor pycor) 0 0)]
    ]
    [
      set color violet
      set sensor-radius static-radius
    ]
  ]
  ask sensors with [sink = true] [
    set shape "flag"
    set color yellow
    set size 2
    set mobile false
  ]
  ask events [
    set color red
    set size 2
  ]

  ;============================= Initialize Protocol Variables =============================
  ;; some protocols need to consider the sink as a special node, others do not consider sinks
  ask sensors [
    ;; Time of Last Encounter is used by Gradient (FRESH), PRoPHET, and DoNothing protocols
    set TLE n-values (n-of-mobile-sensors + n-of-static-sensors + 1) [0]
    ;; [Binary] Spray & Wait OR 2000000000000 & Wait
    set nMC 0
    ;; PRoPHET protocol
    if Routing = "PRoPHET" [
      set P n-values (n-of-mobile-sensors + n-of-static-sensors + 1) [0] ;; P(A,x)
      set P replace-item who P 1  ;; P(A,A) = 1
    ]
    if Routing = "ExtractingFriendshipRelations" [
      set P n-values (n-of-mobile-sensors) [0] ;; [F(A,x)]
      set CS_d n-values (n-of-mobile-sensors) [0] ;; [CS_d(A,d)]
      set I_d n-values (n-of-mobile-sensors) [0] ;; [I_d(A,d)]
    ]
    if Routing = "SpreadingToStrongOrWeakTies" [
      set Ti-list [] ;; The List of Encounter Frequencies
      set CSTi [] ;; The List of Strong Ties
      set CWTi [] ;; The List of Weak Ties
      set strong-estimate []; The list of strong ties estimate
      set weak-estimate   []  ;The list of weak ties estimate
      set strong-ties 0       ;The estimated number of strong ties
      set weak-ties 0         ;The estimated number of weak ties
    ]
    if Routing = "DoNothing" [
      set P n-values (n-of-mobile-sensors) [0] ;; [F(A,x)]
      set CS_d n-values (n-of-mobile-sensors) [0] ;; [CS_d(A,d)]
      set I_d n-values (n-of-mobile-sensors) [0] ;; [I_d(A,d)]
    ]
  ]

  ;; these kind of walks are incompatible with Song's model, so disable it!!!
  if (type-of-walk = "Brownian motion (Wiener)") or (type-of-walk = "correlated directions")
    [ set preferential-return false ]

  ;; If training time is set, let's the model warm up before starting the simulation
  reset-ticks ; initialize tick counter
  set training false
  if Training-Time > 0 [
    set training true
    repeat Training-Time [go]
    set training false
    ;; reset the overlapping so that % of area covered start from 0%
    ask patches [ set overlapping 0 ]
    ;; reset sensor state
    ask sensors with [ mobile ] [
      ; send sensors back to home so that they reflect the initial spatial distribution
      set current-x home-x
      set current-y home-y
      setxy home-x home-y
      set displacement 0
      set jump-size 0
      set waiting-time 0
      ; set time of last visit to 0 because ticks starts back from 0
      set visited-locations map [(list (item 0 ?) (item 1 ?) 0)] visited-locations
    ]
  ]
  clear-all-plots ; so the initial state is the state of the model after warm up
  reset-ticks
end 
;########################################################################################################################################################
;########################################################################################################################################################
;########################################################################################################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;; Deployment Functions ;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to deploy [place-in which number]
  if (place-in = "none") [ set number 0 ] ;; do not deploy sinks or events
  if (which = "event" and number > 1 and event-location != "random") [ error "Only random deployment is supported in multi event environment" ]
  if (which = "sink" and number > 1 and sink-location != "random") [ error "Only random deployment is supported in multi sink environment" ]
  if (event-location = "diagonal") xor (sink-location = "diagonal") [ error "Both event and sink must be on diagonal" ]
  repeat number [
    let x 0
    let y 0
    if place-in = "random"[
      set x random-xcor
      set y random-ycor
    ]
    if place-in = "center"[
      if (any? sensors with [sink = true and xcor = 0 and ycor = 0]) or (any? events with [xcor = 0 and ycor = 0]) [ error "You cannot have both event and sink in the center" ]
      set x 0
      set y 0
    ]
    if place-in = "corner"[
      let corners (list max-pxcor max-pycor min-pxcor min-pycor)
      ;; since there are only 4 corners it's highly probable an overlap of target and sink, so we check to avoid it!
      set x one-of corners
      set y one-of corners
      while [(any? sensors with [sink = true and xcor = x and ycor = y]) or (any? events with [xcor = x and ycor = y])]
      [
        set x one-of corners
        set y one-of corners
      ]
    ]
    if place-in = "diagonal" [
      let diag 0
      ifelse diag-dist > 1 [
        set diag (sqrt ((world-width - 1) ^ 2 + (world-height - 1) ^ 2))
      ][ error "Please specify a greater distance" ]
      ifelse (diag - diag-dist) > 0 [
        ;; we put sink and target equally distant from corners
        let cosine (world-width - 1) * (diag - diag-dist) / (2 * diag)
        let sine (world-height - 1) * (diag - diag-dist) / (2 * diag)
        if which = "event" [
          set x min-pxcor + cosine
          set y min-pycor + sine
        ]
        if which = "sink" [
          set x max-pxcor - cosine
          set y max-pycor - sine
        ]
      ][ error "Distance between sink and event is higher than length of diagonal" ]
    ]

    if which = "event" [create-events 1 [setxy x y]]
    if which = "sink" [
      create-sensors 1 [
        set current-x x
        set current-y y
        set sink true
      ]
    ]
  ];;END of repeat number
end 

to deploy-sensors [ sensor-distribution n-of-sensors ]
  if n-of-sensors = 0 [stop]
  if sensor-distribution = "lattice" [
    if (n-of-sensors < 4) [ error "You need at least 4 sensors!" ]
    let x-side floor sqrt(n-of-sensors)
    let x-increment (max-pxcor - min-pxcor) / (x-side - 1)
    let y-side ceiling (n-of-sensors / x-side)
    let y-increment (max-pycor - min-pycor) / (y-side - 1)
    let x min-pxcor
    let y min-pycor
    repeat y-side [
      repeat x-side [
        create-sensors 1 [
          set current-x x
          set current-y y
        ]
        set x x + x-increment
      ]
      set y y + y-increment
      set x min-pxcor
    ]
  ] ;; END "lattice"
  if sensor-distribution = "uniform" [
    create-sensors n-of-sensors [
      set current-x random-xcor
      set current-y random-ycor
    ]
  ]
  if sensor-distribution = "exponential" or sensor-distribution = "normal" or sensor-distribution = "power-law" [
    let next-one FALSE
    create-sensors n-of-sensors [
      let hypotenuse 0
      while [ not next-one ]
      [
        set heading random-float 360
        if sensor-distribution = "exponential" [
          ;; we want 1 - e^(-lambda * world-width / 2) = 0.95; we pass the mean mu = 1 / lambda
          set hypotenuse exponential (- 0.5 * world-width / ln (1 - 0.95))
        ]
        if sensor-distribution = "normal" [
          ;; We want 2*sigma = world-width / 2 => P(X <= world-width / 2) = 0.954499736104
          ;; 3*sigma => 0.997300203937
          set hypotenuse Rayleigh (0.5 * world-width / 2)
        ]
        if sensor-distribution = "power-law" [
          ;; P (X > x) = x^(1-alpha) = 0.05 => x = (1/0.05)^(1/(alpha-1)) = world-width / 2
          ;; => alpha = 1 + 1 / (log (world-width / 2) (1 / 0.05))
          set hypotenuse Levy ( 1 + 1 / (log (world-width / 2) (1 / 0.05)) )
        ]
        if not (patch-at-heading-and-distance heading hypotenuse = nobody) [ set next-one true ]
      ]
      set current-x hypotenuse * dx
      set current-y hypotenuse * dy
      set next-one FALSE
    ]
  ]
end 
;########################################################################################################################################################
;########################################################################################################################################################
;########################################################################################################################################################

to go
  ;; if simulation is running ripristinate overlapping layer
  if (current-showing >= 0) [ display-overlapping-layer ]

  if (Routing = "Gradient") and (sink-location != "none") [
    ;; TLE of sinks is always the most recent (that is, the time elapsed from the last encounter is 0)
    ask sensors with [sink = true] [set TLE replace-item who TLE ticks]
  ]

  ;profiler:start
  ask sensors with [ mobile ] [ move ]
  ;profiler:stop
  ;; WATCH OUT!!! Do NOT use WITH-hack here because we need sensors to check in a random (not synchronized) order!!!
  ;profiler:start
  ask sensors [
    if not training [
      if (event-location != "none") and (not found) and (any? events in-radius sensor-radius) [
        set nMC L ;; if sensor found an event then spray L copies
        set found true
        set color red
      ]
      if (sink-location != "none") and (found) and (any? sensors with [sink = true] in-radius sensor-radius) [
        set nMC 0
        set reported true
        ;; if we reach the sink simulation should stop, so exit from ask sensors
        stop
      ]
    ]
    ;; Forward events towards the sink accordingly to the chosen algorithm
    ; let route-date out of if not training to let sensors update their delivery probabilities
    ;profiler:start
    route-data
    ;profiler:stop
  ] ;;END ask sensors

  if Routing = "ExtractingFriendshipRelations" [
    set CDT (ticks ^ (1 / 3)) ;update Cumulatyve Dynamic Threshold
  ]

  ;; if the event has been reported to a sink then stop the simulation
  if reported [stop]
  tick
end 
;########################################################################################################################################################
;########################################################################################################################################################
;########################################################################################################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;; Move Functions ;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to move
  if (jump-size <= 0) ;and (waiting-time <= 0) ;a new flight must be chosen (note that waiting-time will be always = 0 if wait-time is false)
  [
    let sensor-index who
    ;; we MUST do all this here because only the destination of the travel count for displacement,
    ;; location frequencies and preferential return!!!
    set displacement ( displacement + sqrt ((home-x - current-x) ^ 2 + (home-y - current-y) ^ 2) )
    set sum-of-frequencies sum-of-frequencies + 1
    ;; the following code is doing the same operation as
    ;; ask patch-here [ set frequencies replace-item sensor-index frequencies (item sensor-index frequencies + 1) ]
    ;; but more efficiently
    let patch-index (p-index pxcor pycor)
    ;let  v-loc array:from-list visited-locations
    let location-index position patch-index (n-values (length visited-locations) [ item 0 (item ? visited-locations)] )
    ;let location-index position patch-index (n-values (array:length v-loc) [ item 0 (array:item v-loc ?)] )
    ifelse location-index != false
      [
        set visited-locations replace-item location-index visited-locations (list patch-index (item 1 (item location-index visited-locations) + 1) ticks)
        ;array:set v-loc location-index (list patch-index (item 1 (array:item v-loc location-index) + 1))
        ;set visited-locations array:to-list v-loc
      ]
      [ set visited-locations lput (list patch-index 1 ticks) visited-locations ]
    ifelse preferential-return
    [
      ;; avoid use of (count patches with [ item sensor-index frequencies > 0 ] ) ^ (- gamma) for better performance
      ;set Pnew ro * (S ^ (- gamma))
      ifelse ( random-float 1 < ro * (S ^ (- gamma)) )
      [ ;; if true Explore
        let newLocation patch-here
        while [
          (newLocation = nobody) or
          ;( member? (patch-at-heading-and-distance heading jump-size) ([self] of patches with [ item sensor-index frequencies > 0 ]) )
          member? ( p-index ([pxcor] of newLocation) ([pycor] of newLocation) )  (n-values (length visited-locations) [ item 0 (item ? visited-locations)] )
        ][ ;; we want a new location!!! That is, a patch not visited before
          set heading random-float 360
          set jump-size FlightLength
          set newLocation patch-at-heading-and-distance heading jump-size
        ]
        set S S + 1
      ]
      [ ;; else do a Return jump
        ifelse (random-float 1 < lambda)
        [
          do-frequency-return
        ][
          do-recency-return
        ]
      ]
    ][
      while [ ((patch-at-heading-and-distance heading jump-size) = nobody) or (jump-size = 0)]
      [
        ifelse type-of-walk != "correlated directions"
          [ set heading random-float 360 ]
          [ rt random-normal 0 stdev-angle ]
        set jump-size FlightLength
      ]
    ]

    if wait-time [ set waiting-time round Levy-cutoff beta (1 / cutoff-time) ]
  ] ;; END if jump-size <= 0

  ;; sensors go back "home" at regular time steps (e.g., 24h-48h-72h)
  if back-to-home [
    if (ticks mod back-time = 0) [
      let x (home-x - current-x)
      let y (home-y - current-y)
      if (x != 0 and y != 0) [
        set heading atan x y ;; set sensor heading towards home
        set jump-size sqrt(x ^ 2 + y ^ 2)
      ]
    ]
  ]

  ifelse (waiting-time > 0)
  [ set waiting-time waiting-time - 1 ]
  [
    ;; sensors moves at a fixed speed V <= 1 step/tick
    ifelse(jump-size < 1) [
      set current-x current-x + jump-size * dx ; dx and dy are like cos and sin
      set current-y current-y + jump-size * dy
      set jump-size 0
    ][
      set current-x current-x + dx
      set current-y current-y + dy
      set jump-size jump-size - 1
    ]
    ;; update sensor position
    setxy current-x current-y

    foreach [self] of patches in-radius sensor-radius [
      ask ? [
        set overlapping overlapping + 1
        ;; we visited the patch so color it!
        set pcolor scale-color blue overlapping 1 150
      ]
    ]
  ];;END of if waiting-time > 0
end  ;; END of movend

to do-frequency-return
  let throw random sum-of-frequencies
  let flag false
  let partial-sum 0
  let x 0
  let y 0

  foreach visited-locations
  [
    if not flag [
      set partial-sum partial-sum + item 1 ?
      if (partial-sum > throw) [
        set flag true
        set x (px-index (item 0 ?))
        set y (py-index (item 0 ?))
      ]
    ]
  ]
  facexy x y  ;; set heading towards the new location
  set jump-size sqrt((current-x - x) ^ 2 + (current-y - y) ^ 2)
end 

to do-recency-return
  ; obtain the value of the quantile function from a zipfian distribution and round it to closest integer. This is the recency rank selcted.
  let k (round (Levy nu) - 1) ; because xmin = 1
  while [ k >= length visited-locations] ; TODO: fix the way to select recent location because we loop many times when length of visited-locations is small
  [
    set k (round (Levy nu) - 1)
  ]
  ; order locations according to the visiting time
  let location (item k (sort-by [item 2 ?1 > item 2 ?2] visited-locations))
  let x (px-index (item 0 location))
  let y (py-index (item 0 location))
  facexy x y  ;; set heading towards the new location
  set jump-size sqrt((current-x - x) ^ 2 + (current-y - y) ^ 2)
end 

to-report FlightLength
  if (type-of-walk = "Brownian motion (Wiener)") or (type-of-walk = "correlated directions") [ report 1 ]
  if type-of-walk = "exponential" [ report exponential (1 / lambda) ]
  if type-of-walk = "Rayleigh flight" [ report Rayleigh stdev ]
  if type-of-walk = "Cauchy flight" [ report Cauchy ]
  if type-of-walk = "Levy flight" [ report Levy alpha ]
  if type-of-walk = "Levy with Exp cutoff" [ report Levy-cutoff alpha (1 / cutoff-length) ]
end 
;########################################################################################################################################################
;########################################################################################################################################################
;########################################################################################################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Routing Functions ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;

to route-data
  if Routing = "Epidemic" and found [ ;; if event not found we don't need to update data
    ask other sensors in-radius sensor-radius with [not found] [updateData]
    stop
  ]

  if Routing = "Probabilistic Flooding" and found [
    ask other sensors in-radius sensor-radius with [not found] [ if ((random-float 1) < delta) [updateData] ]
    stop
  ]

  if Routing = "Gradient" and found [
    if any? sensors with [sink = True] in-radius sensor-radius [
      set TLE replace-item who TLE ticks ;; update Time of Last [sink] Encounter
    ]
    let myTLE item who TLE
    ;; we look for sensors which encountered sink more recently
    ask other sensors in-radius sensor-radius with [ (item who TLE) > myTLE ] [updateData]
    stop
  ]

  if Routing = "PRoPHET" [ ;; we must update deliver predictability even if we didn't find the event
    let sensors-in-radius [who] of other sensors in-radius sensor-radius ;; remove myself because P(A,A) == 1
    ;; 1 - When two nodes A and B meet, the first thing they do is to update the delivery predictability for each other
    if (empty? sensors-in-radius) [stop]
    foreach sensors-in-radius [
      let Pold (item ? P)
      ifelse ( Pold < 0.1 ) ;; if delivery predictability is less than 1% than it's like if it never encountered that node
      ;; If node B has not met node A for a long time or has never met node B, such that P_(A,B) < P_first_threshold, then P_(A,B) should be set to P_encounter_first.
      ;; P_encounter_first SHOULD be set to 0.5 unless the node has extra information obtained other than through PRoPHET about the likelihood of future encounters.
        [
          set P replace-item ? P 0.5
          set TLE replace-item ? TLE ticks
        ]
        ;; P_(A,B) = P_(A,B)_old + ( 1 - delta - P_(A,B)_old ) * P_encounter
        [
          set P replace-item ? P (Pold + (1 - 0.01 - Pold) * Pinit)
          set TLE replace-item ? TLE ticks
        ]
    ]
    let other-sensors [who] of other sensors
    ;; 2 - The predictabilities for all other destinations must be 'aged'.
    ;if (not empty? sensors-in-radius) [ ;; age probabilities of other sensors only if a sensor is met, because update of delivery predictability id done only in that case
    foreach other-sensors [
      let Pold (item ? P)
      let tle-x (item ? TLE)
      ;; If a pair of nodes do not encounter each other during an interval, they are less likely to be good forwarders of bundles to each other,
      ;; thus the delivery predictability values must age.
      ;; The delivery predictabilities are aged before being passed to an encountered node so that they reflect the time
      ;; that has passed since the node had its last encounter with any other node.
      ;; P_(A,B) = P_(A,B)_old * gamma^K
      ;; where 0 <= gamma <= 1 is the aging constant, and K is the number of time units that have elapsed since the last time the metric was aged.
      if not ((member? ? sensors-in-radius) or (Pold = 0)) [
        set P replace-item ? P ( Pold * aging ^ (ticks - tle-x) )
        set TLE replace-item ? TLE ticks
      ]
    ]
    ;; 3 - Predictabilities are exchanged between A and B and the 'transitive' property of predictability is used to update the predictability of destinations C
    ;;     for which B has a P(B,C) value on the assumption that A is likely to meet B again:
    foreach sensors-in-radius [
      set nmessages nmessages + 1 ;; we consider 1 message each time we exchange predictabilities
      ;; P_(A,C) = MAX( P_(A,C)_old, P_(A,B) * P_(B,C)_recv * beta )
      ;; where 0 <= beta <= 1 is a scaling constant that controls how large an impact the transitivity should have on the delivery predictability.
      let Pab (item ? P)
      let tle-x (item ? TLE)
      let current-sensor ?
      foreach other-sensors [ ;; skip P(A,A)
        let Pac_old (item ? P)
        let Pbc (item ? [P] of turtle current-sensor)
        let Pac_new (Pab * Pbc * trans)
        if not ((Pac_new < Pac_old) or (member? ? sensors-in-radius) or (Pbc = 0)) [ ;; skip P(B,B) == 1
          set P replace-item ? P Pac_new
          set TLE replace-item ? TLE tle-x
        ]
      ]
    ]

    if found [
      let PA_sink item 0 (item (n-of-mobile-sensors + n-of-static-sensors ) P)
      ;; implements GRTR forwarding strategy described here: http://tools.ietf.org/html/draft-irtf-dtnrg-prophet-10#section-3.6
      ask other sensors in-radius sensor-radius with [ (item (n-of-mobile-sensors + n-of-static-sensors) P) > PA_sink and not found] [updateData]
    ]
    stop
  ] ;; END "PRoPHET"

  if Routing = "Spray & Wait" [
    if (nMC <= 1) [stop] ;; wait phase
    if found [ ;; not needed actually, I keep it just for safety reasons if I mess with go procedure
      foreach ( [self] of other sensors in-radius sensor-radius with [not found] )
      [
        if (nMC > 1) [ ask ? [updateData] ] ;; spray phase
        set nMC nMC - 1
      ]
    ]
    stop
  ] ;; END Spray & Wait

  if Routing = "Binary Spray & Wait" [
    if (nMC <= 1) [stop] ;; wait phase
    if found [
      foreach ( [self] of other sensors in-radius sensor-radius with [not found] )
      [
        if (nMC > 1) [ ;; spray phase
          ask ? [
            updateData
            set nMC int ([nMC] of myself / 2)
          ]
        ]
        set nMC nMC - int (nMC / 2)
      ]
    ]
    stop
  ] ;; END "Binary Spray & Wait"

;; ############ Peer-to-Peer Communication Protocols ##########
  if Routing = "ExtractingFriendshipRelations" [
    let sensors-in-radius [who] of other sensors in-radius sensor-radius

    ;; 1 - When two nodes A and B meet, the first thing they do is to update encounter frequency
    foreach sensors-in-radius [
      ifelse (? < n-of-mobile-sensors) [
        let fk (item ? P) ;; encounter frequency at time t-1
        set P replace-item ? P (fk + 1) ;; F(A,B) = F(A,B)_old + 1
      ]
      [
        error "Out of P range!"
      ]
    ]

    if (ticks mod aging) = 0 [
      ;; 2 - Update the Cumulative Score of Device d
      ;; we recompute Dunbar's number every so often, so reset it
      set D 0
      ;; compute AVG(F) Average of all encounter frequencies within given interval
      let avgF (mean P)
      if avgF = 0 [ set avgF 1 ]
      let pos 0
      foreach P [
        if (? > 0) [
          let CS_d_old (item pos CS_d)
          let I_d_old (item pos I_d)
          set I_d replace-item pos I_d (I_d_old + 1)
          ;; T_i is the length in second (in our case ticks, but it doesn't really matter!) of the interval
          ;; let T_i aging
          ;; let F_d ? ;(item pos P)
          let CS_d_new ( CS_d_old + (? / avgF) * aging / T_d )
          set CS_d replace-item pos CS_d CS_d_new
        ]
        ;; 3 - compute Dunbar's number estimate
        ;; Pinit is the Static Baseline threshold, in the paper is named alpha
        ;; CDT is the Cumulative Dynamic Threshold, in the paper is named beta
        ;; I_d is the total intervals the devic d were present
        ;; if true then it is a "Familiar" device
        if item pos CS_d > (Pinit + CDT * (item pos I_d)) [ set D (D + 1) ]
        set pos (pos + 1)
      ]
      ;; reset frequencies
      set P n-values (n-of-mobile-sensors) [0]
    ]
    stop
  ]  ;; ExtractingFreiendshipRelations

  if Routing = "SpreadingToStrongOrWeakTies" [
    let sensors-in-radius [who] of other sensors in-radius sensor-radius
    ;show sensors-in-radius
    ; when a sensor meet another this encounter is reported in Ti-list
    foreach sensors-in-radius
    [
       ;Multiply a sensor id by 0.0003 in order to use the form (f.id)
       ;the integer part (f) represents the frequancy while the float represents sensor id
       ;for example, 2.0001 means sensor1's frequency is 2.
       let temp1 ?
       let  temp2 (? * 0.0001)
       ; If the history of encounters is empty THEN
       ; we put the currrent encounter into Ti-list and make its frequency = 1
       ifelse(empty? Ti-list) [set Ti-list fput (temp2 + 1) Ti-list]
       [
         foreach Ti-list
         [
           ; we check whether the current encounter is in the history of encounters
           ; if YES we just update the encounter frequency by 1
           ; if NO this means a new item will be inserted into the history of encounters and make its frequency = 1
           ifelse(member? temp1 map[round((? - int ?) * 10000)] Ti-list)
           [ set Ti-list replace-item position ? Ti-list Ti-list (? + 1) ]
           [ set Ti-list fput (temp2 + 1) Ti-list]
         ] ; foreach Ti-list
       ] ;ifelse 1
    ] ; foreach sensor-in-radius
     ; sort Ti-list according to the frequencies
     set Ti-list sort Ti-list
     ;Extracting the number of strong ties
     set strong-ties 0 ;reset the strong-ties buffer
     set strong-ties 0 ;reset the weak-ties buffer
     ; Get the length of 20% (strong ties) in Ti list
     set strong-ties length (sublist Ti-list round(length Ti-list * 0.8) round( length Ti-list))
     ; Get the length of 80% (weak ties) in Ti list
     set weak-ties length (sublist Ti-list 0 round(length Ti-list * 0.8))

     ; NOW this part invovlves CSTi and CWTi
     ; putting the strong and weak ties in their lists when they are empty
     ;Extracting the strong and weak ties according to 80/20 rule
     ; 1- Strong Ties
     foreach map [round((? - int ?) * 10000)] sublist Ti-list (length Ti-list * 0.8) ( length Ti-list)
     [
       ; if this strong tie sensor is not in CSTi
       if(not member? ? CSTi)
       [
         if(member? ? CWTi)
         [ ; remove it from CWTi
           set CWTi remove-item position ? CWTi CWTi
         ]
         ;add it to CSTi
         set CSTi fput ? CSTi
       ]
     ] ; Foreach - For Adding Strong Ties
;
;     ; 2- Weak Ties
     foreach map [round((? - int ?) * 10000)] sublist Ti-list 0 (length Ti-list * 0.8)
     [
       if(not member? ? CWTi)
       [
         set CWTi fput ? CWTi
       ]
     ] ; Foreach - For Adding Weak Ties

     ;; Spreading Phase
     foreach sensors-in-radius
     [
       if(DataSpreadingTo = "Strong Ties")
       [
         if(not found and member? ? CSTi)
         [
           updatedata       ; Transfer the event from sensor ? to those which are in CWTi
         ]
       ]

       if(DataSpreadingTo = "Weak Ties")
       [
         if(not found and member? ? CWTi)
         [
           updatedata        ; Transfer the event from sensor ? to those which are in CWTi
         ]
       ]
     ]
     stop
  ] ;"SpreadingToStrongOrWeakTies"
end  ;; END Route-Data

to updateData
  ;; increase the number of messages exchanged
  set nmessages nmessages + 1
  ;print (word "from " [who] of myself " to " who " at time " ticks)
  set color red
  set found true
end 
;########################################################################################################################################################
;########################################################################################################################################################
;########################################################################################################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Quantile Functions ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;

to-report exponential [ rlambda ]
  ;; step length is chosen from an exponential distribution, with mean (1/lambda) = rlambda
  report random-exponential rlambda
end 

to-report Rayleigh [ sigma ]
  ;; uses a normal distribution with standard deviation (sigma) = sigma and mean (mu) = 0
  report abs random-normal 0 sigma
end 

to-report Cauchy
  ;; quantile function (inverse cdf) of the Cauchy distribution is:
  ;;               x0 + gamma * tan[pi* (p - 1/2)]
  ;; x0 is the location parameter, specifyng the location of the peak of the distribution
  ;; gamma is the scale parameter and is sometimes called the probable error
  ;; when x0 = 0 and gamma = 1 is called standard Cauchy distribution
  report abs tan((random-float 1 - 0.5) * 180 )
end 

to-report Levy [ scaling ]
  ;; length of flight is given by:
  ;;               x = xmin * (1 - r)^(-1/(alpha-1))
  ;; r is a random uniformly distributed real number
  ;; xmin is the lower bound to the power-law behaviour. We assume xmin = 1
  ifelse(scaling > 1)
  [ report (1 - random-float 1) ^ (1 / (1 - scaling)) ]
  [ error "power law distribution must have exponent greater than 1" ]
end 

to-report Levy-cutoff [ scaling lambd ]
  if(scaling < 1) [ error "power law distribution must have exponent greater than 1" ]
  ;; length of fly is choosed according to power law with exponential cutoff
  ;;               x^(-alpha) * e^(-x*lambda)
  ;; where alpha = scaling and lambda = cutoff (in `Understanding individual human mobility patterns - Nature-2008`, lambda = 1/k)
  ;; For the case of the power law with cutoff there is no closed-form expression for quantile,
  ;; but one can generate an exponentially distributed random number using the formula
  ;;               x = xmin − 1/lambda ln(1 − r)
  ;; where r is uniformly distributed and then accept or reject it with probability p or 1 − p respectively,
  ;; where
  ;;               p = (x/xmin)^(-alpha).
  ;; Repeating the process until a number is accepted then gives an x with the appropriate distribution.
  ;;
  ;; This algorithm is a port of randht.py (Python, by Joel Ornstein) showed here: http://tuvalu.santafe.edu/~aaronc/powerlaws/
  let x (list)
  let y (list)
  let xmin 1.
  let n 1 ;; number of samples to return
  let mu (1. / lambd) ; try to avoid recomputing it when q < 0
  ;repeat (10 * n) [ set y lput (xmin - mu * ln(1 - random-float 1)) y ]
  let samples n-values (10 * n) [?] ;; this is a list [ 0 1 2 3 ... 10*n-1 ]
  loop [
    set y (list)
    ;repeat (10 * n) [ set y lput (xmin + random-exponential mu) y]
    repeat (10 * n) [ set y lput (xmin - mu * ln(1 - random-float 1)) y ]
    let ytemp (list)
    foreach samples [
      ;if ( random-float 1 < ((item ? y) / xmin ) ^ (- scaling) ) [ set ytemp lput (item ? y) ytemp ]
      if ( random-float 1 < (item ? y) ^ (- scaling) ) [ set ytemp lput (item ? y) ytemp ] ; do not divide by xmin because xmin = 1
    ]
    ;;set y ytemp
    ;;set x sentence x y ;; concatenates lists
    set x sentence x ytemp ;; no point of setting y when it is not used later because we either return or overwrite it
    let q (length x - n)
    if (q = 0) [ report item 0 x ]
    if (q > 0) [
      let r n-values (length x) [?] ;; this is a list [ 0 1 2 3 ... length(x)-1 ]
      let perm shuffle r
      let xtemp (list)
      foreach r [ if (not member? ? (sublist perm 0 q)) [ set xtemp lput (item ? x) xtemp ] ]
      ;set x xtemp
      ;report item 0 x
      report item 0 xtemp ; no need to reallocate if we are going to return
    ]
    ;if (q < 0) [ we do not need to check, if we didn't get out of the loop this condition is always true, so I moved it at the beginning of loop
    ;  set y (list)
      ;repeat (10 * n) [ set y lput (xmin + random-exponential mu) y]
    ;  repeat (10 * n) [ set y lput (xmin - mu * ln(random-float 1)) y ] ; 1 - random-float 1 = random-float 1
    ;]
  ] ;; END loop
end 
;########################################################################################################################################################
;########################################################################################################################################################
;########################################################################################################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Monitor & Reporters ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;

to-report p-index [px py]
  ;; this is a linear transform from (x,y) matrix style coordinates to [0,(max-pxcor-min-pxcor)] ==>  y * (world-width) + x
  ;; since min-pxcor <= 0 and min-pycor <= 0 we can easily traslate pxcor and pycor to be >= 0
  report (py - min-pycor) * (world-width) + px - min-pxcor
end 

to-report px-index [pindex]
  report ifelse-value (pindex > 0) [ pindex mod world-width + min-pxcor ] [min-pxcor]
end 

to-report py-index [pindex]
  report ifelse-value (pindex > 0) [ int ( pindex / world-width ) + min-pxcor] [min-pycor]
end 

to-report fraction-of-covered-area
  report (count patches with [ overlapping > 0 ] ) / (count patches)
end 

to-report fraction-of-acknowledged-nodes
  report (count sensors with [found]) / (count sensors)
end 

to-report sensor-density
  ;; unit square here is 10x10 patches
  report (n-of-mobile-sensors + n-of-static-sensors) / (world-width * world-height) * 100
end 

to-report mean-S
  ;;n-values count sensors [ count patches with [item ? frequencies > 0] ]
  report mean ([S] of sensors with [mobile])
end 

to-report mean-D
  let Ds [D] of sensors with [mobile and D > 0]
  report ifelse-value (not empty? Ds) [mean Ds][0]
end 

to-report MSD
  report ifelse-value (ticks > 0) [mean [(displacement / sum-of-frequencies) ^ 2] of sensors with [mobile]][0]
end 

to export-patches-own-variables
  if file-exists? (word "frequencies-run-" behaviorspace-run-number ".csv") [ stop ]
  file-open (word "frequencies-run-" behaviorspace-run-number ".csv")
  let nsensors count sensors with [ mobile = true ]
  ;; write the header
  file-type "\"pxcor\",\"pycor\",\"overlapping\""
  let sensor-index 0
  repeat nsensors [
    file-type (word ",\"f" sensor-index "\"")
    set sensor-index sensor-index + 1
  ]
  file-type "\n"
  file-flush
  ;; write values
  ask patches [
    file-type (word pxcor "," pycor "," overlapping)
    set sensor-index 0
    let patch-index (p-index pxcor pycor)
    let freq 0
    repeat nsensors [
      ask sensor sensor-index [
        let visited-patches (n-values (length visited-locations) [ item 0 (item ? visited-locations)] )
        let location-index position patch-index visited-patches
        ifelse (location-index != false) [set freq (item 1 (item location-index visited-locations))][set freq 0]
      ]
      file-type (word "," freq)
      set sensor-index sensor-index + 1
    ]
    file-type"\n"
  ]
  file-close
end 
;########################################################################################################################################################
;########################################################################################################################################################
;########################################################################################################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Display Functions ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;

to display-overlapping-layer
  ask patches [ set pcolor scale-color blue overlapping 1 1000 ]
  set current-showing -1
end 

to display-sensor-layer [ sensor-index ]
  if (sensor-index >= 0) and (sensor-index < n-of-mobile-sensors) [
    ;; update current index
    set current-showing sensor-index
    ;; change color with sensor frequency
    ask patches [ set pcolor 0 ]
    ask sensor sensor-index [
      foreach visited-locations
      [
        ask patch (px-index (item 0 ?)) (py-index (item 0 ?)) [ set pcolor scale-color lime (item 1 ?) 1 20 ]
      ]
    ]
  ]
end