; This file is in the public domain and is ungoverned by the Cypherpunks' Anti-License.
; Do with it as you will.

; The utgraph concept uses the n(n-1)/2 cells in an upper triangular (ut)
; matrix to represent the adjacencies between the n vertices of a graph.
; An embedded (in R^d) graph is this "adjacency ut" cons'ed to
; a list of n d'tuples, representing the coordinates of the n vertices
; in d dimensions.

; utgraph 2.0, an attempt to improve the efficiency of the utgraph concept.

; for utgraph 1.0 see http://geocities.com/n8chz/utgraph.txt and http://geocities.com/n8chz/tweak.txt

; The goal is to find a "low energy" embedding for a given adjacency-ut,
; where low energy is taken to mean a large negative correlation between distance and adjacency (nominally 0 or 1).

; Version 1.0 randomly swapped randomly selected coordinate pairs of randomly selected vertex pairs.
; Version 2.0 attempts a more elegant approach, based roughly on genetic algorithms.
; The assumption is made that embeddings that are "unusually low energy" are "nearly similar"
; (or else "nearly inverted" in the Euclidean sense.)
; In hopes of finding other low energy configurations, we construct an "intermediate embedding"
; which is the "halfway" state if we were to "morph" one low energy embedding
; into another, using a method that varies angle and radius
; (literally, the "r" and "theta" of polar coordinates in the two dimensional case)
; at a constant rate.

; The unit (hyper)sphere is used as a boundary enclosing all vertices, to avoid problems of scale creep.

; Convergence seems very rapid compared to 1.0, but optimality may not be assured.
; Testing continues.  Hopefully a useful graph embedding tool will result.

; Coded in Gnu Common Lisp 2.2.2

; The following few routines set up a basic library of vector arithmetic.

(defun square (number) (* number number))

(defun vector-length (vector)
 (sqrt (apply #'+ (mapcar #'square vector))))

(defun unit-vector (vector)
 (let ((len (vector-length vector)))
  (if (plusp len) (mapcar #'(lambda (x) (/ x len)) vector))))

(defun vector-sum (&rest vectors)
 (apply #'mapcar (cons #'+ vectors)))

(defun average (&rest numbers)
 (if (cdr numbers) (/ (apply #'+ numbers) (length numbers))))

(defun set-vector-length (vector lenth)
 (mapcar #'(lambda (x) (* lenth x)) (unit-vector vector)))

(defun polar-midpoint (&rest vectors)
 (set-vector-length
  (apply
   #'vector-sum
   (mapcar #'unit-vector vectors))
  (apply
   #'average
   (mapcar #'vector-length vectors))))

(defun random-subbox-vector (&optional (ndims 3)) ; uniformly distributed random vector inside [-1,1]^ndims
 (let (result)
  (dotimes (i ndims result)
   (setf result (cons (- (random 2.0) 1.0) result)))))

(defun random-subunit-vector (&optional (ndims 3)) ; uniformly distributed random vector inside unit (hyper)sphere
 (do
  ((result (random-subbox-vector ndims) (random-subbox-vector ndims)))
  ((<= (vector-length result) 1) result)))

; generate a random list of {1,0} of length n(n-1)/2
; this is essentially the upper triangular of an adjacency matrix
; whose underlying graph is non-directed and is a unigraph (univalent graph)

; no provisions have been made for disconnectedness,
; so more than one graph is a likely outcome

(defun random-adjacency-ut (n)
 (let ((result nil))
  (dotimes (i (/ (* n (1- n)) 2) result)
   (setf result (cons (random 2) result)))))

; Generate a random embedding of vertices in d dimensional space.
; This is modified slightly from 1.0.
; This is because 2.0 embeds in the interiors
; of spheres rather than cubes.

(defun random-embed (n &optional (d 3))
 (let ((result nil))
  (dotimes (i n result)
   (setf result (cons (random-subunit-vector d) result)))))
                                                  
; The next few functions are carried over from utgraph 1.0

; The #'random-graph function in utgraph.lsp
; generates a random embedded graph in ut
; representation.

; The car of the resulting list is of length n(n-1)/2
; and gives the adjacencies between the n vertices

; The cdr of the result is a list of the d dimensional
; coordinates of the n vertices.
; Each of these is given as a list of d numbers.

(defun random-graph (n &optional (d 3))
 (cons (random-adjacency-ut n) (random-embed n d)))

; The 'utpair' is a list
; of the C(n,2) pairings of
; items in 'nodelist'

; list C(length(nodelist),2) distinct pairings of nodelist
; this results in a utmatrix not entirely unlike the
; familiar "mileage" (i.e. distance) tables
; on an "obsolete" technology called printed maps

; The order in which the pairings are presented is as follows:

; n-1,n
; n-2,n
; n-2,n-1
; n-3,n
; ...
; 1,n
; ...
; 1,3
; 1,2

(defun utpair (nodelist)
 (let
  ((result nil))
  (do
   ((i nodelist (cdr i)))
   ((null (cdr i)) result)
   (do
    ((j (cdr i) (cdr j)))
    ((null j) result)
    (setf result (cons (cons (car i) (car j)) result))))))

; give the square of the distance between two vectors

(defun dist (v1 v2)
 (let
  ((s 0))
  (do
   ((w1 v1 (cdr w1))
    (w2 v2 (cdr w2)))
   ((null w1) (sqrt s))
   (setf s (+ (expt (- (car w1) (car w2)) 2) s)))))

; correlation function for paired data
; it assumes the data come as a list of pairs (i.e. cons cells)

(defun corr (pairs)
 (let
  (x
   y
   (sx 0)
   (sy 0)
   (sxy 0)
   (sx2 0)
   (sy2 0))
  (do
   ((rest pairs (cdr rest))
    (n 0 (1+ n)))
   ((null rest)
    (/
     (- (* n sxy) (* sx sy))
      (sqrt
       (*
        (- (* n sx2) (* sx sx))
        (- (* n sy2) (* sy sy))))))
   (setf
    x (caar rest)
    y (cdar rest)
    sx (+ sx x)
    sy (+ sy y)
    sxy (+ sxy (* x y))
    sx2 (+ sx2 (* x x))
    sy2 (+ sy2 (* y y))))))

; calculate correlation between adjacency and distance.
; for our purposes, strong negative correlations are most desired.

(defun correlate-graph (ut embed)
 (corr
  (mapcar
   #'cons
   ut
   (mapcar
    #'(lambda (x) (dist (car x) (cdr x)))
    (utpair embed)))))

; If graph1 and graph2 are "low energy,"
; the intermediate graph might be even lower energy!
; Or so we hope.

(defun intermediate-embed (embed1 embed2)
 (mapcar #'polar-midpoint embed1 embed2))

(defun flip (embed) ; flip an embedding along x axis, creating a mirror image
 (mapcar #'(lambda (x) (cons (- (car x)) (cdr x))) embed))

; test-method allows one to watch the correlations converge.
; Dramatic improvements seem equally likely attributable
; to intermediate embedding or dumb luck, which is not encouraging.

(defun test-method (n &optional (d 3))
 (do
  (p q (adj (random-adjacency-ut n)) (pool (list (random-embed n d))) (i 1 (1+ i)) (old 0) (new 0))
  (nil) ; Yup it's an endless looop.
  (let ((surplus (- (length pool) 3)))
   (setf pool (sort pool #'<= :key #'(lambda (x) (correlate-graph adj x))))
   (if (equal (car pool) p) (print "intermediate embedding utilized!"))
   (if (equal (car pool) q) (print "intermediate embedding to mirror image utilized!"))
   (setf old new)
   (setf new (correlate-graph adj (car pool)))
   (if (/= old new) (progn (print (cons i new))))
   (setf pool
    (cond
     ((minusp surplus) (setf p nil q nil) (cons (random-embed n d) pool))
     ((zerop surplus)
      (cons
       (setf p (intermediate-embed (car pool) (cadr pool)))
       (cons
        (setf q (intermediate-embed (car pool) (flip (cadr pool))))
        pool)))
     (t
      (setf p nil q nil)
      (nthcdr (1+ surplus) (reverse pool))))))))

(defun old-test-method (n &optional (d 3))
 (do
  (p q (adj (random-adjacency-ut n)) (pool (list (random-embed n d))) (i 1 (1+ i)) (old 0) (new 0))
  (nil) ; Yup it's an endless looop.
  (let ((surplus (- (length pool) 3)))
   (setf pool (sort pool #'<= :key #'(lambda (x) (correlate-graph adj x))))
   (if (equal (car pool) p) (print "intermediate embedding utilized!"))
   (if (equal (car pool) q) (print "intermediate embedding to mirror image utilized!"))
   (setf old new)
   (setf new (correlate-graph adj (car pool)))
   (if (/= old new) (progn (print (cons i new))))
   (setf pool
    (cond
     ((minusp surplus) (setf p nil q nil) (cons (random-embed n d) pool))
     ((zerop surplus)
      (cons
       (setf p (intermediate-embed (car pool) (cadr pool)))
       (cons
        (setf q (intermediate-embed (car pool) (flip (cadr pool))))
        pool)))
     (t
      (setf p nil q nil)
      (nthcdr (1+ surplus) (reverse pool))))))))

; A recursive test of the method...
; This seems more fruitful, though painfully slow.

(defun ranq (ut pool) (sort pool #'<= :key #'(lambda (x) (correlate-graph  ut x)))) ; rank pool by correlation

(defun threesome (level ut d &optional (n (floor (/ (1+ (sqrt (1+ (* (length ut) 8)))) 2))))
 (if
  (zerop level)
  (list
   (random-embed n d)
   (random-embed n d)
   (random-embed n d))
  (list
   (onesome (1- level) ut d n)
   (onesome (1- level) ut d n)
   (onesome (1- level) ut d n))))

(defun onesome (level ut d &optional (n (floor (/ (1+ (sqrt (1+ (* (length ut) 8)))) 2))))
 (let ((3some (ranq ut (threesome level ut d n))))
  (car (ranq ut
   (cons (intermediate-embed (car 3some) (cadr 3some))
    (cons (intermediate-embed (flip (car 3some)) (cadr 3some))
     3some))))))

; It is recommended to run test with level <= 6.

(defun test (level n &optional (d 3))
 (print (multiple-value-list (get-decoded-time)))
 (let (result ut)
  (prog1
   (setf result
    (onesome level (setf ut (random-adjacency-ut n)) d n))
   (print (correlate-graph ut result))
   (print (multiple-value-list (get-decoded-time))))))
  


; Look up "universal choosing"
; Look up "workers' independent news service"