HIV and Polygamy

; 6/30/10 - added simple birth and death rates
; 6/30/10 - added plots
; 8/10/10 - more plots
;         - % on treatment (total; < 200;  200 - 350; 350 - 500)
;         - Incidence (total; due to acute; due to chronic)
;         - % infected (total; acute; chronic)
;         - more inputs
;         - treatment threshold (CD4 count)
;         - entry rate
;         - acute duration
;         - initial percent infected
;         - initial percent acute
;         - proportion on treatment

undirected-link-breed [
  non-spousal-links non-spousal-link
]
undirected-link-breed [
  spousal-links spousal-link
  ]

globals [
  Pyramid
  f-add-prob ; item i of f-add-prob applies when node in state item i of R
  ; item i is a list [a b c] where a is the probability of not adding any links, b of a spousal link, and c of a non-spousal link
  m-add-prob ; similar for males
  P_m ;Transition Probability Matrix
  P_f ;Transition Probability Matrix
  R ;List of possible relationship combinations
  MtoF-ratio
  male-remarry-pos-odds ;odds and ever-married? infected will have a new relationship
  female-remarry-pos-odds
  sero-sorting-odds ;odds one chooses partner with same sero-status (is everything symmetric???)
  ever-married-pos-choose-polygamous-odds   ; odds that ever-married HIV+/- women choose a husband who already has a partner
  ever-married-neg-choose-polygamous-odds
  F+M-divorce-odds ;relative risk of a discordant couple disolving
  F-M+divorce-odds
  relative-coital-freq
  number-of-female-nodes
  male-proportion-spousal 
  female-proportion-spousal
  male-proportion-0
  male-proportion-1
  male-proportion-2
  male-proportion-3
  male-proportion-4
  female-proportion-0
  female-proportion-1
  female-proportion-2
  female-proportion-3
  female-proportion-4
  entry-rate-m ;ignored here
  entry-rate-f
  spousal-change-rate
  non-spousal-change-rate
  t-hat
  t-hat-probs
  m15 m16 m17 m18 m19 m20
  m21 m22 m23 m24 m25 m26 m27 m28 m29 m30
  m31 m32 m33 m34 m35 m36 m37 m38 m39 m40
  m41 m42 m43 m44 m45 m46 m47 m48 m49
  f15 f16 f17 f18 f19 f20
  f21 f22 f23 f24 f25 f26 f27 f28 f29 f30
  f31 f32 f33 f34 f35 f36 f37 f38 f39 f40
  f41 f42 f43 f44 f45 f46 f47 f48 f49

  total-deaths
  curr-component-size ; for recursion
  component-size ; list of component sizes (one per component)
  
  outside-death-rate   
  current-yearly-incidence
  current-yearly-incidence-acute
  current-yearly-incidence-chronic
  current-yearly-incidence-final
  current-yearly-incidence->1-partner
  current-yearly-incidence-1-partner
  current-yearly-incidence-non-spousal
  current-yearly-incidence-spousal
  cum-infections
  cum-infections-acute
  
  run-number
]

to startup
    carefully [
      file-open "default.txt"
      while [not file-at-end?] [
        run file-read-line
      ]
      file-close
    ] [print error-message]
    
  set male-proportion-4 0
  set female-proportion-4 0
end 

turtles-own [
  randnum
  num-partners
  num-non-spousal-partners
  num-spousal-partners
  num-non-spousal-missing
  num-spousal-missing
  num-non-spousal-requests
  num-spousal-requests
  ever-married? ;actually both types of relationships
  age-years
  age-months
  female?
  infected?
  explored?
  
  ; HIV Variables
  stage-3-indicator
  months-since-infected
  months-left-to-live
  infectivity
  HIV-stage
  CD4-count-original
  CD4-count-present
  CD4-at-treatment
  CD4-decline-rate-no-treatment
  on-treatment?
  eligible-for-treatment?
  failed-treatment?
  time-on-treatment 
  HAART-supp-death-rate
]

to go
  age
  deaths-births
  progress-disease
  put-on-treatment
  spread-HIV
  change-relationships
  find-all-components
  plot-statistics
  tick
  ;if ticks = 120 [ stop ]
  ;if current-incidence != current-incidence-acute + current-incidence-chronic [ stop ]
end 

to setup-keep
  setup
  set run-number run-number + 1  
end 

to setup-clear
  clear-drawing
  clear-all-plots
  clear-output
  call-globals
  setup
  set run-number 1
end 

to setup
  clear-patches
  ask patches [set pcolor white] ; turns world white
  auto-plot-on
  set-default-shape turtles "person"
  
;  let tuned-prob-file (word scenario "T.txt")
;  ifelse (file-exists? tuned-prob-file) [
;    file-open tuned-prob-file
;    while [not file-at-end?] [ run file-read-line ]
;    file-close
;  ] [ ;otherwise tune probabilities and write them to file
    clear-turtles  
    reset-ticks
    setup-nodes
    initialize-HIV-population
    put-on-treatment
    ask non-spousal-links [set color black]
    ask spousal-links [set color red]
    ask turtles [set size 10 ]  
    repeat 120 [
      age
      deaths-births
      progress-disease
      put-on-treatment
      spread-HIV
      change-relationships
      tweak-probs
      tick
    ]
;    file-open tuned-prob-file
;    file-print (word "set m-add-prob " m-add-prob)
;    file-print (word "set f-add-prob " f-add-prob)
;    file-close
;  ]
  
  clear-turtles
  reset-ticks
  setup-nodes
  initialize-HIV-population
  put-on-treatment
  ask non-spousal-links [set color black]
  ask spousal-links [set color red]
  ask turtles [set size 10 ]
  set cum-infections 0
  set cum-infections-acute 0
  
  find-all-components
end 

; tweaks the new relationship probabilities to better match the desired distribution
; uses simulation optimization

to tweak-probs
  let m-props (list male-proportion-0 male-proportion-1 male-proportion-2 male-proportion-3)
  let f-props (list female-proportion-0 female-proportion-1 female-proportion-2 female-proportion-3)
  
  let f-curr-probs (n-values 4 [count turtles with [female? and num-partners = ?] / count turtles with [female?]])
  let m-curr-probs (n-values 4 [count turtles with [not female? and num-partners = ?] / count turtles with [not female?]])
  let f-ratio (map [(?2 + 1) / (?1 + 1)] f-props f-curr-probs)
  let m-ratio (map [(?2 + 1) / (?1 + 1)] m-props m-curr-probs)
  let new-m-add-prob []
  let new-f-add-prob []
  (foreach R m-add-prob f-add-prob [
    let partner-num (sum ?1)
    let f-ratio-num (item partner-num f-ratio)
    let m-ratio-num (item partner-num m-ratio)
    let m-prob-1 ((item 1 ?2) * m-ratio-num)
    let m-prob-2 ((item 2 ?2) * m-ratio-num)
    let f-prob-1 ((item 1 ?3) * f-ratio-num)
    let f-prob-2 ((item 2 ?3) * f-ratio-num)
    let m-prob (list (1 - m-prob-1 - m-prob-2) m-prob-1 m-prob-2)
    let f-prob (list (1 - f-prob-1 - f-prob-2) f-prob-1 f-prob-2)
    set new-m-add-prob (lput m-prob new-m-add-prob)
    set new-f-add-prob (lput f-prob new-f-add-prob)
  ])
  set m-add-prob new-m-add-prob
  set f-add-prob new-f-add-prob
end 

to ignore [stuff]
end 

to call-globals
  file-open (word scenario "P.txt")
  ignore file-read-line
  set P_f file-read
  ignore file-read-line
  set P_m file-read
  file-close
  
;  let Pyramid []
;  carefully [
;      file-open "default.txt"
;      while [not file-at-end?] [
;        run file-read-line
;      ]
;      file-close
;    ] [print error-message]
  file-open (word scenario ".txt")
  while [not file-at-end?] [
   run file-read-line
  ]
  file-close

  foreach Pyramid [
    let lower item 0 ?
    let upper item 1 ?
    let m item 2 ?
    let f item 3 ?
    ;show lower show upper show m show f
    (foreach (n-values (upper - lower + 1) [word "m" (? + lower)]) (n-values (upper - lower + 1) [word "f" (? + lower)]) [
        run ((word "set " ?1 word " " (m / (upper - lower + 1))))
        run ((word "set " ?2 word " " (f / (upper - lower + 1))))
    ])
  ]
  
  ; create m-add-prob and f-add-prob here (should probably be done in matlab and in the P files)
  set m-add-prob []
  set f-add-prob []
  (foreach R P_m P_f [ ; R entry is current state
    let new-state-1 (position (list ((item 0 ?1) + 1) (item 1 ?1)) R) ; index into R with extra spouse
    let new-state-2 (position (list (item 0 ?1) ((item 1 ?1) + 1)) R) ; index into R with extra non-spouse
    let m-prob-1 0
    let f-prob-1 0
    let m-prob-2 0
    let f-prob-2 0
    if new-state-1 != false [
      set m-prob-1 (item new-state-1 ?2)
      set f-prob-1 (item new-state-1 ?3)
    ]
    if new-state-2 != false [
      set m-prob-2 (item new-state-2 ?2)
      set f-prob-2 (item new-state-2 ?3)
    ]
    let m-prob (list (1 - m-prob-1 - m-prob-2) m-prob-1 m-prob-2)
    let f-prob (list (1 - f-prob-1 - f-prob-2) f-prob-1 f-prob-2)
    set m-add-prob (lput m-prob m-add-prob)
    set f-add-prob (lput f-prob f-add-prob)
  ])
  
  set current-yearly-incidence 0
  set current-yearly-incidence-acute 0
  set current-yearly-incidence-chronic 0
end 

to turtles-own-default
  set ever-married? (any? my-links)
  
  set infected? false

  ;; HIV Variables
  set months-since-infected 0
  set months-left-to-live 0
  set infectivity 0
  set HIV-stage 0
  set on-treatment? false
  set eligible-for-treatment? false
  set failed-treatment? false
  set time-on-treatment 0
  set HAART-supp-death-rate 0
  
  ;; Set CD4-count-present 
  set CD4-count-original exp ((ln 1116) + ((random-normal 0 1) * .303))
  set CD4-count-present CD4-count-original
  set CD4-decline-rate-no-treatment 0
  set stage-3-indicator 0
end 

to make-male
  setxy (random-xcor * .95) (random-ycor * .95)
  set color blue
  set female? false
end 

to make-female
  setxy (random-xcor * .95) (random-ycor * .95)
  set color pink
  set female? true
end 

to setup-nodes  
  crt number-of-male-nodes [
    make-male
  ]
  ask turtles [
      set num-partners 0
      set num-non-spousal-partners 0
      set num-spousal-partners 0
  ]
  let last-cut-off 0
  let current-cut-off round(male-proportion-0 * number-of-male-nodes)
  (foreach (list (male-proportion-1) (male-proportion-2) (male-proportion-3) (male-proportion-4)) ([1 2 3 4])  [
      set last-cut-off current-cut-off
      set current-cut-off (last-cut-off + round(number-of-male-nodes * (?1 / (1 - male-proportion-0)) * (1 - male-proportion-spousal - male-proportion-0)))  
      ;show (word last-cut-off " " current-cut-off)   
      ask turtles with [who + 1 > last-cut-off and who + 1 <= current-cut-off] [
        set num-partners ?2
        set num-non-spousal-partners ?2
        set num-spousal-partners 0
      ]
  ])
  (foreach (list (male-proportion-1) (male-proportion-2) (male-proportion-3) (male-proportion-4)) ([0 1 2 3])  [
      set last-cut-off current-cut-off
      set current-cut-off (last-cut-off + round(number-of-male-nodes * (?1 / (1 - male-proportion-0)) * male-proportion-spousal))  
      ;show (word last-cut-off " " current-cut-off)     
      ask turtles with [who + 1 > last-cut-off and who + 1 <= current-cut-off] [
        set num-partners ?2 + 1
        set num-non-spousal-partners ?2
        set num-spousal-partners 1
      ]
  ])

  ; create female turtles to match number of links of male turtles
  let number-of-male-links (count turtles with [not female? and num-partners = 1]) + 2 * (count turtles with [not female? and num-partners = 2]) + 3 * (count turtles with [not female? and num-partners = 3]) + 4 * (count turtles with [not female? and num-partners = 4])
  set number-of-female-nodes round(number-of-male-links / ((female-proportion-1 + 2 * female-proportion-2 + 3 * female-proportion-3 + 4 * female-proportion-4)))
  ; define female-proportion-spousal
  set female-proportion-spousal ((number-of-male-nodes * male-proportion-spousal) / number-of-female-nodes)
  crt number-of-female-nodes [
    make-female
  ]
  ask turtles with [female?] [
      set num-partners 0
      set num-non-spousal-partners 0
      set num-spousal-partners 0
  ]
  set last-cut-off 0
  set current-cut-off round(female-proportion-0 * number-of-female-nodes)
  (foreach (list (female-proportion-1) (female-proportion-2) (female-proportion-3) (female-proportion-4)) ([1 2 3 4])  [
      set last-cut-off current-cut-off
      set current-cut-off (last-cut-off + round(number-of-female-nodes * (?1 / (1 - female-proportion-0)) * (1 - female-proportion-spousal - female-proportion-0)))  
      ;show (word last-cut-off " " current-cut-off)   
      ask turtles with [who + 1 - number-of-male-nodes > last-cut-off and who + 1 - number-of-male-nodes <= current-cut-off] [
        set num-partners ?2
        set num-non-spousal-partners ?2
        set num-spousal-partners 0
      ]
  ])
  (foreach (list (female-proportion-1) (female-proportion-2) (female-proportion-3) (female-proportion-4)) ([0 1 2 3])  [
      set last-cut-off current-cut-off
      set current-cut-off (last-cut-off + round(number-of-female-nodes * (?1 / (1 - female-proportion-0)) * female-proportion-spousal))  
      ;show (word last-cut-off " " current-cut-off)     
      ask turtles with [who + 1 - number-of-male-nodes > last-cut-off and who + 1 - number-of-male-nodes <= current-cut-off] [
        set num-partners ?2 + 1
        set num-non-spousal-partners ?2
        set num-spousal-partners 1
      ]
  ])  
  ask turtles [
    set num-non-spousal-missing num-non-spousal-partners
    set num-spousal-missing num-spousal-partners
    set age-years -1
    set age-months -1
  ]
  
  let number-of-female-links (count turtles with [female? and num-partners = 1]) + 2 * (count turtles with [female? and num-partners = 2]) + 3 * (count turtles with [female? and num-partners = 3]) + 4 * (count turtles with [female? and num-partners = 4])
  ;show (word number-of-female-links " " number-of-male-links)
  
  ask turtles with [num-non-spousal-missing > 0 or num-spousal-missing > 0] [
    while [num-spousal-missing > 0 and (any? other turtles with [(female? = not ([female?] of myself)) and num-spousal-missing > 0 and (link-neighbor? myself = false)])] [
      ask one-of turtles with [(female? = not ([female?] of myself)) and num-spousal-missing > 0] [
        create-spousal-link-with myself [
          ask both-ends [set num-spousal-missing num-spousal-missing - 1]
        ]
      ]
    ]
    while [num-non-spousal-missing > 0 and (any? other turtles with [(female? = not ([female?] of myself)) and num-non-spousal-missing > 0 and (link-neighbor? myself = false)])] [
      ask one-of turtles with [(female? = not ([female?] of myself)) and num-non-spousal-missing > 0] [
        create-non-spousal-link-with myself [
          ask both-ends [set num-non-spousal-missing num-non-spousal-missing - 1]
        ]
      ]
    ]   
  ] 
  
  ask turtles with [num-non-spousal-missing > 0] [
    set num-partners num-partners - num-non-spousal-missing
    set num-non-spousal-partners num-non-spousal-partners - num-non-spousal-missing
    set num-spousal-missing 0
  ]
  ask turtles with [num-spousal-missing > 0] [
    set num-partners num-partners - num-spousal-missing
    set num-spousal-partners num-spousal-partners - num-spousal-missing
    set num-non-spousal-missing 0
  ]
  
  (foreach (n-values 35 [runresult (word "m" (? + 15))]) (n-values 35 [? + 15]) [
      ask n-of (min (list (round(number-of-male-nodes * ?1)) (count turtles with [age-years = -1 and not female?]))) turtles with [age-years = -1 and not female?] [
        set age-months (?2 * 12 + (random 12))
        set age-years (age-months / 12)
      ]
  ])
  (foreach (n-values 35 [runresult (word "f" (? + 15))]) (n-values 35 [? + 15]) [
      ask n-of (min (list (round(number-of-female-nodes * ?1)) (count turtles with [age-years = -1 and female?]))) turtles with [age-years = -1 and female?] [
        set age-months (?2 * 12 + (random 12))
        set age-years (age-months / 12)
      ]
  ])

  ask non-spousal-links [set color black]
  ask spousal-links [set color red]
  
  
  ;; Initialize turtles with default turtles-own values
  ask turtles [
    turtles-own-default
  ]
end 

to change-relationships
  ; setup book keeping
  ask turtles [
    set num-spousal-requests false
    set num-non-spousal-requests false
  ]
  
  ;deleting links
  ask non-spousal-links [
    let multiplier 1
    if (any? both-ends with [infected? and female?]) and (any? both-ends with [not infected? and not female?]) [
      set multiplier F+M-divorce-odds
    ]
    if (any? both-ends with [not infected? and female?]) and (any? both-ends with [infected? and not female?]) [
      set multiplier F-M+divorce-odds
    ]
    if ((random-float 1) < non-spousal-change-rate * multiplier) [
        die
    ]
  ]
  ask spousal-links [
    let multiplier 1
    if (any? both-ends with [infected? and female?]) and (any? both-ends with [not infected? and not female?]) [
      set multiplier F+M-divorce-odds
    ]
    if (any? both-ends with [not infected? and female?]) and (any? both-ends with [infected? and not female?]) [
      set multiplier F-M+divorce-odds
    ]
    if ((random-float 1) < spousal-change-rate * multiplier) [
        die
    ]
  ]
  
  ; creating requests for new links
  ask turtles [
    let rownum (position (list num-spousal-partners num-non-spousal-partners) R)
    let pmf []
    ifelse female? [
      set pmf (item rownum f-add-prob)
    ] [
      set pmf (item rownum m-add-prob)
    ]
    if ever-married? [
      let multiplier 1
      if infected? [ifelse female? [set multiplier female-remarry-pos-odds] [set multiplier male-remarry-pos-odds]]
      let p0 (item 0 pmf)
      set p0 (p0 / ((1 - p0) * multiplier + p0))
      let p1 0
      if p0 < 1 [set p1 ((1 - p0) * ((item 1 pmf) / (1 - (item 0 pmf))))]
      set pmf (list p0 p1 (1 - p0 - p1))
    ]
    let change (random-discrete pmf)
    if change = 1 [ set num-spousal-requests true ]
    if change = 2 [ set num-non-spousal-requests true ]
  ]
  
  ; create links
  ask turtles with [ num-spousal-requests or num-non-spousal-requests ] [
    if num-spousal-requests [
      set num-spousal-requests false
      let eligible turtles with [ (female? = not [female?] of myself) ; opposite gender
          and num-spousal-requests and (link-neighbor? myself = false)]
      if any? eligible [
        if female? and ever-married? [set eligible (ever-married-filter eligible)]
        set eligible (sero-sorting-filter eligible)
        ask one-of eligible [
          create-spousal-link-with myself
          set num-spousal-requests false
        ]
      ]
    ]
    if num-non-spousal-requests [
      set num-non-spousal-requests false
      let eligible turtles with [ (female? = not [female?] of myself) ; opposite gender
          and num-non-spousal-requests and (link-neighbor? myself = false)]
      if any? eligible [
        if female? and ever-married? [set eligible (ever-married-filter eligible)]
        set eligible (sero-sorting-filter eligible)
        ask one-of eligible [
          create-non-spousal-link-with myself
          set num-non-spousal-requests false
        ]
      ]
    ]
  ]
  
  ; update book keeping
  ask turtles [
    set num-partners (count my-links)
    set num-spousal-partners (count my-spousal-links)
    set num-non-spousal-partners (count my-non-spousal-links)
    set ever-married? (ever-married? or num-partners >= 1)
  ]
  
  ask non-spousal-links [set color black]
  ask spousal-links [set color red]
end 

to initialize-HIV-population
  ; This is the acute infection population
  let initial-percent-acute (item 0 t-hat-probs) * (acute-duration / 12)
  set t-hat-probs replace-item 0 t-hat-probs ((item 0 t-hat-probs) - initial-percent-acute)
  let adjustment-factor sum t-hat-probs
  
  (foreach (n-values acute-duration [?]) [
        ask n-of round ((initial-percent-acute / acute-duration) * ((initial-percent-infected) * (count turtles))) turtles [
          set months-since-infected ?1
          become-infected
          set months-left-to-live (random-weibull 10.9522 2.25) + (acute-duration - ?1)
          set stage-3-indicator (months-left-to-live - acute-duration + ?1) * .1
          set CD4-decline-rate-no-treatment (CD4-count-present * .75) / (months-left-to-live - acute-duration + ?1)
          set CD4-count-present CD4-count-present - (?1 * ((CD4-count-original * .25) / acute-duration))
        ]
  ])
  ; Now, we'll infect those in the chronic/final stages of HIV by...
  ; Determining how long ago when they seroconverted
  ; Determining how much longer they have to live
  ; Then determining the rate of CD4 decline and starting CD4 value by 
  ; how much longer they are scheduled to live 


  ask n-of ((initial-percent-infected * count turtles) - count turtles with [infected?]) turtles with [not infected?] [
    become-infected
    let lower 0
    let upper 0
    set randnum random-float adjustment-factor
    (foreach t-hat-probs t-hat [
      set lower upper
      set upper upper + ?1
      if randnum < upper and randnum >= lower [
        set months-since-infected ?2 * 12
        set months-left-to-live (12 * 10.5922 * ((((?2 / 10.5922) ^ (2.25)) - ln (1 - random-float 1)) ^ (1 / 2.25))) - (?2 * 12)
        set CD4-decline-rate-no-treatment (CD4-count-original * .75) / (months-since-infected + months-left-to-live)
        set CD4-count-present (CD4-count-original * .75) - (CD4-decline-rate-no-treatment * months-since-infected) 
        set stage-3-indicator .1 * (months-left-to-live + (?2 * 12))
      ]
    ])      
  ]
end 

to age
  ask turtles [
    set age-months age-months + 1
    set age-years age-months / 12
  ]
end 

to deaths-births
  let num-deaths 0
  ask turtles [
    ;; Note: deaths are not age/gender specific
    ; (not any? link-neighbors) and ; 
    if ((random-float 1 < outside-death-rate + HAART-supp-death-rate) or (CD4-count-present <= 0) or (age-years >= 50)) [
      ask my-links [
        die
      ]
      set num-deaths (num-deaths + 1)
      set age-months 15 * 12
      set age-years 15
      turtles-own-default
    ]
  ]
  ;; Remember to reset num-partners, etc. after links disappear (this needs to be in separate block
  ;; since (ask turtles) performs code for each turtle and not in unison
  ask turtles [
    set num-partners (count my-links)
    set num-spousal-partners (count my-spousal-links)
    set num-non-spousal-partners (count my-non-spousal-links)
  ]
  
  set total-deaths (total-deaths + num-deaths)
  set-current-plot "births and deaths"
  plot num-deaths
end 

to progress-disease
  ask turtles with [infected?] [
    determine-HIV-stage
    determine-CD4-count
    determine-infectivity
    determine-supp-death-rates
    set months-since-infected months-since-infected + 1
    set months-left-to-live months-left-to-live - 1
  ]
end 

to determine-HIV-stage
  if (months-since-infected < acute-duration) [
    set HIV-stage 1
  ]
  if (months-since-infected >= acute-duration) [
    if months-left-to-live < stage-3-indicator [ set HIV-stage 3 ]
    if months-left-to-live >= stage-3-indicator [ set HIV-stage 2 ]
  ]
end 

to determine-infectivity
  ifelse (on-treatment?) 
  [ 
    set infectivity chronic-infectivity * treatment-infectivity-ratio
  ]
  [ 
    if (HIV-stage = 1) [ set infectivity chronic-infectivity * acute-infectivity-ratio ]
    if (HIV-stage = 2) [ set infectivity chronic-infectivity ]
    if (HIV-stage = 3) [ set infectivity chronic-infectivity * final-infectivity-ratio ]
  ]
  if (num-partners > 1) [
    set infectivity infectivity * (item (num-partners - 2) relative-coital-freq)
  ]
  if (not female?) [
    set infectivity infectivity * MtoF-ratio
  ]
end  

to determine-supp-death-rates
  ask turtles with [infected?] [
    ifelse (on-treatment?) [
      if CD4-at-treatment <= 25 [ set HAART-supp-death-rate .175 / 12 ]
      if CD4-at-treatment > 25 and CD4-at-treatment <= 50 [ set HAART-supp-death-rate .121 / 12 ]
      if CD4-at-treatment > 50 and CD4-at-treatment <= 100 [ set HAART-supp-death-rate .738 / 12 ]
      if CD4-at-treatment > 100 and CD4-at-treatment <= 200 [ set HAART-supp-death-rate .483 / 12 ]
      if CD4-at-treatment > 200 and CD4-at-treatment <= 250 [ set HAART-supp-death-rate .05 / 12 ]
      if CD4-at-treatment > 250 and CD4-at-treatment <= 350 [ set HAART-supp-death-rate .05 / 12 ]
      if CD4-at-treatment > 350 and CD4-at-treatment <= 500 [ set HAART-supp-death-rate .03 / 12 ]
      if CD4-at-treatment > 500 [ set HAART-supp-death-rate .015 / 12 ]
    ]
    [
      set HAART-supp-death-rate 0
    ]
  ]
end 

to determine-CD4-count
  ifelse (on-treatment?) [
    ifelse (failed-treatment?) 
    [ set CD4-count-present CD4-count-present - CD4-decline-rate-no-treatment ]
    [ set CD4-count-present CD4-count-present ]
  ]
  [ 
    ifelse (HIV-stage = 1)
      [ set CD4-count-present CD4-count-present - ((CD4-count-original * .25) / acute-duration) ]
      [ set CD4-count-present CD4-count-present - CD4-decline-rate-no-treatment ]
  ]
end 

to spread-HIV
  if ((ticks mod 12) = 0) [
    set current-yearly-incidence-acute 0
    set current-yearly-incidence-chronic 0
    set current-yearly-incidence-final 0
    set current-yearly-incidence->1-partner 0    
    set current-yearly-incidence-1-partner 0
    set current-yearly-incidence-non-spousal 0
    set current-yearly-incidence-spousal 0
    set current-yearly-incidence 0
  ]
  
  ask turtles with [infected? = true] [
    ask link-neighbors with [infected? = false] [
      if random-float 1 < ([infectivity] of myself) [
        become-infected
        set months-left-to-live (random-weibull 10.9522 2.25) + acute-duration
        set stage-3-indicator (months-left-to-live - acute-duration) * .1
        set CD4-decline-rate-no-treatment (CD4-count-present * .75) / (months-left-to-live - acute-duration)
        set cum-infections cum-infections + 1
        if (([HIV-stage] of myself) = 1) [
          set current-yearly-incidence-acute current-yearly-incidence-acute + 1
          set cum-infections-acute cum-infections-acute + 1
          ]
        if (([HIV-stage] of myself) = 2) [
          set current-yearly-incidence-chronic current-yearly-incidence-chronic + 1
          ]
        if (([HIV-stage] of myself) = 3) [
          set current-yearly-incidence-final current-yearly-incidence-final + 1
        ]
        if (([num-partners] of myself) > 1) [
          set current-yearly-incidence->1-partner current-yearly-incidence->1-partner + 1
          ]
        if (([num-partners] of myself) = 1) [
          set current-yearly-incidence-1-partner current-yearly-incidence-1-partner + 1
          ]
        if (non-spousal-link-neighbor? myself) [
          set current-yearly-incidence-non-spousal current-yearly-incidence-non-spousal + 1 
        ]      
        if (spousal-link-neighbor? myself) [
          set current-yearly-incidence-spousal current-yearly-incidence-spousal + 1  
        ]
      ]
    ]
  ]
  
  ;; Calculate Incidence
  set current-yearly-incidence current-yearly-incidence + (count turtles with [infected? and months-since-infected = 0])
end 

to become-infected
  set infected? true
  if female? = true [set color pink - 2]
  if female? = false [set color blue - 2]
end 

to put-on-treatment
  ; tick time-on-treatment variable for those on treatment
  ask turtles with [on-treatment?] [
    set time-on-treatment time-on-treatment + 1
    ; determine if the treatment has failed
    ;if (not failed-treatment?) [
     ; if (random-float 1 < treatment-failure-rate) [
     ;   set failed-treatment? true
     ;   set months-left-to-live CD4-count-present / CD4-decline-rate-no-treatment
     ; ]
    ;]
  ]
  
  ask turtles with [infected? and months-left-to-live > 0 and (CD4-count-present < treatment-threshold) and (not failed-treatment?)] [
    set eligible-for-treatment? true
  ]
  
  ;; a constant (maximum) fraction of the population will be on treatment
  if (any? turtles with [eligible-for-treatment?]) [
  while [(count turtles with [on-treatment?] / count turtles with [eligible-for-treatment?]) < (proportion-on-treatment)] [
    ask one-of turtles with [eligible-for-treatment? and (not on-treatment?)] [
      set on-treatment? true
      set CD4-at-treatment CD4-count-present
    ]
  ]
  ]
end 

to make-prettier
  ask links [set color white]
  ask turtles [set size 3]
  repeat 10
  [
    layout-spring turtles links .8 (world-width / (sqrt count turtles)) 4
  ]
end 

to set-plot-pen [plot-pen keep?]
  if keep? [
  create-temporary-plot-pen word plot-pen run-number
  set-plot-pen-color black
  ]
  if (not keep?) [
    set-current-plot-pen plot-pen
    if ticks = 0 [ plot-pen-reset ]
  ]
end 

to plot-statistics
  set-current-plot "Women / Men"
  set-current-plot-pen "overall"
  plot count turtles with [female?] / count turtles with [not female?]
  set-current-plot-pen "prevalence"
  plot count turtles with [female? and infected?] / count turtles with [not female? and infected?]
  
  set-current-plot "Prevalence"
  set-plot-pen "All-HIV-Positive" true
  ;set-current-plot-pen "All-HIV-Positive"
  if (not draw-all-HIV-positives?) [ plot-pen-up ]
  if any? turtles with [infected?] [
    plot count turtles with [infected?] / count turtles
  ]
  set-plot-pen "Acute" false
  if (not draw-acute?) [ plot-pen-up ]
  ifelse any? turtles with [HIV-stage = 1 and infected?] [
    plot count turtles with [HIV-stage = 1 and infected?] / count turtles
  ]
  [ plot 0 ]
  set-plot-pen "Chronic" false
  if (not draw-chronic?) [ plot-pen-up ]
  ifelse any? turtles with [HIV-stage = 2 and infected?] [
    plot count turtles with [HIV-stage = 2 and infected?] / count turtles
  ]
  [ plot 0 ]
  set-plot-pen "Final" false
  if (not draw-final?) [ plot-pen-up ]
  ifelse any? turtles with [HIV-stage = 3 and infected?] [
    plot count turtles with [HIV-stage = 3 and infected?] / count turtles
  ]
  [ plot 0 ]
  
  if (ticks mod 12 = 0) [
    set-current-plot "Incidence"
    set-plot-y-range 0 precision ((current-yearly-incidence / count turtles) + .05) 2
    set-plot-pen "Incidence" false
    plot (current-yearly-incidence / count turtles)
  
    set-current-plot "Proportion of Incidence Due to HIV Stages"
    set-plot-pen "Acute" false
    ifelse current-yearly-incidence != 0
    [plot current-yearly-incidence-acute / current-yearly-incidence]
    [plot 0]
    set-plot-pen "Chronic" false
    ifelse current-yearly-incidence != 0
    [plot current-yearly-incidence-chronic / current-yearly-incidence]
    [plot 0]
    set-plot-pen "Final" false
    ifelse current-yearly-incidence != 0
    [plot current-yearly-incidence-final / current-yearly-incidence]
    [plot 0]
    
    set-current-plot "Proporition of Incidence Due to Types of Partnerships"
    set-plot-pen "Non-Spousal" false
    ifelse current-yearly-incidence != 0
    [plot current-yearly-incidence-non-spousal / current-yearly-incidence]
    [plot 0]
    set-plot-pen "Spousal" false
    ifelse current-yearly-incidence != 0
    [plot current-yearly-incidence-spousal / current-yearly-incidence]
    [plot 0]

    set-current-plot "Proportion of Incidence Due to Concurrent Partnerships"
    set-plot-pen ">1 Partner" false
    ifelse current-yearly-incidence != 0
    [plot current-yearly-incidence->1-partner / current-yearly-incidence]
    [plot 0]
    set-plot-pen "1 Partner" false
    ifelse current-yearly-incidence != 0
    [plot current-yearly-incidence-1-partner / current-yearly-incidence]
    [plot 0]
    
  ]

  set-current-plot "Male Partnership Distribution"
  foreach n-values 4 [?] [
    set-plot-pen (word ? " Partners") false
    plot count turtles with [(not female?) and (count my-links = ?)] / count turtles with [(not female?)]
  ]
  set-current-plot "Female Partnership Distribution"
  foreach n-values 4 [?] [
    set-plot-pen (word ? " Partners") false
    plot count turtles with [(female?) and (count my-links = ?)] / count turtles with [(female?)]
  ]
  
  set-current-plot "Histogram of Years Left to Live"
  histogram [months-left-to-live / 12] of turtles with [infected?]
  set-current-plot "Histogram of Years Since Infected"
  histogram [months-since-infected / 12] of turtles with [infected?]
  
  set-current-plot "Mean CD4 Count at Treatment Initiation"
  set-plot-pen "CD4-count" false
  ifelse not any? turtles with [infected? and on-treatment?] 
  [ plot 0 ]
  [ plot mean [CD4-at-treatment] of turtles with [infected? and on-treatment?] ]
end 

to-report random-weibull [a b]
 report (a * ((- ln (1 - random-float 1)) ^ (1 / b)) * 12)
end 

to-report random-discrete [pmf] ;returns discrete random variable with given pmf (0-based)
  let itemnum 0
  let lower 0
  let upper 0
  let randval random-float 1
  foreach pmf [
    set upper upper + ?
    if (randval < upper and randval >= lower) [
      report itemnum
    ]
    set lower upper
    set itemnum (itemnum + 1)
  ]
end 

to-report ever-married-filter [eligible]
  ;we assume eligible is a non-empty agentset
  ;guaranteed to return a non-empty subset of eligible
  let withpartner eligible with [num-partners >= 1]
  let nopartner eligible with [num-partners = 0]
  let prob 0
  ifelse infected?
     [ set prob ((count nopartner) / (ever-married-pos-choose-polygamous-odds * (count withpartner) + (count nopartner)))]
     [ set prob ((count nopartner) / (ever-married-neg-choose-polygamous-odds * (count withpartner) + (count nopartner)))]
  ifelse ((random-float 1)  >= prob) [report withpartner] [report nopartner]
end 

to-report sero-sorting-filter [eligible]
  ;we assume eligible is a non-empty agentset
  ;guaranteed to return a non-empty subset of eligible
  let same eligible with [infected? = [infected?] of myself]
  let diff eligible with [infected? = not [infected?] of myself]
  let prob ((count diff) / (sero-sorting-odds * (count same) + (count diff)))
  ifelse ((random-float 1)  >= prob) [report same] [report diff]
end 

; based on Giant Component code from model library
;; to find all the connected components in the network, their sizes and starting turtles

to find-all-components
  set component-size []
  ask turtles [ set explored? false ]
  ;; keep exploring till all turtles get explored
  loop [
    ;; pick a node that has not yet been explored
    let start one-of turtles with [ not explored? ]
    if start = nobody [ stop ]
    ;; reset the number of turtles found to 0
    ;; this variable is updated each time we explore an
    ;; unexplored node.
    set curr-component-size 0
    ask start [ explore ]
    set component-size (lput curr-component-size component-size)
  ]
end 

;; Finds all turtles reachable from this node

to explore ;; node procedure
  if explored? [ stop ]
  set explored? true
  set curr-component-size curr-component-size + 1
  ask link-neighbors [ explore ]
end