; some helper routines, look for the word "hypergraph" to get started

(defun range (n)
  (reverse (let ((x nil)) (dotimes (i n x) (setq x (cons i x))))))

(defun row-major-aref (array index)
  (aref (make-array (array-total-size array)
		    :displaced-to array
		    :element-type (array-element-type array))
	index))

(defun setf-row-major-aref (array index value)
  (setf (aref (make-array (array-total-size array)
			  :displaced-to array
			  :element-type (array-element-type array))
	      index) value))

(defsetf row-major-aref setf-row-major-aref)

; returns a copy of the array (the elements are eq)
(defun copy-array (a)
  (let ((new-a (make-array (array-dimensions a))))
    (dotimes (index (array-total-size a) new-a)
	     (setf (row-major-aref new-a index) 
		   (row-major-aref a index)))))

; conses e to each of the lists
(defun cons-all (e lists)
  (mapcar #'(lambda (l) (cons e l)) lists))

; (pos expr test) gives a list of the positions at which objects
; satisfying test appear in expr.

; depth-first search

; wenn eine Liste expr den Test erfüllt, soll dann noch weitergesucht
; werden?
(defun pos (expr test)
  (cond ((funcall test expr) (list nil))
	((null expr) nil)
	((listp expr) (append (cons-all 1 (pos (car expr) test))
			      (mapcar #'(lambda (l)
					  (when l (rplaca l (1+ (car l)))))
				      (pos (cdr expr) test))))
	(T nil)))

; hier ebendiese Variante:
(defun pos-a (expr test)
  (cond ((null expr) nil)
	((listp expr) (append (if (funcall test expr) (list nil))
			      (cons-all 1 (pos-a (car expr) test))
			      (mapcar #'(lambda (l)
					  (when l (rplaca l (1+ (car l)))))
				      (pos-a-cdr (cdr expr) test))))
	(T (if (funcall test expr) (list nil)))))

(defun pos-a-cdr (expr test)
  (cond ((null expr) nil)
	((listp expr) (append (cons-all 1 (pos-a (car expr) test))
			      (mapcar #'(lambda (l)
					  (when l (rplaca l (1+ (car l)))))
				      (pos-a-cdr (cdr expr) test))))
	(T (error "pos-a-cdr takes only lists"))))

(defun substitute-in-level (new old lst level &rest args)
  (if (listp lst)
      (if (= 0 level) 
	  (apply #'substitute (append (list new old lst) args))
	(mapcar #'(lambda(sublst) 
		    (apply #'substitute-in-level 
			   (append (list new old sublst (1- level)) args)))
		lst))
    lst))	      
 
; eg: g:hypergraph(vertices,edges);
;     print(matchings(g));
;     print(matchings(remove-edge(g,edge)));
;     kill(g);
; the second line should store the matchings of g somewhere, so that
; little has to be done in the third line. The last line removes
; everything, except of the pointers kept by the labels.
;
; Question: should "remove-edge" be destructive or not?
; Answer: rather not...
;  
; Thus, "remove-edge" returns a fresh hypergraph, with all the
; properties of the original one preserved.
;
; Philosophy: functions that return a new hypergraph are
; non-destructive, and try to preserve as many properties as possible,
; functions that return a property of the hypergraph store the value
; on the plist of the hypergraph
;
; delayed evaluation:
; if a function modifies the hypergraph, it does not need to
; re-adjust all the properties immediately. Rather, it may return a
; closure.

(defun make-hypergraph ()
  (list '(%hypergraph) (make-hash-table)))

; sets the property and returns value
(defun prop-set (hgraph prop value)
  (setf (gethash prop (cadr hgraph)) value))

(defun prop-get (hgraph prop)
  (gethash prop (cadr hgraph)))

; calls fun (prop value) for every prop of hgraph and returns nothing
(defun map-props (fun hgraph)
  (maphash fun (cadr hgraph)))

(defun print-props (hgraph)
  (map-props #'(lambda (prop value) (print (list prop value)))
	     hgraph))		

; force evaluation
(defun force (value)
  (if (functionp value) (funcall value) value))

; delay evaluation
; for some reason this results in a larger closure. why?
(defmacro delay (hgraph prop value) 
  `(lambda () (prop-set ,hgraph ,prop ,value)))

;********************************************************************
; OPERATIONS
;********************************************************************

;********************************************************************
; hypergraph, remove-edges
;********************************************************************

; it would be very useful if the vertices were numbered 0 to n. Maybe
; a vertex should be of the form

; (properties), 

; an edge of the form

; (indices properties), 

; where indices is a list of vertex-indices. Properties should be
; user-extensible, but I'm not sure how. 

; Doing so, rename-vertex becomes unnecessary...  Not quite: removing
; a vertex makes a renumbering necessary, at least for those functions
; that depend on the fact that all vertices from 0 to n exist,
; perfect-matchings for example. 

; However, hypergraph should take the arguments just the way it does
; now

(defstruct vertex (name))
(defstruct edge vertices (weight 1) key)

; when edge is a list of (names of) vertices, the function returns the
; corresponding list of vertex-indices
(defun vertex-indices (vertices edge)
  (mapcar #'(lambda (vtx) (position vtx vertices :test #'equal)) edge))

; vertices should be an mlist and edges a mlist of mlists
; should I enable also other possibilities?
; The question really is: do I want to have a conversion at in and
; output, or do I want to work with maxima objects throughout?

(defun create-edge (vertices edge-vertices edge-weight edge-key)
  (make-edge :vertices (vertex-indices vertices edge-vertices)
	     :weight edge-weight
	     :key edge-key))

(defun hypergraph (vertices edges)
; sollte Argumente prüfen
  (let ((hgraph (make-hypergraph)))
    (prop-set hgraph 'vertices 
	      (mapcar #'(lambda (v) (make-vertex :name v))
		      vertices))
    (prop-set hgraph 'edges
	      (let ((key -1))
		(mapcar #'(lambda (e)
			    (create-edge vertices 
					 (car e) 
					 (if (= (length e) 1) 1
					   (cadr e))
					 (incf key)))
			edges)))
    hgraph))

; from here on, the key to an edge is a list of vertex-indices!
; _____________________________________________________________

; what exactly is the key to an edge? Do I want to allow for multiple
; weighted edges? The main difference is the underlying graph, I
; think... Also the result of perfect-matchings has to be changed.  If
; yes, I have to be very careful when identifying an edge. Bigger: how
; to specify a subset of edges then? I want to be able to specify an
; edge without knowing its weight. A sophisticated way out would be
; the following: If there are multiple edges with differing weights,
; each connecting the given set of vertices, the one with weight 1
; wins, if present, otherwise error. This implies that I have to check
; whether multiple edges are present.
;
; Maybe the easiest way would be to keep a unique index with each
; edge. 
; I could then have the following functions:
; UnderlyingGraph making the weight of all edges equal to 1
; UnderlyingSimpleGraph replacing multiple weighted edges by a single
; edge with weight 1
; UnderlyingWeightedGraph adding up the weights of multiple edges,
; just as the adjacency matrix does
; plus all the undirected versions...
;
; Is this worth it? For example: computing weight of the matchings of
; a graph would take much longer if done on the original graph than on
; the underlying weighted graph
; Who knows...
; 
; is an edge with weight 0 removed? (this is important for the concept
; of components...) I might also have a global switch controlling this.
; Sometimes it is necessary that a loop counts twice, isn't it?

; for edge-structures only
(defun member-edge (edge edges)
  (member (edge-key edge) edges :key 'edge-key))

; also if edge is not an edge-structure, edge being of the form
; ((vertex-indices)) or ((vertex-indices) weight)
(defun find-edge (edge edges)
  (if (edge-p edge) (car (member-edge edge edges))
    (let ((res (remove-if (if (= (length edge) 1)
			      #'(lambda (e)
				  (not (equal (car edge)
					      (edge-vertices e))))
			    #'(lambda (e)
				(not (equal edge
					    (list (edge-vertices e)
						  (edge-weight e))))))
			  edges)))
      (if (= (length res) 1) (car res)
	(find '1 res :key 'edge-weight)))))

(defun find-edges (edges all-edges)
  (mapcar #'(lambda (e) (cond ((find-edge e all-edges))
			      (T (error "Not an edge"))))
	  edges))			       

; **************************** remove-edges ***************************

(defun remove-edges-adjacency-matrix (edges value)
  (let ((new-a (copy-array value)))
    (dolist (edge edges new-a)
	    (decf (apply #'aref (cons new-a (edge-vertices edge)))
		  (edge-weight edge)))))

(defun remove-edges-edges (edges value)
  (set-difference value edges :key #'edge-key))

(defun remove-edges-perfect-matchings (edges value)
  (do ((edg edges (cdr edg))
       (res value (remove (car edg) res :count 1 :test #'member-edge)))
      ((or (null edg)
	   (null res)) res)))

(defun remove-edges-components (edges value new-edges)
  (do ((affected-vertices (remove-duplicates 
			   (mapcar #'(lambda (e) (car (edge-vertices e)))
				   edges))
			  (if (cdr cmp)
			      (set-difference affected-vertices (car cmp))))
       (cmp value 
	    (cdr cmp))
       (res nil
	    (append (if (some #'(lambda (v) (member v (car cmp)))
			      affected-vertices)
			(components (car cmp) new-edges)
		      (list (car cmp))) 
		    res)))
      ((null cmp) res)))

; maybe it would be better to delay the computation of all
; properties... Mind though, that some properties need more
; information than the others. It is probably not a good idea to pass
; the whole old hgraph
(defun remove-edges (hgraph edges)
  (let ((s-edges (find-edges edges (G-edges hgraph)))
	(newhgraph (make-hypergraph)))
    (map-props 
     #'(lambda (prop value)
	 (block nil
		(prop-set newhgraph prop
			  (case prop
				('adjacency-matrix 
				 (remove-edges-adjacency-matrix 
				  s-edges 
				  (force value)))
				('vertices value)
				('edges
				 (remove-edges-edges 
				  s-edges
				  (force value)))
				('components
; for some reason this results in a larger closure. why?
; 				   (delay newhgraph prop 
; 					  (remove-edges-components 
; 					   s-edges
; 					   (force value)
; 					   (G-edges newhgraph))))
				 `(lambda () 
				    (prop-set ',newhgraph ',prop
					      (remove-edges-components 
					       ',s-edges
					       (force ',value)
					       (G-edges	',newhgraph)))))
				('perfect-matchings
				 `(lambda ()
				    (prop-set ',newhgraph ',prop
					      (remove-edges-perfect-matchings
					       ',s-edges 
					       (force ',value)))))
				(T (return))))))
     hgraph)
    newhgraph)))

;********************************************************************
; PROPERTIES
;********************************************************************

;********************************************************************
;  edges, vertices, adjacency-matrix, perfect-matchings, components
;********************************************************************

(defun G-edges (hgraph) (prop-get hgraph 'edges))

(defun G-edge-vertices (hgraph) (mapcar #'edge-vertices (G-edges hgraph)))

(defun G-vertices (hgraph) (prop-get hgraph 'vertices))

; what to do with non-uniform hypergraphs? At least I should test...
(defun G-adjacency-matrix (hgraph)
  (multiple-value-bind 
   (value presentP) (prop-get hgraph 'adjacency-matrix)
   (if presentP (force value)
     (let* ((edges (G-edges hgraph))
	    (dim (length (G-vertices hgraph)))
	    (res (make-array (list dim dim) :initial-element 0)))
       (dolist (e edges (prop-set hgraph 'adjacency-matrix res))
	       (incf (apply #'aref res (edge-vertices e)) (edge-weight e)))))))

(defun G-perfect-matchings (hgraph)
  (multiple-value-bind 
   (value presentP) (prop-get hgraph 'perfect-matchings)
   (if presentP (force value)
     (prop-set hgraph 'perfect-matchings 
	       (perfect-matchings (range (length (G-vertices hgraph)))
				  (G-edges hgraph))))))

(defun G-components (hgraph)
  (multiple-value-bind 
   (value presentP) (prop-get hgraph 'components)
   (if presentP (force value)
     (prop-set hgraph 'components
	       (components (range (length (G-vertices hgraph)))
			   (G-edges hgraph))))))

;********************************************************************
;  helper routines:
;
; These routines should not depend on the fact that all the vertices
; from 0 to n are present
;********************************************************************

; the edges which are not incident with the edge e, ie., the vertex
; induced subgraph G\e
(defun non-incident (e edges)
  (remove-if #'(lambda (f)  
		 (some #'(lambda (v) (member v (edge-vertices f)
					     :test #'equal)) 
		       (edge-vertices e)))
	     edges))

(defun find-isolated-vertex (vertices edges)
  (set-difference vertices (reduce #'append 
				   (mapcar 'edge-vertices edges))))

(defun neighbourhood+ (vertex edges)
  (remove-duplicates
   (reduce #'append 
	   (mapcar #'cdr (remove-if #'(lambda (e) 
					(not (equal vertex 
						    (car e))))
				    (mapcar #'edge-vertices edges))))))

(defun neighbourhood (vertex edges)
  (remove-duplicates
   (reduce #'append 
	   (remove-if #'(lambda (e) (not (member vertex e)))
		      (mapcar #'edge-vertices edges)))))

; Finds all the vertices in the same (weakly connected) component as
; v, when the graph is given by edges. Visited is a list of vertices
; which shall not be visited. Could be improved by keeping also an
; array containing nil for the unvisited, T for the visited vertices.
(defun component (edges v &optional visited)
  (let ((res (cons v visited)))
    (dolist (vtx (neighbourhood v edges) res)
	    (if (not (member vtx res))
		(setq res (component edges vtx res))))))

(defun components (vertices edges)
  (do* ((cmp nil 
	     (component edges (car vtx) nil))
	(vtx vertices
	     (set-difference vtx cmp))
	(res nil (cons cmp res)))
       ((null vtx) res)))

; a matching is a collection of edges that are mutually non-incident
; take an edge e, construct the matchings of the vertex induced
; subgraph G\e, and adjoin e to each of those. Then do the same for
; the edge induced graph G\e, and so on. If G has no matching, don't
; try for the edge induced graph G\e
(defun perfect-matchings (vertices edges)
  (cond ((null vertices) '(nil))
	((find-isolated-vertex vertices edges) nil)
	(T (mapcon #'(lambda (edgesrest) 
		       (cons-all (car edgesrest) 
				 (perfect-matchings 
				  (set-difference vertices (edge-vertices 
							    (car edgesrest)))
				  (non-incident (car edgesrest)
						(cdr edgesrest)))))
		   edges))))

;********************************************************************
; GRAPHS
;********************************************************************

;********************************************************************
;  COMPLETE GRAPH *** complete-graph (n)
;********************************************************************

(defun complete-graph (n)
  (hypergraph
   (range (1+ n))
  
   (do ((i n (1- i))
	(result nil (append (do ((j n (1- j))
				 (edgelist nil (cons (list (list i j))
						     edgelist)))
				((= i j) edgelist))
			    result)))
       ((= -1 i) result))))


;********************************************************************
; AZTEC DIAMOND 3D *** 3aztec (n) *** n ungerade
;********************************************************************
;
; Kanten lassen sich durch Tripel (x,y,z) mit 
; |x|+|y|+|z| = 2k+1 <= 2m+1 = n, x oder y oder z gerade
; indizieren.  Zwei Kanten (a,b,c) und (x,y,z) sind inzident, wenn
; |a-x|+|b-y|+|c-z|=2 und |k-l|=2 nur für gerade k.
;
; das hilft mir leider nichts...
;
; die Kante (x,y,z) ist zu folgenden Knoten (mit i,j aus {-1, +1}) inzident:
; (x+i,y+j,z) wenn z ungerade
; (x+i,y,z+j) wenn y ungerade
; (x,y+i,z+j) wenn x ungerade
;
(defun edge-aux (e)
  (cond ((oddp (car e))
	 (list (list (car e) (1+ (cadr e)) (1+ (caddr e)))
	       (list (car e) (1+ (cadr e)) (1- (caddr e)))
	       (list (car e) (1- (cadr e)) (1+ (caddr e)))
	       (list (car e) (1- (cadr e)) (1- (caddr e)))))
	((oddp (cadr e))
	 (list (list (1+ (car e)) (cadr e) (1+ (caddr e)))
	       (list (1+ (car e)) (cadr e) (1- (caddr e)))
	       (list (1- (car e)) (cadr e) (1+ (caddr e)))
	       (list (1- (car e)) (cadr e) (1- (caddr e)))))
	(T ; (oddp (caddr e))
	 (list (list (1+ (car e)) (1+ (cadr e)) (caddr e))
	       (list (1- (car e)) (1+ (cadr e)) (caddr e))
	       (list (1+ (car e)) (1- (cadr e)) (caddr e))
	       (list (1- (car e)) (1- (cadr e)) (caddr e))))))

; noch einfacher:
;
; die Knoten sind durch Tripel (2x+1,2y+1,2z+1) mit 
; |2x+1|+|2y+1|+|2z+1| <= 2m+3 = n+2
; indiziert.
;
; erzeuge die Kanten und verwandle sie in Knoten-Quadrupel?
;
; der 3aztec hat 4n^3+12n^2+14n+6 Kanten

(defun 3aztec (n)
  (let ((vertices (do ((x (- n) 
			  (+ x 2))
		       (r nil
			  (append (do ((y (+ (abs x) (- n) -1)
					  (+ y 2))
				       (r nil
					  (append (do ((z (+ (abs x) (abs y) (- n) -2)
							  (+ z 2))
						       (r nil (cons (list x y z)
								    r)))
						      ((> z (- n -2 (abs x) (abs y))) r))
						  r)))
				      ((> y (- n -1 (abs x))) r))
				  r)))
		      ((> x n) r)))
        (edges (mapcar #'(lambda (e) (list (edge-aux e)))
		       (do ((x (- n)
			       (1+ x))
			    (r nil
			       (append (do ((y (+ (abs x) (- n))
					; x ungerade -> y gerade
					; x gerade -> y beliebig
					       (if (oddp x) (+ y 2) (1+ y)))
					    (r nil
					       (append (do ((z (+ (abs x) 
								  (abs y) 
								  (- n))
							       (+ z 2))
							    (r nil 
							       (cons 
								(list x y z)
								r)))
							   ((> z (- n 
								    (abs x) 
								    (abs y))) 
							    r))
						       r)))
					   ((> y (- n (abs x))) r))
				       r)))
			   ((> x n) r)))))
	(hypergraph vertices edges)))