HOTnet-maps

extensions [ gis ]
globals [ cities-dataset
          countries-dataset
          citiesVectorFeature ]
breed [ cities city ]
;; the number of hops from a fixed center of the tree
cities-own [ nhop ]

;;;;;;;;;;;;;;;;;;;;;;;;
;;; Setup Procedures ;;;
;;;;;;;;;;;;;;;;;;;;;;;;

to setup
  clear-all
  ; Note that setting the coordinate system here is optional, as
  ; long as all of your datasets use the same coordinate system.
  gis:load-coordinate-system (word "data/" projection ".prj")
  ; Load all of our datasets
  set cities-dataset gis:load-dataset "data/cities/cities.shp"
  ;set cities-dataset gis:load-dataset "data/cities.shp"
  set countries-dataset gis:load-dataset "data/countries.shp"
  ;; Set the world envelope to the union of all of our dataset's envelopes
  gis:set-world-envelope (gis:envelope-union-of (gis:envelope-of cities-dataset)
                                                (gis:envelope-of countries-dataset))
  ;; Display countries
  gis:set-drawing-color white
  gis:draw countries-dataset 1
   
  ;; list of cities represented by VectorFeature ordered by country
  ;set citiesVectorFeature sort-by [gis:property-value ?1 "CNTRY_NAME" > gis:property-value ?2 "CNTRY_NAME"] gis:feature-list-of cities-dataset
  ;set citiesVectorFeature sort-by [gis:property-value ?1 "COUNTRY" > gis:property-value ?2 "COUNTRY"] gis:feature-list-of cities-dataset
  ;; list of cities represented by VectorFeature in a random order
  set citiesVectorFeature shuffle gis:feature-list-of cities-dataset
  
  set-default-shape cities "circle" 
  ;; create first 2 cities connected by a backbone
  create-cities 1 
  [ 
    let location gis:location-of (first (first (gis:vertex-lists-of item 0 citiesVectorFeature)))
    setxy item 0 location item 1 location
    set size 1
    ;; set color of the city proportionally with population: darker = bigger
    set color scale-color red (gis:property-value item 0 citiesVectorFeature "POP_RANK") 1 7 
    ;set color scale-color red (gis:property-value item 0 citiesVectorFeature "POPULATION") 5000000 1000 
    set nhop 0
  ]
  create-cities 1 
  [ 
    let location gis:location-of (first (first (gis:vertex-lists-of item 1 citiesVectorFeature)))
    setxy item 0 location item 1 location
    set size 1
    ;; set color of the city proportionally with population: darker = bigger
    set color scale-color red (gis:property-value item 1 citiesVectorFeature "POP_RANK") 1 7 
    ;set color scale-color red (gis:property-value item 0 citiesVectorFeature "POPULATION") 5000000 1000 
    create-link-with turtle 0 [ set color green ]
    set nhop 1
  ]
  
  reset-ticks
end 

;;;;;;;;;;;;;;;;;;;;;;;
;;; Main Procedures ;;;
;;;;;;;;;;;;;;;;;;;;;;;

to go
  if (ticks + 2) = length citiesVectorFeature [ stop ]
  ;; new edge is green, old edges are gray
  ask links [ set color gray ]
  ;; index of the list cities
  let current-city item (ticks + 2) citiesVectorFeature
  ; a feature in a point dataset may have multiple points, so we
  ; have a list of lists of points, which is why we need to use
  ; first twice here
  let location gis:location-of (first (first (gis:vertex-lists-of current-city)))
  ; location will be an empty list if the point lies outside the
  ; bounds of the current NetLogo world, as defined by our current
  ; coordinate transformation
  if not empty? location 
  [ 
    ;; The behavior of the model depends crucially on the value of alfa:
    ;; if alfa is less than a particular constant depending on the shape of the region, 
    ;; then Euclidean distances are not important, and the resulting network is easily seen to be a star,
    ;; the ultimate in degree concentration, and, depending on how you look at it, the exact opposite, or absurd extreme, of a power law.
    ;; If alfa grows at least as fast as sqrt(n), where n is the final number of points, then Euclidean distance becomes too important, 
    ;; and the resulting graph is a dynamic version of the Euclidean minimum spanning tree, in which high degrees do occur, 
    ;; but with exponentially vanishing probability.
    ;; If, however, alfa is anywhere in between — is larger than a certain constant, but grows slower than sqrt(n) if at all —
    ;; then, almost certainly, the degrees obey a power law.
    let x item 0 location
    let y item 1 location
    let partner nobody
    ;; Node i attaches itself to the node j that minimizes the weighted sum of the two objectives:
    ;; alfa * dij + hj
    ;; where dij is the /normalized/ Euclidean distance, and hj is some measure of the “centrality” of node j, such as 
    ;; (a) the average number of hops from other nodes; 
    ;; (b) the maximum number of hops from another node; 
    ;; (c) the number of hops from a fixed center of the tree;
    ;; all three measures result in similar power laws.
    ;; We use metric (b). To compute it we choose tourtle 0 as the center of our network. Then the maximum number of hops from a node is:
    ;; number of hops of the node from the center + maximum number of hops a node is from the center.
    ;; we must check the case when nodes with maximum number of hop are children of current node to don't overstimate hj by an excess of max_nhop - 1
    ;; Optionally there is a preference to attach to bigger city
    set partner min-one-of cities
    [ 
      alfa * 
      sqrt 
      ( 
        ;;( (x - min-pxcor + 0.5) / (max-pxcor - min-pxcor) - (xcor - min-pxcor + 0.5) / (max-pxcor - min-pxcor) ) ^ 2 + 
        ( (x - xcor) / (max-pxcor - min-pxcor) ) ^ 2 +
        ;;( (y - min-pycor + 0.5) / (max-pycor - min-pycor) - (ycor - min-pycor + 0.5) / (max-pycor - min-pycor) ) ^ 2 
        ( (y - ycor) / (max-pycor - min-pycor) ) ^ 2
      ) 
      + hj self
      + ifelse-value (city_size_pref? and gis:property-value current-city "POP_RANK" > 0)
      ;+ ifelse-value (city_size_pref? and gis:property-value current-city "POPULATION" > 0)
          ;; the greater is the city the less is the added value
          ;; 7 less than 50K
          ;; 6 50K < people < 100K
          ;; 5 100K < people < 250K
          ;; 4 250K < people < 500K
          ;; 3 500K < people < 1M
          ;; 2 1M < people < 5M
          ;; 1 greater than 5M
          [ gis:property-value current-city "POP_RANK" ] [0]
    ]
  
    if partner != nobody
    [ 
      create-cities 1 
      [ 
        setxy x y
        set size 1
        ;; set color of the city proportionally with population: darker = bigger
        set color scale-color red (gis:property-value current-city "POP_RANK") 1 7 
        ;set color scale-color red (gis:property-value current-city "POPULATION") 5000000 1000 
        create-link-with partner [ set color green ]
        set nhop 1 + [ nhop ] of partner
      ]
    ]
  ] ;; END if not empty? location
 
  tick
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Compute hj heuristic ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to-report hj [node]
    let max_nhop max [nhop] of cities
    ;; if all cities at max_nhop are children of current-city then decrease hj
    while [ all? cities with [nhop = max_nhop] [ is-child node max_nhop] ]
    [ 
      set max_nhop max_nhop - 1
    ]
    report max_nhop + [nhop] of node
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Find if city at max_ops is child of current-city ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to-report is-child [root max-nhop]
  let child false
  if root != city 0 ;; all cities are children of the root of the tree, so don't check it!
  [
    ask root
    [
      ask link-neighbors with [nhop > [nhop] of root]
      [
        ifelse nhop = max-nhop
          [ set child true ]
          [ set child is-child self max-nhop] ;; use of recursion to traverse the tree
      ]
    ]
  ]
  report child
end 

;;;;;;;;;;;;;;;;;;;;
;;; Compute s(g) ;;;
;;;;;;;;;;;;;;;;;;;;

to-report log-likelihood
  let s 0
  ;; for each link compute di*dj and sum it to s
  ask links 
  [ 
    set s s + [ count link-neighbors ] of end1 * [ count link-neighbors ] of end2 
  ]
  report s
end 
;;;;;;;;;;;;;;;;;;;;;;;;
;;; Compute S-metric ;;;
;;;;;;;;;;;;;;;;;;;;;;;; 

to-report relative-log-likelihood
  let smax 0
  let counter 0
  let di 0
  let child 0

  ;; D = { d1, d2, d3, ..., dn }, d1 >= d2 >= d3 >= ... >= dn
  let degree-sequence sort-by > [ count  link-neighbors ] of turtles
  set di item 0 degree-sequence
  set degree-sequence remove-item 0 degree-sequence
  foreach degree-sequence
  [
    set smax smax + di * ?
    set counter counter + 1
    if di = counter ;; we have iterated through all di's childs; if di = 0 select the highest degree.
    [
      set counter 1 ;; count the parent if it's not the root
      set di item child degree-sequence ;; select child; child = 0 is the root.
      set child child + 1
    ]
  ]

  report log-likelihood / smax
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Save Nodes Degrees to file ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to save-node-degree-to-file
  if file-exists? "NodeDegrees.txt" [ file-delete "NodeDegrees.txt" ]
  
  file-open "NodesDegrees.txt"
  
  ;; save in descending orders
  ;; D = { d1, d2, d3, ..., dn }, d1 >= d2 >= d3 >= ... >= dn
  foreach sort-by > [ count link-neighbors ] of turtles
  [
    file-print ?
  ]
  file-close
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Export Graph Connectivity to txt ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to export-graph
  if file-exists? "HOTGraph.txt" [ file-delete "HOTGraph.txt" ]
  
  file-open "HOTGraph.txt"
  
  ;; write each linked couple of tourtles and their degree
  ask links 
  [ 
    file-type [who] of end1 ;; writes without blank spaces
    file-write [who] of end2 ;; write a space value space
    file-print "" ;; write carriage return
  ]
  file-close
end 

;;;;;;;;;;;;;;
;;; Layout ;;;
;;;;;;;;;;;;;;
;; resize-nodes, change back and forth from size based on degree to a size of 1

to resize-nodes
  ifelse all? turtles [size = 1]
  [
    ;; a node is a circle with diameter determined by
    ;; the SIZE variable; using SQRT makes the circle's
    ;; area proportional to its degree
    ask turtles [ set size sqrt count link-neighbors ]
  ]
  [
    ask turtles [ set size 1 ]
  ]
end 


; Copyright 2013 Tomasini Marcello.
; See Info tab for full copyright and license.