marriage-age-gap-simulation.lsp

#!/usr/bin/newlisp -n

;;; marriage-age-gap-simulation.lsp Copyright (c) 2016 TheMost
;;;
;;; Permission to use, copy, modify, and/or distribute this software for
;;; any purpose with or without fee is hereby granted, provided that the
;;; above copyright notice and this permission notice appear in all
;;; copies.
;;;
;;; THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
;;; WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
;;; WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
;;; AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
;;; DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA
;;; OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
;;; TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
;;; PERFORMANCE OF THIS SOFTWARE.
;;;
;;; to run:
;;; Windows: newlisp.exe marriage-age-gap-simulation.lsp
;;; Mac OSX: newlisp ./marriage-age-gap-simulation.lsp
;;; Linux:   newlisp ./marriage-age-gap-simulation.lsp
;;;
;;; The polygamy problem: If some men have more than one wife, other men will
;;; have to go without!  Is this problem real?  Or does it depend on false
;;; assumptions?
;;;
;;; Assuming an age gap between the age at which men and women marry, how many
;;; men can have 2 wives?
;;;
;;; Assuming every man has a certain number of wives, what age gap is required?
;;; Also, how many times will a woman marry, compared to a man?
;;;
;;; First model:
;;;   * everyone lives their full lifespan
;;;   * every is sufficiently healthy and attractive to mate
;;;   * no abortions, miscarriages, murders
;;;   * marriage is for life, no divorce
;;; Future models: include hypergamy, SMV, more realistic death
;;; model.

(define def define)
(define defm define-macro)

(def total-years 2000)        ; number of years to run the simulation
(def max-population 1000000L) ; end simulation when this number of people is reached
(def num-wives 2)             ; number of wives per man
(def average-number-of-children-per-woman 8)
(def male-age-of-death       72)
(def female-age-of-death     79)
(def male-age-of-majority    18)
(def female-age-of-majority  16)
(def female-age-of-menopause 48)
;(def female-age-of-widowhood 60)
(def female-age-of-widowhood 34) ; women don't remarry after this age
(def birth-odds (div average-number-of-children-per-woman
                  (sub female-age-of-menopause female-age-of-majority)))

;; source for information about sex ratio by fathers age:
;; https://www.psychologytoday.com/blog/the-scientific-fundamentalist/201104/why-are-older-parents-more-likely-have-daughters
;; by Satoshi Kanazawa, posted April 18, 2011
(def male-birth-ratio-by-fathers-age
  (list
     (list 20 0.5329)
     (list 25 0.5121)
     (list 30 0.4840)
     (list 35 0.5224)
     (list 40 0.5256)
     (list male-age-of-death 0.3592)))

(def male-birth-ratio
  (lambda (x)
    (let (a (- yr ((males (husband-of x)) 0))) ; a is the father's age
      (catch
        (dolist (b male-birth-ratio-by-fathers-age)
          (if (< a (b 0))
            (throw (b 1))))))))

(def init-years female-age-of-widowhood) ; number of years to generate people ex nihilo

(def max-32bit-int 4294967295) ; in binary, 0b11111111111111111111111111111111
(def fd            (open "/dev/urandom" "read"))
(seed (time-of-day)) ; initialize random number generator
(def /dev/urandom
  (if (= ostype "Windows")
      random ; fallback to builtin RNG on Windows
    (lambda ()
      (read fd buf 4)
      (div (first (unpack "lu" buf)) max-32bit-int))))

(def males (array (+ 1 max-population)))   ; (birthyear (current-wives) nwives ...)
(def females (array (+ 1 max-population))) ; (birthyear current-husband nhusbands ...)
(def oldest-man 0)
(def youngest-male 0)
(def youngest-man 0)  ; a man can marry and sire children
(def oldest-woman 0)
(def youngest-female 0)
(def youngest-woman 0) ; a woman can marry and bear children

(def wives-of   (lambda (x) ((males x) 1)))
(def husband-of (lambda (x) ((females x) 1)))
(def married?   (lambda (x) (husband-of x)))

;; yearly stats
(def mdeaths 0) ; male deaths
(def fdeaths 0) ; female deaths
(def mbirths 0) ; male births
(def fbirths 0) ; female births
(def mages 0)   ; used to calculate average male age at marriage
(def fages 0)   ; used to calculate average female age at marriage
(def nmarriages 0) ; number of marriages

;; advances indexes into the arrays of people
;; this let's us convert some iteration into simple arithmetic
;; and thereby avoids memory accesses and cache misses/flushes
(def advance-age-status
  (lambda (yr)
    ;; youngest-male and youngest-female are updated in do-births
    ;; oldest-man and oldest-woman are updated in do-deaths
    (let (y (- yr male-age-of-majority))    (while (> y ((males youngest-man) 0))          (++ youngest-man)))
    (let (y (- yr female-age-of-majority))  (while (> y ((females youngest-woman) 0))      (++ youngest-woman)))
    (let (y (- yr female-age-of-menopause)) (while (> y ((females oldest-fertile) 0))      (++ oldest-fertile)))
    (let (y (- yr female-age-of-widowhood)) (while (> y ((females oldest-marriageable) 0)) (++ oldest-marriageable)))
  ))

(def initialize-population
  (lambda ()
    (dotimes (yr init-years)
      (setf (males yr)   (list yr '()))
      (setf (females yr) (list yr nil)))
    (setf youngest-male   (- init-years 1)
          youngest-female (- init-years 1)
          youngest-man        0
          youngest-woman      0
          oldest-man          0
          oldest-woman        0
          oldest-fertile      0
          oldest-marriageable 0)))

(def (do-births yr)
  (let (i oldest-fertile yf youngest-female ym youngest-male)
    (while (<= i youngest-woman)
      (when (and (married? i) (< (/dev/urandom) birth-odds))
        (when (>= (max youngest-female youngest-male) max-population)
          (println "Ran out of memory. Stopping now.")
          (exit))
        (if (< (/dev/urandom) (male-birth-ratio i))
          (begin (setf   (males (++ youngest-male))   (list yr '())))
          (begin (setf (females (++ youngest-female)) (list yr nil)))))
      (++ i))
    (set 'fbirths (- youngest-female yf))
    (set 'mbirths (- youngest-male ym))))

(def (do-deaths yr)
  (let (af (- yr female-age-of-death)
        am (- yr male-age-of-death)
        om oldest-man
        ow oldest-woman)
    (while (> af ((females oldest-woman) 0))
      (when (married? oldest-woman)
        (setf ((males (husband-of oldest-woman)) 1)
              (clean (fn (x) (= x oldest-woman)) (wives-of (husband-of oldest-woman)))))
      (++ oldest-woman))
    (while (> am ((males oldest-man) 0))
      (dolist (wife (wives-of oldest-man))
        (setf ((females wife) 1) nil))
      (++ oldest-man))
    (++ fdeaths (- oldest-woman ow))
    (++ mdeaths (- oldest-man om))))

;; a man will only marry one woman per year
(def (do-marriages yr)
  (let (i oldest-man k youngest-woman)
    (while (and (<= i youngest-man) (>= k oldest-marriageable))
      (when (< (length ((males i) 1)) num-wives) ; man doesn't have his full complement of wives
        (while (and ((females k) 1) (>= k oldest-marriageable))
          (-- k))
        (unless (married? k) ; found a single woman, marry her off.
          (setf ((females k) 1) i)
          (push k ((males i) 1))
          (++ fages (- yr ((females k) 0)))
          (++ mages (- yr ((males i) 0)))
          (++ nmarriages)))
      (++ i))))

(def (print-year yr)
  (let (num-men (- youngest-male oldest-man -1) num-women (- youngest-female oldest-woman -1))
    (println (format
      "Year %d. population %d (%.2f) %d births (%.2f) %d deaths (%.2f) %d marriages (m: %d f: %d)"
      yr
      (+ num-men num-women) (div num-men (add num-men num-women))
      (+ mbirths fbirths) (div mbirths (add mbirths fbirths))
      (+ mdeaths fdeaths) (div mdeaths (add mdeaths fdeaths))
      nmarriages (round (div mages nmarriages)) (round (div fages nmarriages)))))
  (set 'mdeaths 0 'fdeaths 0 'mbirths 0 'fbirths 0 'mages 0 'fages 0 'nmarriages 0))

(initialize-population)
(def (do-year yr)
  (when (and (> youngest-male oldest-man) (> youngest-female oldest-woman)) ; while there are living people
      (advance-age-status yr)
      (do-births yr)
      (do-deaths yr)
      (do-marriages yr)
      (print-year yr)))

(for (yr (+ 1 init-years) (+ init-years total-years))
  (do-year yr)))

(exit)