rc 0.1.0.1 → 0.3.0.0
raw patch · 7 files changed
+238/−101 lines, 7 filesdep +lineardep ~Learningdep ~dde
Dependencies added: linear
Dependency ranges changed: Learning, dde
Files
- README.md +3/−0
- examples/NTC/Main.hs +13/−9
- rc.cabal +13/−8
- rc/RC/Helpers.hs +13/−1
- rc/RC/NTC.hs +104/−83
- rc/RC/NTC/Reservoir.hs +58/−0
- rc/RC/NTC/Types.hs +34/−0
README.md view
@@ -33,6 +33,9 @@ ## Further reading +* Appeltant, L., et al. “Information Processing Using a Single+ Dynamical Node as Complex System.” Nature Communications, vol. 2,+ 2011, p. 468., doi:10.1038/ncomms1476. * Larger, L., et al. “Photonic Information Processing beyond Turing: an Optoelectronic Implementation of Reservoir Computing.” Optics Express, vol. 20, no. 3, 2012, p. 3241., doi:10.1364/oe.20.003241. * Rabinovič, Mihail Izrailevič, et al. Principles of Brain Dynamics: Global State Interactions. The MIT Press, 2012. * Jaeger, H. “Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication.” Science, vol. 304, no. 5667, Feb. 2004, pp. 78–80., doi:10.1126/science.1091277.
examples/NTC/Main.hs view
@@ -5,16 +5,23 @@ import RC.NTC as RC -- Get training data-takePast :: Int -> Matrix Double -> Matrix Double+takePast :: Element a => Int -> Matrix a -> Matrix a takePast horizon xs = xs ?? (All, Take (len - horizon)) where len = cols xs --- Get teacher data-takeFuture :: Int -> Matrix Double -> Matrix Double+-- Get target data+takeFuture :: Element a => Int -> Matrix a -> Matrix a takeFuture horizon = (?? (All, Drop horizon)) -splitAt' :: Int -> Matrix Double -> (Matrix Double, Matrix Double)+-- Manipulate the data matrix: horizontal split+splitAtRatio :: Element a => Double -> Matrix a -> (Matrix a, Matrix a)+splitAtRatio ratio m = splitAt' spl m+ where+ total = cols m+ spl = round $ ratio * fromIntegral total++splitAt' :: Element a => Int -> Matrix a -> (Matrix a, Matrix a) splitAt' i m = (m ?? (All, Take i), m ?? (All, Drop i)) main :: IO ()@@ -22,11 +29,8 @@ -- Load and transpose time series to predict dta <- tr <$> loadMatrix "examples/data/mg.txt" - let splitRatio = 0.50 -- Train on 50% of data- total = cols dta- spl = round $ splitRatio * fromIntegral total- -- Split the data into training and validation data sets- (train, validate) = splitAt' spl dta+ -- Train on 50% of data+ let (train, validate) = splitAtRatio 0.50 dta -- Configure a new NTC network let p = RC.par0 { RC._inputWeightsRange = (0.1, 0.3) }
rc.cabal view
@@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 080028c4b027364671d6876826a008e2c27817e0f1b4394a51c0a890ac59f86e+-- hash: 17efef2c866ea47d85997d6bc72d2cfd2fe68053ff3dfecd15aed57a1004fe33 name: rc-version: 0.1.0.1+version: 0.3.0.0 synopsis: Reservoir Computing, fast RNNs description: Please see the README on Github at <https://github.com/masterdezign/rc#readme> category: Machine Learning@@ -31,15 +31,18 @@ hs-source-dirs: rc build-depends:- Learning+ Learning >=0.1.0 , base >=4.7 && <5- , dde ==0.0.0+ , dde >=0.3.0 , hmatrix >=0.18.0.0+ , linear , random , vector exposed-modules: RC.Helpers RC.NTC+ RC.NTC.Reservoir+ RC.NTC.Types other-modules: Paths_rc default-language: Haskell2010@@ -49,10 +52,11 @@ hs-source-dirs: examples/NTC build-depends:- Learning+ Learning >=0.1.0 , base >=4.7 && <5- , dde ==0.0.0+ , dde >=0.3.0 , hmatrix >=0.18.0.0+ , linear , random , rc , vector@@ -67,10 +71,11 @@ test ghc-options: -threaded -rtsopts -with-rtsopts=-N build-depends:- Learning+ Learning >=0.1.0 , base >=4.7 && <5- , dde ==0.0.0+ , dde >=0.3.0 , hmatrix >=0.18.0.0+ , linear , random , rc , vector
rc/RC/Helpers.hs view
@@ -1,5 +1,7 @@ module RC.Helpers ( addBiases+ , flatten'+ , unflatten' , randList , randMatrix , randSparse@@ -8,7 +10,7 @@ import System.Random import Data.List ( unfoldr )-import qualified Numeric.LinearAlgebra as LA+import Numeric.LinearAlgebra as LA -- | Hard sigmoid hsigmoid :: (Fractional a, Ord a)@@ -21,6 +23,8 @@ f | x < offset = 0.0 | x < width = β * (x - offset) | otherwise = β * (width - offset)+{-# SPECIALISE hsigmoid :: (Double, Double, Double) -> Double -> Double #-}+{-# SPECIALISE hsigmoid :: (Float, Float, Float) -> Float -> Float #-} -- | Prepend a row of ones --@@ -64,3 +68,11 @@ randList n = take n. unfoldr (Just. random) {-# SPECIALISE randList :: Int -> StdGen -> [Float] #-} {-# SPECIALISE randList :: Int -> StdGen -> [Double] #-}++-- | Flatten by columns+flatten' :: Matrix Double -> Vector Double+flatten' = flatten. tr++-- | Reshape using columns instead of rows+unflatten' :: Int -> Vector Double -> Matrix Double+unflatten' nodes = tr. reshape nodes
rc/RC/NTC.hs view
@@ -10,66 +10,28 @@ module RC.NTC ( new , learn+ , learnClassifier , predict , par0 , NTCParameters (..)- , DDEModel.Par (..)- , DDEModel.BandpassFiltering (..) ) where import Numeric.LinearAlgebra import System.Random ( StdGen , split )+import qualified Data.List as List import qualified Data.Vector.Storable as V+import qualified Data.Vector.Storable.Mutable as VM+import System.IO.Unsafe ( unsafePerformIO ) import qualified Learning-import qualified Numeric.DDE as DDE import qualified Numeric.DDE.Model as DDEModel +import RC.NTC.Types+import qualified RC.NTC.Reservoir as Reservoir import qualified RC.Helpers as H --- | Reservoir abstraction.------ Reservoir (a recurrent neural network) is an essential--- component in the NTC framework.------ In this project, we exploit an established analogy between--- spatially extended and delay systems ^1-3. That makes possible--- to employ DDEs as a reservoir substrate. The choice of substrate--- does not affect the `Reservoir` definition below.------ ^1 Arecchi, F. T., et al. “Two-Dimensional Representation of--- a Delayed Dynamical System.” Physical Review A, vol. 45, no. 7,--- Jan. 1992, doi:10.1103/physreva.45.r4225.--- ^2 Appeltant, L., et al. “Information Processing Using a Single--- Dynamical Node as Complex System.” Nature Communications, vol. 2,--- 2011, p. 468., doi:10.1038/ncomms1476.--- ^3 Virtual Chimera States for Delayed-Feedback Systems - Laurent Larger,--- Bogdan Penkovsky, Yuri Maistrenko. Physical Review Letters - 08 / 2013----newtype Reservoir = Reservoir { _transform :: Matrix Double -> Matrix Double }---- | Customizable NTC parameters-data NTCParameters = Par- { _preprocess :: Matrix Double -> Matrix Double- -- ^ Modify data before masking (e.g. compression)- , _inputWeightsRange :: (Double, Double) -- ^ Input weights (mask) range- , _inputWeightsGenerator :: StdGen -> (Int, Int) -> (Double, Double) -> Matrix Double- , _postprocess :: Matrix Double -> Matrix Double- -- ^ Modify data before training or prediction (e.g. add biases)- , _reservoirModel :: DDEModel.Par- }---- | NTC network structure-data NTC = NTC- { _inputWeights :: Matrix Double- , _reservoir :: Reservoir- , _outputWeights :: Maybe (Matrix Double)- -- ^ Trainable part of NTC- , _par :: NTCParameters- }---- | Creates an untrained NTC network+-- | An untrained NTC network new :: StdGen -> NTCParameters@@ -80,7 +42,7 @@ let iwgen = _inputWeightsGenerator par iw = iwgen g (nodes, ind) (_inputWeightsRange par) ntc = NTC { _inputWeights = iw- , _reservoir = genReservoir (_reservoirModel par)+ , _reservoir = _reservoirModel par , _outputWeights = Nothing , _par = par }@@ -93,10 +55,7 @@ , _inputWeightsGenerator = H.randMatrix , _postprocess = H.addBiases -- Usually `id` will work , _inputWeightsRange = undefined -- To be manually set, e.g. (-1, 1)- , _reservoirModel = DDEModel.RC { DDEModel._filt = filt'- , DDEModel._rho = 3.25- , DDEModel._fnl = H.hsigmoid (1.09375, 1.5, 0.0)- }+ , _reservoirModel = Reservoir.genReservoir p } where filt' = DDEModel.BandpassFiltering {@@ -104,35 +63,10 @@ , DDEModel._theta = recip 0.34375 } --- | Substrate-specific low-level reservoir implementation-genReservoir :: DDEModel.Par -> Reservoir-genReservoir par@DDEModel.RC {- DDEModel._filt = DDEModel.BandpassFiltering { DDEModel._tau = tau }- } = Reservoir _r- where- _r sample = unflatten response- where- flatten' = flatten. tr- unflatten = tr. reshape nodes-- oversampling = 1 :: Int -- No oversampling- detuning = 1.0 :: Double -- Delay detuning factor, 1 = no detuning- nodes = rows sample- delaySamples = round $ detuning * fromIntegral (oversampling * nodes)-- -- Matrix to timetrace- trace1 = flatten' sample-- -- Duplicate the last element (DDE.integHeun2_2D consumes one extra input)- trace = trace1 V.++ V.singleton (V.last trace1)-- -- Empirically chosen integration time step:- -- twice faster than the system response time tau- hStep = tau / 2-- !(_, !response) = DDE.integHeun2_2D delaySamples hStep (DDEModel.rhs par) (DDE.Input trace)--genReservoir _ = error "Unsupported DDE model"+ p = DDEModel.RC { DDEModel._filt = filt'+ , DDEModel._rho = 3.25+ , DDEModel._fnl = H.hsigmoid (1.09375, 1.5, 0.0)+ } -- | Nonlinear transformation performed by an NTC network forwardPass :: NTC -- ^ NTC network@@ -141,12 +75,10 @@ forwardPass NTC { _par = Par { _preprocess = prep, _postprocess = post } , _inputWeights = iw , _reservoir = Reservoir res- } sample =+ } !sample = let pipeline = post. res. (iw <>). prep in pipeline sample --- TODO: introduce an explicit `learnClassifier` function- -- | NTC training: learn the readout weights offline learn :: NTC@@ -159,11 +91,100 @@ -> Either String NTC learn ntc forgetPts inp out = ntc' where- state' = (forwardPass ntc inp) ?? (All, Drop forgetPts)+ state' = forwardPass ntc inp ?? (All, Drop forgetPts) teacher' = out ?? (All, Drop forgetPts) ntc' = case Learning.learn' state' teacher' of Nothing -> Left "Cannot create a readout matrix" w -> Right $ ntc { _outputWeights = w }++_concatForwardPass :: NTC+ -> Int+ -- ^ Number of samples in the list (for+ -- memory allocation)+ -> [Matrix Double]+ -- ^ List of samples+ -> Matrix Double+_concatForwardPass ntc m (state0_:samples') = unsafePerformIO $ do+ -- Detect the postprocessing output dimension `nodes`.+ -- Alternatively, use Static shapes from hmatrix.+ let state0 = forwardPass ntc state0_+ nodes = rows state0+ state0' = H.flatten' state0++ -- Allocate memory for a mutable vector+ v <- VM.new (nodes * m)++ -- Copy the first computed state0'+ copy state0' v 0 0 nodes++ let foldA _ [] = return ()+ foldA f ((y, i):ys) = do+ let v0 = f y+ copy v0 v 0 (i * nodes) ((i + 1) * nodes)+ foldA f ys++ -- NB: Does compiler know (H.flatten'. H.unflatten') == id?+ -- If not, consider refactoring genReservoir and forwardPass+ foldA (H.flatten'. forwardPass ntc) $ zip samples' [1..m - 1]++ processed' <- V.unsafeFreeze v+ let processed = H.unflatten' nodes processed'+ return processed+ where+ copy !v0 !v !i0 !i !k+ | i < k = do+ VM.unsafeWrite v i (v0 V.! i0)+ copy v0 v (i0 + 1) (i + 1) k+ | otherwise = return ()++-- | NTC training specific to classification task.+-- The readout weights are learned offline.+--+-- Alternatively, one could use `learn` function. However,+-- to make sure the training samples do not mix in the reservoir RNN,+-- a significant padding of zeros would be needed. Although that+-- way directly corresponds to a physical experiment, that would be not+-- memory-efficient on a computer.+learnClassifier+ :: NTC+ -- ^ NTC network+ -> [Int]+ -- ^ Target labels+ -> Int+ -- ^ Number of inputs (for memory allocation)+ -> [Matrix Double]+ -- ^ Training inputs+ -> Either String (Learning.Classifier Int)+learnClassifier ntc labels m samples =+ case Learning.learn' state teacher' of+ Nothing -> Left "Cannot create a readout matrix"+ Just w -> let clf = Learning.winnerTakesAll w klasses. forwardPass ntc+ in Right (Learning.Classifier clf)+ where+ klasses = fromList. List.sort. List.nub $ labels+ klassesNo = V.length klasses+ teacher' = concatHor. map (flip (Learning.teacher klassesNo) 1) $ labels++ -- Alternatively, a streaming interface might be a solution+ -- https://hackage.haskell.org/package/pipes-4.3.7/docs/Pipes-Tutorial.html+ state = _concatForwardPass ntc m samples++-- | Horizontal matrix concatenation, alternative to fromBlocks+concatHor :: Element a => [Matrix a] -> Matrix a+concatHor ms@(m:_) = foldr (|||) (zeroCols (rows m)) ms+-- concatHor = concatMapHor id+{-# SPECIALIZE concatHor :: [Matrix Double] -> Matrix Double #-}++concatMapHor+ :: Element b =>+ (Matrix a -> Matrix b) -> [Matrix a] -> Matrix b+concatMapHor f ms@(m:_) = foldr (\a b -> f a ||| b) (zeroCols (rows m)) ms+{-# SPECIALIZE concatMapHor :: (Matrix Double -> Matrix Double) -> [Matrix Double] -> Matrix Double #-}++-- | Matrix with zero elements+zeroCols :: V.Storable a => Int -> Matrix a+zeroCols rows' = (rows'><0) []+{-# SPECIALIZE zeroCols :: Int -> Matrix Double #-} -- | Run prediction using a "clean" (uninitialized) reservoir and then -- forget the reservoir's state.
+ rc/RC/NTC/Reservoir.hs view
@@ -0,0 +1,58 @@+-- |+-- = Single-node reservoir+--+-- In this project, we exploit an established analogy between+-- spatially extended and delay systems (refs. 1-3). That makes possible+-- to employ DDEs as a reservoir substrate.+--+-- 1. Arecchi, F. T., et al. “Two-Dimensional Representation of+-- a Delayed Dynamical System.” Physical Review A, vol. 45, no. 7,+-- Jan. 1992, doi:10.1103/physreva.45.r4225.+-- 2. Appeltant, L., et al. “Information Processing Using a Single+-- Dynamical Node as Complex System.” Nature Communications, vol. 2,+-- 2011, p. 468., doi:10.1038/ncomms1476.+-- 3. Virtual Chimera States for Delayed-Feedback Systems - Laurent Larger,+-- Bogdan Penkovsky, Yuri Maistrenko. Physical Review Letters - 08 / 2013.++{-# LANGUAGE BangPatterns #-}+module RC.NTC.Reservoir+ ( Reservoir (..)+ , genReservoir+ ) where++import Numeric.LinearAlgebra+import qualified Data.Vector.Storable as V+import qualified Linear.V2 as V2+import qualified Numeric.DDE as DDE+import qualified Numeric.DDE.Model as DDEModel++import RC.NTC.Types+import qualified RC.Helpers as H+++-- | Substrate-specific low-level reservoir implementation+genReservoir :: DDEModel.RC -> Reservoir+genReservoir par@DDEModel.RC {+ DDEModel._filt = DDEModel.BandpassFiltering { DDEModel._tau = tau }+ } = Reservoir _r+ where+ _r sample = H.unflatten' nodes responseX+ where+ oversampling = 1 :: Int -- No oversampling+ detuning = 1.0 :: Double -- Delay detuning factor, 1 = no detuning+ nodes = rows sample+ delaySamples = round $ detuning * fromIntegral (oversampling * nodes)++ -- Matrix to timetrace+ trace1 = H.flatten' sample++ -- Duplicate the last element (DDE.integHeun2_2D consumes one extra input)+ trace = trace1 V.++ V.singleton (V.last trace1)++ -- Empirically chosen integration time step:+ -- twice faster than the system response time tau+ hStep = tau / 2++ (_, response) = DDE.integHeun2_2D [delaySamples] hStep (DDEModel.bandpassRhs par) (DDE.Input trace)+ responseX = V.map (\(V2.V2 x _) -> x) response+
+ rc/RC/NTC/Types.hs view
@@ -0,0 +1,34 @@+module RC.NTC.Types+ ( Reservoir (..)+ , NTCParameters (..)+ , NTC (..)+ ) where++import System.Random ( StdGen )+import Numeric.LinearAlgebra ( Matrix )++-- | Reservoir abstraction.+--+-- Reservoir (a recurrent neural network) is an essential+-- component in the NTC framework.+newtype Reservoir = Reservoir { _transform :: Matrix Double -> Matrix Double }++-- | Customizable NTC parameters+data NTCParameters = Par+ { _preprocess :: Matrix Double -> Matrix Double+ -- ^ Modify data before masking (e.g. compression)+ , _inputWeightsRange :: (Double, Double) -- ^ Input weights (mask) range+ , _inputWeightsGenerator :: StdGen -> (Int, Int) -> (Double, Double) -> Matrix Double+ , _postprocess :: Matrix Double -> Matrix Double+ -- ^ Modify data before training or prediction (e.g. add biases)+ , _reservoirModel :: Reservoir+ }++-- | NTC network structure+data NTC = NTC+ { _inputWeights :: Matrix Double+ , _reservoir :: Reservoir+ , _outputWeights :: Maybe (Matrix Double)+ -- ^ Trainable part of NTC+ , _par :: NTCParameters+ }