\begin{code}
{-# OPTIONS -fglasgow-exts #-}

module Main
where

import IO
import List
import Maybe
import Data.Set (Set)
import qualified Data.Set as Set
import System
import Monad (foldM,msum)
\end{code}
Permutations are written as lists of integers on the set $E_k$ where 
the $i$-th entry determines where to map $i$. For example, \verb+[0,2,3,1,5]+ symbolizes
a permutation $\pi$ on $E_5$ with $\pi = (1\,2\,3)$. The function \verb+allperm+ generates all 
permutations on the set $E_h$ where \verb+allperm'+ is a helper function which works on any set of integers.
\begin{code}
type Perm = [Int]
type SPerm = Set Perm

allperm' :: [Int] -> [Perm]
allperm' [x] = [[x]]
allperm' xs = concat [ map (\ys -> x:ys) (allperm' (xs \\ [x])) | x <- xs ]

allperm :: Int -> [Perm]
allperm h = allperm' $ [0..h-1]
\end{code}

In the program we need all subgroups of $S_h$ at one point. To generate all these groups we need to multiply
two permutations $a,b: E_h \to E_h$ and this can be done by \verb+multperm+. If $c = \verb+multperm+\ a\ b$
then $c$ is a permutation with $c(x) = b(a(x))$ for all $x \in E_h$. Given two sets of permutations 
$p_a, p_b \subseteq S_h$ we can form the set
of all binary products \[ 
p_{ab} := \verb+prodperm+\ p_a\ p_b = \{ \verb+multperm+\ a\ b \mid  a \in p_a, b \in p_b \}.
\]
\begin{code}
multperm :: Perm -> Perm -> Perm
multperm a b = map (b !!) a

prodperm :: SPerm -> SPerm -> SPerm
prodperm pa pb = Set.fromList . concat $ map (\b -> map (\a -> multperm a b) (Set.toList pa)) (Set.toList pb) 
\end{code}

We generate a group $G = [A]$ from a set $A \subseteq S_h$ recursively.
Let $A_0 := A$ and $A_{n+1} := A \cup \left(\verb+prodperm+\ A\ A_n\right)$ for $n \geq 0$.
Since $G \leq S_h$ is a finite group we get $A_{n+1} = A_n$ for some $n \in \IN$ and this process
terminates. The function \verb+closure'+ implements the recursion and
\verb+closureS+ is the interface for sets of permutations and 
\verb+closureL+ is the interface for lists of permutations. 
\begin{code}
closure' :: SPerm -> SPerm -> SPerm
closure' p0 ps = if ps == px then ps else closure' p0 px
  where
    px = Set.union ps (prodperm p0 ps) 

closureS :: SPerm -> SPerm
closureS ps = closure' ps ps

closureL :: [Perm] -> [Perm]
closureL = Set.toAscList . closureS . Set.fromList
\end{code}

\begin{comment}
\begin{code}
gengroups' :: Set SPerm -> Set SPerm -> Set SPerm
gengroups' gs0 gs = if Set.null gs0 then gs else gsn
  where
    (g0,gr) = Set.deleteFindMin gs0
    gsn = gengroups' gr $! Set.fold addclosure gs gs
    addclosure :: SPerm -> Set SPerm -> Set SPerm
    addclosure s gss = if Set.size (s `Set.intersection` g0) > 1 
      then gss else Set.insert (closureS $ Set.union s g0) gss

primes :: [Int]
primes = sieve [2..]
  where sieve (p:xs) = p : sieve [x | x<-xs, x `mod` p /= 0]

primeordercyclicgroups :: [Perm] -> [SPerm]
primeordercyclicgroups ps = nub . filter (\s -> length (primefactors s) == 1) 
                                . map (closureS . Set.singleton) $ ps
  where
    primefactors s = filter (`divides` Set.size s) . takeWhile (<= Set.size s) $ primes 
    p `divides` n = n `mod` p == 0 

lift1SPerm :: SPerm -> SPerm
lift1SPerm = Set.map (\p -> p ++ [length(p)])

allgroups :: Int -> Set SPerm
allgroups h = if h == 1 then Set.singleton e else Set.insert e . gengroups' cycgr 
                . Set.map lift1SPerm . allgroups $ h-1
  where
    e = closureS . Set.singleton $ [0..h-1]
    ps = allperm h
    psnotfix = filter (\p -> p !! (h-1) /= (h-1)) ps
    cycgr = Set.fromList $ primeordercyclicgroups psnotfix
\end{code}
\end{comment}

The following code generates all subgroups of $S_h$. The algorithm is not very efficient 
but works. Starting with $P_0 := \{ X_1,\dots,X_l \} := \{ [\{x\}] \mid x \in S_h \}$ we generate
$P_{n} := P_{n-1} \cup \{ [X_n \cup X] \mid X \in P_{n-1} \}$ for $n \in \{1,\dots,l\}$.
\begin{code} % this is not the actual code used to generate the groups; it could be used but is very slow 
allgroups' :: [SPerm] -> Set SPerm -> Set SPerm
allgroups' []     pss = pss
allgroups' (x:xs) pss = allgroups' xs $! Set.union pss pn
  where
    pn = Set.map (closureS . Set.union x) . Set.filter (not . Set.isSubsetOf x) $ pss

allgroups :: Int -> Set SPerm
allgroups h = allgroups' ps pss
  where
    ps = nub . map (closureS . Set.singleton) $ allperm h
    pss = Set.fromList $! ps
\end{code}

The following functions deal with partitions representing the equivalence relations $\varepsilon$ used in the
description of coherent relations in the terms $\delta_\varepsilon$.
Since we can permute the rows of relations without affecting the corresponding clone we sort 
each partition such that the longest equivalence class comes first. The function
\verb+alllinpart+ generates all such partitions on $E_h$.

The functions \verb+vecinpart+ and \verb+vecinmultipart+ test if a given vector $v$ belongs to 
$\delta_{\verb+pa+}$ and $\bigcup_{\verb+pa+ \in \verb+pas+} \delta_{\verb+pa+}$, respectively.
\begin{code}
lencomp :: [a] -> [a] -> Ordering
lencomp a b = compare (length a) (length b)

longest :: [ [a] ] -> [a]
longest = maximumBy lencomp

type Partition = [[Int]] 

alllinpart' :: [Int] -> [Partition]
alllinpart' [i] = [ [[i]] ]
alllinpart' (i:is) = concat [ [[i]:(c0:cs),(i:c0):cs] | (c0:cs) <- alllinpart' is ]

alllinpart :: Int -> [Partition]
alllinpart h = sort 
               . filter (\p ->  (reverse . sortBy lencomp . reverse $ p) == p  )
               . alllinpart' $ [0..h-1]

vecinpart :: Partition -> [Int] -> Bool
vecinpart [] _ = False
vecinpart pa v = all (\(i1:is) -> and [ (v !! i1) == (v !! i2) | i2 <- is ] ) pa

vecinmultipart :: [Partition] -> [Int] -> Bool
vecinmultipart pas v = any (\pa -> vecinpart pa v) pas
\end{code}

The function \verb+filtergroupprevpart gs pa+ gets all groups \verb+g+ from \verb+gs+ such that 
all permutations $\pi \in \verb+g+$ preserve the partition \verb+pa+, 
i.e., $\delta_{\verb+pa+}\B{\pi} = \delta_{\verb+pa+}$.
Additionally \verb+map Set.toAscList+ enforces that the identity permutation $[0,1,\dots,h-1]$ 
is the first element in the list presentation of a group.
\begin{code}
applypermtopart :: Perm -> Partition -> Partition
applypermtopart pi pa = sort . map (sort . map (\x -> pi !! x)) $ pa

piprevpart :: Perm -> Partition -> Bool
piprevpart pi pa = (pa == (applypermtopart pi pa))

groupprevpart :: SPerm -> Partition -> Bool
groupprevpart g pa = all (\pi -> piprevpart pi pa) $ Set.toList g

filtergroupprevpart :: Set SPerm -> Partition -> [[Perm]]
filtergroupprevpart gs pa = sortBy lencomp . map Set.toAscList . Set.toList . Set.filter (\g -> groupprevpart g pa) $ gs
\end{code}

The function \verb+applyperm+ takes a permutation \verb+p+ and applies it to a vector \verb+x+.
\begin{code}
applyperm :: Perm -> [Int] -> [Int]
applyperm p x = map (x !!) p

applygroup :: [Perm] -> [Int] -> [ [Int] ]
applygroup g x = [ applyperm p x | p <- g ]

listmodgroup :: [Perm] -> [ [Int] ] -> [ [Int] ]
listmodgroup _ [] = []
listmodgroup _ [a] = [a]
listmodgroup g (a:as) = a : (listmodgroup g (as \\ (applygroup g a)))
\end{code}

We need to have all tuples from $E_k^h \setminus \iota_k^h$ which are generated by \verb+allareflvec+.
The tuples in $\iota_k^h$ are handled differently as the Theorem of Haddad and Rosenberg suggests.
\begin{code}
allareflvec' :: [Int] -> Int -> [ [Int] ]
allareflvec' _  0 = [ [] ] 
allareflvec' [] _ = [    ] 
allareflvec' xs h = concat [ map (x:) (allareflvec' (xs \\ [x]) (h-1))| x <- xs ]

allareflvec :: Int -> Int -> [ [Int] ]
allareflvec k h = sort $ allareflvec' [0..k-1] h
\end{code}

The definition of coherent relations includes the use of relational homomorphisms from 
$\varrho$ to $M(\varrho)$. Candidates for these relational homomorphisms are all 
surjective functions $\varphi$ from $E_k$ to $E_h$ since there shall be an $v \in \varrho$
with $\varphi(v) = \eta_h$. The function \verb+allsurjfunc+ generates all of them.
\begin{code}
allfunc' :: [Int] -> [Int] -> [ [Int] ]
allfunc' [_] ys = [ [y] | y <- ys ]
allfunc' (_:xs) ys = [ y:r | y <- ys, r <- allfunc' xs ys ]

allsurjfunc' :: Int -> Int -> [ [Int] ]
allsurjfunc' k h = filter (\f -> h == length (nub f)) (allfunc' [0..k-1] [0..h-1])

transfunc :: Int -> [Int] -> (Int -> Int)
transfunc _ fi x = fi !! x

allsurjfunc :: Int -> Int -> [ Int -> Int ]
allsurjfunc k h = map (transfunc k) ( allsurjfunc' k h )
\end{code}

The function \verb+testhom+ checks if a given function \verb+f+ maps the tuple \verb+r+
into the model of $\varrho$ and additionally returns information if the tuple is not in $\iota_h^h$.
If the last bit is true and \verb+f+ is a relational homomorphism, then there is some \verb+f'+ with
$\verb+f'+(\verb+r+) = \eta_h$. The relation $\varrho = \verb+sigma+ \cup \delta(\varrho)$
is encoded in $\verb+gsigma+ = G_\sigma$ and the function \verb+patest+ representing
$\delta(\varrho)$. By the way we represent relations and permutations 
we have $\verb+gsigma+ = M(\varrho) \setminus \iota_k^h$. Additionally \verb+patest+ does not
change when switching to the model of $\varrho$. 

The function \verb+ishomwitheta+ then checks if \verb+f+ is a relational homomorphism such that
some $\verb+r+ \in \varrho$ exists with $\verb+f+(\verb+r+) = \eta_h$. \verb+existrelhom+ tests 
all surjective functions if they are relational homomorphisms, and \verb+iscoherent+ does this test for
all subsets of $\sigma(\varrho)$, i.e., checking if $\varrho$ is coherent, since $\varrho$ is constructed such
that the condition on $G_\sigma$ is fullfilled. 
\begin{code}
testhom :: [Perm] -> ([Int] -> Bool) -> (Int -> Int) -> [Int] -> (Bool, Maybe [Int])
testhom gsigma patest f r = case ( fr `elem` msigma, patest fr) of
    (True,_) -> (True, Just r) 
    (_,True) -> (True, Nothing)
    (_,_)    -> (False, Nothing)
  where
    fr = map f r
    msigma = gsigma -- model of sigma equals G_sigma because of representation of gsigma

ishomwitheta :: [ [Int] ] -> [Perm] -> ([Int] -> Bool) -> (Int -> Int) -> Maybe [ [Int] ]
ishomwitheta sigma gsigma patest f = case (and inMrho, catMaybes mappedtoeta) of
    (True, []) -> Nothing
    (True, ls) -> Just ls
    (False, _) -> Nothing
  where
    (inMrho, mappedtoeta) = unzip . map (testhom gsigma patest f) $ sigma

existrelhom :: [ Int -> Int ] -> [ [Int] ] -> ([Int] -> Bool) -> [Perm] -> Maybe ([ [Int] ])
existrelhom ksurjh sigma patest gsigma = msum . map (ishomwitheta sigma gsigma patest) $ ksurjh

iscoherent :: [ Int -> Int ] -> [ [Int] ] -> ([Int] -> Bool) -> [Perm] -> Bool
iscoherent _ [] _ _ = True
iscoherent _ [_] _ _ = True
iscoherent ksurjh sigma patest g = value mrs
  where
    mrs = existrelhom ksurjh sigma patest g -- vectors mapped into Model(sigma)
    value :: Maybe ([ [Int] ]) -> Bool
    value Nothing = False
    value (Just rs) = iscoherent ksurjh (sigma \\ rs) patest g
\end{code}

The type \verb+RelInfo+ stores all kinds of information and is used and updated in the recursive search.
\verb+rel+ stores only the tuples such that 
\[ \verb+fullrel+ = \sigma(\varrho) = \bigcup_{\pi \in \verb+gsigma+} \verb+rel+\B{\pi} \]
with $\verb+gsigma+ = G_\sigma$. This is used in the function \verb+iscoherent+ to reduce the number of checks at that 
point.

The entry \verb+pa+ takes the partition corresponding to the equivalence
relation $\varepsilon$ for areflexive and quasi-diagonal relations. The function \verb+patest+ is the test
used in \verb+iscoherent+ and \verb+parep+ is the entry which will be displayed in the list of relations representing
$\delta(\varrho)$.

The function \verb+isgood+ is used to decide if the given relation is non-trivial, i.e., check that 
$\pPOL_k \varrho \neq \wP_k$.

The entries \verb+ksurjh+, \verb+allbij+, \verb+k+ and \verb+h+ are parameters of the relation and
added to this record to avoid recalculation of them at every step, i.e., reducing computations.

Let $G'$ be the set of permutations that preserve $\delta(\varrho)$. Then $\verb+gsigma+ \leq G'$ and there
are cosets of $\verb+gsigma+$ in $G'$. The list \verb+coreps+ contains the coset representatives of the left
cosets of \verb+gsigma+ in $G'$; see also \verb+cosetreps+. It is useful to use \verb+coreps+ instead of $G'$
since $\varrho\B{\pi'} = \varrho\B{\pi}$ if $\pi' \in \pi \verb+gsigma+$, i.e., if $pi'$ and $\pi$ are in the same
coset of \verb+gsigma+.

Let $\verb+allbij+ =: [ f_1,\dots,f_l ]$ and $\verb+coreps+ =: [ \pi_1, \dots, \pi_p ]$.
The entry \verb+xclass+ is a list of lists of relations, where
\begin{itemize}
  \item $\verb+xclass+ = [ L_1,\dots,L_l ]$, indexed by the functions in \verb+allbij+;
  \item $L_i = [ r_{i1}, \dots, r_{ip} ]$, where $i \in \{1,\dots,l\}$;
  \item 
    $r_{ij} = \{ f_i(v) \mid v \in \verb+fullrel+\B{\pi_j} \}$, i.e., the relation generated
    from $\verb+fullrel+ = \sigma(\varrho)$ by applying the permutation $\pi_j \in \verb+coreps+$ 
    and then mapping by the bijection $f_i \in \verb+allbij+$.
\end{itemize}
Let $r,r' \in L_i$ for some $i \in \{1,\dots,l\}$. Then $\pPOL_k r = \pPOL_k r'$ since we can write 
\[ r' = f_i(\verb+fullrel+\B{\pi'}) = f_i(\verb+fullrel+)\B{\pi'} = 
f_i(\verb+fullrel+)\B{\pi\pi^{-1}\pi'} = r\B{\pi^{-1}\pi'}. \]
Thus every $L_i$ describes exactly one maximal partial clone.
Since the first element in \verb+coreps+ is the identity permutation (see \verb+filtergroupprevpart+)
the first element in $L_i$ is $f_i(\verb+fullrel+)$. Therefore \verb+quasiclass+ generates the quasi-relation-class
$\qrelclass(\varrho)$. The relation-class \verb+relclass+ yields a list of the minimal representatives for a given
maximal clone. These coincide if and only if the corresponding maximal clones are equal. Thus we can count 
the maximal partial clones with a given minimal relation by counting the different relations in \verb+relclass+.  
\begin{code}
data RelInfo a = RelInfo {
  rel :: [[Int]],
  fullrel :: [[Int]],
  gsigma :: [Perm],
  pa :: [[Int]],
  patest :: [Int] -> Bool,
  parep :: [a], 
  coreps :: [Perm],
  xclass :: [[ [[Int]] ]],
  ksurjh :: [ Int -> Int ],
  allbij :: [ Int -> Int ],
  k :: Int,
  h :: Int,
  isgood :: RelInfo a -> Bool
}

quasiclass :: RelInfo a -> [ [[Int]] ]
quasiclass = map head . xclass

relclass :: RelInfo a -> [ [[Int]] ]
relclass = map minimum . xclass

checkQuasiMin :: RelInfo a -> Bool
checkQuasiMin ri = all ( >= (fullrel ri)) (quasiclass ri)

checkMin :: RelInfo a -> Bool
checkMin ri = all ( >= (fullrel ri)) (relclass ri)
\end{code}

The function \verb+preInitRelInfo+ initializes a \verb+RelInfo+ structure for an empty relation with 
given \verb+gsigma+ and  
\begin{code}
cosetreps :: [Perm] -> [Perm] -> [Perm]
cosetreps g subg = nub . map (\pi -> minimum . map (multperm pi) $ subg) $ g  

preInitRelInfo :: Int -> Int -> [ Int -> Int ] -> [ Int -> Int ] -> (b,[Perm],[Perm],[ [Int] ]) -> RelInfo a
preInitRelInfo k' h' allbij' allksurjh' (_,gofpa',g',_) = RelInfo { 
    rel = [], fullrel = [],
    gsigma = g',  
    coreps = coreps',
    xclass = [ [ [] | pi <- coreps' ] | f <- allbij' ],
    ksurjh = allksurjh', allbij = allbij', k = k', h = h',
    pa = [], patest = \_ -> True, parep = [], isgood = \_ -> True
    }
  where
    coreps' = cosetreps gofpa' g'
\end{code}

The representative of the relation class and the size of the relation class are used
for the output of coherent relations with \verb+computeRelRep+ and \verb+printGetSize+. 
\verb+printGetSize+ only prints output and returns non-zero class size, if the relation
is really a coherent relation (i.e., \verb+isgood+ is true) and the relation is a minimal
relation. The check \verb+iscoherent+ done before, does not check these criteria.

The output 
consists of a tuple $(\verb+length+, \verb+parep+, \verb+rel+)$ where
\begin{itemize}
  \item \verb+length+ is the size of relation class;
  \item \verb+parep+ is 
    \begin{itemize}
      \item \verb+[]+ if the relation is areflexive,
      \item the partition \verb+pa+ as a list representing $\varepsilon$ if the relation
        $\varrho = \sigma \cup \delta_\varepsilon$
        is quasi-diagonal,
      \item \verb+iota_k^h+ if the relation is totally symmetric and totally reflexive,
      \item \verb+rho_1+ or \verb+rho_2+ if $\linrel_1 \subseteq \varrho$ or $\linrel_2 \subseteq \varrho$
        for the relation $\varrho$;
    \end{itemize}
  \item \verb+rel+ are the tuples of $\sigma(\varrho)$ for the relation class representative $\varrho$.
\end{itemize}
The function \verb+checkCoherent+
lifts the \verb+iscoherent+ function to the \verb+RelInfo+ type.
\begin{code}
checkCoherent :: RelInfo a -> Bool
checkCoherent ri = iscoherent (ksurjh ri) (rel ri) (patest ri) (gsigma ri)

computeRelRep :: Show a => RelInfo a -> (Int, [a], [[Int]])
computeRelRep ri = (length rc, parep ri, head rc)
  where
    rc = map head . group . sort . relclass $ ri

printGetSize :: Show a => RelInfo a -> IO Integer
printGetSize ri = 
    if isgood ri ri && checkMin ri
    then (print $! relRep) >> return sizeRep
    else return 0
  where
    relRep = computeRelRep ri
    sizeRep = toInteger $ (\(s,_,_) -> s) relRep
\end{code}

The function \verb+insertNew+ inserts an element \verb+v+ into an ordered list if it is not included yet.

The function \verb+relInfoAddVec v ri+ updates the structure \verb+ri+ by adding \verb+v+ to \verb+rel+,
\verb+fullrel+ and \verb+xclass+. The term $\verb+gv+ = \{ v\B{\pi} \mid \pi \in \verb+gsigma+ \}$ is used
for $\verb+fr+$ and $\verb+xclass+$. The new \verb+xclass+ is generated from the old one by
by replacing $r_{ij}$ by $r_{ij} \cup \{ w\B{\pi} \mid \pi \in \verb+gsigma+ \}$ where $w := f_i(v)\B{\pi_j}$.
\begin{code}
insertNew :: Ord a => a -> [a] -> [a]
insertNew v []     = [v]
insertNew v (r:rs) = case compare v r of
  LT -> v:r:rs
  EQ -> r:rs
  GT -> r:(insertNew v rs)

relInfoAddVec :: [Int] -> RelInfo a -> RelInfo a
relInfoAddVec v ri = ri { rel = rn, fullrel = fr, xclass = xc }
  where
    rn = insertNew v (rel ri)
    fr = foldr insertNew (fullrel ri) gv
    xc = map yc $ zip (xclass ri) (allbij ri)
    yc (c,f) = map (\(d,pi) -> foldr (insertNew . map f . applyperm pi) d gv) $ zip c (coreps ri)
    gv = map (\pi -> applyperm pi v) (gsigma ri)
\end{code}
        
\verb+coherentSubsetsWith+ is the central function traversing the tree of coherent relations for given 
initialisation with the root at $\sigma(\varrho) =  \emptyset$.
If the current relation is not quasi-minimal then the rest of this subtree can be omitted. Thus the function returns
$0$. The sub-nodes are created in \verb+aav'+ and then checked to be coherent. If the relation is not coherent for some
\verb+v+ then we can not add \verb+v+ anywhere in this subtree. Therefore \verb+aavn+ are the tuples from
\verb+aav+ which are allowed in the subtrees. Since we only want to generate a given relation once
we allow in the subtree to \verb+v+ only tuples greater than \verb+v+, see \verb+aavngtv+ in \verb+recurse+.

The \verb+foldM+ line traverses all pairs $(\verb+v+,\verb+rj+)$ and sums up the number of maximal partial clones.
\begin{code}
coherentSubsetsWith :: Show a => RelInfo a -> [ [Int] ] -> IO Integer
coherentSubsetsWith ri aav = 
    if not (checkQuasiMin ri) 
    then return 0 
    else do
      size <- printGetSize ri
      let aav' = map (\v -> (v,relInfoAddVec v ri)) (aav \\ (fullrel ri))
          vrjs = filter (\(_,rj) -> checkCoherent rj) $ aav'
          aavn = map fst vrjs
      foldM (recurse aavn) size vrjs

recurse :: Show a => [[Int]] -> Integer -> ([Int], RelInfo a) -> IO Integer
recurse aavn s (v,rj) = coherentSubsetsWith rj aavngtv >>= strictAddToS
  where
    strictAddToS x = return $! (s + x)
    aavngtv = filter (>v) aavn
\end{code}

The function \verb+isEkh+ is used for areflexive and quasi-diagonal relations to check if $\varrho = E_k^h$.
There are the following cases
\begin{itemize}
  \item $h = 1$; then $\varrho = E_k^h$ if $|\sigma(\varrho)| = k$;
  \item $h = 2$; then $\varrho = E_k^h$ if $\delta(\varrho) = \delta_{\{0,1\}}\A{2} = \iota_k^2$ and
    $|\sigma(\varrho)| = \binom{k}{2} = k \cdot (k-1)$;
  \item for every $h > 2$ we have $\varrho \neq E_k^h$.
\end{itemize}
A relation is a valid candidate for a maximal partial clone if $\varrho \neq E_k^h$ and $\sigma(\varrho) \neq \emptyset$
(see \verb+isgood+ definition in \verb+findcoherent+). 

\verb+findcoherent+ takes a tuple \verb+i+ consisting of the partition \verb+pa'+, the group preserving \verb+pa'+,
the group \verb+gsigma'+ and a list of tuples \verb+av+.
The tuples in \verb+av+ represent all tuples in $\sigma(E_k^h)$,
but reduced with respect to \verb+gsigma'+. Let 
$\{ V_1,\dots,V_n \} := \{ \{ v\B{\pi} \mid \pi \in \verb+gsigma'+ \} \mid v \in \sigma(E_k^h) \}$.
Then $\verb+av+ = \{v_1,\dots,v_n \}$ where $v_i \in V_i$ for all $i \in \{1,\dots,n\}$. This is done in
\verb+getpgvs+ or more specifically in \verb+listmodgroup+. 

\verb+getpgvs+ takes all groups preserving \verb+pa+ (denoted by \verb+ga0+) and emits all tuples 
for \verb+findcoherent+. \verb+printcoherent+ just sums the number of different maximal partial clones 
over all different partitions for a given $h$
and returns this sum.
\begin{code}
isEkh :: Int -> Int -> Partition -> [ [Int] ] -> Bool
isEkh k 1 _  rel         = length (rel) == k
isEkh k 2 pa rel         = pa == [[0,1]] && length (rel) == k*(k-1)
isEkh _ h _  _   | h > 2 = False

findcoherent :: Int -> Int -> [ Int -> Int ] -> [ Int -> Int ] -> (Partition,[Perm],[Perm],[ [Int] ]) -> IO Integer
findcoherent k' h' allbij' allksurjh' i@(pa',gofpa',gsigma',av) = do
  let ri0 = (preInitRelInfo k' h' allbij' allksurjh' i) {
        pa = pa', patest = vecinpart pa', parep = pa', 
        isgood = \ri -> not (isEkh (k ri) (h ri) (pa ri) (fullrel ri)) && fullrel ri /= []
        }
  coherentSubsetsWith ri0 av >>= (\x -> return $! x)

getpgvs :: Set SPerm -> [[Int]] -> Partition -> [ (Partition,[Perm],[Perm],[ [Int] ]) ] 
getpgvs gs aav pa = map (\g -> (pa,longest ga0,g,listmodgroup g aav)) ga0 
  where
    ga0 = filtergroupprevpart gs pa

printcoherent :: Int -> Int -> IO (Integer)
printcoherent k h = do
  let pts = []:((alllinpart h) \\ [[[x] | x <- [0..h-1]]])
      allbij = allsurjfunc k k
      allksurjh = allsurjfunc k h
      gs = allgroups h
      aav = allareflvec k h
      pgvs = concatMap (getpgvs gs aav) pts

  mapM (findcoherent k h allbij allksurjh) pgvs >>= (return . sum)
\end{code}

The totally reflexive, totally symmetric relation $\varrho$
is not $E_k^h$ if $|\sigma(\varrho)| \neq |\sigma(E_k^h)| = k\cdot k-1 \cdot \ldots \cdot k-h+1$, the number of $h$-tuples
with pairwise different entries. Thus \verb+isgood+ is defined that way. Any tuple not in $\sigma(E_k^h)$ is in
$\iota_k^h$ and thus in $\varrho$. Therefore \verb+patest+ is always true. 

By the Theorem of Haddad and Rosenberg $\verb+gsigma+ = \verb+gofpa+ = S_h$.
\begin{code}
prod :: [Int] -> Int
prod xs = foldl (*) 1 xs

findIota :: Int -> Int -> [ Int -> Int ] -> [ Int -> Int ] -> (a,[Perm],[Perm],[ [Int] ]) -> IO Integer
findIota k' h' allbij' allksurjh' i@(_,_,_,av) = do
  let ri0 = (preInitRelInfo k' h' allbij' allksurjh' i) {
        pa = [], patest = (\_ -> True), 
        parep = "iota_" ++ (show k') ++ "^" ++ (show h'),
        isgood = \ri -> length (fullrel ri) /= prod [k'-h'+1..k']
        }
  coherentSubsetsWith ri0 av >>= (\x -> return $! x)

printIotakh :: Int -> Int -> IO (Integer)
printIotakh k h = do
  let g = allperm h
      aav = allareflvec k h                
      pg0 = ([],g,g,listmodgroup g aav) 
      allbij = allsurjfunc k k
      allksurjh = allsurjfunc k h
  findIota k h allbij allksurjh pg0 
\end{code}

The quaternary relations with $\linrel_i \subseteq \varrho$ for some $i \in \{1,2\}$ are generated next
where \verb+rhoiparam+ gives the group $G_\sigma$, the set of partitions of $\linrel_i$, and the output string. 
Since any such relation (which is generated in this program) is a valid coherent relation, we have \verb+isgood+ always
be true.
\begin{code}
rhoiparam :: Int -> ([Perm],[Partition], String)
rhoiparam 1 = (closureL [[2,0,3,1],[3,1,2,0]], [[[0,1],[2,3]],[[0,2],[1,3]]]              , "rho_1")
rhoiparam 2 = (allperm 4                     , [[[0,1],[2,3]],[[0,2],[1,3]],[[0,3],[1,2]]], "rho_2")

printRhoi :: Int -> Int -> IO (Integer)
printRhoi k i = do
  let h = 4
      (g, pas, delta) = rhoiparam i
      aav = allareflvec k h                
      pg0 = (pas,g,g,listmodgroup g aav) 
      allbij = allsurjfunc k k
      allksurjh = allsurjfunc k h
  findRhoi k h allbij allksurjh pg0 delta

findRhoi :: Int -> Int -> [ Int -> Int ] -> [ Int -> Int ] -> ([Partition],[Perm],[Perm],[ [Int] ]) -> String -> IO Integer
findRhoi k' h' allbij' allksurjh' i@(pas',_,_,av) parep' = do
  let ri0 = (preInitRelInfo k' h' allbij' allksurjh' i) {
        pa = [], patest = vecinmultipart pas', parep = parep', 
        isgood = \_ -> True
        }
  coherentSubsetsWith ri0 av >>= (\x -> return $! x)
\end{code}

Print the list of coherent relations in the following order:
\begin{itemize}
  \item the areflexive and quasi-diagonal relations,
  \item the totally symmetric, totally reflexive relations with $h \geq 3$,
  \item the quaternary relations with $\linrel_i \subseteq \varrho$ for some $i \in \{1,2\}$,
  \item $P_k \cup C_\notdef$.
\end{itemize} 
In the end the numbers are printed grouped as above, and a grand total.

The $k$ for which to generate the list of coherent relations is given as a command line parameter.
\begin{code}
runPrint :: Int -> IO ()
runPrint k = do
  sizes <- mapM (printcoherent k) [1..k]
  sizes2 <- mapM (printIotakh k) [3..k]
  sizes3 <- mapM (printRhoi k) [1,2]
  putStrLn ("(1,P_" ++ (show k) ++ " \\cup C_\\infty)")
  sizes4 <- return [1] -- P_k \cup C_\infty
  print (sizes, sizes2, sizes3, sizes4, sum (sizes ++ sizes2 ++ sizes3 ++ sizes4))

main = do
  x <- getArgs
  let k :: Int = read (head x)
  runPrint k
  return 0
\end{code}