Plots related to the Collatz function

collatz.rkt

#lang racket

;; http://www2.bakersfieldcollege.edu/resperic/ca/cellular_automata_dynamics.pdf

(require graph racket/system 2htdp/image)

(struct memf [f h domain])

(define (new-memf f [dom #f])
  (memf f (make-hash) dom))

(define (mem-size f)
  (length (hash-keys (memf-h f))))


(define (saturate f)
  (define rg
    (match (memf-domain f)
      [(? number? n) (range n)]
      [(list a ...) a]
      [(? pair? p) (range (car p) (cdr p))]
      [e (error "Invalid domain" e)]))
  (for ([k rg])(value f k)))
    

(define (orbit f k)
  (stream-cons k (orbit f (value f k))))

(define (orbit/list f k)
  (define (has-repetitions? l)
    (< (set-count (list->set l)) (length l)))
  (let loop ([o '()] [s (orbit f k)])
    (if (has-repetitions? o)
        (reverse o)
        (loop (cons (stream-first s) o) (stream-rest s)))))
               

(define (value f k)
  (let* ([h (memf-h f)]
         [upd (λ ()
                (let ([v ((memf-f f) k)])
                  (hash-set! h k v)
                  v))])
    (hash-ref h k upd)))
  


(define (iterate f g . ns)
  (define (it-once n)
    (let ([ni (f n)])
      (unless
          (hash-has-key? g n)
        (hash-set! g n ni)
        (iterate f g ni))))
  (for ([n ns]) (it-once n)))

(define (orbit1 g n)
  (let iter ([n n] [l '()])
    (if (member n l)
        (reverse l)
        (iter (hash-ref g n) (cons n l)))))


(define (show-graph f [tool "dot"])
  (define (new-temp ext)
    (make-temporary-file (string-append "~a." ext) #f))
  (saturate f)
  (let* ([dot-cmd (format "/usr/local/bin/~a" tool)]
         [dot-path (new-temp "dot")]
         [png-path (new-temp "png")]
         [g (memf-h f)]
         [coll-graph
          (unweighted-graph/directed
           (for/list ([(k v) g]) (list k v)))]
         [out
          (open-output-file dot-path #:exists 'replace)])
    (graphviz coll-graph #:output out)
    (close-output-port out)
    (with-output-to-file
        png-path
      (λ () (system (format "~a ~a -Nfontsize=10 -Nfixedsize=true -Nshape=circle -Tpng" dot-cmd dot-path)))
      #:exists 'truncate)
    (let ([bm (bitmap/file png-path)])
      (delete-file png-path)
      (delete-file dot-path)
      bm)))

(define (show-orbit g n)
  (local-require plot)
  (let ([o (orbit g n)])
    (plot
     (lines (map vector (range (length o)) o)))))

(define (show-times f n)
  (local-require plot)
  (let ([lg (λ (x) (log x 10))])
    (plot
     (points
      (for/list ([k (range 1 n)])
        (vector k (lg (length (orbit/list f k)))))
      #:sym 'dot))))



(define (pdf f t #:samples [samples 1000] xs)
  (local-require racket/random)
  (define fq (make-hash))
  (for ([n (in-range samples)])
    (let* ([x0 (random-ref xs)]
           [x (last (take (orbit f x0) t))])
      (hash-update! fq x add1 1)))
  (for/hash ([(x f) fq]) (values x (/ (exact->inexact f) samples))))
  



(define (take s k)
  (for/list ([x s] [k (in-range k)]) x))


(define (entropy f t)
  (define pd (pdf f t (memf-domain f)))
  (define (plogp x)
    (let ([p (hash-ref pd x 0)])
      (if (zero? p) 0 (* p (log p 2)))))
  (- (for/sum ([k (memf-domain f)]) (plogp k))))
      

(define-syntax define/mem
  (syntax-rules ()
    [(_ (id x rg) body)
      (define id
       (new-memf
        (λ (x) body) rg))]
    [(_ (id x) body)
      (define id
       (new-memf
        (λ (x) body)))]
    [(_ id proc)
     (define id (new-memf proc))]))
    
        

(define/mem (f17 x (range 17))
  (modulo (add1 (sqr x)) 17))

(define/mem (f37 x (range 37))
  (modulo (add1 (sqr x)) 37))

(define/mem (f221 x (range 221))
  (modulo (add1 (sqr x)) 221))


(define/mem (collatz n)
  (if
   (zero? (modulo n 2))
   (/ n 2)
   (+ (* 3 n) 1)))

(define/mem (g28 x (range 28))
  (modulo (+ (* x x x) (* 2 x x) (* 3 x) 2) 28))

Generalise to any function.
Better display of graphs