diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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. diff --git a/examples/NTC/Main.hs b/examples/NTC/Main.hs --- a/examples/NTC/Main.hs +++ b/examples/NTC/Main.hs @@ -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) } diff --git a/rc.cabal b/rc.cabal --- a/rc.cabal +++ b/rc.cabal @@ -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 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 diff --git a/rc/RC/Helpers.hs b/rc/RC/Helpers.hs --- a/rc/RC/Helpers.hs +++ b/rc/RC/Helpers.hs @@ -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 diff --git a/rc/RC/NTC.hs b/rc/RC/NTC.hs --- a/rc/RC/NTC.hs +++ b/rc/RC/NTC.hs @@ -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. diff --git a/rc/RC/NTC/Reservoir.hs b/rc/RC/NTC/Reservoir.hs new file mode 100644 --- /dev/null +++ b/rc/RC/NTC/Reservoir.hs @@ -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 + diff --git a/rc/RC/NTC/Types.hs b/rc/RC/NTC/Types.hs new file mode 100644 --- /dev/null +++ b/rc/RC/NTC/Types.hs @@ -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 + }