#!/usr/bin/newlisp

;;; marriage-age-gap-simulation.lsp Copyright (c) 2016 Mycroft Jones
;;;
;;; to run: newlisp.exe 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
;;;   * no abortions, miscarriages, murders
;;;   * marriage is for life, no divorce
;;;   * if father is over 40, female birthrate goes to 0.6408

(define maxwives 2)
(define male-age-of-death 72) ; 72
(define female-age-of-death 79)
(define male-age-of-majority 18)
(define female-age-of-majority 16)
(define female-age-of-menopause 42)
(define female-age-of-widowhood 60)
(define male-birth-ratio (ratio 1 1.05))
(define ratio-after-40 0.6408)
(define average-number-of-children-per-woman 8)
(define birth-odds (div average-number-of-children-per-woman
                        (sub female-age-of-menopause female-age-of-majority)))

(context '/dev/urandom)
(define max-32bit-int 4294967295)
(set 'fd (open "/dev/urandom" "read"))
(define (/dev/urandom:/dev/urandom)
  (read fd buf 4)
  (div (max 1 (first (unpack "lu" buf))) max-32bit-int))
(context MAIN)

;(define-macro (nq a b) (push (eval a) (eval b) -1))
(define-macro (nq a b) (eval (list 'push a b -1)))
(define-macro ($# a) (+ -1 (length (eval a))))

(define males (list))
(define females (list))
(define everyone (list)) ; (m birthyear ((wife-id marriage-year) ...) (child-id ...) deathyear (current-wives))
                         ; (f birthyear ((husband-id marriage-year) ...) (child-id ...) deathyear current-husband)

(define (age x) (- yr (everyone x 1)))

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

(define initialize-population
  (lambda ()
    (dotimes (yr (max male-age-of-majority female-age-of-majority))
      (nq (list 'm yr '() '() nil '()) everyone)
      (nq (+ -1 (length everyone)) males)
      (nq (list 'f yr '() '() nil nil) everyone)
      (nq (+ -1 (length everyone)) females)
    )))

;; married women can give birth between ages of 18 and 44
;; average number of children, 5
;; let's say a 1 in 5 chance of giving birth in any given year

(define (do-births yr)
  (letn (bfemales (filter (fn (x) (and (everyone x 5) (> (age x) female-age-of-majority) (< (age x) female-age-of-menopause) (< (/dev/urandom) birth-odds))) females)
         gfemales (filter (fn (x) (< (/dev/urandom) (if (< (age x) 40) male-birth-ratio ratio-after-40))) bfemales)
         mfemales (difference bfemales gfemales)
         leveryone (+ -1 (length everyone)))
    (dotimes (i (length mfemales)) (push (list 'm yr '() '() nil '()) everyone -1) (push (++ leveryone) males -1))
    (dotimes (i (length gfemales)) (push (list 'f yr '() '() nil nil) everyone -1) (push (++ leveryone) females -1))
    (++ mbirths (length mfemales))
    (++ fbirths (length gfemales))))

(define (do-deaths yr)
  (let (dmales   (filter (fn (y) (> (age y) male-age-of-death)) males)
        dfemales (filter (fn (y) (> (age y) female-age-of-death)) females))
    (set 'males (difference males dmales))
    (set 'females (difference females dfemales))
    (set 'mdeaths (length dmales))
    (set 'fdeaths (length dfemales))
    (dolist (m dmales)
      (setf (everyone m 4) yr)
      (dolist (w (everyone m 5)) (setf (everyone w 5) nil))
    )
    (dolist (f dfemales)
      (setf (everyone f 4) yr)
      (let (h (everyone f 5))
        (when h
          (setf (everyone h 5)
              (clean (fn (x) (= x f)) (everyone h 5)))) ; free up husband to take another wife
      )
    )
  )
)

(define (do-marriages yr)
  (let (smales (filter (fn (x) (and (< (length (everyone x 5)) maxwives) (> (age x) male-age-of-majority))) males)
        sfemales (clean (fn (x) (or (everyone x 5) (< (age x) female-age-of-majority) (> (age x) female-age-of-widowhood))) females))
  (catch (dolist (m smales)
    (when (empty? sfemales) (throw))
    (let (f (pop sfemales))
      (push f (everyone m 5))
      (setf (everyone f 5) m)
      (++ fages (age f))
      (++ mages (age m))
      (++ nmarriages)
  )))))

(define (print-year yr)
  (println (format "Year %d. population %d (m: %d f: %d) %d births (m: %d f: %d) %d deaths (m: %d f: %d) %d marriages (m: %d f: %d)"
    yr
    (+ (length males) (length females)) (length males) (length females)
    (+ mbirths fbirths) mbirths fbirths
    (+ mdeaths fdeaths) 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)
(define (do-year yr)
      (when (and (< 0 (length males)) (< 0 (length females)))
      (do-births yr)
      (do-deaths yr)
      (do-marriages yr)
      (print-year yr)))

(for (yr male-age-of-majority 1000)
     (do-year yr))

(exit)