3
\$\begingroup\$

For a while now I've been meaning to post some of my Haskell code here so that someone can tell me what parts of the language/base library I've been completely overlooking. This is the first thing I've brought into a "working/finished" state that isn't a PE problem, and unfortunately it is a tad big.

It is rather mathematical in nature, and contains multiple nontrivial alglorithms which could almost certainly be improved... but I'm hoping to focus less on the algorithms and more on the minutae. I.e. how is my usage of Haskell?

What does it do? Some symmetry analysis for materials research. If you run it, it will print out all of the possible supercell shapes containing 8 16 atoms which are unique under symmetry for a diamond cubic crystal.

Last minute Hawthorne-effect edits aside, the majority of this code was written with the expectation that I would be the only person to ever look at it. I can only hope that at least one soul is brave enough to continue after reading that statement...

Directory structure

$ tree -P '*.hs'
.
├── Main.hs
└── My
    ├── Common.hs
    ├── GroupTheory.hs
    ├── IntegerRref.hs
    └── Matrix.hs

Main.hs

Mostly functions specific to the problem. Here you'll see from all the commented code that my workflow for debugging largely centers around modifying main and recompiling.

{-# OPTIONS_GHC -fno-ignore-asserts #-}
{-# LANGUAGE BangPatterns #-}
import qualified Math.NumberTheory.Primes.Factorisation as Factor -- package arithmoi

import qualified Data.Set as Set -- package containers
import Data.Set (Set)

import qualified Data.List as List
import Control.Exception
import Text.Printf
import Debug.Trace

import My.Matrix
import My.GroupTheory
import My.IntegerRref
import My.Common(decorate)

-------------------------------------
-- Factorization

-- positive, ordered tuples (a,b) such that `a*b == x`
factorPairs :: Integer -> [(Integer, Integer)]
factorPairs x = decorate (x `div`) $ Set.toList $ Factor.divisors x

-- ordered n-tuplets of positive factors `fs` such that `product fs == x`
factorTuplets :: Int -> Integer -> [[Integer]]
factorTuplets n x = (assert (x>0)) $ -- true for expected inputs in this program
    case n `compare` 1 of
        LT -> error "factorTuplets: n < 1"
        EQ -> [[x]]
        GT -> concatMap forPair (factorPairs x)
    where forPair (a,b) = map (a:) $ factorTuplets (pred n) b

-------------------------------------
-- Groups and sets

universeForDiagonal :: Vec -> [Mat]
universeForDiagonal [] = [[]]
universeForDiagonal (a:bcs) = assert (a>0) $ do
    let left = a:(fmap (const 0) bcs)
    inner <- universeForDiagonal bcs
    top <- sequence $ fmap (\b -> [0..b-1]) bcs
    return $ (prependCol left.prependRow top) inner

universeForVolume :: Integer -> [Mat]
universeForVolume vol = concatMap universeForDiagonal (factorTuplets 3 vol)

-------------------------------------
-- Supercell symmetry group

scGenerators :: [Mat]
scGenerators =
    -- Twofold rotations
    [[0,1,0],[1,0,0],[-1,-1,-1]]:
    [[-1,0,0],[0,-1,0],[1,1,1]]:
    [[1,0,0],[0,0,1],[-1,-1,-1]]:
    -- Threefold rotation
    [[0,0,1],[1,0,0],[0,1,0]]:
    []

scMul :: Mat -> Mat -> Mat
scMul = mulMatMat

-- sc group matrices are written to operate on a matrix whose columns are the sc vecs,
-- but we have them in columns;
scAction :: Mat -> Mat -> Mat
scAction g x = integerRref $ mulMatMat x (List.transpose g)

scGroup :: Set Mat
scGroup = generateGroup scMul scGenerators

-------------------------------------
-- Random nonsense

assertEq :: (Eq a, Show a) => a -> a -> b -> b
assertEq expected actual
    | expected /= actual = error (printf "\n Expected: %s \n Actual: %s" (show expected) (show actual))
    | otherwise = id

tests :: ()
tests = id
    .(assertEq 155 $ length $ universeForVolume 8)
    .(assertEq [[-10,5],[2,3]] $ opAddMultiple 3 1 0 $ [[-16,-4],[2,3]])
    .(assertEq 77 $ length $ equivalenceClasses scAction (Set.toList scGroup) (universeForVolume 20))
    $()

main :: IO ()
main = do
    let !_ = tests
--  mapM_ print $ factorTuplets 3 10
--  mapM_ print $ upperTriangleFromDiags [[1,2,3,4],[5,6,7],[8,9],[1]]
--  mapM_ print $ concat $ universeForDiagonal [1,2,3]
--  mapM_ print $ opAddMultiple (-2) 2 1 $ [[1,4,7,6],[2,3,1,6],[7,7,7,2],[1,2,1,1]]
--  mapM_ print $ fmap (`getDiag` [[1,4,7,6],[2,3,1,6],[7,7,7,2],[1,2,1,1]]) $ [-3..3]
--  mapM_ print $ fmap (`getDiag` [[1,4,7,6],[2,3,1,6],[7,7,7,2],[1,2,1,1]]) $ [-3..3]
--  mapM_ print $ integerRref [[1,1,1],[5,2,2],[3,3,4]]
--  mapM_ print $ integerRref [[5,2,2],[1,1,1],[3,3,4]]
--  print $ Set.size scGroup
--  mapM_ print $ equivalenceClasses scAction (Set.toList scGroup) (universeForVolume 77)
--  mapM_ print $ integerRref [[0,1,0],[0,0,1],[4,0,0]]
    mapM_ print $ map Set.findMin $ equivalenceClasses scAction (Set.toList scGroup) (universeForVolume 8)

My/GroupTheory.hs

A small number of algorithms related to group theory.

module My.GroupTheory (
    generateGroup,
    equivalenceClasses,
    checkedEquivalenceClasses,
    ) where

import qualified Data.Set as Set -- package containers
import Data.Set (Set)

setUnionAll :: Ord a => [Set a] -> Set a
setUnionAll = foldl Set.union Set.empty

-- Generate all elements of a finite group from a generating subset.
generateGroup :: Ord g => (g -> g -> g) -> [g] -> Set g
generateGroup _   [] = error "generateGroup: empty group"
generateGroup mul generators =
    loop (Set.fromList generators) Set.empty
        where
        loop recent output
            | Set.null recent = output
            | otherwise = loop new (output `Set.union` new)
                where new =
                    (`Set.difference` output) $ Set.fromList $
                    fmap (\[a,b] -> a `mul` b) $
                    sequence [Set.toList recent, generators]

equivalenceClasses :: Ord x => (g -> x -> x) -> [g] -> [x] -> [Set x]
equivalenceClasses _ [] _  = error "equivalenceClasses: empty group"
equivalenceClasses action group xs =
    loop (Set.fromList xs)
    where loop remaining
        | Set.null remaining = []
        | otherwise = newClass:loop (remaining `Set.difference` newClass)
        where newClass = Set.fromList $ fmap (`action` (Set.findMin remaining)) group

checkedEquivalenceClasses :: Ord x => (g -> x -> x) -> [g] -> [x] -> [Set x]
checkedEquivalenceClasses action group xs = validate classes where
    classes = equivalenceClasses action group xs
    validate = validateClosed.validateDisjoint

    -- each element of xs should appear once and only once in the classes
    validateDisjoint = case totalOutCount `compare` uniqueOutCount of
        LT -> error "checkedEquivalenceClasses: internal error"
        GT -> error "checkedEquivalenceClasses: classes not disjoint"
        EQ -> id

    -- the union of the classes must equal xs
    validateClosed
        | (not . Set.null) (uniqueIn `Set.difference` uniqueOut)
            = error "checkedEquivalenceClasses: internal error in equivalenceClasses"
        | (not . Set.null) (uniqueOut `Set.difference` uniqueIn)
            = error "checkedEquivalenceClasses: group action not closed on xs"
        | otherwise = id

    uniqueIn = Set.fromList xs
    uniqueOut = setUnionAll classes
    uniqueOutCount = Set.size uniqueOut
    totalOutCount = sum (map Set.size classes)

-- findInverses :: Set g -> Map g g
-- findIdentity :: Set g -> g
-- validateClosure :: Set g -> Set g       -- O(n^2)
-- validateAssociativity :: Set g -> Set g   -- O(n^3)
-- validateGroup :: Set g -> Set g
-- validateAction :: Set g -> [x] -> ()

My/Matrix.hs

Defines operations on vectors and matrices, "implemented" as little more than lists.

module My.Matrix (
    Mat, Vec,
    prependCol, prependRow,
    deleteCol,
    overRow, overLowerRight,
    getCol, getDiag,
    upperTriangleFromDiags,
    height, width,
    innerProd, mulMatMat, mulMatVec,
    ) where

import qualified Data.List as List
import Control.Exception

import My.Common(listSet,deleteAt,zipWithExact)

type Mat = [[Integer]]
type Vec = [Integer]


prependCol :: Vec -> Mat -> Mat
prependRow :: Vec -> Mat -> Mat
prependCol col mat = zipWith (:) col mat
prependRow = (:)

deleteCol :: Int -> Mat -> Mat
deleteCol i = fmap (deleteAt i)

-- think "over (_!! i)", if _!! were some sort of Lens for lists
overRow :: Int -> (Vec -> Vec) -> Mat -> Mat
overRow i f mat = listSet i (f (mat!!i)) mat

-- applies a function to the lower right submatrix excluding `nDropped`
--  rows and columns.
overLowerRight :: Int -> (Mat -> Mat) -> Mat -> Mat
overLowerRight nDropped f mat =
    let (top, bottom) = List.splitAt nDropped mat in
    let (botLs, botRs) = unzip $ map (List.splitAt nDropped) bottom in
    top ++ (zipWith (++) botLs (f botRs))

getCol :: Int -> Mat -> Vec
getCol i = fmap (!!i)

getDiag :: Int -> Mat -> Vec
getDiag n rows = case n `compare` 0 of
    LT -> getDiag 0 (drop (-n) rows)
    GT -> getDiag 0 (fmap (drop n) rows)
    EQ -> fmap (\(i,row) -> row!!i) $ take w (zip [0,1..] rows)
        where w = (length.head) rows

upperTriangleFromDiags :: [Vec] -> Mat
upperTriangleFromDiags [] = []
upperTriangleFromDiags diags = assert checks $ topRow:otherRows where
    checks = length diags == length (head diags)
    topRow = map head diags
    otherRows = map (0:) $ upperTriangleFromDiags (map tail (init diags))

height :: Mat -> Int
height = length

width :: Mat -> Int
-- We do not explicitly store a width, so none can be determined from
-- a zero row matrix.  That said, I never plan to use one, so better
-- safe than sorry:
width [] = error "width: null matrix"
width mat = (length.head) mat

innerProd :: Vec -> Vec -> Integer
innerProd a b = sum $ zipWithExact (*) a b

mulMatMat :: Mat -> Mat -> Mat
mulMatMat a b
    | (width a) /= (height b) = error "mulMatMat: dimension"
    | otherwise = [[innerProd row col | col<-List.transpose b] | row <- a]

mulMatVec :: Mat -> Vec -> Vec
mulMatVec m v
    | (width m) /= (length v) = error "mulMatVec: dimension"
    | otherwise = map (`innerProd` v) m

My/IntegerRref.hs

The largest piece of code, implementing an algorithm based on various handwritten proofs. It is very closely related to the Hermite normal form of a matrix.

module My.IntegerRref (
    opAddRow, opSubRow, opNegate2, opAddMultiple,
    integerRref,
    validateIrref,
    ) where

import qualified Data.List as List
import Control.Exception

import My.Matrix
import My.Common(listSet,pDiv,compose)

---------------------------------------------------
-- This module implements an analogue to Reduced Row Echelon Form where
-- the only primitive symmetry operations are the following:
--   * adding one row into another, different row.
--   * subtracting one row from another, different row.
--
-- These operations preserve the determinant of a matrix.

opSubRow :: Int -> Int -> Mat -> Mat
opSubRow = opAddMultiple (-1)
opAddRow :: Int -> Int -> Mat -> Mat
opAddRow = opAddMultiple 1

-- We can't negate one row, but we can negate two:
opNegate2 :: Int -> Int -> Mat -> Mat
opNegate2 i1 i2 = overRow i1 (fmap negate).overRow i2 (fmap negate)

opAddMultiple :: Integer -> Int -> Int -> Mat -> Mat
opAddMultiple b src dest mat
    | dest == src = error "rowAdd: src == dest (not volume-conserving)"
    | otherwise = listSet dest newRow mat where
        newRow = zipWith (+) (mat!!dest) (fmap (b*) (mat!!src))

---------------------------------------------------

-- For a square, invertible, integer matrix, produces the unique matrix that is:
--  * Reachable from the original by a finite sequence of operations consisting
--    of adding an integer multiple of one row into another (different) row.
--  * Is upper triangular.
--  * Is entirely nonnegative with the SINGLE possible exception of the
--    lower-right most element.
--  * For each column, the absolute value of the element on the main diagonal
--    is strictly greater than all other values in the column.
integerRref :: Mat -> Mat
integerRref [] = []
integerRref mat = validateIrref $ compose operations $ partialRref mat
    where operations = do
        rPivot <- [0..(length mat)-1]
        rMod <- [0..rPivot-1]
        return $ (reduceRowModuloRow rPivot rPivot) rMod

-- This produces the correct values in the lower triangular part of the matrix
-- I.e. the main diagonal is positive with the possible exception of the last
--  element, and off-diagonals below the main diagonal are zero.
partialRref :: Mat -> Mat
partialRref = rec where
    rec [] = []
    rec mat = (fixRest.fixLeft) mat where
        -- change column into [g,g,g..], and then into [g,0,0..]
        fixLeft = compose $
            (reduceColumnToConstant 0):
            map (opSubRow 0) [1..(length mat)-1]
        -- recurse
        fixRest = overLowerRight 1 rec

-- reduces a column [a,b,c..] into [g,g,g..], where g = gcd (a,b,c..).
--  by performing primitive row ops.  As long as the original matrix
--  is larger than 1x1, g will be non-negative.
reduceColumnToConstant :: Int -> Mat -> Mat
reduceColumnToConstant _ [] = error "reduceColumnToConstant: empty mat" -- no columns!
reduceColumnToConstant _ [row] = [row] -- do NOT call ensureColumnNonNegative
reduceColumnToConstant c mat'@(_:_) = (verify.compose operations) mat' where
    verify mat = assert (all (==mat!!0!!c) (getCol c mat)) mat

    -- The column is made positive to simplify working with the gcd.
    -- Then two passes are made where each row is reduced with the top row by
    --   adding multiples of one to another until they share a value in column c.
    -- The first pass puts the final correct value in the first and last rows;
    --   the second pass propagates it to the other rows.
    operations = (ensureColumnNonNegative c):(onePass ++ onePass) where
        onePass = map (reduceTwoRows 0) [1..(length mat')-1]

    reduceTwoRows r1 r2 = verify.loop where
        -- like a bizarre variant of euclid's algorithm, where we're performing
        --  operations on entire rows, and we're not allowed to swap their order
        verify mat = assert (mat!!r1!!c == mat!!r2!!c) mat
        loop mat = case (mat!!r1!!c, mat!!r2!!c) of
            (0, 0) -> mat
            (0, _) -> (opAddRow r2) r1 mat
            (_, 0) -> (opAddRow r1) r2 mat
            (g, h) -> case (g `compare` h) of
                -- use lesser to reduce larger
                LT -> loop $ (reduceRowModuloRow c r1) r2 mat
                GT -> loop $ (reduceRowModuloRow c r2) r1 mat
                EQ -> mat

ensureColumnNonNegative :: Int -> Mat -> Mat
ensureColumnNonNegative c = verify.loop where
    verify mat = assert (all (0<=) (getCol c mat)) mat
    loop mat = case (positives, negatives, zeros) of

        -- Column is non-negative
        (_,  [],   _) -> mat

        -- Column is hopeless
        ([], [_], []) -> error "ensureColumnPositive: only one row and value is negative"

        -- At least one negative and one other: we can negate them
        ([], n1:n2:_, _) -> loop $ opNegate2 n1 n2 mat
        ([], [n1], z1:_) -> loop $ opNegate2 n1 z1 mat

        -- At least one positive p:
        --    use it to bring the negatives into the range 0..p-1
        (p:_, ns@(_:_), _) -> compose (fmap fixRow ns) mat
            where fixRow = reduceRowModuloRow c p

        where
        negatives = List.findIndices ( <0) column
        zeros     = List.findIndices (==0) column
        positives = List.findIndices ( >0) column
        column = getCol c mat

-- add a multiple of row `pivotRow` into row `modRow` that reduces
--  the element in column `col` into the range [0..|pivot| - 1].
-- The pivot must be nonzero.
reduceRowModuloRow :: Int -> Int -> Int -> Mat -> Mat
reduceRowModuloRow col pivotRow modRow mat =
    let moderand = mat!!modRow!!col in
    let pivot = mat!!pivotRow!!col in
    let d = moderand `pDiv` pivot in
    opAddMultiple (-d) pivotRow modRow mat

-- note: this assumes the matrix is square and invertible, so that
--       the pivots MUST be on the main diagonal
-- warning: this inspects every element, causing any deferred lazy
--          computations to occur conditionally based on whether
--          or not assertions are enabled.  I imagine there is the
--          potential for this to hide nasty memory bugs...
validateIrref :: Mat -> Mat
validateIrref mat = checkDims $ checkDiag $ checkZeros $ checkReduced $ mat
    where
    checkDims = assert $ all ((==length mat).length) mat

    -- the final pivot is the only element allowed to be negative in the
    --  entire matrix
    checkDiag = assert $ (all (0<) (init pivots)) && (0 /= last pivots)

    checkZeros = assert $ all test $ zip [0,1..] mat
        where test (i,row) =
            i == (length $ takeWhile (0==) $ row)

    checkReduced = assert $ all test $ List.zip3 [0,1..] columns pivots
        where test (i,col,pivot) =
            all (\x -> 0 <= x && x < abs pivot) $ (take i col)

    columns = List.transpose mat
    pivots = getDiag 0 mat

My/Common.hs

Various helper functions that I have deemed "missing" from the Haskell base library. (Your job is to tell me why they aren't!)

module My.Common (
    deleteAt, listSet,
    zipWithExact,
    decorate,
    compose,
    indices,
    windows,
    traceWith,
    pDivMod,pDiv,pMod,
    ) where

import Debug.Trace
import qualified Data.List as List

deleteAt :: Int -> [a] -> [a]
deleteAt _ [] = error "deleteAt: i >= length"
deleteAt i xs@(_:xt)
    | i < 0     = error "deleteAt: i < 0"
    | i == 0    = xt
    | otherwise = deleteAt (i-1) xs

listSet :: Int -> a -> [a] -> [a]
listSet _ _   []     = error "listSet: index (or empty list)"
listSet 0 new (_:xs) = new:xs
listSet n new (x:xs) = x:listSet (n-1) new xs

zipWithExact :: (a -> b -> c) -> [a] -> [b] -> [c]
zipWithExact f as' bs' = iter as' bs' where
    iter []     []     = []
    iter []     (_:_)  = error "zipWithExact: first list ended early"
    iter (_:_)  []     = error "zipWithExact: second list ended early"
    iter (a:as) (b:bs) = (f a b):iter as bs

decorate :: Functor f => (a -> b) -> f a -> f (a,b)
decorate f = fmap (\x -> (x, f x))

-- compose [f,g,h..] x ==  (..h.g.f) x
compose :: [(a -> a)] -> a -> a
compose = foldl (flip (.)) id

-- equivalent to [0..(length xs-1)] but possibly less painful?
indices :: [a] -> [Int]
indices xs = map fst $ zip [0..] xs

-- overlapping windows
windows :: Int -> [a] -> [[a]]
windows n xs = takeWhile ((n==).length) $ fmap (take n) $ List.tails xs

-- threads a value through a function, printing the output of the
--  function and returning the original value
traceWith :: (a -> String) -> a -> a
traceWith f x = trace (f x) x

-- HACK; Debug.Trace apparently should have this but mine doesn't.
-- (my base library is probably still the one from the Canonical repos >_>)
traceShowId :: Show a => a -> a
traceShowId = traceWith show

-- another variant of divMod and quotRem satisfying the property
-- that `pMod` is never negative, even if the divisor is.
pDivMod :: Integral a => a -> a -> (a,a)
pDivMod a b = (d, a-b*d)
    where d = (signum b) * (a `div` (abs b))

pDiv :: Integral a => a -> a -> a
pMod :: Integral a => a -> a -> a
pDiv a b = fst $ pDivMod a b
pMod a b = snd $ pDivMod a b

Specific points of concern:

  • In IntegerRef you'll see the use of a helper function called compose to string long chains of operations together into one function that performs them in sequence. I understand that this is precisely the sort of problem that Monads are meant to solve, but I'm not sure how one could help. (the only ones I really understand are [] and Maybe)
  • Something I mention in the notes above validateIrref (in IntegerRref.hs); a typical debugging strategy of mine is to verify post-conditions of functions through expensive checks which are only enabled in debug mode. But in Haskell, with lazy evaluation, that might cause the "strictness" of a function to depend on whether assertions are enabled! Seems troubling...
  • Better ways to write more things in a pointfree style. Or conversely, places in my code where I used a pointfree style to the detriment of readability.
  • Any "anti-idioms" I use that I should be aware of
\$\endgroup\$
  • \$\begingroup\$ Why not use the aptly named matrix library? 1) it's much faster, 2) you shouldn't need to worry whether their functions are correct so can debug less. \$\endgroup\$ – Michael Klein Jun 9 '16 at 15:05
  • \$\begingroup\$ @MichaelKlein a valid question. Obvious true answer aside (which is "I didn't look for one"), I think part of it is that most of my proofs were structured recursively and so lists seemed a natural fit. The matrix-related type aliases and methods are largely things I pulled out as an afterthought when it began to bother me that code working on columns tended to look so different from code working on rows. \$\endgroup\$ – Exp HP Jun 10 '16 at 0:15
2
\$\begingroup\$

Here's a little to start with:

Safe.Exact implements zipWithExact.

Lens implements listSet n as set (ix n), and overRow n as over (ix n). These don't error on being out of range, instead doing nothing. listSet is a bad name because list is also a verb and set is also a noun.

overLowerRight n is over (foldr (.) id (replicate n _tail) . each . foldr (.) id (replicate n _tail) . each).

generateGroup should be called generateAbelianMagma or generateSemigroup, because you aren't generating the neutral element or inverses, and are using either commutativity or associativity to only append the original generators to any new elements, and only to the right.

generateSemigroup :: Ord g => (g -> g -> g) -> [g] -> Set g
generateSemigroup op generators = foldr foo Set.empty generators where
  foo :: g -> Set g -> Set g
  foo x set = if S.member x set
    then set
    else foldr foo (insert x set) $ (`op` x) <$> generators
\$\endgroup\$
  • \$\begingroup\$ Hmm... Perhaps I should have called it generateFiniteGroup? I think that ought to at least address the inverses and identity problem because then each generator has the identity in its cyclic subgroup. I'm still trying to digest your new definition, but it's impressive! \$\endgroup\$ – Exp HP Jun 6 '16 at 17:54
  • \$\begingroup\$ Ah, that explains why you error out on an empty generator list - the neutral element isn't generated then. And for infinite groups, this doesn't halt anyway, I guess. \$\endgroup\$ – Gurkenglas Jun 7 '16 at 21:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.