diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,3 @@ +# Changelog for HNumeric + +## Unreleased changes diff --git a/HNumeric.cabal b/HNumeric.cabal new file mode 100644 --- /dev/null +++ b/HNumeric.cabal @@ -0,0 +1,67 @@ +-- This file has been generated from package.yaml by hpack version 0.28.2. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: 0a4c04fe65db3e146424c0f3450ab7549b8db8ec8c21cb1d08bc26575227481f + +name: HNumeric +version: 0.2.0.0 +synopsis: Haskell Numeric Library with pure functionality, R & MATLAB Syntax. +description: Please see the README on GitHub at +category: HNum, library, Numeric, LinearAlgebra, Statistics, bsd3 +homepage: https://github.com/Axect/HNumeric#readme +bug-reports: https://github.com/Axect/HNumeric/issues +author: Tae Geun Kim +maintainer: edeftg@gmail.com +copyright: 2018 Tae Geun Kim +license: BSD3 +license-file: LICENSE +build-type: Simple +cabal-version: >= 1.10 +extra-source-files: + ChangeLog.md + README.md + +source-repository head + type: git + location: https://github.com/Axect/HNumeric + +library + exposed-modules: + HNum.Stats + HNum.Vector + other-modules: + Paths_HNumeric + hs-source-dirs: + src + build-depends: + base >=4.7 && <5 + , normaldistribution + default-language: Haskell2010 + +executable HNumeric-exe + main-is: Main.hs + other-modules: + Paths_HNumeric + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + HNumeric + , base >=4.7 && <5 + , normaldistribution + default-language: Haskell2010 + +test-suite HNumeric-test + type: exitcode-stdio-1.0 + main-is: Spec.hs + other-modules: + Paths_HNumeric + hs-source-dirs: + test + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + HNumeric + , base >=4.7 && <5 + , normaldistribution + default-language: Haskell2010 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Tae Geun Kim (c) 2018 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Author name here nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,193 @@ +# HNumeric + +![Travis](https://travis-ci.org/Axect/HNumeric.svg?branch=master) + +## Packages + +* HNum.Vector : Contain vector, matrix, linear algebra +* HNum.Stats : Contain statistical functions + +## Installation + +### 1. Native Use + +You can use this package just change `app/Main.hs` +Then, just type next command + +```bash +git clone https://github.com/Axect/HNumeric +``` + +### 2. Cabal Install + +First, install prerequisite. + +```sh +cabal install normaldistribution +``` + +Second, install tarball & unpack. + +```sh +export VERSION=0.2.0.0 + +wget -O HNumeric-${VERSION}.tar.gz https://github.com/Axect/HNumeric/blob/master/dist/HNumeric-${VERSION}.tar.gz\?raw\=true + +tar -xvf HNumeric-${VERSION}.tar.gz +``` + +Finally, build & install + +```sh +cabal build +cabal install +``` + +Finish! + +### 3. Import to Stack project + +If you use this package to your own project, then you should change `stack.yaml` and `package.yaml` + +#### 1) Change `stack.yaml` + +```yaml +# In stack.yaml +extra-deps: + - git: https://github.com/Axect/HNumeric.git + commit: [Latest Commit] + - normaldistribution-1.1.0,3 +``` + +* Replace [Latest Commit] with latest commit in [HNumeric Commit](https://github.com/Axect/HNumeric/commits/master) + +#### 2) Change `package.yaml` + +```yaml +# In package.yaml +dependecies: +- base +- HNumeric +- normaldistribution +``` + +Then enjoy! + +## Usage + +### Import Module + +* HNum.Vector +* HNum.Stats + +### Basic Vector Usage + +```haskell +-- HNumeric-0.2.0.0 Documentation + +let a = vector [1,2,3] -- Vector declaration +let b = vector [4,5,6] + +-- Print Vector +print a + +-- You can (+1) by fmap (Vector is functor) +(+1) <$> a + +-- Or MATLAB-like operator (.+, .-, .*, ./, .^) +a .+ 1 + +-- You can make list from vector +toList a -- [1, 2, 3] + +-- You can make vector from list +fromList [1,2,3] -- Vector [1,2,3] + +-- You can add (subtract, multiply, divide) vectors +a + b -- Vector [5,7,9] + +-- Also dot product is here. +a .*. b -- 1*4 + 2*5 + 3*6 = 32 + +-- Declare Matrix (Syntactic Sugar) +let c = matrix [[1,2],[3,4]] + +-- or Declare using R Syntax +let d = Matrix {val = Vector [5,6,7,8], row = 2, col = 2, byRow = True} + +-- Determinant +det c + +-- Inverse +inv c + +-- Transpose +transpose c + +-- Matrix ops with Constant (+, -, *, /, ^) +c .+ 1 -- Matrix [[2,3],[4,5]] + +-- Matrix ops with Matrix (+, -) +c + c -- Matrix [[2,4],[6,8]] + +-- Matrix Multiplication +c %*% d + +-- Matrix - Inverse Multiplication +c %/% d + +-- Vector Concatenate +hcat a b -- Vector [1,2,3,4,5,6] +vcat a b -- Matrix [[1,2,3],[4,5,6]] + +-- Matrix Concatenate +hcat c d -- Matrix [[1,2,5,6],[3,4,7,8]] +vcat c d -- Matrix [[1,2],[3,4],[5,6],[7,8]] + +-- Insert Vector to Matrix +Vector [1, 2] .: c -- Matrix [[1,2],[1,2],[3,4]] +``` + +### Basic Stats Usage + +```haskell +-- Sample Vector (import Vector) +v = Vector [1..10] +w = Vector [10, 9 .. 1] + +-- Mean +mean v + +-- Var +var v + +-- Std +std v + +-- Cov Matrix +cov v w + +-- Only Cov +cov' v w + +-- Linear Fit +(intercept, slope) = lm v w -- (11.0, -1.0) -- (Intercept, Slope) + +-- Linear Fit function +lineFit (intercept, slope) (Vector [1 .. 20]) + +-- RSS +rss v w + +-- RSE +rse v w +``` + +### TODO + +* ~~Effective Matrix Multiplication~~ +* Write Vector to CSV +* ~~Haddock~~ +* DataFrame using Map +* ~~Fix Matrix Implementation~~ +* ~~Numeric Class Define~~ diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/app/Main.hs b/app/Main.hs new file mode 100644 --- /dev/null +++ b/app/Main.hs @@ -0,0 +1,12 @@ +module Main where + +import HNum.Vector +import HNum.Stats +import Data.Random.Normal + +main :: IO () +main = do + let a = Vector [1, 3, 4] + print a + print $ map ($ a) [mean, var, std] + diff --git a/src/HNum/Stats.hs b/src/HNum/Stats.hs new file mode 100644 --- /dev/null +++ b/src/HNum/Stats.hs @@ -0,0 +1,48 @@ +module HNum.Stats where + +import HNum.Vector + +type Coeff a = (a, a) + +-- | Expectation Value +mean :: Fractional a => Vector a -> a +mean v = sum v / fromIntegral (length v) + +-- | Covariance (Single-Valued) +cov' :: Floating a => Vector a -> Vector a -> a +cov' x y + | length x <= 1 || length y <= 1 = error "Samples are not enough" + | length x /= length y = error "Length is not same" + | otherwise = ((x .- mean x) .*. (y .- mean y)) / fromIntegral (length x - 1) + +-- | Variance +var :: Floating a => Vector a -> a +var v = cov' v v + +-- | Standard Deviation +std :: Floating a => Vector a -> a +std = sqrt . var + +-- | Covariance Matrix +cov :: Floating a => Vector a -> Vector a -> Matrix a +cov x y = matrix [[var x, cov' x y], [cov' y x, var y]] + +-- | Least Square Method - (Intercept, Slope) +lm :: Floating a => Vector a -> Vector a -> Coeff a +lm x y = (my - b1 * mx, b1) + where + mx = mean x + my = mean y + b1 = (x .- mx) .*. (y .- my) / ((x .- mx) .*. (x .- mx)) + +-- | Line Fitting with (Intercept, Slope) & Range of x +lineFit :: Floating a => Coeff a -> Vector a -> Vector a +lineFit (n, m) x = x .* m .+ n + +-- | Residual Sum of Squares +rss :: Floating a => Vector a -> Vector a -> a +rss x y = sum ((y - lineFit (lm x y) x) .^ 2) + +-- | Relative Standard Error +rse :: Floating a => Vector a -> Vector a -> a +rse x y = sqrt (1 / fromIntegral (length x - 2) * rss x y) diff --git a/src/HNum/Vector.hs b/src/HNum/Vector.hs new file mode 100644 --- /dev/null +++ b/src/HNum/Vector.hs @@ -0,0 +1,403 @@ +{-| +Module : HNumeric.Vector +Description : Haskell Vector & Matrix & Linear Algebra Library to do machine learning +CopyRight : (c) Tae Geun Kim, 2018 +License : GPL-3 +Maintainer : edeftg@gmail.com +Stability : Experimental +-} +module HNum.Vector where + +import Data.Functor ( ) +import Control.Applicative ( ) + +--------------------------------------------------- +-- Vector +-------------------------------------------------- +-- Type Section +newtype Vector a = Vector [a] deriving (Show, Eq) + +vector :: [a] -> Vector a +vector = Vector + +-- Instance Section +instance Functor Vector where + fmap f (Vector x) = Vector (fmap f x) + +instance Applicative Vector where + pure a = Vector [] + Vector fs <*> Vector xs = Vector (zipWith ($) fs xs) + +instance (Num a) => Num (Vector a) where + negate v = negate <$> v + (+) v1 v2 = (+) <$> v1 <*> v2 + (*) v1 v2 = (*) <$> v1 <*> v2 + fromInteger n = fromInteger <$> Vector [n] + signum v = signum <$> v + abs v = abs <$> v + +instance (Fractional a) => Fractional (Vector a) where + recip v = recip <$> v + (/) v1 v2 = (*) <$> v1 <*> recip v2 + fromRational n = fromRational <$> Vector [n] + +instance (Floating a) => Floating (Vector a) where + pi = Vector [pi] + exp v = exp <$> v + log v = log <$> v + sqrt v = sqrt <$> v + sin v = sin <$> v + cos v = cos <$> v + tan v = tan <$> v + asin v = asin <$> v + acos v = acos <$> v + atan v = atan <$> v + sinh v = sinh <$> v + cosh v = cosh <$> v + tanh v = tanh <$> v + asinh v = asinh <$> v + acosh v = acosh <$> v + atanh v = atanh <$> v + +instance Foldable Vector where + foldr _ z (Vector []) = z + foldr f z (Vector xs) = foldr f z xs + + foldl _ z (Vector []) = z + foldl f z (Vector xs) = foldl f z xs + +-- Class for Vector with List +class List m where + toList :: m a -> [a] + fromList :: [a] -> m a + +instance List Vector where + toList (Vector xs) = xs + fromList = Vector + +--------------------------------------------------- +-- Matrix +--------------------------------------------------- + +-- |Definition of Matrix +data Matrix a = Matrix {val :: Vector a, row :: Int, col :: Int, byRow :: Bool} deriving (Eq) + +-- |matrix is syntactic sugar to create Matrix +matrix :: [[a]] -> Matrix a +matrix = formMat + +-- |Matrices is necessary class for Matrix +class Matrices m where + matForm :: m a -> [[a]] + formMat :: [[a]] -> m a + +instance Matrices Matrix where + matForm (Matrix (Vector v) r c b) + | r*c /= length v = error "Matrix Dimension mismatch!" + | b = ctake c v + | otherwise = dtake c v + where ctake :: Int -> [a] -> [[a]] + ctake _ [] = [] + ctake n m = take n m : ctake n (drop n m) + dtake :: Int -> [a] -> [[a]] + dtake _ [] = [] + dtake n m = [ptake n m r | r <- [0..(length m `div` n - 1)]] + ptake n v r = [v !! x | x <- idx v, x `mod` (length v `div` n) == r] + idx v = take (length v) [0..] + formMat [] = Matrix (Vector []) 0 0 True + formMat xs = Matrix (Vector (concat xs)) (length xs) (length (head xs)) True + +instance Show a => Show (Matrix a) where + show m = "Matrix " ++ show (matForm m) + +instance Functor Matrix where + fmap f mat = mat { val = fmap f (val mat) } + +instance Applicative Matrix where + pure a = matrix [] + mf <*> mx = mx { val = val mf <*> val mx } + +instance Num a => Num (Matrix a) where + negate m = negate <$> m + (+) m n = (+) <$> m <*> n + (*) m n = (*) <$> m <*> n + fromInteger a = fromInteger <$> matrix [[a]] + signum m = signum <$> m + abs m = abs <$> m + +instance Fractional a => Fractional (Matrix a) where + recip m = recip <$> m + (/) m n = (*) <$> m <*> recip n + fromRational n = fromRational <$> matrix [[n]] + +instance (Floating a) => Floating (Matrix a) where + pi = matrix [[pi]] + exp v = exp <$> v + log v = log <$> v + sqrt v = sqrt <$> v + sin v = sin <$> v + cos v = cos <$> v + tan v = tan <$> v + asin v = asin <$> v + acos v = acos <$> v + atan v = atan <$> v + sinh v = sinh <$> v + cosh v = cosh <$> v + tanh v = tanh <$> v + asinh v = asinh <$> v + acosh v = acosh <$> v + atanh v = atanh <$> v + +instance Foldable Matrix where + foldr _ z (Matrix (Vector []) _ _ _) = z + foldr f z (Matrix (Vector xs) _ _ _) = foldr f z xs + + foldl _ z (Matrix (Vector []) _ _ _) = z + foldl f z (Matrix (Vector xs) _ _ _) = foldl f z xs + +--------------------------------------------------- +-- Operation +--------------------------------------------------- +{-| + (.) is an operation Vector(or Matrix) with Constant. + Dot means position of Vector. + Example: a .* 2 = twice whole elements of a + a .*. b = Dot product +-} +class Functor f => VecOps f where + (.+) :: Num a => f a -> a -> f a + (.-) :: Num a => f a -> a -> f a + (.*) :: Num a => f a -> a -> f a + (./) :: Fractional a => f a -> a -> f a + (.^) :: Floating a => f a -> a -> f a + (.*.) :: Num a => f a -> f a -> a + norm :: Floating a => f a -> a + +{- + MatOps is just additional operations for Matrices. +-} +class Functor f => MatOps f where + (%*%) :: Num a => f a -> f a -> f a + (%/%) :: Fractional a => f a -> f a -> f a + det :: Fractional a => f a -> a + inv :: Fractional a => f a -> f a + transpose :: f a -> f a + +instance VecOps Vector where + v .+ n = (+ n) <$> v + v .- n = (+ negate n) <$> v + v .* n = (* n) <$> v + v ./ n = (/ n) <$> v + v .^ n = (** n) <$> v + v .*. w = sum $ v * w + norm v = sqrt $ v .*. v + +instance VecOps Matrix where + v .+ n = (+ n) <$> v + v .- n = (+ negate n) <$> v + v .* n = (* n) <$> v + v ./ n = (/ n) <$> v + v .^ n = (** n) <$> v + v .*. w = sum $ v * w + norm v = sqrt $ v .*. v + +instance MatOps Matrix where + m %*% n | col m /= row n = error "Can't Multiply - Dimension mismatch!" + | otherwise = matrix $ matForm m %-*-% matForm n + m %/% n = m %*% inv n + det m | col m /= row m = error "Can't calculate determinant of non-square matrix" + | otherwise = detMat (matForm m) + inv m | col m /= row m = error "Can't calculate inverse of non-square matrix" + | otherwise = (matrix . invMat . matForm) m + transpose m = m {row = col m, col = row m, byRow = not (byRow m)} + +-- |Block Partitioning +bp :: Int -> Matrix a -> Matrix a +bp n m = matrix $ bpMat n (matForm m) + +--------------------------------------------------- +-- Concatenate +--------------------------------------------------- +class Functor f => Concatable f where + hcat :: f a -> f a -> f a + vcat :: f a -> f a -> Matrix a + +instance Concatable Vector where + hcat v w = fromList (toList v ++ toList w) + vcat v w = matrix (toList v : [toList w]) + +instance Concatable Matrix where + hcat m n | row m == row n = matrix (zipWith (++) mf nf) + | otherwise = error "Can't concatenate matrices horizontally which have different row" + where mf = matForm m + nf = matForm n + vcat m n | col m == col n = m {val = hcat (val m) (val n), row = row m + row n} + | otherwise = error "Can't concatenate matrices vertically which have different col" + +-- |(.:) inserts vector to head of matrix. +(.:) :: Vector a -> Matrix a -> Matrix a +v .: m | length v == col m = matrix (toList v : matForm m) + | otherwise = error "Can't insert length(Vector) /= col(Matrix)" + +--------------------------------------------------- +-- Backend Functions (Do not Understand) +--------------------------------------------------- +-- Transpose +transposeMat :: [[a]] -> [[a]] +transposeMat m = map (\l -> map (!! l) m) [0 .. (length (head m) - 1)] + +-- indexMat +indexMat :: [[a]] -> [[(Int, Int)]] +indexMat m@(xs : xss) = do + i <- [0 .. (length m - 1)] + [zip (replicate (length xs) i) [0 .. (length xs - 1)]] + +-- dropAtMat +dropAtMat :: Int -> Int -> [[a]] -> [[a]] +dropAtMat i j mat = map (dropAt j) $ dropAt i mat + +-- postSplitAt +postSplitAt (x, y) = x ++ tail y + +-- dropAt +dropAt :: Int -> [a] -> [a] +dropAt i = postSplitAt . splitAt i + +-- dropAtMat' : drop nth array +dropAtMat' :: Int -> [[a]] -> [[a]] +dropAtMat' n mat | n /= (length mat - 1) = dropAt n mat + | otherwise = take n mat + +-- Block Partitioning +bpMat :: Int -> [[a]] -> [[a]] +bpMat _ [] = [] +bpMat n m | n == 1 = (map (take sl) . take sl) m + | n == 2 = (map (drop sl) . take sl) m + | n == 3 = (map (take sl) . drop sl) m + | n == 4 = (map (drop sl) . drop sl) m + | otherwise = error "Please input 1 ~ 4" + where + l = length m + sl = (floor . sqrt . fromIntegral) l + +-- Matrix + Matrix +(%-+-%) :: Num a => [[a]] -> [[a]] -> [[a]] +m %-+-% [] = m +[] %-+-% m = m +[[]] %-+-% m = m +m %-+-% [[]] = m +m %-+-% n = zipWith (zipWith (+)) m n + +negMap :: Num a => [[a]] -> [[a]] +negMap = map (map negate) + +-- Matrix - Matrix +(%---%) :: Num a => [[a]] -> [[a]] -> [[a]] +m %---% [] = m +[] %---% m = map (map negate) m +[[]] %---% m = map (map negate) m +m %---% [[]] = m +m %---% n = zipWith (zipWith (-)) m n + +-- Matrix Multiplication +(%-*-%) :: Num a => [[a]] -> [[a]] -> [[a]] +_ %-*-% [] = [] +[] %-*-% _ = [] +_ %-*-% [[]] = [[]] +[[] ] %-*-% _ = [[]] +[[x]] %-*-% [[y]] = [[x * y]] +m %-*-% n = zipWith (++) a11 a12 ++ zipWith (++) a21 a22 + where + (m11, n11) = (bpMat 1 m, bpMat 1 n) + (m12, n12) = (bpMat 2 m, bpMat 2 n) + (m21, n21) = (bpMat 3 m, bpMat 3 n) + (m22, n22) = (bpMat 4 m, bpMat 4 n) + a11 = (m11 %-*-% n11) %-+-% (m12 %-*-% n21) + a12 = (m11 %-*-% n12) %-+-% (m12 %-*-% n22) + a21 = (m21 %-*-% n11) %-+-% (m22 %-*-% n21) + a22 = (m21 %-*-% n12) %-+-% (m22 %-*-% n22) + +zerosVec :: Int -> [Int] +zerosVec n = take n [0, 0 ..] + +eyeMat :: Int -> [[Int]] +eyeMat n = [ basisVec x n | x <- [0 .. (n - 1)] ] + +-- Position -> Length +basisVec :: Int -> Int -> [Int] +basisVec n m = zerosVec n ++ [1] ++ zerosVec (m - n - 1) + +permMat :: Int -> Int -> [[a]] -> [[Int]] +permMat i j m + | i < j + = take i idx + ++ [idx !! j] + ++ take (j - i - 1) (drop (i + 1) idx) + ++ [idx !! i] + ++ drop (j + 1) idx + | otherwise + = permMat j i m + where idx = eyeMat (length m) + +whichMax :: Ord a => [a] -> Int +whichMax v = whichMax' v 0 m + where + m = maximum v + whichMax' :: Ord a => [a] -> Int -> a -> Int + whichMax' (x : xs) n m' = if x == m' then n else whichMax' xs (n + 1) m' + +colMat :: [[a]] -> Int -> [a] +colMat m n = map (!! n) m + +colMaxIdx :: Ord a => [[a]] -> Int -> Int +colMaxIdx m n = whichMax $ colMat m n + +-- | Another Block Partitioning +bpMat' :: Int -> [[a]] -> [[a]] +bpMat' _ [] = [] +bpMat' _ [x] = [x] +bpMat' n m | n == 1 = (map (take l) . take l) m + | n == 2 = (map (drop 1) . take l) m + | n == 3 = (map (take l) . drop 1) m + | n == 4 = (map (drop 1) . drop 1) m + | n == 0 = (map (drop 1 . take l) . drop 1 . take l) m + where l = length m - 1 + +-- | Determinant for Double List - Order ~ 4^n +detMat :: Fractional a => [[a]] -> a +detMat [[x]] = x +detMat m + | l == 2 = detMat m11 * detMat m22 - detMat m12 * detMat m21 + | otherwise = (detMat m11 * detMat m22 - detMat m12 * detMat m21) / detMat m00 + where + l = length m + m11 = bpMat' 1 m + m12 = bpMat' 2 m + m21 = bpMat' 3 m + m22 = bpMat' 4 m + m00 = bpMat' 0 m + +-- | Inverse for Double List - Order ~ n * 2^n +invMat :: Fractional a => [[a]] -> [[a]] +invMat [] = [] +invMat [[] ] = [[]] +invMat [[x]] = [[x]] +invMat m + | length m == 2 + = map (map (/ detMat m)) + $ zipWith (++) m22 (negMap m12) + ++ zipWith (++) (negMap m21) m11 + | otherwise + = zipWith (++) a11 a12 ++ zipWith (++) a21 a22 + where + m11 = bpMat 1 m + m12 = bpMat 2 m + m21 = bpMat 3 m + m22 = bpMat 4 m + a00 = invMat m11 + s = m22 %---% (m21 %-*-% a00 %-*-% m12) + s00 = invMat s + a11 = a00 %-+-% (a00 %-*-% m12 %-*-% s00 %-*-% m21 %-*-% a00) + a12 = negMap a00 %-*-% m12 %-*-% s00 + a21 = negMap s00 %-*-% m21 %-*-% a00 + a22 = s00 diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,2 @@ +main :: IO () +main = putStrLn "Test suite not yet implemented"