diff --git a/hmatrix-morpheus.cabal b/hmatrix-morpheus.cabal --- a/hmatrix-morpheus.cabal +++ b/hmatrix-morpheus.cabal @@ -1,5 +1,5 @@ name: hmatrix-morpheus -version: 0.1.0.2 +version: 0.1.1.0 synopsis: Low-level machine learning auxiliary functions. description: Purely functional interface to morpheus based on hmatrix. @@ -23,6 +23,7 @@ , src/simple_blas.h , src/matrix_reduce.h , src/activation.h + , src/statistics.h flag openblas description: Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries. @@ -34,6 +35,8 @@ exposed-modules: Numeric.Morpheus , Numeric.Morpheus.Activation , Numeric.Morpheus.MatrixReduce + , Numeric.Morpheus.Statistics + other-modules: Numeric.Morpheus.Utils build-depends: base >= 4.7 && < 5 , hmatrix >= 0.17.0.1 default-language: Haskell2010 @@ -42,6 +45,7 @@ C-sources: src/simple_blas.c , src/matrix_reduce.c , src/activation.c + , src/statistics.c cc-options: -O3 -Wall -std=c99 if arch(x86_64) cc-options: -msse2 @@ -76,6 +80,7 @@ , Test.HUnit.Plus , Numeric.Morpheus.ActivationTest , Numeric.Morpheus.MatrixReduceTest + , Numeric.Morpheus.StatisticsTest build-depends: base , hmatrix-morpheus , hmatrix >= 0.17.0.1 @@ -101,15 +106,6 @@ default-extensions: BangPatterns if os(OSX) extra-lib-dirs: /usr/lib - -executable hmatrix-morpheus-example - hs-source-dirs: hmatrix-morpheus/examples - main-is: Main.hs - build-depends: base - , hmatrix-morpheus - , hmatrix >= 0.17.0.1 - ghc-options: -threaded -rtsopts -with-rtsopts=-N - default-language: Haskell2010 source-repository head type: git diff --git a/hmatrix-morpheus/CHANGELOG.md b/hmatrix-morpheus/CHANGELOG.md --- a/hmatrix-morpheus/CHANGELOG.md +++ b/hmatrix-morpheus/CHANGELOG.md @@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html). +## [0.1.1.0] - 2017-08-27 +### Added +- Numeric.Morpheus.Statistics module - Statistics Functions. + ## [0.1.0.0] - 2017-08-26 ### Added - Numeric.Morpheus.Activation module - activation functions for neural network implementations; diff --git a/hmatrix-morpheus/examples/Main.hs b/hmatrix-morpheus/examples/Main.hs deleted file mode 100644 --- a/hmatrix-morpheus/examples/Main.hs +++ /dev/null @@ -1,41 +0,0 @@ -module Main where - -import Numeric.LinearAlgebra -import Numeric.Morpheus - - -a = matrix 5 [ - 71, 11, 3, -9, -7, - 21, -7, -2, 23, 11, - -11, 32, 53, -49, 37, - 1, -24, 78, 90, 17 - ] - - -main :: IO () -main = do - putStrLn "Numeric.Morpheus.MatrixReduce functions: " - - putStr "sum of elements of rows: " - print $ rowSum a - - putStr "product of elements of rows: " - print $ rowPredicate (*) a - - putStr "max values of columns: " - print $ fst $ columnMaxIndex a - - putStr "indices of max values of columns: " - print $ snd $ columnMaxIndex a - - - putStrLn "Numeric.Morpheus.Activation functions: " - - putStrLn "\nSigmoid:" - disp 3 (sigmoid a) - - putStrLn "ReLu:" - disp 3 (relu a) - - putStrLn "ReLu gradient:" - disp 3 (reluGradient a) diff --git a/hmatrix-morpheus/src/Numeric/Morpheus.hs b/hmatrix-morpheus/src/Numeric/Morpheus.hs --- a/hmatrix-morpheus/src/Numeric/Morpheus.hs +++ b/hmatrix-morpheus/src/Numeric/Morpheus.hs @@ -43,7 +43,7 @@ main :: IO () main = do - putStrLn "Numeric.Morpheus.MatrixReduce functions: " + putStrLn "\nNumeric.Morpheus.MatrixReduce functions: " putStr "sum of elements of rows: " print $ rowSum a @@ -58,7 +58,7 @@ print $ snd $ columnMaxIndex a - putStrLn "Numeric.Morpheus.Activation functions: " + putStrLn "\n\nNumeric.Morpheus.Activation functions: " putStrLn "\nSigmoid:" disp 3 (sigmoid a) @@ -68,6 +68,16 @@ putStrLn "ReLu gradient:" disp 3 (reluGradient a) + + + putStrLn "\n\nNumeric.Morpheus.Statistics functions: " + + putStrLn "\nColumn Means:" + let colmeans = columnMean a + print colmeans + + putStrLn "Column Standard Deviations: " + print (columnStddev_m colmeans a) @ -} @@ -75,9 +85,11 @@ ( module Numeric.Morpheus.Activation , module Numeric.Morpheus.MatrixReduce + , module Numeric.Morpheus.Statistics ) where import Numeric.Morpheus.Activation import Numeric.Morpheus.MatrixReduce +import Numeric.Morpheus.Statistics diff --git a/hmatrix-morpheus/src/Numeric/Morpheus/MatrixReduce.hs b/hmatrix-morpheus/src/Numeric/Morpheus/MatrixReduce.hs --- a/hmatrix-morpheus/src/Numeric/Morpheus/MatrixReduce.hs +++ b/hmatrix-morpheus/src/Numeric/Morpheus/MatrixReduce.hs @@ -21,16 +21,14 @@ where +import Numeric.Morpheus.Utils(morpheusLayout) + import Numeric.LinearAlgebra import Numeric.LinearAlgebra.Devel import System.IO.Unsafe(unsafePerformIO) import Foreign import Foreign.C.Types import Foreign.Ptr(Ptr) - - -morpheusLayout :: CInt -> CInt -> CInt -morpheusLayout xCol cols = if xCol == 1 || cols == 1 then 21 else 22 type Predicate = R -> R -> R diff --git a/hmatrix-morpheus/src/Numeric/Morpheus/Statistics.hs b/hmatrix-morpheus/src/Numeric/Morpheus/Statistics.hs new file mode 100644 --- /dev/null +++ b/hmatrix-morpheus/src/Numeric/Morpheus/Statistics.hs @@ -0,0 +1,145 @@ +{-| +Module: Numeric.Morpheus.Statistics +Description: Statistics Functions +Copyright: (c) Alexander Ignatyev, 2017 +License: BSD-3 +Stability: experimental +Portability: POSIX +-} + +module Numeric.Morpheus.Statistics +( + mean + , stddev_m + , stddev + , columnMean + , rowMean + , columnStddev_m + , rowStddev_m +) + +where + +import Numeric.Morpheus.Utils(morpheusLayout) + +import Numeric.LinearAlgebra +import Numeric.LinearAlgebra.Devel +import System.IO.Unsafe(unsafePerformIO) +import Foreign.C.Types +import Foreign.Ptr(Ptr) +import Foreign.Storable + +{- morpheus_mean -} +foreign import ccall unsafe "morpheus_mean" + c_morpheus_mean :: CInt -> Ptr R -> IO R + +-- | Calculates mean (average). +mean :: Vector R -> R +mean x = unsafePerformIO $ do + apply x id c_morpheus_mean + + +{- morpheus_stddev_m -} +foreign import ccall unsafe "morpheus_stddev" + c_morpheus_stddev_m :: R -> CInt -> Ptr R -> IO R + +-- | Calculates sample standard deviation using given mean value. +stddev_m :: R -> Vector R -> R +stddev_m m x = unsafePerformIO $ do + apply x id (c_morpheus_stddev_m m) + + +{- morpheus_stddev -} +foreign import ccall unsafe "morpheus_stddev" + c_morpheus_stddev :: CInt -> Ptr R -> IO R + +-- | Calculates sample standard deviation. +stddev :: Vector R -> R +stddev x = unsafePerformIO $ do + apply x id c_morpheus_stddev + + +{- morpheus_column_mean -} +foreign import ccall unsafe "morpheus_column_mean" + c_morpheus_column_mean :: CInt -> CInt -> CInt -> Ptr R -> Ptr R -> IO () + + +call_morpheus_column_mean :: CInt -> CInt -> CInt -> CInt -> Ptr R + -> CInt -> Ptr R + -> IO () +call_morpheus_column_mean rows cols xRow xCol mat _ means = do + let layout = morpheusLayout xCol cols + c_morpheus_column_mean layout rows cols mat means + + +-- | Calculates mean (average) for every column. +columnMean :: Matrix R -> Vector R +columnMean m = unsafePerformIO $ do + v <- createVector (cols m) + apply m (apply v id) call_morpheus_column_mean + return v + + +{- morpheus_row_mean -} +foreign import ccall unsafe "morpheus_row_mean" + c_morpheus_row_mean :: CInt -> CInt -> CInt -> Ptr R -> Ptr R -> IO () + + +call_morpheus_row_mean :: CInt -> CInt -> CInt -> CInt -> Ptr R + -> CInt -> Ptr R + -> IO () +call_morpheus_row_mean rows cols xRow xCol mat _ means = do + let layout = morpheusLayout xCol cols + c_morpheus_row_mean layout rows cols mat means + + +-- | Calculates mean (average) for every row. +rowMean :: Matrix R -> Vector R +rowMean m = unsafePerformIO $ do + v <- createVector (rows m) + apply m (apply v id) call_morpheus_row_mean + return v + + +{- morpheus_column_stddev_m -} +foreign import ccall unsafe "morpheus_column_stddev_m" + c_morpheus_column_stddev_m :: CInt -> CInt -> CInt -> Ptr R -> Ptr R -> Ptr R -> IO () + + +call_morpheus_column_stddev_m :: CInt -> Ptr R + -> CInt -> CInt -> CInt -> CInt -> Ptr R + -> CInt -> Ptr R + -> IO () +call_morpheus_column_stddev_m _ means rows cols xRow xCol mat _ stddevs = do + let layout = morpheusLayout xCol cols + c_morpheus_column_stddev_m layout rows cols means mat stddevs + + +-- | Calculates sample standard deviation using given mean value for every column. +columnStddev_m :: Vector R -> Matrix R -> Vector R +columnStddev_m means m = unsafePerformIO $ do + v <- createVector (cols m) + apply means (apply m (apply v id)) call_morpheus_column_stddev_m + return v + + +{- morpheus_row_stddev_m -} +foreign import ccall unsafe "morpheus_row_stddev_m" + c_morpheus_row_stddev_m :: CInt -> CInt -> CInt -> Ptr R -> Ptr R -> Ptr R -> IO () + + +call_morpheus_row_stddev_m :: CInt -> Ptr R + -> CInt -> CInt -> CInt -> CInt -> Ptr R + -> CInt -> Ptr R + -> IO () +call_morpheus_row_stddev_m _ means rows cols xRow xCol mat _ stddevs = do + let layout = morpheusLayout xCol cols + c_morpheus_row_stddev_m layout rows cols means mat stddevs + + +-- | Calculates sample standard deviation using given mean value for every row. +rowStddev_m :: Vector R -> Matrix R -> Vector R +rowStddev_m means m = unsafePerformIO $ do + v <- createVector (rows m) + apply means (apply m (apply v id)) call_morpheus_row_stddev_m + return v diff --git a/hmatrix-morpheus/src/Numeric/Morpheus/Utils.hs b/hmatrix-morpheus/src/Numeric/Morpheus/Utils.hs new file mode 100644 --- /dev/null +++ b/hmatrix-morpheus/src/Numeric/Morpheus/Utils.hs @@ -0,0 +1,20 @@ +{-| +Module: Numeric.Morpheus.Utils +Description: Utility Functions +Copyright: (c) Alexander Ignatyev, 2017 +License: BSD-3 +Stability: experimental +Portability: POSIX +-} + +module Numeric.Morpheus.Utils +( + morpheusLayout +) + +where + +import Foreign.C.Types + +morpheusLayout :: CInt -> CInt -> CInt +morpheusLayout xCol cols = if xCol == 1 || cols == 1 then 21 else 22 diff --git a/hmatrix-morpheus/test/Main.hs b/hmatrix-morpheus/test/Main.hs --- a/hmatrix-morpheus/test/Main.hs +++ b/hmatrix-morpheus/test/Main.hs @@ -1,10 +1,12 @@ import Test.Framework (defaultMain, testGroup) -import qualified Numeric.Morpheus.MatrixReduceTest as MatrixReduce import qualified Numeric.Morpheus.ActivationTest as Activation +import qualified Numeric.Morpheus.MatrixReduceTest as MatrixReduce +import qualified Numeric.Morpheus.StatisticsTest as Statistics main = defaultMain tests -tests = [ testGroup "Morpheus.MatrixReduce" MatrixReduce.tests - , testGroup "Morpheus.Activation" Activation.tests +tests = [ testGroup "Morpheus.Activation" Activation.tests + , testGroup "Morpheus.MatrixReduce" MatrixReduce.tests + , testGroup "Morpheus.Statistics" Statistics.tests ] diff --git a/hmatrix-morpheus/test/Numeric/Morpheus/StatisticsTest.hs b/hmatrix-morpheus/test/Numeric/Morpheus/StatisticsTest.hs new file mode 100644 --- /dev/null +++ b/hmatrix-morpheus/test/Numeric/Morpheus/StatisticsTest.hs @@ -0,0 +1,53 @@ +module Numeric.Morpheus.StatisticsTest +( + tests +) + +where + +import Test.Framework (testGroup) +import Test.Framework.Providers.HUnit +import Test.HUnit +import Test.HUnit.Approx +import Test.HUnit.Plus + +import Numeric.LinearAlgebra +import Numeric.Morpheus.Statistics + + +x = fromList [35,32,213,54,81,92,63,70] +meanValue = 80 +stddevValue = 57.60456 + +a = matrix 4 [35,32,213,54, + 81,92,63,70, + 77,14,91,20] +aT = tr a + +columnMeansRM = fromList [64.33333, 46, 122.33333, 48] +rowMeansRM = fromList [83.5, 76.5, 50.5] +columnMeansCM = rowMeansRM +rowMeansCM = columnMeansRM + +columnStddevsRM = fromList [25.48202, 40.84116, 79.75797, 25.53429] +rowStddevsRM = fromList [86.88114, 12.71482, 39.17908] +columnStddevsCM = rowStddevsRM +rowStddevsCM = columnStddevsRM + + +tests = [ testGroup "Statistics" [ + testCase "mean" $ assertApproxEqual "" 1e-5 meanValue (mean x) + , testCase "stddev_m" $ assertApproxEqual "" 1e-5 stddevValue (stddev_m meanValue x) + , testCase "stddev" $ assertApproxEqual "" 1e-5 stddevValue (stddev x) + ] + , testGroup "Matrix Statistics" [ + testCase "columnMean - row major" $ assertVector "" 1e-5 columnMeansRM (columnMean a) + , testCase "columnMean - column major" $ assertVector "" 1e-5 columnMeansCM (columnMean aT) + , testCase "rowMean - row major" $ assertVector "" 1e-5 rowMeansRM (rowMean a) + , testCase "rowMean - column major" $ assertVector "" 1e-5 rowMeansCM (rowMean aT) + , testCase "columnStddev_m - row major" $ assertVector "" 1e-5 columnStddevsRM (columnStddev_m columnMeansRM a) + , testCase "columnStddev_m - column major" $ assertVector "" 1e-5 columnStddevsCM (columnStddev_m columnMeansCM aT) + , testCase "rowStddev_m - row major" $ assertVector "" 1e-5 rowStddevsRM (rowStddev_m rowMeansRM a) + , testCase "rowStddev_m - column major" $ assertVector "" 1e-5 rowStddevsCM (rowStddev_m rowMeansCM aT) + ] + ] diff --git a/src/statistics.c b/src/statistics.c new file mode 100644 --- /dev/null +++ b/src/statistics.c @@ -0,0 +1,120 @@ +#include "statistics.h" +#include "matrix_reduce.h" + +#include +#include +#include + +double morpheus_mean(int n, const double *x) { + double m = 0; + for (int i = 0; i < n; ++i) { + m += x[i]; + } + return m/n; +} + +double morpheus_stddev_m(double mean, int n, const double *x) { + if (n < 2) return 0; + + double s = 0; + for (int i = 0; i < n; ++i) { + double v = x[i] - mean; + s += v*v; + } + return sqrt(s/(n-1)); +} + +double morpheus_stddev(int n, const double *x) { + if (n < 2) return 0; + + double m = morpheus_mean(n, x); + return morpheus_stddev_m(m, n, x); +} + +void morpheus_column_mean(morpheus_layout_e layout, + int nrows, int ncols, + const double *x, + double *means) { + morpheus_column_sum(layout, nrows, ncols, x, means); + for (int i = 0; i < ncols; ++i) { + means[i] /= nrows; + } +} + +void morpheus_row_mean(morpheus_layout_e layout, + int nrows, int ncols, + const double *x, + double *means) { + morpheus_row_sum(layout, nrows, ncols, x, means); + for (int i = 0; i < nrows; ++i) { + means[i] /= ncols; + } +} + +void morpheus_column_stddev_m(morpheus_layout_e layout, + int nrows, int ncols, + const double *means, + const double *x, + double *s) { + if (layout == morpheus_row_major) { + for (int j = 0; j < ncols; ++j) { + s[j] = 0; + } + for (int i = 0; i < nrows; ++i) { + const double *row = x + i*ncols; + for (int j = 0; j < ncols; ++j) { + double v = row[j] - means[j]; + s[j] += v*v; + } + } + } else if (layout == morpheus_col_major) { + for (int j = 0; j < ncols; ++j) { + const double *column = &x[j*nrows]; + s[j] = 0; + for (int i = 0; i < nrows; ++i) { + double v = column[i] - means[j]; + s[j] += v*v; + } + } + } else { + perror("unkwnown matrix layout"); + assert(0); + } + for (int j = 0; j < ncols; ++j) { + s[j] = sqrt(s[j]/(nrows-1)); + } +} + +void morpheus_row_stddev_m(morpheus_layout_e layout, + int nrows, int ncols, + const double *means, + const double *x, + double *s) { + if (layout == morpheus_row_major) { + for (int i = 0; i < nrows; ++i) { + const double *row = x + i*ncols; + s[i] = 0; + for (int j = 0; j < ncols; ++j) { + double v = row[j] - means[i]; + s[i] += v*v; + } + } + } else if (layout == morpheus_col_major) { + for (int i = 0; i < nrows; ++i) { + s[i] = 0; + } + for (int j = 0; j < ncols; ++j) { + const double *column = x + j*nrows; + for (int i = 0; i < nrows; ++i) { + double v = column[i] - means[i]; + s[i] += v*v; + } + } + } else { + perror("unkwnown matrix layout"); + assert(0); + } + for (int i = 0; i < nrows; ++i) { + s[i] = sqrt(s[i]/(ncols-1)); + } +} diff --git a/src/statistics.h b/src/statistics.h new file mode 100644 --- /dev/null +++ b/src/statistics.h @@ -0,0 +1,55 @@ +/*! \file statistics.h + \brief Statistics functions +*/ + +#ifndef MORHEPUS_STATISTICS_H +#define MORHEPUS_STATISTICS_H + +#include "simple_blas.h" + +/*! \defgroup statistics statistics + \brief Statistics functions + \{ +*/ + +/*! Calculates mean (average). */ +double morpheus_mean(int n, const double *x); + +/*! Calculates sample standard deviation using given mean value. */ +double morpheus_stddev_m(double mean, int n, const double *x); + +/*! Calculates sample standard deviation.*/ +double morpheus_stddev(int n, const double *x); + +/*! Calculates mean (average) for every column. */ +void morpheus_column_mean(morpheus_layout_e layout, + int nrows, int ncols, + const double *x, /*!< matrix of nrows x ncols */ + double *means /*!< output vector of mean values of size ncols */ + ); + +/*! Calculates mean (average) for every row. */ +void morpheus_row_mean(morpheus_layout_e layout, + int nrows, int ncols, + const double *x, /*!< matrix of nrows x ncols */ + double *means /*!< output vector of mean values of size nrows */ + ); + +/*! Calculates sample standard deviation using given mean value for every column. */ +void morpheus_column_stddev_m(morpheus_layout_e layout, + int nrows, int ncols, + const double *means, /*!< vector of mean values of size ncols */ + const double *x, /*!< matrix of nrows x ncols */ + double *s /*!< output vector of stddev values of size ncols */ + ); + +/*! Calculates sample standard deviation using given mean value for every row. */ +void morpheus_row_stddev_m(morpheus_layout_e layout, + int nrows, int ncols, + const double *means, /*!< vector of mean values of size nrows */ + const double *x, /*!< matrix of nrows x nrows */ + double *s /*!< output vector of stddev values of size nrows */ + ); +/*! \} */ + +#endif /* MORHEPUS_STATISTICS_H */