hmatrix-morpheus (empty) → 0.1.0.0
raw patch · 23 files changed
+1698/−0 lines, 23 filesdep +HUnitdep +MonadRandomdep +basesetup-changed
Dependencies added: HUnit, MonadRandom, base, criterion, hmatrix, hmatrix-morpheus, test-framework, test-framework-hunit, test-framework-quickcheck2
Files
- LICENSE +29/−0
- Setup.hs +2/−0
- hmatrix-morpheus.cabal +112/−0
- hmatrix-morpheus/CHANGELOG.md +12/−0
- hmatrix-morpheus/README.md +64/−0
- hmatrix-morpheus/bench/Main.hs +34/−0
- hmatrix-morpheus/bench/Naive.hs +19/−0
- hmatrix-morpheus/examples/Main.hs +41/−0
- hmatrix-morpheus/src/Numeric/Morpheus.hs +19/−0
- hmatrix-morpheus/src/Numeric/Morpheus/Activation.hs +151/−0
- hmatrix-morpheus/src/Numeric/Morpheus/MatrixReduce.hs +223/−0
- hmatrix-morpheus/test/Main.hs +10/−0
- hmatrix-morpheus/test/Numeric/Morpheus/ActivationTest.hs +38/−0
- hmatrix-morpheus/test/Numeric/Morpheus/MatrixReduceTest.hs +90/−0
- hmatrix-morpheus/test/Test/HUnit/Approx.hs +85/−0
- hmatrix-morpheus/test/Test/HUnit/Plus.hs +44/−0
- src/activation.c +41/−0
- src/activation.h +33/−0
- src/matrix_reduce.c +267/−0
- src/matrix_reduce.h +78/−0
- src/simple_blas.c +129/−0
- src/simple_blas.h +123/−0
- src/types.h +54/−0
+ LICENSE view
@@ -0,0 +1,29 @@+BSD 3-Clause License++Copyright (c) 2017, Alexander Ignatyev+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 the copyright holder nor the names of its+ 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 HOLDER 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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ hmatrix-morpheus.cabal view
@@ -0,0 +1,112 @@+name: hmatrix-morpheus+version: 0.1.0.0+synopsis: Low-level machine learning auxiliary functions.+description:+ Purely functional interface to morpheus based on hmatrix.+ Morpheus library contains a bunch of cache line aware numerical algorithms+ suitable for using as a low-level primitives to build machine learning applications.+homepage: https://github.com/alexander-ignatyev/morpheus+license: BSD3+license-file: LICENSE+author: Alexander Ignatyev+maintainer: ignatyev.alexander@gmail.com+copyright: Alexander Ignatyev+category: math+build-type: Simple+extra-source-files: hmatrix-morpheus/README.md hmatrix-morpheus/CHANGELOG.md+cabal-version: >=1.10+extra-source-files: src/types.h+ , src/simple_blas.h+ , src/matrix_reduce.h+ , src/activation.h++flag openblas+ description: Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries.+ default: False+ manual: True++library+ hs-source-dirs: hmatrix-morpheus/src+ exposed-modules: Numeric.Morpheus+ , Numeric.Morpheus.Activation+ , Numeric.Morpheus.MatrixReduce+ build-depends: base >= 4.7 && < 5+ , hmatrix+ default-language: Haskell2010+ default-extensions: ForeignFunctionInterface+ include-dirs: src+ C-sources: src/simple_blas.c+ , src/matrix_reduce.c+ , src/activation.c+ cc-options: -O3 -Wall -std=c99+ if arch(x86_64)+ cc-options: -msse2+ if arch(i386)+ cc-options: -msse2++ if os(OSX)+ frameworks: Accelerate+ cc-options: -DACCELERATE+ extra-lib-dirs: /opt/local/lib/+ include-dirs: /opt/local/include/+ extra-lib-dirs: /usr/local/lib/+ include-dirs: /usr/local/include/+ if arch(i386)+ cc-options: -arch i386++ if os(linux)+ if flag(openblas)+ extra-lib-dirs: /usr/lib/openblas/lib+ extra-libraries: openblas+ else+ extra-libraries: blas lapack++ if arch(x86_64)+ cc-options: -fPIC++test-suite hmatrix-morpheus-test+ type: exitcode-stdio-1.0+ hs-source-dirs: hmatrix-morpheus/test+ main-is: Main.hs+ other-modules: Test.HUnit.Approx+ , Test.HUnit.Plus+ , Numeric.Morpheus.ActivationTest+ , Numeric.Morpheus.MatrixReduceTest+ build-depends: base+ , hmatrix-morpheus+ , hmatrix+ , MonadRandom >= 0.4.2.3+ , test-framework >= 0.8.1.1+ , test-framework-hunit >= 0.3.0.2+ , test-framework-quickcheck2 >= 0.3.0.3+ , HUnit >= 1.3.1.1+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ default-language: Haskell2010++benchmark hmatrix-morpheus-bench+ type: exitcode-stdio-1.0+ hs-source-dirs: hmatrix-morpheus/bench+ main-is: Main.hs+ other-modules: Naive+ build-depends: base+ , hmatrix-morpheus+ , hmatrix+ , criterion+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ default-language: Haskell2010+ 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+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ default-language: Haskell2010++source-repository head+ type: git+ location: https://github.com/alexander-ignatyev/morpheus
+ hmatrix-morpheus/CHANGELOG.md view
@@ -0,0 +1,12 @@+# Changelog for hmatrix-morpheus++All notable changes to this project will be documented in this file.++The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/)+and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html).+++## [0.1.0.0] - 2017-08-26+### Added+- Numeric.Morpheus.Activation module - activation functions for neural network implementations;+- Numeric.Morpheus.MatrixReduce module - Matrix reduce by column/row operations.
+ hmatrix-morpheus/README.md view
@@ -0,0 +1,64 @@+# hmatrix-morpheus++Purely functional interface to morpheus based on hmatrix. Morpheus library contains a bunch of cache line aware numerical algorithms suitable for using as a low-level primitives to build machine learning applications.++Naive implementations, which do not take into account cache lines, would suffer order of magnitude performance degradation if tried to traverse a row major matrix in column order.++The library functions take into account CPU cache and work well for both row major and column major matrices, as you can see from the benchmark below which measure time for finding indices of max elements for every column and every row in comparison with naive approach:++```+| Benchmark | Morheus | Naive |++-------------------------------------+-----------+----------++| max index in columns - row major | 0.9077 ms | 12.93 ms |++-------------------------------------+-----------+----------++| max index in columns - column major | 0.8778 ms | 1.490 ms |++-------------------------------------+-----------+----------++| max index in rows - row major | 0.8622 ms | 1.504 ms |++-------------------------------------+-----------+----------++| max index in rows - column major | 0.8913 ms | 12.80 ms |++-------------------------------------+-----------+----------++```++## Examples++```haskell+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)+```
+ hmatrix-morpheus/bench/Main.hs view
@@ -0,0 +1,34 @@+module Main where++import Criterion.Main (defaultMain, bgroup, bench, nf, nfIO, whnf)+import qualified Numeric.LinearAlgebra as LA+import qualified Numeric.Morpheus as M+import qualified Naive as N++main :: IO ()+main = do+ !ar <- LA.rand 1000 1000+ let !ac = LA.tr ar++ defaultMain [+ bgroup "Activation" [+ bench "sigmoid" $ nf M.sigmoid ar+ , bench "sigmoid gradient" $ nf M.sigmoidGradient ar+ ]+ , bgroup "Matrix Reduce" [+ bench "sum by columns - row major" $ nf M.columnSum ar+ , bench "sum by columns - column major" $ nf M.columnSum ac+ , bench "sum by rows - row major" $ nf M.rowSum ar+ , bench "sum by rows - column major" $ nf M.rowSum ac+ , bench "max index in columns - row major" $ nf M.columnMaxIndex ar+ , bench "max index in columns - column major" $ nf M.columnMaxIndex ac+ , bench "max index in rows - row major" $ nf M.rowMaxIndex ar+ , bench "max index in rows - column major" $ nf M.rowMaxIndex ac+ ]+ , bgroup "Matrix Reduce (naive)" [+ bench "max index in columns - row major" $ nf (N.reduceByColumnsV (fromIntegral . LA.maxIndex) ) ar+ , bench "max index in columns - column major" $ nf (N.reduceByColumnsV (fromIntegral . LA.maxIndex) ) ac+ , bench "max index in rows - row major" $ nf (N.reduceByRowsV (fromIntegral . LA.maxIndex) ) ar+ , bench "max index in rows - column major" $ nf (N.reduceByRowsV (fromIntegral . LA.maxIndex) ) ac+ ]+ ]
+ hmatrix-morpheus/bench/Naive.hs view
@@ -0,0 +1,19 @@+module Naive+(+ reduceByRowsV+ , reduceByColumnsV+)++where++import qualified Numeric.LinearAlgebra as LA++type Matrix = LA.Matrix LA.R+type Vector = LA.Vector LA.R+type R = LA.R++reduceByRowsV :: (Vector -> R) -> Matrix -> Vector+reduceByRowsV f = LA.vector . map f . LA.toRows++reduceByColumnsV :: (Vector -> R) -> Matrix -> Vector+reduceByColumnsV f = LA.vector . map f . LA.toColumns
+ hmatrix-morpheus/examples/Main.hs view
@@ -0,0 +1,41 @@+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)
+ hmatrix-morpheus/src/Numeric/Morpheus.hs view
@@ -0,0 +1,19 @@+{-|+Module: Numeric.Morpheus+Description: Low-level machine learning auxiliary functions+Copyright: (c) Alexander Ignatyev, 2017+License: BSD-3+Stability: experimental+Portability: POSIX+-}++module Numeric.Morpheus+(+ module Numeric.Morpheus.Activation+ , module Numeric.Morpheus.MatrixReduce+)++where++import Numeric.Morpheus.Activation+import Numeric.Morpheus.MatrixReduce
+ hmatrix-morpheus/src/Numeric/Morpheus/Activation.hs view
@@ -0,0 +1,151 @@+{-|+Module: Numeric.Morpheus.Activation+Description: Activation Functions+Copyright: (c) Alexander Ignatyev, 2017+License: BSD-3+Stability: experimental+Portability: POSIX+-}++module Numeric.Morpheus.Activation+(+ sigmoid+ , sigmoidGradient+ , relu+ , reluGradient+ , tanh_+ , tanhGradient+)++where+++import Numeric.LinearAlgebra+import Numeric.LinearAlgebra.Devel+import System.IO.Unsafe(unsafePerformIO)+import Foreign.C.Types+import Foreign.Ptr(Ptr)+import Foreign.Storable+++createMatrixOfShape :: Storable a => Matrix a -> IO (Matrix a)+createMatrixOfShape m = createMatrix (orderOf m) (rows m) (cols m)+++{- morpheus_sigmoid -}+foreign import ccall unsafe "morpheus_sigmoid"+ c_morpheus_sigmoid :: CInt -> Ptr Double -> Ptr Double -> IO ()+++call_morpheus_sigmoid :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> IO ()+call_morpheus_sigmoid rows cols _ _ xPtr _ _ _ _ yPtr = do+ c_morpheus_sigmoid (rows*cols) xPtr yPtr+++-- | Calculates sigmoid.+sigmoid :: Matrix R -> Matrix R+sigmoid m = unsafePerformIO $ do+ y <- createMatrixOfShape m+ apply m (apply y id) call_morpheus_sigmoid+ return y+++{- morpheus_sigmoid_gradient -}+foreign import ccall unsafe "morpheus_sigmoid_gradient"+ c_morpheus_sigmoid_gradient :: CInt -> Ptr Double -> Ptr Double -> IO ()+++call_morpheus_sigmoid_gradient :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> IO ()+call_morpheus_sigmoid_gradient rows cols _ _ xPtr _ _ _ _ yPtr = do+ c_morpheus_sigmoid_gradient (rows*cols) xPtr yPtr+++-- | Calculates derivates of sigmoid.+sigmoidGradient :: Matrix R -> Matrix R+sigmoidGradient m = unsafePerformIO $ do+ y <- createMatrixOfShape m+ apply m (apply y id) call_morpheus_sigmoid_gradient+ return y+++{- morpheus_relu -}+foreign import ccall unsafe "morpheus_relu"+ c_morpheus_relu :: CInt -> Ptr Double -> Ptr Double -> IO ()+++call_morpheus_relu :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> IO ()+call_morpheus_relu rows cols _ _ xPtr _ _ _ _ yPtr = do+ c_morpheus_relu (rows*cols) xPtr yPtr+++-- | Calculates ReLu.+relu :: Matrix R -> Matrix R+relu m = unsafePerformIO $ do+ y <- createMatrixOfShape m+ apply m (apply y id) call_morpheus_relu+ return y+++{- morpheus_relu_gradient -}+foreign import ccall unsafe "morpheus_relu_gradient"+ c_morpheus_relu_gradient :: CInt -> Ptr Double -> Ptr Double -> IO ()+++call_morpheus_relu_gradient :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> IO ()+call_morpheus_relu_gradient rows cols _ _ xPtr _ _ _ _ yPtr = do+ c_morpheus_relu_gradient (rows*cols) xPtr yPtr+++-- | Calculates derivates of ReLu.+reluGradient :: Matrix R -> Matrix R+reluGradient m = unsafePerformIO $ do+ y <- createMatrixOfShape m+ apply m (apply y id) call_morpheus_relu_gradient+ return y++{- morpheus_tanh -}+foreign import ccall unsafe "morpheus_tanh"+ c_morpheus_tanh :: CInt -> Ptr Double -> Ptr Double -> IO ()+++call_morpheus_tanh :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> IO ()+call_morpheus_tanh rows cols _ _ xPtr _ _ _ _ yPtr = do+ c_morpheus_tanh (rows*cols) xPtr yPtr+++-- | Calculates tanh.+tanh_ :: Matrix R -> Matrix R+tanh_ m = unsafePerformIO $ do+ y <- createMatrixOfShape m+ apply m (apply y id) call_morpheus_tanh+ return y+++{- morpheus_tanh_gradient -}+foreign import ccall unsafe "morpheus_tanh_gradient"+ c_morpheus_tanh_gradient :: CInt -> Ptr Double -> Ptr Double -> IO ()+++call_morpheus_tanh_gradient :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> IO ()+call_morpheus_tanh_gradient rows cols _ _ xPtr _ _ _ _ yPtr = do+ c_morpheus_tanh_gradient (rows*cols) xPtr yPtr+++-- | Calculates derivates of tanh.+tanhGradient :: Matrix R -> Matrix R+tanhGradient m = unsafePerformIO $ do+ y <- createMatrixOfShape m+ apply m (apply y id) call_morpheus_tanh_gradient+ return y
+ hmatrix-morpheus/src/Numeric/Morpheus/MatrixReduce.hs view
@@ -0,0 +1,223 @@+{-|+Module: Numeric.Morpheus.MatrixReduce+Description: Matrix reduce by column/row operations.+Copyright: (c) Alexander Ignatyev, 2017+License: BSD-3+Stability: experimental+Portability: POSIX+-}++module Numeric.Morpheus.MatrixReduce+(+ columnPredicate+ , rowPredicate+ , columnSum+ , rowSum+ , columnMaxIndex+ , columnMinIndex+ , rowMaxIndex+ , rowMinIndex+)++where++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+foreign import ccall "wrapper" mkPredicate :: Predicate -> IO (FunPtr Predicate)+++{- morpheus_column_predicate -}+foreign import ccall safe "morpheus_column_predicate"+ c_morpheus_column_predicate :: FunPtr Predicate -> CInt -> CInt -> CInt+ -> Ptr R -> Ptr R -> IO ()+++call_morpheus_column_predicate :: FunPtr Predicate+ -> CInt -> CInt -> CInt -> CInt -> Ptr R+ -> CInt -> Ptr R+ -> IO ()+call_morpheus_column_predicate f rows cols xRow xCol matPtr _ vecPtr = do+ let layout = morpheusLayout xCol cols+ c_morpheus_column_predicate f layout rows cols matPtr vecPtr+++-- | Scan every column of the given matrix.+-- Predicate takes and accumulator and next value of the column, returns new accumulator.+-- Returns accumulator values for every column.+columnPredicate :: (R -> R -> R) -> Matrix R -> Vector R+columnPredicate f m = unsafePerformIO $ do+ v <- createVector (cols m)+ fpred <- mkPredicate f+ apply m (apply v id) (call_morpheus_column_predicate fpred)+ return v+++{- morpheus_row_predicate -}+foreign import ccall safe "morpheus_row_predicate"+ c_morpheus_row_predicate :: FunPtr Predicate -> CInt -> CInt -> CInt+ -> Ptr R -> Ptr R -> IO ()+++call_morpheus_row_predicate :: FunPtr Predicate+ -> CInt -> CInt -> CInt -> CInt -> Ptr R+ -> CInt -> Ptr R+ -> IO ()+call_morpheus_row_predicate f rows cols xRow xCol matPtr _ vecPtr = do+ let layout = morpheusLayout xCol cols+ c_morpheus_row_predicate f layout rows cols matPtr vecPtr+++-- | Scan every row of the given matrix.+-- Predicate takes and accumulator and next value of the row, returns new accumulator.+-- Returns accumulator values for every row.+rowPredicate :: (R -> R -> R) -> Matrix R -> Vector R+rowPredicate f m = unsafePerformIO $ do+ v <- createVector (rows m)+ fpred <- mkPredicate f+ apply m (apply v id) (call_morpheus_row_predicate fpred)+ return v+++{- morpheus_column_sum -}+foreign import ccall unsafe "morpheus_column_sum"+ c_morpheus_column_sum :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> IO ()+++call_morpheus_column_sum :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> Ptr Double+ -> IO ()+call_morpheus_column_sum rows cols xRow xCol matPtr _ vecPtr = do+ let layout = morpheusLayout xCol cols+ c_morpheus_column_sum layout rows cols matPtr vecPtr+++-- | Calculates sums of elements of every column of the given matrix+columnSum :: Matrix Double -> Vector Double+columnSum m = unsafePerformIO $ do+ v <- createVector (cols m)+ apply m (apply v id) call_morpheus_column_sum+ return v+++{- morpheus_row_sum -}+foreign import ccall unsafe "morpheus_row_sum"+ c_morpheus_row_sum :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> IO ()+++call_morpheus_row_sum :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> Ptr Double+ -> IO ()+call_morpheus_row_sum rows cols xRow xCol matPtr _ vecPtr = do+ let layout = morpheusLayout xCol cols+ c_morpheus_row_sum layout rows cols matPtr vecPtr+++-- | Calculates sums of elements of every row of the given matrix+rowSum :: Matrix Double -> Vector Double+rowSum m = unsafePerformIO $ do+ v <- createVector (rows m)+ apply m (apply v id) call_morpheus_row_sum+ return v+++{- morpheus_column_max_index -}+foreign import ccall unsafe "morpheus_column_max_index"+ c_morpheus_column_max_index :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> Ptr CInt -> IO ()+++call_morpheus_column_max_index :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> Ptr Double+ -> CInt -> Ptr CInt+ -> IO ()+call_morpheus_column_max_index rows cols xRow xCol matPtr _ vecPtr _ idxPtr = do+ let layout = morpheusLayout xCol cols+ c_morpheus_column_max_index layout rows cols matPtr vecPtr idxPtr+++-- | Finds maximum values and their indices of every column of the given matrix+columnMaxIndex :: Matrix Double -> (Vector R, Vector I)+columnMaxIndex m = unsafePerformIO $ do+ v <- createVector (cols m)+ i <- createVector (cols m)+ apply m (apply v (apply i id)) call_morpheus_column_max_index+ return (v, i)+++{- morpheus_column_min_index -}+foreign import ccall unsafe "morpheus_column_min_index"+ c_morpheus_column_min_index :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> Ptr CInt -> IO ()+++call_morpheus_column_min_index :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> Ptr Double+ -> CInt -> Ptr CInt+ -> IO ()+call_morpheus_column_min_index rows cols xRow xCol matPtr _ vecPtr _ idxPtr = do+ let layout = morpheusLayout xCol cols+ c_morpheus_column_min_index layout rows cols matPtr vecPtr idxPtr+++-- | Finds minimum values and their indices of every column of the given matrix+columnMinIndex :: Matrix Double -> (Vector R, Vector I)+columnMinIndex m = unsafePerformIO $ do+ v <- createVector (cols m)+ i <- createVector (cols m)+ apply m (apply v (apply i id)) call_morpheus_column_min_index+ return (v, i)+++{- morpheus_row_max_index -}+foreign import ccall unsafe "morpheus_row_max_index"+ c_morpheus_row_max_index :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> Ptr CInt -> IO ()+++call_morpheus_row_max_index :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> Ptr Double+ -> CInt -> Ptr CInt+ -> IO ()+call_morpheus_row_max_index rows cols xRow xCol matPtr _ vecPtr _ idxPtr = do+ let layout = morpheusLayout xCol cols+ c_morpheus_row_max_index layout rows cols matPtr vecPtr idxPtr+++-- | Finds maximum values and their indices of every row of the given matrix+rowMaxIndex :: Matrix Double -> (Vector R, Vector I)+rowMaxIndex m = unsafePerformIO $ do+ v <- createVector (rows m)+ i <- createVector (rows m)+ apply m (apply v (apply i id)) call_morpheus_row_max_index+ return (v, i)+++{- morpheus_row_min_index -}+foreign import ccall unsafe "morpheus_row_min_index"+ c_morpheus_row_min_index :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> Ptr CInt -> IO ()+++call_morpheus_row_min_index :: CInt -> CInt -> CInt -> CInt -> Ptr Double+ -> CInt -> Ptr Double+ -> CInt -> Ptr CInt+ -> IO ()+call_morpheus_row_min_index rows cols xRow xCol matPtr _ vecPtr _ idxPtr = do+ let layout = morpheusLayout xCol cols+ c_morpheus_row_min_index layout rows cols matPtr vecPtr idxPtr+++-- | Finds minimum values and their indices of every row of the given matrix+rowMinIndex :: Matrix Double -> (Vector R, Vector I)+rowMinIndex m = unsafePerformIO $ do+ v <- createVector (rows m)+ i <- createVector (rows m)+ apply m (apply v (apply i id)) call_morpheus_row_min_index+ return (v, i)
+ hmatrix-morpheus/test/Main.hs view
@@ -0,0 +1,10 @@+import Test.Framework (defaultMain, testGroup)++import qualified Numeric.Morpheus.MatrixReduceTest as MatrixReduce+import qualified Numeric.Morpheus.ActivationTest as Activation++main = defaultMain tests++tests = [ testGroup "Morpheus.MatrixReduce" MatrixReduce.tests+ , testGroup "Morpheus.Activation" Activation.tests+ ]
+ hmatrix-morpheus/test/Numeric/Morpheus/ActivationTest.hs view
@@ -0,0 +1,38 @@+module Numeric.Morpheus.ActivationTest+(+ tests+)++where++import Test.Framework (testGroup)+import Test.Framework.Providers.HUnit+import Test.HUnit+import Test.HUnit.Plus++import Numeric.LinearAlgebra+import Numeric.Morpheus.Activation+++x = row [10e10, 0, -10e10]+x_sigmoid = row [1, 0.5, 0]+x_sigmoid_gradient = row [0, 0.25, 0]+x_relu = row [10e10, 0, 0]+x_relu_gradient = row [1, 0, 0]+x_tanh = row [1, 0, -1]+x_tanh_gradient = row [0, 1, 0]+++tests = [ testGroup "Sigmoid" [+ testCase "Sigmoid" $ assertMatrix "" 1e-10 x_sigmoid (sigmoid x)+ , testCase "Sigmoid Gradient" $ assertMatrix "" 1e-10 x_sigmoid_gradient (sigmoidGradient x)+ ]+ , testGroup "ReLu" [+ testCase "ReLu" $ assertMatrix "" 1e-10 x_relu (relu x)+ , testCase "ReLu Gradient" $ assertMatrix "" 1e-10 x_relu_gradient (reluGradient x)+ ]+ , testGroup "Tanh" [+ testCase "Tanh" $ assertMatrix "" 1e-10 x_tanh (tanh_ x)+ , testCase "Tanh Gradient" $ assertMatrix "" 1e-10 x_tanh_gradient (tanhGradient x)+ ]+ ]
+ hmatrix-morpheus/test/Numeric/Morpheus/MatrixReduceTest.hs view
@@ -0,0 +1,90 @@+module Numeric.Morpheus.MatrixReduceTest+(+ tests+)++where++import Test.Framework (testGroup)+import Test.Framework.Providers.HUnit+import Test.HUnit+import Test.HUnit.Plus++import Numeric.LinearAlgebra+import Numeric.Morpheus.MatrixReduce+++xRM = matrix 5 [+ 71, 11, 3, 9, 7,+ 21, 7, 2, 23, 11,+ 11, 32, 53, 49, 37,+ 1, 24, 78, 90, 17+ ]++xCM = tr xRM++columnSumRM = vector [104, 74, 136, 171, 72]+columnSumCM = vector [101, 64,182, 210]++rowSumRM = vector [101, 64,182, 210]+rowSumCM = vector [104, 74, 136, 171, 72]++columnMaxValuesRM = fromList [71, 32, 78, 90, 37]+columnMaxIndicesRM = fromList [0, 2, 3, 3, 2]+columnMaxValuesCM = fromList [71, 23, 53, 90]+columnMaxIndicesCM = fromList [0, 3, 2, 3]++columnMinValuesRM = fromList [1, 7, 2, 9, 7]+columnMinIndicesRM = fromList [3, 1, 1, 0, 0]+columnMinValuesCM = fromList [3, 2, 11, 1]+columnMinIndicesCM = fromList [2, 2, 0, 0]++rowMaxValuesRM = fromList [71, 23, 53, 90]+rowMaxIndicesRM = fromList [0, 3, 2, 3]+rowMaxValuesCM = fromList [71, 32, 78, 90, 37]+rowMaxIndicesCM = fromList [0, 2, 3, 3, 2]++rowMinValuesRM = fromList [3, 2, 11, 1]+rowMinIndicesRM = fromList [2, 2, 0, 0]+rowMinValuesCM = fromList [1, 7, 2, 9, 7]+rowMinIndicesCM = fromList [3, 1, 1, 0, 0]+++assertValuesAndIndices (expectedR, expectedI) (actualR, actualI) = do+ assertVector "values" 1e-10 expectedR actualR+ assertVectorI "indices" expectedI actualI+++tests = [testGroup "Column Sum" [+ testCase "row major" $ assertVector "" 1e-10 columnSumRM (columnSum xRM)+ , testCase "column major" $ assertVector "" 1e-10 columnSumCM (columnSum xCM)+ ]+ , testGroup "Row Sum" [+ testCase "row major" $ assertVector "" 1e-10 rowSumRM (rowSum xRM)+ , testCase "column major" $ assertVector "" 1e-10 rowSumCM (rowSum xCM)+ ]+ , testGroup "Column Max Index" [+ testCase "row major" $ assertValuesAndIndices (columnMaxValuesRM, columnMaxIndicesRM) (columnMaxIndex xRM)+ , testCase "column major" $ assertValuesAndIndices (columnMaxValuesCM, columnMaxIndicesCM) (columnMaxIndex xCM)+ ]+ , testGroup "Column Min Index" [+ testCase "row major" $ assertValuesAndIndices (columnMinValuesRM, columnMinIndicesRM) (columnMinIndex xRM)+ , testCase "column major" $ assertValuesAndIndices (columnMinValuesCM, columnMinIndicesCM) (columnMinIndex xCM)+ ]+ , testGroup "Row Max Index" [+ testCase "row major" $ assertValuesAndIndices (rowMaxValuesRM, rowMaxIndicesRM) (rowMaxIndex xRM)+ , testCase "column major" $ assertValuesAndIndices (rowMaxValuesCM, rowMaxIndicesCM) (rowMaxIndex xCM)+ ]+ , testGroup "Row Min Index" [+ testCase "row major" $ assertValuesAndIndices (rowMinValuesRM, rowMinIndicesRM) (rowMinIndex xRM)+ , testCase "column major" $ assertValuesAndIndices (rowMinValuesCM, rowMinIndicesCM) (rowMinIndex xCM)+ ]+ , testGroup "Column Predicate" [+ testCase "max - row major" $ assertVector "" 1e-10 columnMaxValuesRM (columnPredicate max xRM)+ , testCase "max - col major" $ assertVector "" 1e-10 columnMaxValuesCM (columnPredicate max xCM)+ ]+ , testGroup "Row Predicate" [+ testCase "min - row major" $ assertVector "" 1e-10 rowMinValuesRM (rowPredicate min xRM)+ , testCase "min - col major" $ assertVector "" 1e-10 rowMinValuesCM (rowPredicate min xCM)+ ]+ ]
+ hmatrix-morpheus/test/Test/HUnit/Approx.hs view
@@ -0,0 +1,85 @@+{-# LANGUAGE ImplicitParams, CPP #-}+#if __GLASGOW_HASKELL__ >= 707+{-# LANGUAGE Safe #-} -- Test.HUnit is not Safe in 7.6 and below+#endif+-----------------------------------------------------------------------------+-- |+-- Module : Test.HUnit.Approx+-- Copyright : (C) 2014 Richard Eisenberg+-- License : BSD-style (see LICENSE)+-- Maintainer : Richard Eisenberg (eir@cis.upenn.edu)+-- Stability : intended to be stable+-- Portability : not portable (uses implicit parameters)+--+-- This module exports combinators to allow approximate equality of+-- floating-point values in HUnit tests.+-----------------------------------------------------------------------------++module Test.HUnit.Approx (+ -- * Assertions+ assertApproxEqual, (@~?), (@?~),++ -- * Tests+ (~~?), (~?~)+ ) where++import Test.HUnit+import Control.Monad ( unless )++-- | Asserts that the specified actual value is approximately equal to the+-- expected value. The output message will contain the prefix, the expected+-- value, the actual value, and the maximum margin of error.+-- +-- If the prefix is the empty string (i.e., @\"\"@), then the prefix is omitted+-- and only the expected and actual values are output.+assertApproxEqual :: (Ord a, Num a, Show a)+ => String -- ^ The message prefix+ -> a -- ^ Maximum allowable margin of error+ -> a -- ^ The expected value + -> a -- ^ The actual value+ -> Assertion+assertApproxEqual preface epsilon expected actual =+ unless (abs (actual - expected) <= epsilon) (assertFailure msg)+ where msg = (if null preface then "" else preface ++ "\n") +++ "expected: " ++ show expected ++ "\n but got: " ++ show actual +++ "\n (maximum margin of error: " ++ show epsilon ++ ")"++-- | Asserts that the specified actual value is approximately equal to the+-- expected value (with the expected value on the right-hand side). The margin+-- of error is specified with the implicit parameter @epsilon@.+(@?~) :: (Ord a, Num a, Show a, ?epsilon :: a)+ => a -- ^ The actual value+ -> a -- ^ The expected value+ -> Assertion+x @?~ y = assertApproxEqual "" ?epsilon y x+infix 1 @?~++-- | Asserts that the specified actual value is approximately equal to the+-- expected value (with the expected value on the left-hand side). The margin+-- of error is specified with the implicit parameter @epsilon@.+(@~?) :: (Ord a, Num a, Show a, ?epsilon :: a)+ => a -- ^ The expected value+ -> a -- ^ The actual value+ -> Assertion+x @~? y = assertApproxEqual "" ?epsilon x y+infix 1 @~?++-- | Shorthand for a test case that asserts approximate equality (with the+-- expected value on the left-hand side, and the actual value on the+-- right-hand side).+(~~?) :: (Ord a, Num a, Show a, ?epsilon :: a)+ => a -- ^ The expected value+ -> a -- ^ The actual value+ -> Test+expected ~~? actual = TestCase (expected @~? actual)+infix 1 ~~?++-- | Shorthand for a test case that asserts approximate equality (with the+-- actual value on the left-hand side, and the expected value on the+-- right-hand side).+(~?~) :: (Ord a, Num a, Show a, ?epsilon :: a)+ => a -- ^ The actual value+ -> a -- ^ The expected value + -> Test+actual ~?~ expected = TestCase (actual @?~ expected)+infix 1 ~?~
+ hmatrix-morpheus/test/Test/HUnit/Plus.hs view
@@ -0,0 +1,44 @@+module Test.HUnit.Plus+(+ assertMaybeDouble+ , assertOnFunction+ , assertVector+ , assertVectorI+ , assertMatrix+)+where++import Control.Monad+import Test.HUnit+import Test.HUnit.Approx+import Numeric.LinearAlgebra++assertMaybeDouble :: Maybe Double -> Maybe Double -> Double -> Assertion+assertMaybeDouble Nothing Nothing _ = assertString ""+assertMaybeDouble expected Nothing _ = assertString msg+ where msg = "expected: " ++ show expected ++ "\nbut got: Nothing"+assertMaybeDouble Nothing actual _ = assertString msg+ where msg = "expected: Nothing\nbit got: " ++ show actual+assertMaybeDouble (Just expected) (Just actual) eps = assertApproxEqual "Maybe Double" eps expected actual++assertOnFunction :: (Eq b, Show b) => (a -> b) -> a -> a -> Assertion+assertOnFunction func expected actual = func expected @=? func actual++assertVector :: String -> R -> Vector R -> Vector R -> Assertion+assertVector message eps expected actual =+ let diff = norm_2 (expected-actual)+ msg = message ++ "\nexpected: " ++ show expected ++ "\nbut got: " ++ show actual+ in unless (diff < eps) (assertFailure msg)++assertVectorI :: String -> Vector I -> Vector I -> Assertion+assertVectorI message expected actual =+ let diff = norm_2 (expected-actual)+ msg = message ++ "\nexpected: " ++ show expected ++ "\nbut got: " ++ show actual+ in unless (diff == 0) (assertFailure msg)++assertMatrix :: String -> R -> Matrix R -> Matrix R -> Assertion+assertMatrix message eps expected actual =+ let diff = norm_2 (expected-actual)+ msg = message ++ "\nexpected: " ++ show expected ++ "\nbut got: " ++ show actual+ in unless (diff < eps) (assertFailure msg)+
+ src/activation.c view
@@ -0,0 +1,41 @@+#include "activation.h"++#include <math.h>++void morpheus_sigmoid(int n, const double *x, double *y) {+ for (int i = 0; i < n; ++i) {+ y[i] = 1.0/(1.0+exp(-x[i]));+ }+}++void morpheus_sigmoid_gradient(int n, const double *x, double *y) {+ for (int i = 0; i < n; ++i) {+ const double s = 1.0/(1.0+exp(-x[i]));+ y[i] = s * (1-s);+ }+}++void morpheus_relu(int n, const double *x, double *y) {+ for (int i = 0; i < n; ++i) {+ y[i] = x[i] > 0 ? x[i] : 0;+ }+}++void morpheus_relu_gradient(int n, const double *x, double *y) {+ for (int i = 0; i < n; ++i) {+ y[i] = x[i] > 0 ? 1 : 0;+ }+}++void morpheus_tanh(int n, const double *x, double *y) {+ for (int i = 0; i < n; ++i) {+ y[i] = tanh(x[i]);+ }+}++void morpheus_tanh_gradient(int n, const double *x, double *y) {+ for (int i = 0; i < n; ++i) {+ const double t = tanh(x[i]);+ y[i] = 1 - t*t;+ }+}
+ src/activation.h view
@@ -0,0 +1,33 @@+/*! \file activation.h+ \brief Activation functions+*/++#ifndef MORHEPUS_ACTIVATION_H+#define MORHEPUS_ACTIVATION_H+++/*! \defgroup activation activation+ \brief Activation functions+ \{+*/++/*! Calculates sigmoid. */+void morpheus_sigmoid(int n, const double *x, double *y);++/*! Calculates derivatives of sigmoid. */+void morpheus_sigmoid_gradient(int n, const double *x, double *y);++/*! Calculates ReLu. */+void morpheus_relu(int n, const double *x, double *y);++/*! Calculates derivatives of ReLu. */+void morpheus_relu_gradient(int n, const double *x, double *y);++/*! Calculates tanh. */+void morpheus_tanh(int n, const double *x, double *y);++/*! Calculates derivatives of tanh. */+void morpheus_tanh_gradient(int n, const double *x, double *y);+/*! \} */++#endif /* MORHEPUS_ACTIVATION_H */
+ src/matrix_reduce.c view
@@ -0,0 +1,267 @@+#include "matrix_reduce.h"++#include <math.h>+#include <stdio.h>+#include <assert.h>++void morpheus_column_predicate(double (*pred)(double, double),+ morpheus_layout_e layout,+ int nrows,+ int ncols,+ const double *x,+ double *y) {+ if (layout == morpheus_row_major) {+ /* first row is initial values */+ morpheus_dcopy(ncols, x, y);+ for (int i = 1; i < nrows; ++i) {+ const double *row = x + i*ncols;+ for (int j = 0; j < ncols; ++j) {+ y[j] = pred(y[j], row[j]);+ }+ }+ } else if (layout == morpheus_col_major) {+ for (int j = 0; j < ncols; ++j) {+ const double *column = &x[j*nrows];+ y[j] = column[0];+ for (int i = 1; i < nrows; ++i) {+ y[j] = pred(y[j], column[i]);+ }+ }+ } else {+ perror("unkwnown matrix layout");+ assert(0);+ }+}++void morpheus_row_predicate(double (*pred)(double, double),+ morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y) {+ if (layout == morpheus_row_major) {+ for (int i = 0; i < nrows; ++i) {+ const double *row = x + i*ncols;+ y[i] = row[0];+ for (int j = 1; j < ncols; ++j) {+ y[i] = pred(y[i], row[j]);+ }+ }+ } else if (layout == morpheus_col_major) {+ /* first column is initial values */+ morpheus_dcopy(nrows, x, y);+ for (int j = 1; j < ncols; ++j) {+ const double *column = x + j*nrows;+ for (int i = 0; i < nrows; ++i) {+ y[i] = pred(y[i], column[i]);+ }+ }+ } else {+ perror("unkwnown matrix layout");+ assert(0);+ }+}++void morpheus_column_max_index(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y,+ int *is) {+ if (layout == morpheus_row_major) {+ for (int j = 0; j < ncols; ++j) {+ is[j] = 0;+ y[j] = x[j];+ }+ for (int i = 1; i < nrows; ++i) {+ const double *row = x + i*ncols;+ for (int j = 0; j < ncols; ++j) {+ if (y[j] < row[j]) {+ y[j] = row[j];+ is[j] = i;+ }+ }+ }+ } else if (layout == morpheus_col_major) {+ for (int j = 0; j < ncols; ++j) {+ const double *column = &x[j*nrows];+ y[j] = column[0];+ is[j] = 0;+ for (int i = 1; i < nrows; ++i) {+ if (y[j] < column[i]) {+ y[j] = column[i];+ is[j] = i;+ }+ }+ }+ } else {+ perror("unkwnown matrix layout");+ assert(0);+ }+}+++void morpheus_column_min_index(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y,+ int *is) {+ if (layout == morpheus_row_major) {+ for (int j = 0; j < ncols; ++j) {+ is[j] = 0;+ y[j] = x[j];+ }+ for (int i = 1; i < nrows; ++i) {+ const double *row = x + i*ncols;+ for (int j = 0; j < ncols; ++j) {+ if (y[j] > row[j]) {+ y[j] = row[j];+ is[j] = i;+ }+ }+ }+ } else if (layout == morpheus_col_major) {+ for (int j = 0; j < ncols; ++j) {+ const double *column = &x[j*nrows];+ y[j] = column[0];+ is[j] = 0;+ for (int i = 1; i < nrows; ++i) {+ if (y[j] > column[i]) {+ y[j] = column[i];+ is[j] = i;+ }+ }+ }+ } else {+ perror("unkwnown matrix layout");+ assert(0);+ }+}++void morpheus_row_max_index(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y,+ int *is) {+ if (layout == morpheus_row_major) {+ for (int i = 0; i < nrows; ++i) {+ const double *row = x + i*ncols;+ y[i] = row[0];+ is[i] = 0;+ for (int j = 1; j < ncols; ++j) {+ if (y[i] < row[j]) {+ y[i] = row[j];+ is[i] = j;+ }+ }+ }+ } else if (layout == morpheus_col_major) {+ for (int i = 0; i < nrows; ++i) {+ is[i] = 0;+ y[i] = x[i];+ }+ for (int j = 1; j < ncols; ++j) {+ const double *column = x + j*nrows;+ for (int i = 0; i < nrows; ++i) {+ if (y[i] < column[i]) {+ y[i] = column[i];+ is[i] = j;+ }+ }+ }+ } else {+ perror("unkwnown matrix layout");+ assert(0);+ }+}++void morpheus_row_min_index(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y,+ int *is) {+ if (layout == morpheus_row_major) {+ for (int i = 0; i < nrows; ++i) {+ const double *row = x + i*ncols;+ y[i] = row[0];+ is[i] = 0;+ for (int j = 1; j < ncols; ++j) {+ if (y[i] > row[j]) {+ y[i] = row[j];+ is[i] = j;+ }+ }+ }+ } else if (layout == morpheus_col_major) {+ for (int i = 0; i < nrows; ++i) {+ is[i] = 0;+ y[i] = x[i];+ }+ for (int j = 1; j < ncols; ++j) {+ const double *column = x + j*nrows;+ for (int i = 0; i < nrows; ++i) {+ if (y[i] > column[i]) {+ y[i] = column[i];+ is[i] = j;+ }+ }+ }+ } else {+ perror("unkwnown matrix layout");+ assert(0);+ }+}++void morpheus_column_sum(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y) {+ if (layout == morpheus_row_major) {+ for (int j = 0; j < ncols; ++j) {+ y[j] = 0;+ }+ for (int i = 0; i < nrows; ++i) {+ const double *row = x + i*ncols;+ for (int j = 0; j < ncols; ++j) {+ y[j] += row[j];+ }+ }+ } else if (layout == morpheus_col_major) {+ for (int j = 0; j < ncols; ++j) {+ const double *column = &x[j*nrows];+ y[j] = 0;+ for (int i = 0; i < nrows; ++i) {+ y[j] += column[i];+ }+ }+ } else {+ perror("unkwnown matrix layout");+ assert(0);+ }+}++void morpheus_row_sum(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y) {+ if (layout == morpheus_row_major) {+ for (int i = 0; i < nrows; ++i) {+ const double *row = x + i*ncols;+ y[i] = 0;+ for (int j = 0; j < ncols; ++j) {+ y[i] += row[j];+ }+ }+ } else if (layout == morpheus_col_major) {+ for (int i = 0; i < nrows; ++i) {+ y[i] = 0;+ }+ for (int j = 0; j < ncols; ++j) {+ const double *column = x + j*nrows;+ for (int i = 0; i < nrows; ++i) {+ y[i] += column[i];+ }+ }+ } else {+ perror("unkwnown matrix layout");+ assert(0);+ }+}
+ src/matrix_reduce.h view
@@ -0,0 +1,78 @@+/*! \file matrix_reduce.h+ \brief Matrix reduce by column/row operations.+*/++#ifndef MORPHEPUS_MATRIX_REDUCE_H+#define MORPHEPUS_MATRIX_REDUCE_H++#include "simple_blas.h"++/*! \defgroup matrix_reduce matrix_reduce+ \brief Matrix reduce by column/row operations.+ \{+*/++/*! Scan every column of the matrix x.+ Pred takes and accumulator and next value of the column, returns new accumulator.+ Returns accumulator values for every column.+ */+void morpheus_column_predicate(double (*pred)(double, double),+ morpheus_layout_e layout,+ int nrows,+ int ncols,+ const double *x,+ double *y /*!< vector of size ncols */);++/*! Scan every row of the matrix x.+ Pred takes and accumulator and next value of the row, returns new accumulator.+ Returns accumulator values for every row.+*/+void morpheus_row_predicate(double (*pred)(double, double),+ morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y /*!< vector of size ncols */);+++/*! Finds maximum values and their indices of every column of the matrix x */+void morpheus_column_max_index(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y, /*!< vector of maximum values of size ncols */+ int *is /*!< vector of indices of size ncols */);++/*! Finds minimum values and their indices of every column of the matrix x */+void morpheus_column_min_index(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y, /*!< vector of minimum values of size ncols */+ int *is /*!< vector of indices of size ncols */);++/*! Finds maximum values and their indices of every row of the matrix x */+void morpheus_row_max_index(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y, /*!< vector of maximum values of size nrows */+ int *is /*!< vector of indices of size nrows */);++/*! Finds minimum values and their indices of every row of the matrix x */+void morpheus_row_min_index(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y, /*!< vector of minimum values of size nrows */+ int *is /*!< vector of indices of size nrows */);++/*! Calculates sums of elements of every column of the matrix x */+void morpheus_column_sum(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y /*!< vector of sums of size ncols */);++/*! Calculates sums of elements of every row of the matrix x */+void morpheus_row_sum(morpheus_layout_e layout,+ int nrows, int ncols,+ const double *x,+ double *y /*!< vector of sums of size nrows */);+/*! \} */++#endif /* MORPHEPUS_MATRIX_REDUCE_H */
+ src/simple_blas.c view
@@ -0,0 +1,129 @@+#include "simple_blas.h"++#ifdef ACCELERATE+#include <Accelerate/Accelerate.h>+#else+#include <cblas.h>+#include "lapack.h"+#endif++int to_cblas_transpose(const morpheus_transpose_e trans) {+ switch (trans) {+ case morpheus_no_trans:+ return CblasNoTrans;+ case morpheus_trans:+ return CblasTrans;+ }+ return 0;+}++int to_cblas_layout(const morpheus_layout_e layout) {+ switch (layout) {+ case morpheus_row_major:+ return CblasRowMajor;+ case morpheus_col_major:+ return CblasColMajor;+ }+ return 0;+}++void morpheus_dscal(const int n, const double alpha, double *x) {+ const int inc = 1;+ cblas_dscal(n, alpha, x, inc);+}++void morpheus_dcopy(const int n, const double *x, double *y) {+ const int inc = 1;+ cblas_dcopy(n, x, inc, y, inc);+}++void morpheus_dcopy_scalar(const int n, const double alpha, double *x) {+ cblas_dcopy(n, &alpha, 0, x, 1);+}++void morpheus_dgemv(const morpheus_layout_e layout,+ const morpheus_transpose_e trans,+ const int nrows, const int ncols,+ const double alpha, const double *A, const double *x,+ const double beta, double *y) {+ const int inc = 1;+ const int lda = layout == morpheus_row_major ? ncols : nrows;+ cblas_dgemv(to_cblas_layout(layout), to_cblas_transpose(trans), nrows, ncols,+ alpha, A, lda, x, inc, beta, y, inc);+}++void morpheus_daxpy(const int n, const double alpha,+ const double *x, double *y) {+ const int inc = 1;+ cblas_daxpy(n, alpha, x, inc, y, inc);+}++double morpheus_dnrm2(const int n, const double *x) {+ const int inc = 1;+ return cblas_dnrm2(n, x, inc);+}++double morpheus_ddot(const int n, const double *x, const double *y) {+ const int inc = 1;+ return cblas_ddot(n, x, inc, y, inc);+}++int morpheus_inverse(int n, double *x, int *pivot, double *workspace) {+ workspace = workspace+n;+ int lwork = n*n;++ int rc;+ dgetrf_(&n, &n, x, &n, pivot, &rc);+ if (rc != 0) {+ return rc;+ }+ dgetri_(&n, x, &n, pivot, workspace, &lwork, &rc);+ return rc;+}++void morpheus_dger(const morpheus_layout_e layout,+ const int nrows,+ const int ncols,+ const double alpha,+ const double *x,+ const double *y,+ double *a) {+ const int inc = 1;+ const int lda = layout == morpheus_row_major ? ncols : nrows;+ cblas_dger(to_cblas_layout(layout), nrows, ncols,+ alpha, x, inc, y, inc, a, lda);+}++void morpheus_dgemm(const morpheus_layout_e layout,+ const morpheus_transpose_e trans_a,+ const morpheus_transpose_e trans_b,+ const int m,+ const int n,+ const int k,+ const double alpha,+ const double *a,+ const double *b,+ const double beta,+ double *c) {+ int lda, ldb, ldc;+ if (layout == morpheus_row_major) {+ lda = trans_a == morpheus_no_trans ? k : m;+ ldb = trans_b == morpheus_no_trans ? n : k;+ ldc = n;+ } else {+ lda = trans_a == morpheus_no_trans ? m : k;+ ldb = trans_b == morpheus_no_trans ? k : n;+ ldc = m;+ }+ cblas_dgemm(to_cblas_layout(layout),+ to_cblas_transpose(trans_a), to_cblas_transpose(trans_b),+ m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);+}++void morpheus_identity(const int n, double *a) {+ const int m = n*n;+ const double zero = 0;+ const double one = 1;+ cblas_dcopy(m, &zero, 0, a, 1);+ cblas_dcopy(n, &one, 0, a, n+1);+}
+ src/simple_blas.h view
@@ -0,0 +1,123 @@+/*! \file simple_blas.h+ \brief Simple Wrappers around BLAS functions+*/++#ifndef MORHEPUS_SIMPLE_BLAS_H+#define MORHEPUS_SIMPLE_BLAS_H+++/*! \defgroup simple_blas simple_blas+ \brief Simple Wrappers around BLAS functions+ \{+*/+++/*! Whether to transpose matrix */+typedef enum {+ morpheus_no_trans = 11,+ morpheus_trans+} morpheus_transpose_e;+++/*! Matrix Layout */+typedef enum {+ morpheus_row_major = 21, /*!< C-style layout */+ morpheus_col_major /*!< Fortran-style layout */+} morpheus_layout_e;+++/*! Scales a vector by a constant. */+void morpheus_dscal(const int n, const double alpha, double *x);+++/*! Copies a vector, x, to a vector, y. */+void morpheus_dcopy(const int n, const double *x, double *y);+++/*! Fills vector x with values of scalar alpha. */+void morpheus_dcopy_scalar(const int n, const double alpha, double *x);++++/*! Performs one of the matrix-vector operations.+ y = alpha*A*x + beta*y, or y = alpha*A**T*x + beta*y,++ where alpha and beta are scalars, x and y are vectors and A is an+ nrows by ncols matrix.++ A = nrows x ncols;+ x = trans == no ? ncols : nrows;+ y = trans == no ? nrows : ncols;+*/+void morpheus_dgemv(const morpheus_layout_e layout,+ const morpheus_transpose_e trans, /*!< Specifies whether to transpose matrix A */+ const int nrows,+ const int ncols,+ const double alpha,+ const double *A,+ const double *x,+ const double beta,+ double *y);++/*! Constant times a vector plus a vector.+ y = alpha*x + y;+*/+void morpheus_daxpy(const int n, /*!< number of elements in vectors x and y */+ const double alpha, /*!< scalar alpha */+ const double *x, /*!< vector x */+ double *y /*!< vector y */);+++/*! Returns the euclidean norm of a vector++ sqrt( x'*x )+*/+double morpheus_dnrm2(const int n, const double *x);+++/*! Forms the dot product of two vectors. */+double morpheus_ddot(const int n, const double *x, const double *y);+++/*! Compute the inverse of a given square matrix. */+int morpheus_inverse(int n, /*!< size of the matrix to inverse */+ double *x, /*!< matrix to inverse */+ int *pivot, /*!< vector of size n */+ double *workspace /*!< vector of size of n*n */);++/*!+ A = alpha * x * yT + A+*/+void morpheus_dger(const morpheus_layout_e layout,+ const int nrows,+ const int ncols,+ const double alpha,+ const double *x,+ const double *y,+ double *a);++/*!+ C = alpha * op(A) * op(B) + beta*C;++ op(A) - m x k;+ op(B) - k x n;+ C - m x n.+*/+void morpheus_dgemm(const morpheus_layout_e layout,+ const morpheus_transpose_e trans_a,+ const morpheus_transpose_e trans_b,+ const int m,+ const int n,+ const int k,+ const double alpha,+ const double *a,+ const double *b,+ const double beta,+ double *c);++/*! Creates identity matrix */+void morpheus_identity(const int n, double *a);++/*! \} */++#endif /* MORHEPUS_SIMPLE_BLAS_H */
+ src/types.h view
@@ -0,0 +1,54 @@+/*! \file types.h+ \brief Common types definitions+*/++#ifndef MORHEPUS_TYPES_H+#define MORHEPUS_TYPES_H++#include "simple_blas.h"+#include "regularization.h"++/*! \defgroup types types+ \brief Common types definitions+ \{+*/++/*! \brief Learing Data structure */+typedef struct {+ const double *x; /*!< (num_examples x num_features) matrix */+ const double *y; /*!< (num_examples) vector */+ int num_features; /*!< number of features including bias dimension */+ int num_examples; /*!< number of examples */+ morpheus_layout_e layout; /*!< layout of matrix x */+} morpheus_data_t;+++/*! Type definition of pointer to cost function. */+typedef double (*morpheus_cost_f)(const morpheus_reg_t *reg,+ const morpheus_data_t *data,+ const double *theta,+ double *tmp_buffer);+++/*! Type definition of pointer to gradient function. */+typedef void (*morpheus_gradient_f)(const morpheus_reg_t *reg,+ const morpheus_data_t *data,+ const double *theta,+ double *grad,+ double *tmp_buffer);+++typedef double (*morpheus_f)(const double *x, void *params);++typedef void (*morpheus_df)(const double *x, void *params, double *grad);+++typedef struct {+ morpheus_data_t data;+ morpheus_reg_t reg;+ double *memory_buffer;+} morpheus_params_t;++/*! \} */++#endif /* MORHEPUS_TYPES_H */