MatrixStuff.hs

{-# LANGUAGE NoMonomorphismRestriction #-}
module MatrixStuff where

import Data.List (sortBy, transpose, findIndex, find)
import Data.Ord (comparing)
import Control.Monad (ap)
import Control.Monad.Instances
import Data.Ratio
import Data.Complex

-- ****** general-purpose stuffs

-- | the square root of (-1)
i = 0 :+ 1


-- ****** vector-related stuffs

-- |Inner product of real vectors
dotprod = (sum .) . zipWith (*)

-- |Inner product of complex vectors
dotProdC a b = dotprod a (map conjugate b)

project vector onto = zipWith (*) vector (norm onto)

lenSquared = sum . map ((^2) . abs)

len = sqrt . lenSquared

norm xs = let s = len xs in map (/s) xs

-- ****** matrix-related stuffs

-- |Scale a matrix.
scale scalar = map (map (*scalar))

-- |Multiply matrix*matrix
mmult a tb = let b = transpose tb in
    map (\r -> map (dotprod r) b) a

-- |Add matrices elementwise
madd = zipWith (zipWith (+))
-- |Subtract matrices elementwise
msubtract = zipWith (zipWith (-))

-- | nXn identity matrix
eye m = map (\n -> replicate n 0 ++ [1] ++ replicate (m - n - 1) 0) [0..m - 1]

-- | nXn matrix of zeroes
noughts m = replicate m (replicate m 0)

-- | matrix inverse
inverse xs = 
    let n = length xs in
    map (drop n) $ 
    rref $
    zipWith (++) xs (eye n)

-- | reduced row echelon form
rref xs = map normalize $ 
          reverse $
          scanl1Tails combineUp $
          reverse $ ref xs
              where normalize xs = case (find (/= 0) xs) of
                                     Nothing -> xs
                                     Just x -> map (/x) xs

-- | row echelon form
ref [] = []
ref xs | null (head xs) = xs
ref m = next ref (zeroOutBelow $ sortRows m)


-- ******* quantum stuffs


-- |The (2^n)X(2^n) Hadamard transform, without a normalization factor.
had_unscaled :: (Eq a, Num a, Num t) => a -> [[t]]
had_unscaled 1 = [[1, 1], [1, -1]]
had_unscaled n = 
    let lower = had_unscaled (n - 1) in
    zipWith (++) lower lower ++
    zipWith (++) lower (scale (-1) lower)

-- |The (2^n)X(2^n) Hadamard transform.
had :: (Eq a, Floating a) => a -> [[a]]
had n = scale (2**(-(n/2))) $ had_unscaled n

-- |The nXn quantum fourier transform, numerically.
qft :: RealFloat a => Int -> [[Complex a]]
qft n = 
    qft_base (exp (i * 2 * pi / fromIntegral n)) n

-- |The nXn quantum fourier transform, taking omega as an argument.
qft_base :: Num t => t -> Int -> [[t]]
qft_base omega n = 
     [[omega^a | a <- take n $ [0, r ..]] | r <- [0..n - 1]]

cnot = [[1, 0, 0, 0],
        [0, 1, 0, 0],
        [0, 0, 0, 1],
        [0, 0, 1, 0]]

cswap = [[1, 0, 0, 0, 0, 0, 0, 0],
         [0, 1, 0, 0, 0, 0, 0, 0],
         [0, 0, 1, 0, 0, 0, 0, 0],
         [0, 0, 0, 1, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 0],
         [0, 0, 0, 0, 0, 1, 0, 0],
         [0, 0, 0, 0, 0, 0, 0, 1]]


-- ****** the Approx class

-- | Approx: Can be convenient to use when you're working numerically and are getting results that are very close to integers. Goes inside matrices and the real/imaginary parts of complex numbers.
-- | e.x.: 
-- > [sqrt n * sqrt n | n <- [1..10]] == [1.0,2.0000000000000004,2.9999999999999996,4.0,5.000000000000001,5.999999999999999,7.000000000000001,8.000000000000002,9.0,10.000000000000002]
-- > approx [sqrt n * sqrt n | n <- [1..10]] == [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0]
class Approx a where
    approx :: a -> a

instance Approx Float where
    approx x = 
        let r = approxRational x 0.000001 in
        fromIntegral (numerator r) / fromIntegral (denominator r)


instance Approx Double where
    approx x = 
        let r = approxRational x 0.000001 in
        fromIntegral (numerator r) / fromIntegral (denominator r)



instance Approx a => Approx (Complex a) where
    approx (a :+ b) = approx a :+ approx b

instance Approx a => Approx [a] where
    approx xs = map approx xs


{- | Symbolic algebra, slightly similar on the surface to simple-reflect but meant for simplifying expressions
   | The important bits are the Symbol constructor and the simplify function.
   | Example use:
   
   > let w = Symbol "w"
   > mapM_ print $ map (map simplify) $ qft_base w 

   Results in:

   [1,1,1]
   [1,w,w^2]
   [1,w^2,w^4]
-}

data Expr a = 
    Const a | Symbol String | Plus (Expr a) (Expr a) | Negate (Expr a) | Times (Expr a) (Expr a) | Div (Expr a) (Expr a) | Exp (Expr a) (Expr a) | Fun String (Expr a)
        deriving Eq


parensShow (Const a) = show a
parensShow (Symbol s) = s
parensShow a = "(" ++ show a ++ ")"

instance Show a => Show (Expr a) where
    show (Const a) = show a
    show (Symbol s) = s
    show (Plus a (Negate b)) = (parensShow a) ++ " - " ++ (parensShow b)
    show (Plus a b) = (parensShow a) ++ " + " ++ (parensShow b)
    show (Negate a) = "-" ++ parensShow a
    show (Times a b) = (parensShow a) ++ " * " ++ (parensShow b)
    show (Div a b) = (parensShow a) ++ " / " ++ (parensShow b)
    show (Exp a b) = (parensShow a) ++ "^" ++ (parensShow b)
    show (Fun s a) = s ++ " " ++ parensShow a

instance Num a => Num (Expr a) where
    a + b = a `Plus` b
    a - b = a `Plus` (Negate b)
    a * b = a `Times` b
    negate a = Negate a
    abs a = Fun "abs" a
    signum a = Fun "signum" a
    fromInteger a = Const (fromInteger a)

breakExponent (Exp a b) = (a, b)
breakExponent a = (a, 1)

-- eval and simplify pretty much taken from http://stackoverflow.com/questions/324311/symbolic-simplification-in-haskell-using-recursion
-- though with some extra rules
eval :: (Num a, Eq a) => Expr a -> Expr a
eval (Times (Const a) (Const b)) = Const (a * b)
eval (Times (Const a) b)
    | a == 0 = Const 0
    | a == 1 = b
    | otherwise = (Times (Const a) b)
eval (Times a (Const b))
    | b == 0 = Const 0
    | b == 1 = a
    | otherwise = (Times a (Const b))
eval e@(Times a b) = 
    let (baseA, expA) = breakExponent a
        (baseB, expB) = breakExponent b in
    if baseA == baseB 
    then Exp baseA (eval (expA + expB))
    else e
eval (Plus (Const a) (Const b)) = Const (a + b)
eval (Plus (Const a) (Negate (Const b))) = Const (a - b)
eval (Plus (Negate (Const a)) (Const b)) = Const (-a + b)
eval (Plus (Negate (Const a)) (Negate (Const b))) = Const (-a - b)
eval (Plus (Const a) b)
    | a == 0 = b
    | otherwise = (Plus (Const a) b)
eval (Plus a (Const b))
    | b == 0 = a
    | otherwise = (Plus a (Const b))
eval e = e

simplify :: (Eq a, Num a) => Expr a -> Expr a
simplify (Plus a b) = eval (Plus (simplify a) (simplify b))
simplify (Times a b) = eval (Times (simplify a) (simplify b))
simplify e = e




-- ****** support code
combineUp lower upper
    | all (==0) lower = upper
    | otherwise = let (Just i) = findIndex (/= 0) lower 
                      mult = (upper !! i) / (lower !! i) in
                  zipWith (-) upper $ map (* mult) lower


next f [] = []
next f (x:xs) = 
    x :
    let cols = length $ takeWhile (==0) $ head xs
        (beginSplit, endSplit) = unzip $ map (splitAt cols) xs in
    zipWith (++) beginSplit (f endSplit)

zeroOutBelow [] = []
zeroOutBelow (x:xs) = 
    let h = head x in
    x : map (\(line@(l:ls)) -> 
                 case l of
                   0 -> line
                   _ -> zipWith (\a b -> b - l / h * a) x line) xs


sortPre f = map fst . sortBy (comparing snd) . ap zip (map f)

sortRows = sortPre (length . takeWhile (==0))

scanl1Tails :: (a -> a -> a) -> [a] -> [a]
scanl1Tails f (x:xs) = scanlTails f x xs
scanl1Tails f [] = []

-- |scan through a finite list applying each element to each element past it
-- |(unlike normal scanl which applies each element to the next element)
scanlTails f q ls = q : case ls of
                          [] -> []
                          x:xs -> scanlTails f (f q x) (map (f q) xs)