Segregation model: Model of Connectivity

; ===========
; ## NOTES ##
; ===========

; The model is built partly with code from Crooks et al (2018), Wilensky (1997).
; The A* algorithm applied in NetLogo is taken from http://www.cs.us.es/~fsancho/?e=131

extensions [
  rnd
]

globals [
  percent-unhappy
  prop_neighborhood ;; proportion of agents of type X in a neighborhood // a report variable
  p-valids ;Contain patches which are not barriers.
  Start ;Start position (unhappy turtle)
  Final-Cost ;Store the distance returned by the A* algorithm.
  pot_NBH
  ethnic_share
  similar
  unhappy
  dis ; For dissimilarity index calculation
  indexdissimilarity
]

patches-own [
  NBH_id ;Neighborhoods id.
  father ;Is used by A* to store paths.
  Cost-path ;Is used to store the cost of the path (distance).
  visited? ;A boolean indicating if the patch have been visited before.
  active? ;A boolean indicating if the patch is active in the search for shortest path.
  neigborhood2 ; To be used to calculated segregation index.
]

turtles-own [
  happy? ;boolean variable indicating if the agent is happy with their current living situation (given Schelling, 1971 behavior rules).
  similar-nearby ; similar in neighborhood where agent is located
  myneighbors ; a agent's polygon-neighbor, same as adjacent_polygon of the polygon where the agent is located.
  other-nearby  ; others in neighborhood where agent is located
  total-nearby ; total population of neigborhood where agent is located
  NBH ;Which neighborhoods do the turtle dwell in.
]

; ==========================
; ## SETTING UP THE MODEL ##
; ==========================

to setup
  ca
  set similar []
  set unhappy []

  ; Setting world.
  resize-world 0 (world_size - 1) 0 (world_size  - 1)
  set-patch-size patch_size ; 7 o 36

  ; Identifying neighborhood.
  ; Only works with tens (10, 20, 30, 40, 50 etc...)
  let x sizenbh
  let y sizenbh
  let col 36
  while [y <= max-pycor + 2][
    while [x <= max-pxcor + 2][
      ask patches with [pxcor < x and pxcor >= x - sizenbh and pycor < y and pycor >= y - sizenbh][
        set pcolor col
      ]
      set x x + sizenbh
      set col col + 3
    ]
    set x sizenbh
    set y y + sizenbh
  ]

  ask patches[
    set neigborhood2 [pcolor] of self]

  ; Generating physical features.
  if Environment = "Barrier" or Environment = "Bridges" [
    if Ring-road [
      crt (max-pycor * 3)
      layout-circle turtles (max-pycor * ((random (max-pxcor / 2) + (max-pxcor / 1.5 )) / 100))
      ask turtles [set pcolor brown]
      ask turtles [die]
    ]

    ; Urban form: Monocentric
    if Urban-form  = "Monocentric" [
      let center-region patches with [ pycor = (round(max-pycor / 2)) and pxcor = (round(max-pxcor / 2))]
      ask center-region [sprout 4]
      let i 0
      ask turtles [
        set pcolor brown
      ]

      while [i < (max-pycor )][
        ask one-of turtles [
          fd 1
          set pcolor brown
          fd 1
          set pcolor brown
          fd 1
          set pcolor brown
          fd 1
          set pcolor brown
          ifelse (random 25 <= 1)[
            rt one-of [-90 0 90]
            fd 1
            set pcolor brown
            fd 1
            set pcolor brown
            fd 1
            set pcolor brown
          ][
            set pcolor brown
            fd 1
            set pcolor brown
            fd 1
            set pcolor brown
            fd 1
            set pcolor brown
          ]
          set i i + 1]
      ]
      ask turtles [die]
    ]

    ; Urban form: Polycentric
    if Urban-form = "Polycentric" [
      ask n-of 2 patches with [ (pycor = (round(max-pycor * 0.5 )) and pxcor = (round(max-pxcor * 0.3))) or (pycor = (round(max-pycor * 0.5)) and pxcor = (round(max-pxcor * 0.7))) or (pycor = (round(max-pycor * 0.7)) and pxcor = (round(max-pxcor * 0.50))) or(pycor = (round(max-pycor * 0.3)) and pxcor = (round(max-pxcor * 0.50)))][sprout 2]
      let i 0
      ask turtles [
        set pcolor brown
      ]
      while [i < (max-pycor )][
        ask one-of turtles [
          fd 1
          set pcolor brown
          fd 1
          set pcolor brown
          fd 1
          set pcolor brown
          fd 1
          set pcolor brown
          ifelse (random 25 <= 1)[
            rt one-of [-90 0 90]
            fd 1
            set pcolor brown
            fd 1
            set pcolor brown
            fd 1
            set pcolor brown
          ][
            set pcolor brown
            fd 1
            set pcolor brown
            fd 1
            set pcolor brown
            fd 1
            set pcolor brown
          ]
          set i i + 1 ]
      ]
      ask turtles[die]
    ]
  ]

  ; Bridges.
  if Environment = "Bridges" [
    ask n-of (count patches with [pcolor = brown] * 0.05) patches with [pcolor = brown][
      ask patches in-radius 1 [
        set pcolor 139]
      ]
    ]



  ; Adding agents.
  ask n-of ((count patches with [pcolor != brown]) * ( 1 - (%VacantHH / 100))) patches with [pcolor != brown][sprout 1]
  ask turtles [
    set color red
    set shape "circle"
    set size 0.9
  ]
  ask n-of ((minority-group-size / 100) * count turtles) turtles [set color blue]

  ask patches with [pcolor = 139][set pcolor ([neigborhood2] of self)]

  ; Adding neighborhood.
  ask patches with [pcolor != brown][ set NBH_ID nobody ]
  find-clusters

  ; Giving each neighborhood an ID.
  ask patches with [pcolor != brown] [
    set pcolor black]

  ask turtles[
    set NBH [NBH_ID] of patch-here]

  ; Setting start attributes for the A* algorithm
  ask patches
  [
    set father nobody
    set Cost-path 0
    set visited? false
    set active? false
  ]

  set p-valids patches with [pcolor != brown]
  reset-ticks
end 

; ====================
; ## AGENT BEHAVIOR ##
; ====================

to go
  update-happiness
  if any? turtles with [happy? = false][
    ask one-of turtles with [happy? = false][new-neighborhood]]
  update-globals
  tick
end 

; All agents evaluate their living situation.

to update-happiness
  ask turtles[
    ;Calculating composition of egocentric neighborhood.
    set myneighbors turtles with [NBH = [NBH] of myself]
    set similar-nearby count myneighbors with [color = [color] of myself] ; Number of similar agents in neighborhood
    set other-nearby count myneighbors with [color != [color] of myself]   ; Number of other agents in neighborhood.
    set total-nearby similar-nearby + other-nearby  ; Total agents in neighborhood.
    ;Determine if an agent is happy
    set happy? similar-nearby >= (%-similar-wanted * total-nearby / 100)
  ]
end 

to new-neighborhood
  ; Identifying avaiable neigborhoods.
  ; Listing all potential neighborhoods.
  if Printer [print "-------------------"]
  set Start self
  if Printer [print(word" ""Unhappy agent is: " Start)]

  ; Identifying all vacant households outside the the agent's current neighborhood.
  let pot_NBH1 patches with [not any? turtles-here and pcolor != brown and NBH_id != [NBH] of myself]
  set pot_NBH sort pot_NBH1
  if Printer [print(word" ""My current neighborhood is: " [NBH] of Start)]

  ; 'Ethnic' share of all possible vacant neigborhoods.
  ; Identifying agent's 'ethnicity
  let mycolor [color] of self
  ; Identifying polygon neighborhood
  let pot_NBH2 map [ i -> [NBH_id] of i ] pot_NBH
  ; Computing ethnic share of each neighborhood.
  set ethnic_share map [i -> prop_same_group i mycolor] pot_NBH2

  if Printer [print(word" ""Agent is happy?"  happy?)]
  if Printer [print(word" ""Agent's ethnicity :" mycolor)]

  ; Distance to all neighborhoods with vacant household.
  let distance_NBH map [i -> A* Start i p-valids] pot_NBH

  ; Normalzing distance, makes it weight more even to ethnic share.
  let distance_NBH_normal map [i -> normalize2 i ] distance_NBH

  ; Calculuating utility.
  if Printer [print " Calculating utility"]

  ; Ethnic share
  let utility_share map [i -> i * beta_ethnic] ethnic_share

  ; Household distance
  let utility_distance map [i -> i * beta_distance] distance_NBH_normal

  if Printer [print(word" ""Distance utility :"  utility_distance)]
  if Printer [print(word" ""Share utility :"  utility_share)]
  ; Summing utility
  let added_utility2 []
  (foreach utility_share utility_distance  [ [a b ] -> let suma (a + b ) set added_utility2 lput suma added_utility2 ])

  ; Exp() utility
  let added_utility map [i -> exp i ] added_utility2

  ; Calculating probabilities
  let tot_utility sum added_utility
  let NBH_probabilities map [i -> i / tot_utility] added_utility

  ; Joining final utilities and NBHs' ids
  let final_pool (map list pot_NBH NBH_probabilities)
  if Printer [print(word" ""Each neighborhoods utility :" final_pool)]

  ; Selecting the neighborhood with highest utility and moves.
  let origin_place patch-here
  if Printer [print(word" ""Moving from NBH :"  origin_place)]
  if Printer [ pen-down ]
  let decision select_and_move Start final_pool origin_place
end 

; =============
; ## REPORTs ##
; =============
;To identify neighborhood, made on code from 'Patch Clusters example' model.

to find-clusters
  loop [
    ; Pick a random patch that isn't in a cluster yet
    let seed one-of patches with [NBH_id = nobody]
    ; If we can't find one, then we're done!
    if seed = nobody
    [
      stop
    ]
    ; Otherwise, grow-cluster to find the rest of the cluster
    ask seed
    [ set NBH_id self
      grow-cluster ]
  ]
  display
end 

to grow-cluster
  ; Patch procedure
  ask neighbors4 with [(NBH_id = nobody) and
    (pcolor = [pcolor] of myself)]
  [
    set NBH_id [NBH_id] of myself
    grow-cluster ]
end 

; Calculating fraction of similar neighbors in each neighborhood with a vacant household.

to-report prop_same_group [a mycolor]
    let similar_nearby count turtles with [color = mycolor and NBH = a]
    let total_nearby count turtles with [NBH = a]
    ifelse (total_nearby > 0) [set prop_neighborhood (similar_nearby / total_nearby)] [ set prop_neighborhood 0]
  report prop_neighborhood
end 

to-report normalize2 [ current_distance ]
  let b precision ( current_distance / ( world_size * 3 )) 4
  report b
end 

to-report select_and_move [a b c]
 ask a [
  if ( (length b) > 0 )
    [
      let chosen_NBH first rnd:weighted-one-of-list b [ [p] -> last p ]
      if Printer [print(word"""Moved to NBH = "chosen_NBH)]
      move-to chosen_NBH
      set NBH [NBH_ID] of patch-here
    ]
  ]

  report []
end 

; A* algorithm, used to calculate shortest distance to all vacant households.
; The A* algorithm applied in NetLogo is heavily based on the code from:  http://www.cs.us.es/~fsancho/?e=131

to-report A* [#Start #Goal #valid-map]
  ; Reset everything.
  ask #valid-map with [visited?]
  [
    set father nobody
    set Cost-path 0
    set visited? false
    set active? false
  ]

  ; Active the starting point to begin the searching loop
  ask #Start
  [
    set father self
    set visited? true
    set active? true
  ]

  ; Exists? indicates if the two points can be reached with the A* algorithm.
  let exists? true
  while [not [visited?] of #Goal and exists?]
  [
    ; Only look at 'walk-able' patches which has not been visited before.
    let options #valid-map with [active?]
    ifelse any? options
    [
      ; Chooses active patches with minimal expected cost (heuristic)
      ask min-one-of options [Total-expected-cost #Goal]
      [
        ; Store its real cost (to reach it) to compute the real cost of its children
        let Cost-path-father Cost-path
        ; and deactivate it, because its children will be computed right now
        set active? false
        ; Compute its valid neighbors and look for an extension of the path
        let valid-neighbors neighbors4 with [member? self #valid-map]
        ask valid-neighbors
        [
          let t ifelse-value visited? [ Total-expected-cost #Goal] [2 ^ 20]
          if t > (Cost-path-father + distance myself + Heuristic #Goal)
          [
            ; The current patch becomes the father of its neighbor in the new path.
            set father myself
            set visited? true
            set active? true
            ; Store the distance in the Final-Cost variable
            set Cost-path Cost-path-father + distance father
            set Final-Cost precision Cost-path 3
          ]
        ]
      ]
    ]
    [
      set exists? false ;Can a path be found?
    ]
  ]
  ; If a path can be found (not isolated)
  ifelse exists?
  [
    set Final-Cost (precision [Cost-path] of #Goal 3)
    report Final-Cost
  ]
  [
    ; Otherwise, there is no path: use the Euclidean distance.
    set Final-Cost (precision (distance #Goal) 3) * 3
    report Final-Cost
  ]
end 

to-report Total-expected-cost [#Goal]
  report Cost-path + Heuristic #Goal
end 

; The heuristic is based on Euclidean distance.

to-report Heuristic [#Goal]
  report distance #Goal
end 

to update-globals
  ; Dissimilarity index.
  let tot_red (count turtles with [color = red])
  let tot_blue (count turtles with [color = blue])

  let neighb1 [neigborhood2] of patches
  let neighb remove-duplicates neighb1

  set dis []

  foreach neighb [i -> set dis lput abs((count turtles with [neigborhood2 = i and color = red] / tot_red) - (count turtles with [neigborhood2 = i and color = blue] / tot_blue)) dis]
  set indexdissimilarity sum(dis) / 2
end 

; ================
; ## REFERENCES ##
; ================
; Crooks, A., Malleson, N., Manley, E., & Heppenstall, A. (2018). Agent-based modelling and geographical information systems: a practical primer. SAGE Publications Limited.
; Patch Cluster Example (in NetLogo)
; Sancho Caparrini, F. (2018). A General A* Solver in NetLogo. Retrieved 2020-03-15 from: http://www.cs.us.es/~fsancho/?e=131
; Wilensky, U. (1997). NetLogo Segregation model. http://ccl.northwestern.edu/netlogo/models/Segregation. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL.