Artificial Anasazi revisited

extensions [ gis table ]

breed [ households household ]

globals [ start-year fertility-age fertility-end-age death-age food-requirement
          food-estimate-margin corn-gift-ratio corn-stock-years good-farm-bias
          year pdsi-table population-sequence historical-sequence harvest-list fit-badness desired-fit-length ]

patches-own [ value zone patch-quality base-yield historical-base-yield being-farmed? ]

households-own [ harvest age corn-storage unsatisfied-hunger ]


;;;;;;;;;;;;;;;;;;;;;
;; main setup code ;;
;;;;;;;;;;;;;;;;;;;;;

to setup
  clear-all
  set-default-shape households "person"
  
  gis-map-load "adata/Map.asc"
  set pdsi-table load-zone-file "adata/ZoneAdjPDSI.txt"
  ask patches [ initialise-patch ]

  set start-year           800        ;; initial year
  set food-requirement     (5 * 160)  ;; a person needs 160 kg of corn each year, while a typical household consist of 5 persons
  set fertility-age        17         ;; the minimum age of agents (households) that can reproduce (spin off daughter households)
  set fertility-end-age    36         ;; the maximum age of agents (households) that can reproduce
  set death-age            38         ;; the age after which agents (households) die (dissolve when the matriarch dies)  
  set corn-gift-ratio      0.33       ;; each new household gets this fraction of the corn storage of the parent household
  set corn-stock-years     2          ;; corn can be stored for this many years
  set good-farm-bias       1000       ;; larger numbers mean food yield is taken more strongly into account when choosing a farm
  set food-estimate-margin 1.1        ;; agents want this multiple of food requirements, or they will move to another farm
  set year start-year
  set population-sequence [ ]
  set harvest-list [ ]
  
  set historical-sequence
    [ 14 14 14 14 14 14 14 14 14 14 13 13 13 12 12 12 12 12 12 12 11 11 11 11 11 10 10 10 10 10 9 9 9 9 9 9 9 9 8 8 7 7 7 7 7 7 7 7 7 7
      28 29 29 29 29 29 29 29 29 29 29 30 30 31 31 31 31 32 32 32 32 33 33 33 33 33 37 37 37 37 37 39 39 39 40 40 40 40 41 41 41 42 42
      42 42 42 42 42 42 42 56 58 58 58 58 58 58 58 58 58 58 60 60 60 60 60 60 61 61 61 60 61 61 61 61 59 61 61 61 61 61 62 62 62 62 62
      62 62 62 62 61 63 63 63 63 63 63 63 63 63 66 67 67 67 67 67 67 67 67 67 66 67 67 67 67 67 67 66 66 66 66 69 69 69 69 65 68 68 68
      68 66 67 67 67 66 66 66 66 66 66 66 68 68 68 68 68 68 68 68 68 87 87 87 87 87 87 87 86 86 87 85 86 86 87 87 87 87 88 88 88 87 87
      87 87 87 88 90 90 90 90 89 91 92 92 95 95 95 95 97 97 94 94 94 94 95 95 95 95 95 95 134 139 139 139 139 139 139 142 142 142 142
      143 143 146 146 146 146 151 151 153 151 151 151 151 151 156 164 164 164 164 163 163 163 163 165 165 165 165 164 164 163 164 164
      164 161 162 162 162 162 162 161 164 164 165 165 166 166 166 167 166 166 166 166 159 159 160 160 158 159 160 160 161 162 162 162
      148 150 151 151 151 151 151 149 149 147 148 148 148 148 150 150 150 151 151 152 152 152 153 153 153 116 117 118 118 119 119 120
      121 122 126 124 124 124 126 127 127 129 131 131 133 132 134 134 136 137 133 138 138 139 139 138 140 140 142 142 142 143 143 144
      145 145 147 146 147 147 149 149 150 151 151 145 146 146 147 148 149 149 151 153 154 155 157 159 160 161 163 164 166 167 167 169
      169 171 171 173 170 173 176 176 178 178 180 182 184 184 185 188 189 190 192 191 193 192 194 194 194 195 196 199 200 192 193 194
      196 196 197 200 202 202 204 201 209 208 211 212 212 213 214 215 216 210 208 206 201 196 188 181 176 172 167 159 156 148 146 141
      120 118 114 106 103 95 90 88 83 74 70 68 60 58 56 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
      0 0 0 0 0 0 0 0 0 0 0 0 ]
  
  make-region-label 45 98  "Non-arable uplands" black  ;; region is yellow
  make-region-label 28 47  "General valley" brown
  make-region-label 23 73  "North valley" red
  make-region-label 50 68  "Mid valley" gray
  make-region-label 25 84  "Dunes" green  ;; North and Mid valley dunes
  make-region-label 20 102 "Kinbiko Canyon" pink
  make-region-label 39 11  "Arable uplands" orange
 
  set desired-fit-length (length filter [ ? > 0 ] historical-sequence)
  
  let hist-init-households 14  ;; starting number of households (based on historical data)
  calculate-patch-yields       ;; used to find initial farms
  create-households hist-init-households [
    initialise-household (random 29)  ;; initial ages up to 28
    set corn-storage (n-values (corn-stock-years + 1) [ 600 + random-float 400 ])  ;; corn-storage divided into age categories (0, 1, 2)
    find-best-farm potential-farms
  ]

  reset-ticks
  tick-advance start-year  ;; start counting from the year 800
end 

to initialise-household [ start-age ]
  set age start-age
  set harvest 0
  set size 4
  set color black
end  

to find-best-farm [ farms ]  ;; find the best farm from the set, and leave the valley (i.e. die) if there is none
  ifelse (count farms = 0)
    [ die ]
    [ let existing-farm patch-here
      ; these distance calculations assume that the edges of the world do not wrap around
      let best-farm (min-one-of farms [ distance existing-farm + good-farm-bias / estimated-farm-yield ])
      ask best-farm [ set being-farmed? true ]
      move-to best-farm ]
end 

to-report potential-farms  ;; calculate a set of available farm locations
  report (patches with [ estimated-farm-yield >= food-requirement and not being-farmed? ])
end 

to-report estimated-farm-yield
  ;; just go for something simple: this year's base yield (but we might substitute a better definition later)
  report base-yield
end 


;;;;;;;;;;;;;;;;;;;;
;; main step loop ;;
;;;;;;;;;;;;;;;;;;;;

to go
  calculate-patch-yields
  annual-household-activities
    
  ask households [
    ;; agents with insufficient food for next year try to find a new farm
    if (food-estimate < food-estimate-margin * food-requirement) [
      ask patch-here [ set being-farmed? false ]
      find-best-farm potential-farms
    ]
  ]
  
  let farms potential-farms  ;; avoid recomputing this unnecessarily
  ask households [
    if (age >= fertility-age and age <= fertility-end-age and random-float 1 < fertility and count farms > 0) [ 
      do-reproduce farms
      set farms potential-farms
    ]
  ]
  
  plot-interesting-data
  if (year = 1350 or count households = 0) [ stop ]
  set year year + 1
  tick
end 

to plot-interesting-data
  set-current-plot "Population (households)"
  set-current-plot-pen "Potential"
  plotxy year (count patches with [ mean historical-base-yield >= food-requirement ])
  set-current-plot-pen "Historical"
  plotxy year (item (year - 800) historical-sequence)  ;; index list starting at 0
  set-current-plot-pen "Model"
  plotxy year (count households)
  
  let age-list ([ age ] of households)
  
  set-current-plot "Household ages"
  let xmax (1 + max (fput 0 age-list))
  set-histogram-num-bars xmax
  if (plot-x-max < xmax) [ set-plot-x-range 0 xmax ]  ;; histograms require explicit x-max control
  histogram age-list
  
  set-current-plot "Household harvests"
  let block-size 100
  set xmax (1 + ceiling (max (fput 0 harvest-list) / block-size))
  if (plot-x-max < xmax * block-size) [ set-plot-x-range 0 (xmax * block-size) ]
  set-histogram-num-bars (plot-x-max / block-size)
  histogram harvest-list

  if (length population-sequence < desired-fit-length) [ set population-sequence (lput (count households) population-sequence) ]
  let cut-historical-sequence (n-values (length population-sequence) [ item ? historical-sequence ])
  set fit-badness (sqrt (mean (map [ (?1 - ?2) * (?1 - ?2) ] population-sequence cut-historical-sequence)))
end 


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; code for household actions at each step ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to annual-household-activities  ;; calculate the actual harvest of households, update the corn-storage list, and increment age
  set harvest-list [ ]
  ask households [
    set age age + 1
    if (age > death-age) [
      ask patch-here [ set being-farmed? false ]
      die
    ]

    set harvest (apply-variability harvest-variability) * ([ base-yield ] of patch-here)
    set harvest-list (fput harvest harvest-list)
    
    set corn-storage (lput harvest (but-first corn-storage))  ;; oldest corn is at the front of the list
    set unsatisfied-hunger food-requirement
    set corn-storage (map [ eat-corn ? ] corn-storage)  ;; eat the corn starting from the front of the list (the eat-corn reporter transforms old values to new values)
    if (unsatisfied-hunger > 0) [
      ask patch-here [ set being-farmed? false ]
      die
    ]
  ]
end 

to-report eat-corn [ existing-amount ]  ;; this reporter attempts to subtract unsatisfied-hunger from the input (and as a side-effect also updates unsatisfied-hunger)
  ifelse (existing-amount >= unsatisfied-hunger)
    [ let remaining-amount (existing-amount - unsatisfied-hunger)
      set unsatisfied-hunger 0
      report remaining-amount ]
    [ set unsatisfied-hunger (unsatisfied-hunger - existing-amount)
      report 0 ]
end 

to-report food-estimate  ;; estimate the amount of food available for next year, based on current stocks of corn, and an estimate of the future harvest
  let future-estimate harvest  ;; predict next year's harvest to be the same as this year's harvest
  report (future-estimate + sum (but-first corn-storage))  ;; ignore the corn on the front of the list, which is now too old to eat
end 

to do-reproduce [ farms ]
  let gift-corn (map [ ? * corn-gift-ratio ] corn-storage)
  set corn-storage (map [ ?1 - ?2 ] corn-storage gift-corn)
  
  hatch 1 [ 
    initialise-household 0  ;; initial age of zero
    set corn-storage gift-corn
    find-best-farm farms
  ]
end 


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; code for calculating PDSI etc. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to-report get-pdsi [ z y ]
  ifelse (table:has-key? pdsi-table z)
    [ report item (y - start-year) (table:get pdsi-table z) ]
    [ report 0 ]
end 

to-report thresholded-cubic [ x a b c d lo hi ]
  let res (x * x * x * a + x * x * b + x * c + d)
  if (res < lo) [ set res lo ]
  if (res > hi) [ set res hi]
  report res
end 

to calculate-patch-yields  ;; calculate the crop yield for each patch based on the PDSI data, using a cubic interpolation of published numbers
  ask patches [      
    let pdsi get-pdsi zone year
    let theoretical-yield 0

    if (zone = "North" or zone = "Kinbiko" or (zone = "Mid" and pxcor <= 34)) [
      set theoretical-yield (thresholded-cubic pdsi 4.4167 6.9456 49.583 823.48 617 1153)
    ]
    if (zone = "General" or (zone = "Mid" and pxcor > 34)) [
      set theoretical-yield (thresholded-cubic pdsi 3.65 5.7925 41.65 686.28 514 961)
    ]
    if (zone = "Uplands") [
      set theoretical-yield (thresholded-cubic pdsi 2.9333 4.6599 33.267 548.77 411 769)
    ]
    if (zone = "North Dunes" or zone = "Mid Dunes") [
      set theoretical-yield (thresholded-cubic pdsi 4.5833 7.1871 51.917 858.03 642 1201)
    ]

    set base-yield (theoretical-yield * patch-quality * harvest-adjustment)
    set historical-base-yield (fput base-yield historical-base-yield)  ;; put new base-yield on the front of the list
    if (length historical-base-yield > 1 + corn-stock-years) [ set historical-base-yield (but-last historical-base-yield) ]  ;; truncate to 3 entries, if necessary
  ]
end 


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; code for setting up initial patch data ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to initialise-patch
  set zone (patch-zone-name value)
  if (zone = "?") [ user-message (word "Error in patch data load for " self) ]
  set pcolor (patch-zone-color zone)
  set being-farmed? false
  set historical-base-yield [ ]  ;; this will be expanded to a list of the last 3 yields
  
  ifelse (random-float 1 < usable-farm-fraction)
    [ set patch-quality (apply-variability crop-yield-variability) ]
    [ set patch-quality 0 ]
end 

to-report apply-variability [ var ]
  let res random-normal 1 var  ;; random variation around a mean of 1
  if (res < 0) [ set res 0 ]   ;; because random-normal can return negative values
  report res
end 

to-report load-zone-file [ fname ]  ;; load a file with a list of items for each zone (the zone will be the first element of each list)
  let tbl table:make
  ifelse (file-exists? fname)
    [ file-open fname
      while [ not file-at-end? ] [
        let lst file-read
        table:put tbl (item 0 lst) (but-first lst)  ;; split first item (zone name) from rest of list
      ]
      file-close
    ]
    [ file-error fname ]
  report tbl
end 

to file-error [ fname ]
  user-message (word "Cannot find the file \"" fname "\"")
end 

to gis-map-load [ fname ]
  ifelse (file-exists? fname)
    [ let dataset gis:load-dataset fname
      gis:set-world-envelope (gis:envelope-of dataset)
      gis:apply-raster dataset value
    ]
    [ file-error fname ]
end 

to-report patch-zone-color [ z ]
  ifelse (z = "General") [ report brown ] [
  ifelse (z = "North") [ report red ] [
  ifelse (z = "Mid") [ report gray ] [
  ifelse (z = "Natural") [ report yellow ] [
  ifelse (z = "Uplands") [ report orange ] [
  ifelse (z = "Kinbiko") [ report pink ] [ 
  ifelse (z = "North Dunes") [ report green ] [ 
  ifelse (z = "Mid Dunes") [ report green ] [ 
  ifelse (z = "Empty") [ report white ] [ 
  report magenta ] ] ] ] ] ] ] ] ]  ;; there should be no magenta patches
end   

to-report patch-zone-name [ i ]
  ifelse (i = 0) [ report "General" ] [      ;; General valley floor
  ifelse (i = 10) [ report "North" ] [       ;; North valley floor
  ifelse (i = 15) [ report "North Dunes" ] [ ;; North valley dunes
  ifelse (i = 20) [ report "Mid" ] [         ;; Mid valley floor
  ifelse (i = 25) [ report "Mid Dunes" ] [   ;; Mid valley dunes
  ifelse (i = 30) [ report "Natural" ] [     ;; Natural (non-arable)
  ifelse (i = 40) [ report "Uplands" ] [     ;; Uplands (arable)
  ifelse (i = 50) [ report "Kinbiko" ] [     ;; Kinbiko Canyon
  ifelse (i = 60) [ report "Empty" ] [ 
  report "?" ] ] ] ] ] ] ] ] ]               ;; there should be no "?" patches
end 

to make-region-label [ x y txt clr ]
  ask patch x y [
    set plabel-color clr
    set plabel txt
  ]
end