;; Copyright (C) 2002 Martin Rubey
;; LaBRI, Universite' Bordeaux 1
;;
;; www.mat.univie.ac.at/~rubey
;; rubey@labri.fr
 
;; This program is free software; you can redistribute it and/or
;; modify it under the terms of the GNU General Public License
;; as published by the Free Software Foundation; either version 2
;; of the License, or (at your option) any later version.
 
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
 
; TODO
; FEATURES: NEPS, Disjoint Union, Product, Kneser Graphs, underlying Graphs,
; from_adjacency_matrix
;
; IDEA: the rep of a hypergraph could be
; ((%hypergraph good) props weights) the weights are in sync
; ((%hypergraph) props weights) which means that I don't know.
;
; props is a gensym with a property graph-props, where the properties are stored.
; It might make sense to have also graph-defs, where vertices, edges, dim and
; undirected are stored, i.e.  those properties which have to be present.
;
; CLEANING: better definition of graph_composition

(in-package "MAXIMA")
(macsyma-module graphs)
(displa-def %hypergraph display-graph)

; (append result list2) is a permutation of list1, provided that list2
; is a sublist of list1, modulo permutation.  Elements in list2 that
; do not appear in list1 are ignored.  Hence, disregqrding calls to
; test and key, set-difference (l1 l2) = list-difference (l1 l2) if
; all elements have lower multiplicity in l1 than in l2
(defun list-difference (list1 list2 &key key test)
  (let ((res list1) (tst (if key key #'identity)))
    (dolist (e list2 res)
	    (setq res 
		  (remove (funcall tst e) res :test test :key key
			  :count 1)))))

(defun nlist-difference (list1 list2 &key key test)
  (let ((res list1) (tst (if key key #'identity)))
    (dolist (e list2 res)
	    (setq res 
		  (delete (funcall tst e) res :test test :key key
			  :count 1)))))

(defun single (lst)
  (and (consp lst) (not (cdr lst))))

(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 new-a) new-a)
	     (setf (row-major-aref new-a index) 
		   (row-major-aref a index)))))

(defun array-difference (a1 a2)
  (let ((new-a (make-array (array-dimensions a1))))
        (dotimes (index (array-total-size new-a) new-a)
	     (setf (row-major-aref new-a index) 
		   (meval (list '(mplus) (row-major-aref a1 index)
				(list '(mtimes) -1 
				      (row-major-aref a2 index))))))))

(defun array-add-entry (am ind add)
  (let ((val (apply #'aref am ind)))
    (setf (apply #'aref am ind) (cond ((eql 0 val) add)
;				     ((or (eql 0 val) (null val)) add)
; this would treat nil as 0
				     ((and (numberp val) (numberp add))
				      (+ val add))
				     ((and (listp val) 
					   (equal (car val) '(mplus)))
				      (append val (list add)))
				     (T (list '(mplus) val add))))))

(defun array-subtract-entry (am ind add)
  (let ((val (apply #'aref am ind)))
    (setf (apply #'aref am ind) (cond ((and (numberp val) (numberp add))
				      (- val add))
				     ((and (listp val) 
					   (equal (car val) '(mplus)))
				      (append val (list '(mtimes) -1 add)))
;				     ((or (eql 0 val) (null val))
; this would treat nil as 0
				     ((eql 0 val) (list '(mtimes) -1 add))
				     (T (list '(mplus) val 
					      (list '(mtimes) -1 add)))))))

; (cartesian (list '((a b) (c)) '(((a) (b)) ((c) (d)))))
; produces an error in gcl (bug?!)
(defun cartesian (sets)
 "Returns the Cartesian product of a list of sets
  Example:
       (CARTESIAN4 '((A B) (C D) (E F)))
         ;=> ((B C F) (B C E) (B D F) (B D E) (A C F) (A C E) (A D F) (A D E))
  Non-recursive version.
  "
 (let ((cartesian (list nil)))
   (dolist (set (reverse sets))         ; reverse is only for esthetics
     (let ((new-cartesian nil))
       (dolist (elmt set)
         (dolist (cart cartesian)
           (push (cons elmt cart) new-cartesian)))
       (setf cartesian new-cartesian)))
   cartesian))

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

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

;; ; delete the element specified by index, the first element being 0
;; (defun delete-index (index list)
;;   (if (zerop index) (cdr list)
;;     (progn
;;       (rplacd (nthcdr (1- index) list) (nthcdr (1+ index) list))
;;       list)))

;; ; deletes the (car indices), then (cadr indices) ...
;; (defun delete-indices (indices list)
;;   (do ((ind indices (cdr ind)) 
;;        (res list    (delete-index (car ind) res)))
;;       ((null ind) res)))

; 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 hypergraphp (hgraph)
  (and (listp hgraph) (listp (car hgraph)) (eq (caar hgraph) '%hypergraph)))

(defun display-graph (val res)
  (let ((tmp (dimension-match ($vertices val) 
			      (push-string "hypergraph(" res))))
    (setq tmp (dimension-match ($edges val) (push-string "," tmp)))
    (push-string ")" tmp)))

(defun make-hypergraph ()
  (let ((props (gensym)))
    (setf (get props 'graph-props) nil)
    (list (list '%hypergraph) props)))

; sets the property and returns value
(defun prop-set (hgraph prop value)
  (let ((res (assoc prop (get (cadr hgraph) 'graph-props) 
 		    :test #'eq)))
    (if res (setf (cdr res) value)
      (progn (push (cons prop value) (get (cadr hgraph) 'graph-props))
 	     value))))

; adds and sets the property and returns value
(defun prop-add (hgraph prop value)
  (push (cons prop value) (get (cadr hgraph) 'graph-props))
  value)

; returns (value present), value defaulting to nil, present is nil if
; the property is not present
(defun prop-get (hgraph prop)
  (let* ((graph-props (get (cadr hgraph) 'graph-props))
	 (res (assoc prop graph-props :test #'eq)))
    (values (cdr res) res)))

; calls fun (prop value) for every prop of hgraph and returns nothing
(defun map-props (fun hgraph)
  (let ((graph-props (get (cadr hgraph) 'graph-props)))
    (mapc #'(lambda (p) (funcall fun (car p) (cdr p))) graph-props)))

(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
; I have to use prop-set here! prop will be present, containing the delayed
; value.
(defmacro delay (hgraph prop value) 
  `(lambda () (prop-set ,hgraph ,prop ,value)))

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

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

; index is temporary. It should usually be a bijection vertices <-> [n].
(defstruct vertex name index)
; the key is used to determine edges in subgraphs
; if H is in some way a natural subgraph of G, then the keys of edges in
; H are the same as in G
(defstruct edge vertices (weight 1) (key 0 :type fixnum))
; vertices are vertex-structures now!!!

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; it would be good if the name of a function would indicate the type
; of arguments:

; $xxx take maxima arguments: vertex-names and m-edge-patterns
; m-xxx called by $xxx functions take vertex-names and m-edge-patterns
; or take structures and return maxima-stuff

; there are functions that take vertex indices (vertices)
; there are functions that take vertex names   (n-vertices)
; there functions that take vertex structures  (s-vertices)

; there are functions that take m-edge-patterns (m-edges)  [[names]] or
;                                                          [[names],weight]
; there are functions that take n-edge-patterns (np-edges) ((names) weight) 
; there are functions that take edge-patterns (p-edges)    ((vertices) weight) 
; there are functions that take edge-structures (edges)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 

;; try to replace vertex-indices with vertices ... this shouldn't need more
;; space, since Lisp keeps only the pointer to the vertex

;; one problem arises with functions that take vertices numbered from 0 to n-1,
;; as the adjacency-matrix stuff for example. I need a mapping from vertex to
;; index here! Well, I can store this in the vertex, but this might cause
;; problems when I want to delete or change vertices: Replacing one vertex
;; means to find this vertex in all edges and changing it there, too. If I
;; delete vertex 1 then all vertices have to be replaced...

;; Changing the vertex does not work, since it is used by other graphs too.

;; Hence I need a function vertices-to-indices like
;(defun vertices-to-indices (all-vertices vertices)
;  (mapcar #'(lambda (v) (position v all-vertices :test #'eq)) vertices))

;; it might be good to have a tmp register for each vertex so that I can cache
;; the result. Returns the number of vertices
(defun cache-all-vertex-indices (vertices)
  (let ((i -1))
    (mapcar #'(lambda (v) (setf (vertex-index v) (incf i))) vertices)
    (1+ i)))

(defun cache-some-vertex-indices (all-vertices vertices)
  (mapcar #'(lambda (v) (setf (vertex-index v)
			      (position v all-vertices :test #'eq)))
	  vertices))

(defun vertices-to-indices (vertices undirected)
  (if undirected (sort (mapcar #'vertex-index vertices) #'<)
    (mapcar #'vertex-index vertices)))

; s-vertices is a list of vertex-structures
; n-vertex is a vertex-name 
; returns the vertex-structure corresponding to n-vertex
(defun n-vertex-to-vertex (s-vertices n-vertex)
  (find n-vertex s-vertices :test #'equal :key #'vertex-name))

; s-vertices is a list of vertex-structures
; n-verices is a list of vertex-names
; returns the list of vertices corresponding to n-vertices
(defun n-vertices-to-vertices (s-vertices n-vertices)
  (mapcar #'(lambda (vtx) (n-vertex-to-vertex s-vertices vtx)) n-vertices))

(defun m-edges-to-np-edges (m-edges &optional fun)
  (if (mlistp m-edges)
      (mapcar #'(lambda(e) (if (mlistp (car (margs e)))
			       (cons (margs (car (margs e))) (cdr (margs e)))
			     (merror "The first element of each element of the 
                                     argument to ~M should be a list" fun))) 
	      (margs m-edges))
    (if fun (merror  "Argument to ~M should be a list" fun)
      (merror "Argument to m-edges-to-np-edges should be an mlist"))))

; returns an edge-pattern
(defun m-edges-to-p-edges (s-vertices m-edges)
  (mapcar #'(lambda (m-edge) 
	      (cons (n-vertices-to-vertices s-vertices (car m-edge))
		    (cdr m-edge)))
	  (m-edges-to-np-edges m-edges)))

; returns edge-structures
(defun m-edges-to-edges (hgraph m-edges)
  (p-edges-to-edges (G-edges hgraph) 
		    (m-edges-to-p-edges (G-vertices hgraph) m-edges)))
  
; returns edge-structures
(defun p-edges-to-edges (edges p-edges)
  (mapcar #'(lambda (e) (cond ((find-edge e edges))
			      (T (merror "Not an edge"))))
	  p-edges))
		       
; p-edge is an edge-pattern
; returns the edge-structure
(defun find-edge (p-edge edges)
  (let ((res (remove-if (if (single p-edge)
			    #'(lambda (e)
				(not (equal (car p-edge)
					    (edge-vertices e))))
			  #'(lambda (e)
			      (not (equal p-edge
					  (list (edge-vertices e)
						(edge-weight e))))))
			edges)))
    (if (single res) (car res)
      (find 1 res :key #'edge-weight)))))

; for edge-structures only, cares only about the key
(defun member-edge (s-edge s-edges)
  (member (edge-key s-edge) s-edges :key #'edge-key :test #'=))

; s-edges have to be edge-structures
; returns the edges in edges corresponding via edge-key to those in
; s-edges

; provided that s-edges contains every edge at most once, I could
; remove found edges from all-edges
; here I also assume that both arguments are sorted by key
(defun sorted-p (sequence predicate)
  (every predicate sequence (cdr sequence)))

(defun map-edges-sorted (s-edges all-edges)
  (if (and (sorted-p s-edges #'(lambda (e f) (< (edge-key e)
 						(edge-key f))))
	   (sorted-p all-edges #'(lambda (e f) (< (edge-key e)
 						  (edge-key f)))))
      (map-edges-sorted-aux s-edges all-edges)
    (error "s-edges: ~M or all-edges: ~M are not sorted"
 	   s-edges all-edges)))

; Ivan Boldyrev (boldyrev@uiggm.nsc.ru)
(defun map-edges-sorted-aux (s-edges all-edges)
  (labels ((find-next (key all-edges)
		      (if (= key (edge-key (car all-edges)))
			  (values (car all-edges)
				  (cdr all-edges))
			(find-next key (cdr all-edges))))
	   (find-all (s-edges all-edges)
		     (when s-edges
		       (multiple-value-bind
			   (item all-rest) (find-next (edge-key (car s-edges))
						      all-edges)
			 (cons item (find-all (cdr s-edges) all-rest))))))
    (find-all s-edges all-edges)))

; old versions:
;
;; (defun map-edges-sorted-aux (s-edges all-edges)
;;   (when s-edges
;;     (let ((rest (member (edge-key (car s-edges)) 
;; 			all-edges :key #'edge-key :test #'=)))
;;       (cons (car rest)
;; 	    (map-edges-sorted-aux (cdr s-edges) (cdr rest))))))

;; (defun map-edges-sorted (s-edges all-edges)
;;   (labels ((rec (s-edges all-edges collect)
;; 		(if s-edges
;; 		    (let ((rest 
;; 			   (member (edge-key (car s-edges)) 
;; 				   all-edges :key #'edge-key :test #'=)))
;; 		      (rec (cdr s-edges) 
;; 			   (cdr rest)
;; 			   (cons (car rest) collect)))
;; 		  collect)))
;; 	  (nreverse (rec s-edges all-edges nil))))

; vertices is an mlist, edges an mlist of mlists of type [mlist, weight]
; or of type [mlist] ie. [[[1,2],w1],[[2,3]]]
(defun $hypergraph (m-vertices m-edges &optional undirected)
  (if (mlistp m-vertices)
      (hypergraph (margs m-vertices) 
		  (m-edges-to-np-edges m-edges) undirected)
    (merror "First Argument to hypergraph should be a list of vertex-names")))

(defun $graph_plot (hgraph)
  (if (hypergraphp hgraph)
      (if (eql (G-dim hgraph) 2)
	  (progn
	    (cache-all-vertex-indices (G-vertices hgraph))
	    (with-open-file 
	     (st1 (concatenate 'string 
			       (namestring (car (directory "$TGF")))
			       "/maxout.graph.txt") 
		  :direction :output)
	     (format st1 "graph~%~{~{~D ~}~%~}0 0" 
		     (mapcar #'(lambda (e)  
				 (mapcar #'1+ (vertices-to-indices 
					      (edge-vertices e)
					      (G-undirected hgraph))))
			     (G-edges hgraph))))
	    ($system (concatenate 'string 
				  "pigale -fi "
				  (namestring (car (directory "$TGF")))
				  "/maxout.graph.txt")))
	(merror "I can only plot graphs. Complain to Hubert de
                 Fraysseix <hf@ehess.fr> or Patrice Ossona de Mendez
                 <pom@ehess.fr>")) 
    (merror "The argument to graph_plot has to be a hypergraph")))

; 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!
; edges is a list of ((name1 name2 ...) weight)
(defun hypergraph (n-vertices np-edges &optional undirected)
; sollte Argumente prüfen
  (let ((hgraph (make-hypergraph)))
    (prop-add hgraph 'undirected undirected)
    (prop-add hgraph 'vertices 
	      (mapcar #'(lambda (v) (make-vertex :name v))
		      n-vertices))
    (prop-add hgraph 'edges
	      (let ((key 0)
		    (dim (length (caar np-edges))))
		(prog1 
		    (mapcar 
		     #'(lambda (e)
			 (unless (and dim (= (length (car e)) dim))
			   (setq dim nil))
			 (make-edge :vertices (n-vertices-to-vertices
					       (G-vertices hgraph) (car e))
				    :weight (if (single e) 1 (cadr e))
				    :key (incf key)))
		     np-edges)
		  (prop-add hgraph 'max-edge-key key)
		  (prop-add hgraph 'dim dim))))
    hgraph))

; from here on, the key to an edge is a list of vertices!
; 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?


; **************************** remove-edges ***************************
; new idea: let all sub-operations take the same parameters as the
; operation itself, plus the additional which are needed: eg:

;(defun remove-edges-adjacency-matrix (hgraph edges)
;(defun remove-edges-edges (hgraph edges)
;(defun remove-edges-perfect-matchings (hgraph edges)
;(defun remove-edges-components (hgraph edges 
;                                &optional new-edges
;                                          (remove-edges-edges hgraph
;                                                              edges))
; the last (and only non-trivial) example is rather bad, because in
; this case it probably does make more sense to compute the new hgraph...
; -----------------
; this idea is not too good, because it is nice to have access to these
; functions even without having an hgraph!

(defun remove-edges-adjacency-matrix (vertices edges value undirected)
  (let ((newmatrix (copy-array value)))
    (cache-all-vertex-indices vertices)
    (dolist (edge edges (update-adjacency-matrix newmatrix 
						 edges
						 undirected))
      (array-subtract-entry newmatrix
			    (vertices-to-indices (edge-vertices edge)
						 undirected)
			    (edge-weight edge)))))

(defun remove-edges-edges (edges value)
  (list-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)))

; new-edges are the edges of the new hgraph, i.e., edges(hgraph) \ edges
; only the components which contain one of the affected-vertices need
; to be considered
; if (car cmp) contains one of these, the components of (car cmp) have
; to be recomputed
(defun remove-edges-components (edges value new-edges)
  (do ((affected-vertices (delete-duplicates 
			   (mapcar #'(lambda (e) (car (edge-vertices e)))
				   edges)
			   :test #'eq)
; it is sufficient to mark one vertex of each edge as affected
			  (if (cdr cmp)
			      (nlist-difference affected-vertices (car cmp)
						:test #'eq)))
       (cmp value 
	    (cdr cmp))
       (res nil
	    (if (some #'(lambda (v) (member v (car cmp)))
		      affected-vertices)
		(nconc res (components (car cmp) new-edges))
	      (push (car cmp) res))))
      ((null cmp) res)))

(defun $remove_edges (hgraph m-edges)
  (if (hypergraphp hgraph)
      (remove-edges hgraph (m-edges-to-edges hgraph m-edges))
    (merror "First argument to remove_edges must be a hypergraph")))

; 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. Why???
; edges is a list of edge-structures
(defun remove-edges (hgraph edges)
  (let ((newhgraph (make-hypergraph)))
    (prop-add newhgraph 'undirected   (G-undirected hgraph))
    (prop-add newhgraph 'max-edge-key (G-max-edge-key hgraph))
    (prop-add newhgraph 'vertices     (G-vertices hgraph))
    (prop-add newhgraph 'edges        
	      (remove-edges-edges edges (G-edges hgraph)))
    (map-props 
     #'(lambda (prop value)
	 (block nil
		(prop-add 
		 newhgraph prop
		 (case prop
; needs new-edges, therefore it has to be delayed.
; forcing could be bad, if the result were nil...
; 		       ('dim (cond ((null value) (delay (compute-dim
; 							 (G-edges
; 							  newhgraph))))
; 				   ((numberp value) value)
; 				   (T (let ((d (force value)))
; 					(if (null d)
; 					    (delay newhgraph 'dim
;                                                  (compute-dim (G-edges
; 								 newhgraph)))
; 					  d))))
; or - no forcing - I take this because its simpler
		   ('dim (if (numberp value) value
			   (delay newhgraph prop 
				  (compute-dim (G-edges newhgraph)))))
		   ('adjacency-matrix 
		    (remove-edges-adjacency-matrix (G-vertices hgraph)
						   edges 
						   (force value)
						   (G-undirected hgraph)))

		   ('components
		    (delay newhgraph prop 
			   (remove-edges-components edges
						    (force value)
						    (G-edges newhgraph))))
		   ('perfect-matchings
		    (delay newhgraph prop 
			   (remove-edges-perfect-matchings edges 
							   (force value))))
		   (T (return))))))
     hgraph)
    newhgraph))

; **************************** remove-loops ***************************

(defun $remove_loops (hgraph)
  (if (hypergraphp hgraph) (remove-loops hgraph)
; do I really need equalp here?
    (merror "the argument of remove_loops has to be a hypergraph")))

; I really need support for components here
(defun remove-loops (hgraph &optional (number-of-identical-vertices 1))
  (let ((newhgraph (make-hypergraph))
	removed-loops)
    (prop-add newhgraph 'undirected   (G-undirected hgraph))
    (prop-add newhgraph 'max-edge-key (G-max-edge-key hgraph))
    (prop-add newhgraph 'vertices     (G-vertices hgraph))
    (prop-add newhgraph 'edges        
	      (remove-if #'(lambda (e) 
			     (when (loop-p e number-of-identical-vertices)
			       (push e removed-loops)))
			 (G-edges hgraph)))
    (prop-add newhgraph 'dim 
	      (if (numberp (prop-get hgraph 'dim)) (g-dim hgraph)
		(delay newhgraph 'dim 
		       (compute-dim (G-edges newhgraph)))))
    (map-props 
     #'(lambda (prop value)
	 (block nil
	   (prop-add 
	    newhgraph prop
	    (case prop
	      ('components 		      
	       (if (eql (g-dim hgraph) 
			(1+ number-of-identical-vertices))
		   value
		 (delay newhgraph prop 
			(remove-edges-components removed-loops
						 (force value)
						 (G-edges newhgraph)))))
	      (T (return))))))
     hgraph)
    newhgraph))

; **************************** graph-difference ***************************

(defun graph-difference-adjacency-matrix (hgraph1 hgraph2 newhgraph)
; should I compute the adjacency matrix of hgraph2 in any case? I
; don't think so
  (if (numberp (G-dim newhgraph))
      (multiple-value-bind 
       (value presentP) (prop-get hgraph2 'adjacency-matrix)
       (if presentP (array-difference 
		     (G-adjacency-matrix hgraph1) (force value))
	 (progn
	   (cache-all-vertex-indices (G-vertices hgraph1))
	   (remove-edges-adjacency-matrix (G-vertices hgraph1)
	                                  (G-edges hgraph2)
					  (G-adjacency-matrix hgraph1)
					  (G-undirected newhgraph)))))
    (merror "cannot compute the adjacency matrix of the difference of
             ~M and ~M" hgraph1 hgraph2)))

(defun $graph_difference (hgraph1 hgraph2)
  (if (and (hypergraphp hgraph1) (hypergraphp hgraph2))
; do I really need equalp here?
      (if (and (equalp (G-vertices hgraph1) (G-vertices hgraph2)) 
	       (eq (G-undirected hgraph1) (G-undirected hgraph2)))
	  (graph-difference hgraph1 hgraph2)
	(merror "graph-difference is defined for graphs with identical vertex
            set only"))
    (merror "the arguments of graph-difference have to be hypergraphs")))

; hgraph1 and hgraph2 need to have identical vertex set, the result is
; a graph that has the edges of hgraph1 with positive weight, those of
; hgraph2 with negative weight.
; In a later stage, I might want to allow that the types differ, in
; this case the type would have to be directed.
(defun graph-difference (hgraph1 hgraph2)
  (let ((newhgraph (make-hypergraph)))
; iterate over the properties of hgraph1, if a property is present
; only in hgraph2 - what could I do...
; I would really need to have map-props iterate over the union...
    (prop-add newhgraph 'undirected   (and (G-undirected hgraph1)
					   (G-undirected hgraph2)))
    (prop-add newhgraph 'max-edge-key (+ (G-max-edge-key hgraph1) 
					 (G-max-edge-key hgraph2)))
    (prop-add newhgraph 'vertices     (G-vertices hgraph1))
    (prop-add newhgraph 'edges        
	      (append (G-edges hgraph1)
; make new edges for the edges in hgraph2. Unfortunately, here the key does not
; stay the same - well, the weight doesn't eiter
		      (mapcar #'(lambda (e)
				  (make-edge :vertices (edge-vertices e)
					     :weight (list '(mtimes) -1
							   (edge-weight e))
					     :key (+ (edge-key e) 
						     (G-max-edge-key 
						      hgraph1))))
			      (G-edges hgraph2))))
    (map-props 
     #'(lambda (prop value)
	 (block nil
		(prop-add 
		 newhgraph prop
		 (case prop
		       ('adjacency-matrix 
			(delay newhgraph prop
			       (graph-difference-adjacency-matrix
				hgraph1 hgraph2 newhgraph)))		
		       (T (return))))))
     hgraph1)
    newhgraph))
; just realised that remove-edges is just a variant of
; graph-difference, in that it removes only edges which are really present...
; That is, given graph-difference, I could redo remove-edges by saying
;
; (graph-difference hgraph (hypergraph (g-vertices hgraph) edges))
;
; Not quite: doing so I get each removed edge twice but with opposite
; weights
;
; Maybe I can recycle some thought put into remove-edges...

; **************************** contract-vertices *********************

; now much easier, I only need to replace the vertices in vertex-subset
(defun contract-vertices-replace (newvertex vertices vertex-subset)
  (substitute-if newvertex 
		 #'(lambda (v) (member v vertices :test #'eq))
		 vertex-subset)))

(defun contract-vertices-components (value newvertex vertices)
; replace the vertices in the components by (car vertices) and 
; collect all components containing (car vertices) into one, the
; others remain untouched
  (let (ncmp ocmp)
    (dolist (tcmp (mapcar #'(lambda (cmp) 
			      (contract-vertices-replace newvertex vertices 
							 cmp))
			  value)
		  (cons (delete newvertex ncmp :count (1- (length vertices)) 
				:test #'eq)
			ocmp))
      (if (member newvertex tcmp :test #'eq) 
	  (setq ncmp (nconc ncmp tcmp))
	(push tcmp ocmp)))))

(defun $contract_vertices (hgraph m-vertices)
  (cond ((not (hypergraphp hgraph))
	 (merror "first argument of contract_vertices has to be a hypergraph"))
	((not (mlistp m-vertices))
	 (merror "second argument of contract_vertices has to be a
                  list of vertices"))
	((> (length (margs m-vertices)) 1)
	 (contract-vertices hgraph
			    (n-vertices-to-vertices (G-vertices hgraph)
						    (margs m-vertices))
			    m-vertices))
	(T hgraph)))

; I assume that vertices has at least 2 elements
(defun contract-vertices (hgraph vertices newname)
  (let ((newhgraph (make-hypergraph)))
    (prop-add newhgraph 'undirected   (G-undirected hgraph))
    (prop-add newhgraph 'dim          (G-dim hgraph))
    (prop-add newhgraph 'max-edge-key (G-max-edge-key hgraph))
    (prop-add newhgraph 'vertices     
	      (cons (make-vertex :name newname)
		    (list-difference (G-vertices hgraph) vertices
				     :test #'eq)))
    (prop-add newhgraph 'edges
			(mapcar 
			 #'(lambda (e)
			     (make-edge :vertices 
					(contract-vertices-replace 
					 (car (G-vertices newhgraph))
					 vertices
					 (edge-vertices e))
					:weight (edge-weight e)
					:key (edge-key e)))
			 (G-edges hgraph)))
    (map-props 
     #'(lambda (prop value)
	 (block nil
		(prop-add 
		 newhgraph prop
		 (case prop
		   ('components 
		    (delay newhgraph prop
			   (contract-vertices-components
			    (force value) (car (G-vertices newhgraph))
			    vertices)))
		   (T (return))))))
     hgraph)
    newhgraph))

; **************************** graph-neps ***************************

; **************************** graph-composition ********************

; gets the vertices in newhgraph corresponding to the subgraph corresponding to
; vertex
(defun graph-composition-get-vertices (newhgraph vertex)
  (let (res)
    (maphash #'(lambda (prop v) (push v res))
	     (gethash vertex (prop-get newhgraph 'subgraphs)))
    res))

; gets the vertex in newhgraph corresponding to the vertex in subhgraph
(defun graph-composition-get-vertex (newhgraph subhgraph vertex)
  (gethash vertex (gethash subhgraph (prop-get newhgraph 'subgraphs))))

(defun graph-composition-edges (hgraph hgraphs newhgraph)
  (let ((key 0))
; I really have to copy the vertices, because it is possible that they are
; shared between some hgraphs
    (prog1
	(append 
	 (mapcan #'(lambda (v h)
		     (mapcar #'(lambda (e) ; e is an edge of a subgraph
				 (make-edge 
				  :vertices (mapcar 
					     #'(lambda (w) ; v is a vertex of a subgraph
						 (graph-composition-get-vertex 
						  newhgraph v w))
					     (edge-vertices e))
				  :weight (edge-weight e)
				  :key (incf key)))
			     (G-edges h)))
		 (G-vertices hgraph)
		 hgraphs)
; additionally, for each edge (v1,v2,...) in h there are edges from each vertex
; of hv1 to each vertex of hv2 ...
	 (mapcan #'(lambda (e) ; e is an edge of hgraph
		     (mapcar #'(lambda (vertices)
				 (make-edge
				  :vertices vertices
				  :weight (edge-weight e)
				  :key (incf key)))
			     (cartesian
			      (mapcar 
			       #'(lambda (v) ; v is a vertex of hgraph
					; I need the vertices in newhgraph
					; corresponding to the vertices in h_v 
				   (graph-composition-get-vertices
				    newhgraph v))
			       (edge-vertices e)))))
		 (G-edges hgraph)))
      (prop-add newhgraph 'max-edge-key key))))

(defun $graph_composition (hgraph &rest hgraphs)
  (if (and (hypergraphp hgraph) (every #'hypergraphp hgraphs))
      (let ((l (length hgraphs))
	    (p (length (G-vertices hgraph))))
	(cond ((= l 1)
	       (graph-composition hgraph 
				  (make-list p 
					     :initial-element (car hgraphs))))
	      ((= l p) (graph-composition hgraph hgraphs))
	      (T (merror "The number of arguments to graph_composition must 
                          either be 2 or 1 plus the number of vertices of 
                          the first argument"))))
    (merror "The arguments to graph_composition have to be hypergraphs")))

(defun graph-composition (hgraph hgraphs)
  (let ((newhgraph (make-hypergraph)))
    (prop-add newhgraph 'undirected   (and (G-undirected hgraph)
					   (every #'G-undirected hgraphs)))
    (prop-add newhgraph 'dim          (cond ((null (G-dim hgraph)) 
					     nil)
					    ((zerop (G-dim hgraph))
					     (G-dim (car hgraphs)))
					    ((every #'(lambda (h) (or (eql (G-dim h) 
									   (G-dim hgraph))
								      (and (numberp (G-dim h))
									   (zerop (G-dim h)))))
						    hgraphs)
					     (G-dim hgraph))))
; 'max-edge-key is set in 'edges     
    (prop-add newhgraph 'vertices     
; I make the vertices so that vertex 0 to p1-1 are the vertices of hgraph1, p1
; to p1+p2-1 are the vertices of hgraph2 and so on. Thus if I want to reference
; the second vertex of the third graph I can do (nth p1+p2+2 vertices)
; This means however, that I have to index the vertices of all the subgraphs.
; A different approach would be to set up a table with the vertices of the
; hgraphs as keys and the new vertices as values
	      (let ((subgraphs (make-hash-table :test #'eq))
		    current)
		(prop-add newhgraph 'subgraphs subgraphs)
		(mapcan #'(lambda (v h)
			    (setq current (setf (gethash v subgraphs) 
						(make-hash-table :test #'eq)))
			    (mapcar 
			     #'(lambda (w)
				 (setf (gethash w current)
				       (make-vertex 
					:name (make-mlist (vertex-name v)
							  (vertex-name w)))))
			     (g-vertices h)))
			(g-vertices hgraph)
			hgraphs)))
    (prop-add newhgraph 'edges (graph-composition-edges hgraph hgraphs newhgraph))
;;     (map-props 
;;      #'(lambda (prop value)
;; 	 (block nil
;; 	   (prop-add 
;; 	    newhgraph prop
;; 	    (case prop
;; 	      (T (return))))))
;;      hgraph)
    newhgraph))

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

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

; returns nil if the hgraph is not known to be uniform, i.e., there
; are two edges of different lengths, otherwise the length of each edge. If
; edges is nil, it returns 0.
; edges are edge-structures here
(defun compute-dim (edges)
  (let ((d (length (edge-vertices (car edges)))))
    (when (every #'(lambda (e) (= d (length (edge-vertices e)))) 
		 (cdr edges))
      d)))

(defun G-dim (hgraph) 
  (multiple-value-bind 
   (value presentP) (prop-get hgraph 'dim)
   (if presentP (force value)
     (prop-add hgraph 'dim (compute-dim (G-edges hgraph))))))

(defun $graph_dim (hgraph) 
  (if (hypergraphp hgraph) (G-dim hgraph)
    (merror "The argument to graph_dim has to be a hypergraph")))

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

; there are more than number-of-identical-vertices identical vertices in edge
; this vertex is returned
(defun loop-p (edge &optional (number-of-identical-vertices 1))
  (let (count res)
    (mapc #'(lambda (v)
	      (setq res (assoc v count))
	      (if res 
		  (if (= (cdr res) number-of-identical-vertices)
		      (return-from loop-p v)
		    (incf (cdr res)))
		(setq count (acons v 1 count))))
	  (edge-vertices edge))
    nil))

(defun G-max-edge-key (hgraph) (prop-get hgraph 'max-edge-key))

; edge is an edge-structure, returns the m-edge, i.e. [[names], weight]
(defun m-edge (edge)
  (make-mlist-simp 
   (make-mlist-simp-l (mapcar #'vertex-name (edge-vertices edge)))
   (edge-weight edge)))

(defun m-edges (edges)
  (make-mlist-simp-l (mapcar #'m-edge edges)))

(defun p-edge (edge)
  (mapcar #'vertex-name (edge-vertices edge)))

(defun p-edges (edges)
  (mapcar #'p-edge edges))

(defun $edges (hgraph) 
  (if (hypergraphp hgraph) (m-edges (G-edges hgraph))
    (merror "The argument to edges has to be a hypergraph")))

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

; returns a list of lists of vertices
(defun edges-vertices (edges) (mapcar #'edge-vertices edges))

; vertices are vertex-structure, returns the vertex names
(defun m-vertices (s-vertices)
  (make-mlist-simp-l (mapcar #'vertex-name s-vertices)))

(defun $vertices (hgraph) 
  (if (hypergraphp hgraph) (m-vertices (G-vertices hgraph))
    (merror "The argument to vertices has to be a hypergraph")))

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

(defun $adjacency_matrix (hgraph)
  (if (hypergraphp hgraph)
      (if (numberp (G-dim hgraph))
	  (if (= 2 (G-dim hgraph))
	      (let ((m (G-adjacency-matrix hgraph)))
		(apply #'$genmatrix m
		       (nconc (make-list 2 
					 :initial-element 
					 (1- (array-dimension m 0)))
			      (make-list 2 :initial-element 0))))
	    (G-adjacency-matrix hgraph))
	(merror "a non-uniform hgraph has no adjacency-matrix"))
    (merror "The argument to adjacency_matrix has to be a hypergraph")))

; edges is a list of lists of vertices here, mat the adjacency matrix.
; assumes that the vertices are properly indexed
(defun update-adjacency-matrix (mat edges undirected)
  (let ((ind (remove-duplicates (mapcar #'(lambda (e)
					    (vertices-to-indices e undirected))
; in the undirected case, I need to keep only the sorted lists...
					edges))))
    (if undirected
; 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 (val)
	  (dolist (edge ind mat)
	    (setq val nil)
; val needs to be initialised to nil at every iteration!
	    (fpermute edge
		      #'(lambda () 
			  (if val 
; the first call is the identity, this value needs only to be meval'd
			      (setf (apply #'aref mat edge) val)
			    (setq val (meval (apply #'aref mat edge))))))))
      (dolist (edge ind mat) (setf (apply #'aref mat edge) 
				   (meval (apply #'aref mat edge)))))))

; 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 (cache-all-vertex-indices (G-vertices hgraph)))
	    (res (make-array (make-list (G-dim hgraph)
					:initial-element dim) 
			     :initial-element 0)))
       (dolist (e edges)
	       (array-add-entry res (vertices-to-indices (edge-vertices e) 
							 (G-undirected hgraph))
				(edge-weight e)))
       (update-adjacency-matrix res 
				(edges-vertices edges)
				(G-undirected hgraph))
       (prop-add hgraph 'adjacency-matrix res)))))

(defun $spanning_trees (hgraph)
  (if (hypergraphp hgraph)
      (make-mlist-simp-l (mapcar #'m-edges (G-spanning-trees hgraph)))
    (merror "The argument to spanning_trees has to be a hypergraph")))

; this is only for graphs right now
(defun G-spanning-trees (hgraph)
  (multiple-value-bind 
   (value presentP) (prop-get hgraph 'spanning-trees)
   (if presentP (force value)
     (prop-add hgraph 'spanning-trees
	       (labels ((rec (hgraph)
; t=t\e + e*t|e
; a graph with one vertex only has one spanning tree with no edges
			     (cond ((single (G-vertices hgraph))
				    (list nil))              ; only loops
; I would like to remove the loops here...
				   ((not (single (G-components hgraph)))
				    nil)                     ; disconnected

				   ((and (numberp (g-dim hgraph))
					 (zerop (- (length (G-vertices hgraph))
						   (* (1- (g-dim hgraph))
						      (length (G-edges
							       hgraph)))
						   1)))
				    (list (G-edges hgraph))) ; tree
; find a non-loop
				   (T (let ((e (find-if 
						(lisp::complement #'loop-p)
						(G-edges hgraph))))
; contract - delete and replace the edges via the edge-keys
; with the original ones. neither contracting nor removing may change
; the edge-key to let this work.
					(append 
					 (rec (remove-edges hgraph (list e)))
					 (cons-all 
					  e (rec (remove-loops
						  (contract-vertices 
						   hgraph 
						   (edge-vertices e)
						   nil))))))))))
		 (mapcar #'(lambda (tree)
; map-edges expects to have both arguments sorted by key...
			     (map-edges-sorted tree (G-edges hgraph)))
			 (rec (remove-loops hgraph))))))))

(defun $perfect_matchings (hgraph)
  (if (hypergraphp hgraph)
      (make-mlist-simp-l (mapcar #'(lambda (edges) (m-edges edges))
				 (G-perfect-matchings hgraph)))
    (merror "The argument to perfect_matchings has to be a hypergraph")))


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

(defun $components (hgraph)
  (if (hypergraphp hgraph)
      (make-mlist-simp-l (mapcar #'m-vertices (G-components hgraph)))
    (merror "The argument to components has to be a hypergraph")))

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

; vertex is the name of a vertex here
;(defun $degree (hgraph n-vertex)
;  (if (hypergraphp hgraph)
;      (G-degree hgraph (n-vertex-to-vertex (G-vertices hgraph) n-vertex))
;    (merror "The first argument to degree has to be a hypergraph")))

; this is the total degree of vertex, i.e. for undirected graphs twice
; the degree...
; (defun adjacency-degree (m vertex)
;   (let* ((dim (array-dimension m 0))
; 	 (tot (expt dim (1- (array-rank m)))))
;     (do ((iter 1 (* iter dim))
; 	 (size tot 
; 	       (/ size dim))
; 	 (ofst (* vertex tot) 
; 	       (/ ofst dim))
; 	 (res  nil
; 	       (dotimes (i iter res)
; 			(setq res
; 			      (append res 
; 				      (coerce (make-array 
; 					       size
; 					       :displaced-to m
; 					       :displaced-index-offset 
; 					       (+ (* i size dim) ofst))
; 					      'list))))))
; 	((> iter tot) (cons '(mplus) res)))))
; unfortunately, this is not correct in the case of undirected
; hypergraphs, eg for three-graphs, the degree of non-loops is
; 2 instead of 1, the degree of two-loops is 2 and the degree of a
; triple-loop is 1 instead of 3

; correct way would be: count triple loops (i i i) 
;                       count double loops (i i j) (i j j) with i < j
;                       count single elems (i j k) with i < j < k
; this seems to get rather complicated...
; even for graphs I would have to treat loops specially

(defun $degree_sequence (hgraph)
  (if (hypergraphp hgraph)
      (make-mlist-l (coerce (G-degree-sequence hgraph) 'list))
    (merror "The argument to degree_sequence has to be a hypergraph")))

(defun G-degree-sequence (hgraph)
  (multiple-value-bind 
   (value presentP) (prop-get hgraph 'degree-sequence)
   (if presentP (force value)
     (let ((degs (make-array (list (length (G-vertices hgraph))) 
			     :initial-element 0)))
       (cache-all-vertex-indices (g-vertices hgraph))
       (mapc #'(lambda (e) 
		 (mapcar #'(lambda (v)
			     (array-add-entry degs (list (vertex-index v))
					      (edge-weight e)))
			 (edge-vertices e)))
	     (G-edges hgraph))
       (prop-add hgraph 'degree-sequence degs)))))

(defun $Tutte (hgraph x y)
  (if (and (hypergraphp hgraph) (eql 2 (g-dim hgraph)))
      (g-tutte hgraph x y)
    (merror "The argument to Tutte has to be a graph")))

(defun G-Tutte (hgraph x y)
  (multiple-value-bind 
   (value presentP) (prop-get hgraph 'Tutte)
   (if presentP (force value)
     (prop-add 
      hgraph 'Tutte
      (labels ((rec (hgraph)
		    (cond ((null (G-edges hgraph)) 
; no edges
			   1)
			  ((single (G-vertices hgraph))
; only loops
			   (list '(MEXPT SIMP) y (length (G-edges hgraph))))
			  (T (let ((e (find-if (lisp::complement #'loop-p) 
					       (G-edges hgraph))))
; find a non-loop
			       (if e
				   (let ((h (remove-edges hgraph (list e))))
; contract - delete 
				     (if (< (length (G-components hgraph))
					    (length (G-components h)))
; isthmus		       
					 (list '(mtimes) x (rec h))
; otherwise
				       (list '(mplus) 
					     (rec h)
					     (rec (contract-vertices 
						   h (edge-vertices e) nil)))))
; only loops
				 (list '(MEXPT SIMP) y 
				       (length (G-edges hgraph)))))))))
	(rec hgraph))))))


;********************************************************************
;  helper routines
;********************************************************************

; 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 #'eq)) 
		       (edge-vertices e)))
	     edges))

; returns all isolated vertices
(defun find-isolated-vertex (vertices edges)
  (list-difference vertices (reduce #'append 
				    (mapcar #'edge-vertices edges))
		   :test #'eq))

(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)
  (delete-duplicates
   (reduce #'append 
	   (delete-if #'(lambda (e) (not (member vertex e :test #'eq)))
		      (mapcar #'edge-vertices edges)))
   :test #'eq))

; 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)
	    (unless (member vtx res :test #'eq)
	      (setq res (component edges vtx res))))))

(defun components (vertices edges)
  (do ((res nil)
       (vtx vertices
	    (list-difference vtx (car res) :test #'eq)))
      ((null vtx) res)
      (push (component edges (car 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 
				  (list-difference vertices (edge-vertices 
							     (car edgesrest)))
				  (non-incident (car edgesrest)
						(cdr edgesrest)))))
		   edges))))

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

;********************************************************************
;  path, cycle, wheel, fan
;********************************************************************

(defun $path (n)
  (hypergraph (range n)
	      (let (res) 
		(dotimes (i (1- n) res)
		  (push (list (list i (1+ i))) res))) T))

(defun $cycle (n)
  (hypergraph (range n)
	      (let (res) 
		(dotimes (i (1- n) (cons (list (list 0 (1- n))) res))
		  (push (list (list i (1+ i))) res))) T))

(defun complete-product (hgraph1 hgraph2)
  (graph-composition ($path 2) (list hgraph1 hgraph2)))

(defun $complete_product (hgraph1 hgraph2)
  (if (and (hypergraphp hgraph1) (hypergraphp hgraph2))
      (complete-product hgraph1 hgraph2)
    (merror "the arguments of graph-difference have to be hypergraphs")))

(defun $wheel (n)
  ($complete_product ($complete_graph 1) ($cycle (1- n))))

(defun $fan (n)
  ($complete_product ($complete_graph 1) ($path (1- n))))

;********************************************************************
;  COMPLETE GRAPH *** complete-graph (n)
;********************************************************************
; weights is a function of two arguments ranging from 0 to n
(defun $complete_graph (n &optional weights)
  (threshold-graph (make-list n :initial-element (1- n)) weights))

;********************************************************************
;  THRESHOLD GRAPH *** threshold-graph (partition)
;********************************************************************

(defun $threshold_graph (partition &optional weights)
  (threshold-graph (margs partition) weights))

(defun threshold-graph (partition &optional weights)
  (hypergraph 
   (range (length partition))
   (do ((i 0 (1+ i))
	(p partition (cdr p))
	(res nil (do ((j (1+ i) (1+ j)))
		     ((> j (car p)) res)
		   (push (if weights
			     (list (list i j)
				   (mfuncall weights i j))
			   (list (list i j)))
			 res))))
       ((null p) res))
   T))

;********************************************************************
; 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
;
; 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)
  (labels ((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)))))))
    (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 T))))