(in-package "MAXIMA")
(macsyma-module graphs)

(DEFMACRO MAKE-MLIST (&rest ARGS) `(LIST '(MLIST) . ,ARGS))
(DEFMACRO MAKE-MLIST-L (LLIST) `(CONS '(MLIST) ,LLIST))
(DEFMACRO MAKE-MLIST-SIMP (&rest ARGS) `(LIST '(MLIST SIMP) . ,ARGS))
(DEFMACRO MAKE-MLIST-SIMP-L (LLIST) `(CONS '(MLIST SIMP) ,LLIST))

(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 (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
(defmacro delay (hgraph prop value) 
  `(lambda () (prop-set ,hgraph ,prop ,value)))

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

;********************************************************************
; hypergraph, rename-vertex, rename-vertices, remove-edge
;********************************************************************

; 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)

; vertices is a list of vertex-structures here
(defun vertex-index (vertices vertex)
  (position vertex vertices :test 'equal :key 'vertex-name))

; when edge is a list of (names of) vertices, the function returns the
; corresponding list of vertex-indices
; vertices is a list of vertex-structures here
(defun vertex-indices (vertices edge)
  (mapcar #'(lambda (vtx) (vertex-index vertices vtx)) edge))

; vertices is a list of vertex-structures here; this would be the only place
; where vertex-names would be more appropriate
(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 medges-to-raw-lists (edges &optional fun)
  (if (mlistp edges)
      (mapcar #'(lambda(e) (if (mlistp (car (margs e)))
			       (cons (margs (car (margs e))) (cdr (margs e)))
			     (error "The first element of each element of the 
                                     argument to ~M should be a list"))) 
	      (margs edges))
    (if fun (error  "Argument to ~M should be a list" fun)
      (error "Argument to medges-to-raw-lists should be an mlist"))))

; vertices is a list of vertex-structures here
(defun medges-to-edge-patterns (vertices edges)
  (mapcar #'(lambda (medg) 
	      (cons (vertex-indices vertices (car medg))
		    (cdr medg)))
	  (medges-to-raw-lists edges)))

; vertices is a mlist, edges a mlist of mlists of type [mlist, weight] or of type [mlist]
; ie. [[[1,2],w1],[[2,3]]]
(defun $hypergraph (vertices edges &rest type)
  (if (mlistp vertices)
      (hypergraph (margs vertices) (medges-to-raw-lists edges) type)
    (error "First Argument to hypergraph should be a list")))

; to make defining a graph more straightforward, hypergraph should accept some switches:
; undirected
; unweighted
; ...

; however, should an undirected graph store all its edges in both directions?
; rather not. The question is whether to use sets, or to use lists and order
; them or to assume they come in ordered. I will assume that they come in
; ordered. In any case, they should be stored ordered!
(defun hypergraph (vertices edges &optional type)
; sollte Argumente prüfen
  (let ((hgraph (make-hypergraph)))
    (prop-set hgraph 'type type)
    (prop-set hgraph 'vertices 
	      (mapcar #'(lambda (v) (make-vertex :name v))
		      vertices))
    (prop-set hgraph 'edges
	      (let ((key -1))
		(mapcar #'(lambda (e)
			    (create-edge (G-vertices hgraph)
					 (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!
; ie. the vertex-names are not used anymore
; _____________________________________________________________

; 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)
; returns the edge-structure
(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 graph-type)
  (let ((new-a (copy-array value)))
    (dolist (edge edges (update-adjacency-matrix new-a edges graph-type))
      (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)))

(defun $remove_edges (hgraph edges)
  (remove-edges hgraph (medges-to-edge-patterns (G-vertices hgraph) edges)))

; 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
		       ('type value)
		       ('adjacency-matrix 
			(remove-edges-adjacency-matrix s-edges 
						       (force value)
						       (G-type hgraph)))
		       ('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))))
;; but delay is compiled, while the following is not!
;; 				 `(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-type (hgraph) (prop-get hgraph 'type))

;; this should probably give the vertex names, not the indices
(defun medge (hgraph edge)
  (make-mlist-simp 
   (make-mlist-simp-l (mapcar #'(lambda (v) 
				  (vertex-name (elt (G-vertices hgraph) v)))
			      (edge-vertices edge)))
   (edge-weight edge)))
(defun medges (hgraph edges)
  (make-mlist-simp-l (mapcar #'(lambda (e) (medge hgraph e)) edges)))

(defun $edges (hgraph) 
  (medges hgraph (G-edges hgraph)))

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

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

; doesn't care whether vertices are vertex-structures or not
; is this for components?
(defun mvertices (hgraph vertices)
  (make-mlist-simp-l 
   (mapcar #'(lambda (v) (if (vertex-p v) 
			     (vertex-name v) 
			   (vertex-name (elt (G-vertices hgraph) v))))
	   vertices)))

(defun $vertices (hgraph) (mvertices hgraph (G-vertices hgraph)))

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

(defun $adjacency_matrix (hgraph)
  (G-adjacency-matrix hgraph))

; probably not too good to use this here, continuation takes no argument, list
; contains the current permutation
(defun permute (list continuation)
   (if (cdr list)
       (loop for point on list
             do
             (rotatef (first list) (first point))
             (permute (rest list) continuation)
             (rotatef (first list) (first point)))
       (funcall continuation)))

; edges is a list of edge-structures here, value the adjacency matrix
(defun update-adjacency-matrix (value edges graph-type)
  (if (member '$undirected graph-type)
; for each permutation of the edge-vertices, the corresponding entry must be
; updated - there can only be entries where there is at least one edge...
      (let ((res value)
	    val edge)
	(dolist (edg edges res)
	  (psetq val nil 
; val needs to be initialised to nil at every iteration!
		 edge (edge-vertices edg))
	  (permute edge #'(lambda () (if val 
; the first call is the identity, this value is ok
					 (setf (apply 'aref res edge) val)
				       (setq val (apply 'aref res edge)))))))
    value))

; what to do with non-uniform hypergraphs? At least I should test...  in case
; the graph is undirected, I only keep the upper triangle up to date. (or,
; rather the part where vertex-indices are ordered) This should save some time,
; shouldn't it? NO it doesn't! if only one value is changed, I still would have
; to update all!
(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 (make-list (length (edge-vertices (car edges)))
					:initial-element dim) 
			     :initial-element 0)))
       (dolist (e edges (prop-set hgraph 'adjacency-matrix 
				  (update-adjacency-matrix 
				   (dolist (e edges res)
				     (rplacd (apply 'aref res (edge-vertices e)) nil)
				     (push '(mplus) (apply 'aref res (edge-vertices e))))
				   edges 
				   (prop-get hgraph 'type))))
	 (push (edge-weight e) (apply 'aref res (edge-vertices e))))))))

(defun $perfect_matchings (hgraph)
  (make-mlist-simp-l (mapcar #'(lambda (edges) (medges hgraph edges))
			     (G-perfect-matchings hgraph))))

(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 $components (hgraph)
  (make-mlist-simp-l (mapcar #'(lambda (vertices) (mvertices hgraph vertices))
			     (G-components 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))))))

; vertex is the name of a vertex here
(defun $degree (hgraph vertex)
  (G-degree hgraph (vertex-index (G-vertices hgraph) vertex)))

(defun compute-degree (hgraph vertex)
  (multiple-value-bind 
      (value presentP) (prop-get hgraph 'adjacency-matrix)
    (if presentP (let* ((m (force value))
			(dim (array-dimension m 0)))
		   (cons '(mplus) 
			 (coerce (make-array (list dim) 
					     :displaced-to m
					     :displaced-index-offset 
					     (* vertex dim))
				 'list)))
; otherwise use the list of edges
      (cons '(mplus) (mapcar 'edge-weight 
			     (remove-if #'(lambda (e) (not (member vertex e)))
					(G-edges hgraph)
					:key 'edge-vertices))))))

; this is really the total degree of the vertex
; vertex is the index of vertex here
(defun G-degree (hgraph vertex)
  (multiple-value-bind 
      (value presentP) (prop-get hgraph 'degree-sequence)
    (let ((degs (when presentP (force value))))
      (cond ((and degs (aref degs vertex)))
	    (T (unless degs
		 (setq degs (make-array (list (length (G-vertices hgraph))) 
					:initial-element nil)))
	       (prop-set hgraph 'degree-sequence degs)
	       (setf (aref degs vertex) 
		     (compute-degree hgraph vertex)))))))

(defun $degree_sequence (hgraph)
  (make-mlist-l (coerce (G-degree-sequence hgraph) 'list)))

;; a little overhead could be saved here...
(defun G-degree-sequence (hgraph)
  (dotimes (v (length (G-vertices hgraph)) 
	      (prop-get hgraph 'degree-sequence))
    (G-degree hgraph v)))


;********************************************************************
;  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))
   (list '$undirected)))


;********************************************************************
; 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 $three_aztec (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 (list '$undirected))))