raw-feldspar 0.1 → 0.2

raw patch · 16 files changed

+1131/−108 lines, 16 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• README.md +14/−10
• examples/FFT.hs +119/−60
• examples/FFTW.hs +139/−0
• examples/FFT_bench.hs +71/−0
• examples/Tut2_ExpressionsAndTypes.hs +1/−1
• examples/Tut5_Matrices.hs +1/−1
• examples/Tut8_SequentialVectors.hs +56/−0
• raw-feldspar.cabal +7/−4
• src/Feldspar/Data/Buffered.hs +21/−7
• src/Feldspar/Data/Queue.hs +142/−0
• src/Feldspar/Data/Vector.hs +210/−9
• src/Feldspar/Optimize.hs +3/−1
• src/Feldspar/Processing/Filters.hs +109/−0
• src/Feldspar/Run/Marshal.hs +40/−0
• tests/Examples.hs +19/−5
• tests/Semantics.hs +179/−10

@@ -12,7 +12,7 @@  The installation can be sped up a bit (and the size of the installation reduced) by adding a flag to `language-c-quote` (a dependency of RAW-Feldspar): -    cabal install --constraint="language-c-quote -full-haskell-antiquotes"+    cabal install --constraint="language-c-quote -full-haskell-antiquotes" raw-feldspar  However, this flag should normally only be used when installing in a sandbox that has no other packages depending `language-c-quote`. @@ -20,6 +20,10 @@  The best way to learn how to use RAW-Feldspar at the moment is to look through the [examples](examples/). We suggest going through the files named "TutN_..." in ascending order. The files are well-documented. +The easiest way for users of `cabal` to get access to the examples is to run++    cabal unpack raw-feldspar+ There is also some guidance in the [Haddock documentation](http://hackage.haskell.org/package/raw-feldspar).  The vector library is central to programming in Feldspar. Its general operation is explained in the [Haddock documentation](http://hackage.haskell.org/package/raw-feldspar/docs/Feldspar-Data-Vector.html), and many [examples](examples/) are using vectors.@@ -112,17 +116,17 @@  ### Zeldspar -[Zeldspar](../../../../koengit/zeldspar) is an implementation of the [Ziria DSL](http://dx.doi.org/10.1145/2694344.2694368) for wireless programming on top of RAW-Feldspar.+[Zeldspar](https://github.com/koengit/zeldspar) is an implementation of the [Ziria DSL](http://dx.doi.org/10.1145/2694344.2694368) for wireless programming on top of RAW-Feldspar.  ### raw-feldspar-mcs -[raw-feldspar-mcs](../../../../kmate/raw-feldspar-mcs) extends RAW-Feldspar and Zeldspar with multicore and scratchpad support.+[raw-feldspar-mcs](https://github.com/kmate/raw-feldspar-mcs) extends RAW-Feldspar and Zeldspar with multicore and scratchpad support. -The repository contains many [examples](../../../../kmate/raw-feldspar-mcs/tree/master/examples) written for the [Parallella](http://www.parallella.org) multicore architecture.+The repository contains many [examples](https://github.com/kmate/raw-feldspar-mcs/tree/master/examples) written for the [Parallella](http://www.parallella.org) multicore architecture.  ### feldspar-synch -[feldspar-synch](../../../../emilaxelsson/feldspar-synch) is a library that extends Feldspar with Yampa-style synchronous streams.+[feldspar-synch](https://github.com/emilaxelsson/feldspar-synch) is a library that extends Feldspar with Yampa-style synchronous streams.  It contains a simple polyphonic synthesizer as a demonstration. The synthesizer may serve as a simple example of a complete (toy) application written in RAW-Feldspar. It also demonstrates how to make bindings to an external C library (the ALSA sound library). @@ -132,7 +136,7 @@    * [feldspar-language](http://hackage.haskell.org/package/feldspar-language) -- the language front end   * [feldspar-compiler](http://hackage.haskell.org/package/feldspar-compiler) -- the C-generating back end-  * [feldspar-io](../../../../emilaxelsson/feldspar-io) -- a monadic "IO" layer+  * [feldspar-io](https://github.com/emilaxelsson/feldspar-io) -- a monadic "IO" layer  (`feldspar-io`, which is still at an early stage of development, adds support for writing interactive programs calling external functions, etc.) @@ -146,13 +150,13 @@  RAW-Feldspar has since become a respectable replacement of the previous implementation. RAW-Feldspar typically generates slicker code based on native types and functions. Due to the new design, the user also has more control over array allocations, leading to lower memory usage and fewer array copies. -However, RAW-Feldspar also has some [limitations and lacks some features](../../wiki/Limitations-and-Missing-Features) compared to the previous version. Some features are missing simply because they have not been ported yet; others are missing for more fundamental reasons.+However, RAW-Feldspar also has some [limitations and lacks some features](https://github.com/Feldspar/raw-feldspar/wiki/Limitations-and-Missing-Features) compared to the previous version. Some features are missing simply because they have not been ported yet; others are missing for more fundamental reasons.  ## Limitations and missing features -See [limitations and missing features](../../wiki/Limitations-and-Missing-Features).+See [limitations and missing features](https://github.com/Feldspar/raw-feldspar/wiki/Limitations-and-Missing-Features). -There is also a [list of possible enhancements and fixes](../../wiki/TODOs).+There is also a [list of possible enhancements and fixes](https://github.com/Feldspar/raw-feldspar/wiki/TODOs).  ## Implementation @@ -175,5 +179,5 @@  A prime example of the technique is the [vector library](http://hackage.haskell.org/package/raw-feldspar/docs/Feldspar-Data-Vector.html), which provides high-level vector representations with a rich programming interface. These vectors only exist in the meta-language (i.e. Haskell), and by the time the Feldspar compiler is called, the vectors are already gone and what is left is imperative code with highly optimized loops. We saw an example of this when compiling the `sumSq` example above. -Another example is [feldspar-synch](../../../../emilaxelsson/feldspar-synch), which extends Feldspar with synchronous streams. The whole package is implemented as a shallow extension on top of RAW-Feldspar.+Another example is [feldspar-synch](https://github.com/emilaxelsson/feldspar-synch), which extends Feldspar with synchronous streams. The whole package is implemented as a shallow extension on top of RAW-Feldspar. 

@@ -1,38 +1,18 @@ {-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE GADTs #-} --- Copyright (c) 2013, Emil Axelsson, Peter Jonsson, Anders Persson and---                     Josef Svenningsson--- Copyright (c) 2012, Emil Axelsson, Gergely Dévai, Anders Persson and---                     Josef Svenningsson--- Copyright (c) 2009-2011, ERICSSON AB--- All rights reserved.------ Redistribution and use in source and binary forms, with or without--- modification, are permitted provided that the following conditions are met:+-- | FFT implementation inspired by the paper "Feldspar: Application and+-- implementation": -----     * 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 ERICSSON AB nor the names of its contributors---       may be used to endorse or promote products derived from this software---       without specific prior written permission.+-- <http://publications.lib.chalmers.se/records/fulltext/local_156271.pdf> ----- 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.+-- There are a few differences, partly due to the paper using a different+-- Feldspar implementation. But regardless, the best way to understand the+-- definitions in this file is by reading the paper.  module FFT-  ( fft+  ( tw  -- Exported to allow pre-computation+  , fftCore+  , fft   , ifft   ) where @@ -40,12 +20,18 @@  import Prelude () +import Data.Bool (bool)+ import Feldspar.Run import Feldspar.Data.Vector import Feldspar.Data.Buffered   +----------------------------------------+-- * Helper functions+----------------------------------------+ rotBit :: Data Index -> Data Index -> Data Index rotBit k i = lefts .|. rights   where@@ -57,60 +43,115 @@ riffle :: (Pully pull a, Syntax a) => Data Index -> pull -> Pull a riffle = backPermute . const . rotBit +testBit :: (Bits a, Integral a, PrimType a) => Data a -> Data Index -> Data Bool+testBit a i = i2b (a .&. (1 .<<. i2n i))++-- | @2^n@+twoTo :: (Num a, Bits a, PrimType a) => Data Index -> Data a+twoTo n = 1 .<<. i2n n++flipBit :: (Num a, Bits a, PrimType a) => Data a -> Data Index -> Data a+flipBit i k = i `xor` twoTo k++++----------------------------------------+-- * Bit-reversal permutation+----------------------------------------+ bitRev :: (Manifestable Run vec a, Finite vec, Syntax a)     => Store a+    -> Length  -- ^ Unrolling steps in inner loops (1 means no unrolling)     -> Data Length     -> vec     -> Run (Manifest a)-bitRev st n = loopStore st (1,1,Incl n) $ \i -> return . riffle i+bitRev st u n = loopStore st (1,1,Excl n) $ \i -> return . unroll u . riffle i -testBit :: (Bits a, Num a, PrimType a) => Data a -> Data Index -> Data Bool-testBit a i = a .&. (1 .<<. i2n i) /= 0 -fftCore-    :: ( Manifestable Run vec (Data (Complex a))-       , Finite vec++----------------------------------------+-- * FFT+----------------------------------------++tw :: (Floating a, PrimType a, PrimType (Complex a))+    => Bool  -- ^ Inverse FFT?+    -> Data Index+    -> Data Index+    -> Data (Complex a)+tw inv n k = polar 1 (bool (-2) 2 inv * π * i2n k / i2n n)++twids+    :: ( Pully ts (Data (Complex a))        , RealFloat a        , PrimType a        , PrimType (Complex a)+       , Pully vec (Data (Complex a))        )-    => Store (Data (Complex a))-    -> Bool  -- ^ Inverse?+    => ts+    -> Data Index+    -> Data Index     -> Data Length     -> vec-    -> Run (DManifest (Complex a))-fftCore st inv n = loopStore st (n+1,-1,Incl 1) $ \i -> return . step (i-1)-    -- Note: Cannot loop from n to 0 because 0-1 is `maxBound`, so the loop will-    -- go on forever.+    -> DPull (Complex a)+twids ts n k l vec = Pull l $ \i ->+    let j = (lsbs (i2n k) i) .<<. (n'-1-k')+    in  (testBit i k) ? ((ts!j) * (vec!i)) $ (vec!i)   where-    step k vec = Pull (length vec) ixf-      where-        ixf i = testBit i k ? (twid * (b - a)) $ (a+b)-          where-            k'   = i2n k-            a    = vec ! i-            b    = vec ! (i `xor` k2)-            twid = polar 1 ((if inv then π else -π) * i2n (lsbs k' i) / i2n k2)-            k2   = 1 .<<. k'+    n' = i2n n+    k' = i2n k -fft'-    :: ( Manifestable Run vec (Data (Complex a))+bfly+    :: ( Pully vec (Data (Complex a))+       , RealFloat a+       , PrimType a+       , PrimType (Complex a)+       )+    => Data Index -> vec -> DPull (Complex a)+bfly k as = Pull (length as) $ \i ->+    let a = as ! i+        b = as ! flipBit i k+    in  (testBit i k) ? (b-a) $ (a+b)++-- | Core of the FFT+--+-- It is normally better to use 'fft' or 'ifft' than this functon; however, for+-- doing repeated FFT on vectors of the same size, 'fftCore' can be used to+-- avoid recomputing the twiddle factors and the number of stages.+fftCore+    :: ( Pully ts (Data (Complex a))+       , Manifestable Run vec (Data (Complex a))        , Finite vec        , RealFloat a        , PrimType a        , PrimType (Complex a)        )     => Store (Data (Complex a))-    -> Bool  -- ^ Inverse?+    -> Length       -- ^ Unrolling steps in inner loops (1 means no unrolling)+    -> ts           -- ^ Twiddle factors+    -> Data Length  -- ^ Number of stages     -> vec     -> Run (DManifest (Complex a))-fft' st inv v = do-    n <- shareM (ilog2 (length v) - 1)-    fftCore st inv n v >>= bitRev st n+fftCore st u ts n vec = do+    let step i = return . unroll u . twids ts n i (length vec) . bfly i+    loopStore st ((i2n n :: Data Int32)-1,-1,Incl 0) (step . i2n) vec+      >>= bitRev st u n+      -- `i2n` is used to make the loop index a signed number. Otherwise the+      -- index will wrap to maxBound before the loop test after the final+      -- iteration.+      --+      -- An alternative is to use:+      --+      --     loopStore st (n,-1,Excl 0) (step . subtract 1) vec  -- | Radix-2 Decimation-In-Frequency Fast Fourier Transformation of the given -- complex vector. The given vector must be power-of-two sized, (for example 2, -- 4, 8, 16, 32, etc.) The output is non-normalized.+--+-- The length of the vector must be divisible by the number of unrolling steps.+--+-- The optimal amount of unrolling depends on the target architecture, but a+-- value of 2 might be a reasonable default that gives some performance+-- improvements on many systems and doesn't lead to too much code size increase. fft     :: ( Manifestable Run vec (Data (Complex a))        , Finite vec@@ -119,14 +160,26 @@        , PrimType (Complex a)        )     => Store (Data (Complex a))+    -> Length  -- ^ Unrolling steps in inner loops (1 means no unrolling)     -> vec     -> Run (DManifest (Complex a))-fft st = fft' st False+fft st u vec = do+    n  <- shareM (ilog2 (length vec))+    ts <- manifestFresh $ Pull (twoTo (n-1)) (tw False (twoTo n))+      -- Change `manifestFresh` to `return` to avoid pre-computing twiddle+      -- factors+    fftCore st u ts n vec  -- | Radix-2 Decimation-In-Frequency Inverse Fast Fourier Transformation of the -- given complex vector. The given vector must be power-of-two sized, (for -- example 2, 4, 8, 16, 32, etc.) The output is divided with the input size,--- thus giving 'ifft . fft == id'.+-- thus giving @`ifft` . `fft` == id@.+--+-- The length of the vector must be divisible by the number of unrolling steps.+--+-- The optimal amount of unrolling depends on the target architecture, but a+-- value of 2 might be a reasonable default that gives some performance+-- improvements on many systems and doesn't lead to too much code size increase. ifft     :: ( Manifestable Run vec (Data (Complex a))        , Finite vec@@ -135,9 +188,15 @@        , PrimType (Complex a)        )     => Store (Data (Complex a))+    -> Length  -- ^ Unrolling steps in inner loops (1 means no unrolling)     -> vec     -> Run (DPull (Complex a))-ifft st v = normalize <$> fft' st True v+ifft st u vec = do+    n  <- shareM (ilog2 (length vec))+    ts <- manifestFresh $ Pull (twoTo (n-1)) (tw True (twoTo n))+      -- Change `manifestFresh` to `return` to avoid pre-computing twiddle+      -- factors+    normalize <$> fftCore st u ts n vec   where-    normalize = map (/ (i2n $ length v))+    normalize = map (/ (i2n $ length vec)) 

@@ -0,0 +1,139 @@+-- | Interface to the FFTW routine+--+-- Ubuntu installation of C library:+--+--     sudo apt-get install libfftw3-dev+--+-- Basic usage:+-- <http://www.fftw.org/fftw3_doc/Complex-One_002dDimensional-DFTs.html#Complex-One_002dDimensional-DFTs>+--+-- Linker flags:+--+--     cc file.c -lm -lfftw3++module FFTW where++++import qualified Prelude as P++import qualified Data.Complex as Complex++import qualified Test.QuickCheck as QC+import qualified Test.QuickCheck.Monadic as QC++import Language.Embedded.Expression (varExp)+  -- `varExp` can be used to make named constants in expressions++import Feldspar.Run+import Feldspar.Data.Vector++import DFT+import FFT_bench (printTime_def)++++-- | Wrapper for the FFTW routine+--+-- This wrapper is mainly used for testing. It's not suitable for real code+-- because it constructs a plan every time it's called, and it silently+-- allocates an array for the output.+fftw :: DManifest (Complex Double) -> Run (DManifest (Complex Double))+fftw inp = do+    addInclude "<fftw3.h>"+    out  <- newArr (length inp)+    plan <- newObject "fftw_plan" False+    callProcAssign plan "fftw_plan_dft_1d"+      [ valArg (length inp)+      , iarrArg inp+      , arrArg out+      , valArg (varExp "FFTW_FORWARD"  :: Data Word32)+      , valArg (varExp "FFTW_ESTIMATE" :: Data Word32)+      ]+    callProc "fftw_execute" [objArg plan]+    freezeArr out++++almostEq a b+    =    Complex.magnitude d P.< 1e-7+    P.&& Complex.phase d     P.< 1e-7+  where+    d = abs (a-b)++a ~= b = P.and $ P.zipWith almostEq a b++prop_fft_dft dft' fft' = QC.monadicIO $ do+    n   :: Int              <- QC.pick $ QC.choose (2,5)+    inp :: [Complex Double] <- QC.pick $ QC.vector (2 P.^ n)+    outd <- QC.run $ dft' inp+    outf <- QC.run $ fft' inp+    QC.assert (outd ~= outf)++-- | Compare 'fftw' against 'dft'+testFFTW =+    marshalledM (return . dft) $ \dft'  ->+    marshalledM fftw           $ \fftw' ->+      QC.quickCheck $ prop_fft_dft dft' fftw'+  where+    marshalledM = marshalled' def def+      { externalFlagsPre  = ["-Wno-incompatible-pointer-types"]+      , externalFlagsPost = ["-lm -lfftw3"]+      }++++sizeOf_fftw_complex :: Data Length+sizeOf_fftw_complex = 16+  -- Checked on an x86_64 system++-- | Measure the time for 100 runs of 'fftw' (excluding initialization) for+-- arrays of the given size+benchmark n = do+    addInclude "<stdio.h>"+    addInclude "<string.h>"+    addInclude "<time.h>"+    addInclude "<fftw3.h>"++    addDefinition printTime_def++    inp  <- newObject "fftw_complex" True+    out  <- newObject "fftw_complex" True+    plan <- newObject "fftw_plan" False++    callProcAssign inp "fftw_malloc" [valArg (n*sizeOf_fftw_complex)]+    callProcAssign out "fftw_malloc" [valArg (n*sizeOf_fftw_complex)]+    callProc "memset"+      [ objArg inp+      , valArg (0 :: Data Index)+      , valArg (n*sizeOf_fftw_complex)+      ]++    callProcAssign plan "fftw_plan_dft_1d"+      [ valArg (n :: Data Word32)+      , objArg inp+      , objArg out+      , valArg (varExp "FFTW_FORWARD"  :: Data Word32)+      , valArg (varExp "FFTW_ESTIMATE" :: Data Word32)+          -- Change to `FFTW_MEASURE` to enable tuning+      ]++    start <- newObject "clock_t" False+    end   <- newObject "clock_t" False+    callProcAssign start "clock" []++    for (0,1,Excl 100) $ \(_ :: Data Index) ->+      callProc "fftw_execute" [objArg plan]++    callProcAssign end "clock" []+    callProc "printTime" [objArg start, objArg end]++    callProc "fftw_destroy_plan" [objArg plan]+    callProc "fftw_free" [objArg inp]+    callProc "fftw_free" [objArg out]++runBenchmark n = runCompiled'+    def+    def {externalFlagsPre = ["-O3"], externalFlagsPost = ["-lm","-lfftw3"]}+    (benchmark n)+

@@ -0,0 +1,71 @@+{-# LANGUAGE QuasiQuotes #-}++module FFT_bench where++++import Prelude ()++import Feldspar.Run+import Feldspar.Data.Buffered+import Feldspar.Data.Vector++import FFT++++sizeOf_double_complex :: Data Length+sizeOf_double_complex = 16+  -- Checked on an x86_64 system+  -- TODO Feldspar should have a built-in `sizeof` function++-- | @2^n@+twoTo :: (Num a, Bits a, PrimType a) => Data Index -> Data a+twoTo n = 1 .<<. i2n n++printTime_def = [cedecl|+void printTime(typename clock_t start, typename clock_t end)+{+  printf("CPU time (sec): %f\n", (double)(end-start) / CLOCKS_PER_SEC);+}+|]++-- | Measure the time for 100 runs of 'fftCore' (excluding initialization) for+-- arrays of the given size+benchmark n = do+    addInclude "<stdio.h>"+    addInclude "<string.h>"+    addInclude "<time.h>"++    addDefinition printTime_def++    start <- newObject "clock_t" False+    end   <- newObject "clock_t" False++    st :: Store (Data (Complex Double)) <- newStore n+    inp <- unsafeFreezeStore n st+    callProc "memset"+      [ iarrArg inp+      , valArg (0 :: Data Index)+      , valArg (n*sizeOf_double_complex)+      ]++    n  <- shareM (ilog2 (length inp))+    ts <- manifestFresh $ Pull (twoTo (n-1)) (tw True (twoTo n))+      -- Change `manifestFresh` to `return` to avoid pre-computing twiddle+      -- factors++    callProcAssign start "clock" []++    for (0,1,Excl 100) $ \(_ :: Data Index) ->+      void $ fftCore st 2 ts n inp++    callProcAssign end "clock" []+    callProc "printTime" [objArg start, objArg end]++runBenchmark n = runCompiled'+    def {compilerAssertions = select []}+      -- Note: important to turn off assertions when running the benchmarks+    def {externalFlagsPre = ["-O3"], externalFlagsPost = ["-lm"]}+    (benchmark n)+

@@ -82,7 +82,7 @@ -- A function computing the nth Fibonacci number can be written using `forLoop`:  fib :: Data Word32 -> Data Word32-fib n = fst $ forLoop n (0,1) $ \i (a,b) -> (b,a+b)+fib n = fst $ forLoop n (0,1) $ \_ (a,b) -> (b,a+b)  -- The arguments to `forLoop` are as follows, in order: --

@@ -101,7 +101,7 @@       $ transposePush       $ hideRows (numRows vec)       $ sequens-      $ map (fft st)+      $ map (fft st 1)       $ exposeRows       $ transpose       $ vec

@@ -0,0 +1,56 @@+{-# LANGUAGE PartialTypeSignatures #-}++{-# OPTIONS_GHC -fno-warn-partial-type-signatures #-}++module Tut8_SequentialVectors where++-- This file demonstrates how to use sequential vectors.++-- If you're looking for infinite streams, please see the library:+-- <https://github.com/emilaxelsson/feldspar-synch>++import qualified Prelude++import Feldspar.Run+import Feldspar.Data.Vector++++-- The vector representations `Pull`, `Push` etc. are designed with+-- data-parallel programming in mind. However, it is quite common that elements+-- in a vector depend on each other in a sequential fashion, for example, when+-- filtering signals in the time domain. That is what sequential vectors are+-- for. You can think of sequential vectors as those that can be produced from+-- `scan`-like operations.++++--------------------------------------------------------------------------------++-- The first `n` Fibonacci numbers using `unfold`:++fibSeq :: MonadComp m => Data Length -> DSeq m Index+fibSeq n = unfold n (\(a,b) -> ((b,a+b),a)) (0,1)++fibSeqRun = connectStdIO $ return . (fibSeq :: _ -> _ Run _)++comp_fibSeq = icompile fibSeqRun++run_fibSeq = runCompiled fibSeqRun++++--------------------------------------------------------------------------------++-- Common filters are defined for sequential vectors in the module+-- `Feldspar.Processing.Filters`.++++--------------------------------------------------------------------------------++testAll = do+    comp_fibSeq+    compareCompiled fibSeqRun (putStr "10 0 1 1 2 3 5 8 13 21 34 ") "10\n"+    compareCompiled fibSeqRun (runIO fibSeqRun) "10\n"+

@@ -1,5 +1,5 @@ name:                raw-feldspar-version:             0.1+version:             0.2 synopsis:            Resource-Aware Feldspar description:         An implementation of the Feldspar EDSL with focus on                      resource-awareness.@@ -19,7 +19,8 @@ author:              Emil Axelsson maintainer:          emax@chalmers.se copyright:           Copyright (c) 2016 Anders Persson, Anton Ekblad, Emil Axelsson,-                                        Koen Claessen, Markus Aronsson, Máté Karácsony+                                        Josef Svenningsson, Koen Claessen, Markus Aronsson,+                                        Máté Karácsony                      Copyright (c) 2015 Emil Axelsson homepage:            https://github.com/Feldspar/raw-feldspar bug-reports:         https://github.com/Feldspar/raw-feldspar/issues@@ -55,11 +56,13 @@     Feldspar.Run.Marshal     Feldspar.Run     Feldspar.Data.Array-    Feldspar.Data.Vector     Feldspar.Data.Buffered     Feldspar.Data.Option-    Feldspar.Data.Validated+    Feldspar.Data.Queue     Feldspar.Data.Storable+    Feldspar.Data.Validated+    Feldspar.Data.Vector+    Feldspar.Processing.Filters    other-modules:     Data.Inhabited.TH

@@ -44,6 +44,7 @@ module Feldspar.Data.Buffered   ( Store   , newStore+  , unsafeInplaceStore   , unsafeFreezeStore   , unsafeFreezeStore2   , setStore@@ -76,12 +77,21 @@     , freeBuf   :: Arr a     } --- | Create a new double-buffered 'Store', initialized to a 0x0 matrix+-- | Create a new double-buffered 'Store' -- -- This operation allocates two arrays of the given length. newStore :: (Syntax a, MonadComp m) => Data Length -> m (Store a) newStore l = Store <$> newNamedArr "store" l <*> newNamedArr "store" l +-- | Create a new single-buffered 'Store'+--+-- Using 'unsafeInplaceStore' instead of 'newStore' allows double-buffered+-- algorithms to run inplace.+unsafeInplaceStore :: (Syntax a, MonadComp m) => Data Length -> m (Store a)+unsafeInplaceStore l = do+    arr <- newNamedArr "store" l+    return $ Store arr arr+ -- | Read the contents of a 'Store' without making a copy. This is generally -- only safe if the the 'Store' is not updated as long as the resulting vector -- is alive.@@ -152,15 +162,17 @@     (r,c) = extent2 vec  loopStore-    :: ( Syntax a+    :: ( Integral i+       , PrimType i+       , Syntax a        , Manifestable Run vec1 a        , Finite vec1        , Manifestable Run vec2 a        , Finite vec2        )     => Store a-    -> IxRange (Data Length)-    -> (Data Index -> Manifest a -> Run vec1)+    -> IxRange (Data i)+    -> (Data i -> Manifest a -> Run vec1)     -> vec2     -> Run (Manifest a) loopStore st rng body init = do@@ -175,15 +187,17 @@     unsafeFreezeStore l st  loopStore2-    :: ( Syntax a+    :: ( Integral i+       , PrimType i+       , Syntax a        , Manifestable2 Run vec1 a        , Finite2 vec1        , Manifestable2 Run vec2 a        , Finite2 vec2        )     => Store a-    -> IxRange (Data Length)-    -> (Data Index -> Manifest2 a -> Run vec1)+    -> IxRange (Data i)+    -> (Data i -> Manifest2 a -> Run vec1)     -> vec2     -> Run (Manifest2 a) loopStore2 st rng body init = do

@@ -0,0 +1,142 @@+-- Copyright (c) 2016, Emil Axelsson, Peter Jonsson, Anders Persson and+--                     Josef Svenningsson+-- 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 ERICSSON AB 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.++-- | Indexable FIFO queues++module Feldspar.Data.Queue+  ( Queue (..)+  , initQueueFromBuffer+  , initQueue+  , newQueue+  , initQueueFromBuffer2+  , initQueue2+  , newQueue2+  ) where++++import Prelude ()++import Feldspar+import Feldspar.Data.Vector++++-- | Indexable FIFO queue+data Queue a = Queue+    { indexQ :: forall m   . MonadComp m => Data Index -> m a+    , putQ   :: forall m   . MonadComp m => a -> m ()+    , withQ  :: forall m b . (Syntax b, MonadComp m) => (Pull a -> m b) -> m b+    }++-- Another option would be to represent a queue as its state (the counter and+-- the array), but the above representation leaves room for other+-- implementations.++-- | Create a new cyclic queue using an existing array as buffer. The length of+-- the array determines the queue size.+initQueueFromBuffer :: forall m a . (Syntax a, MonadComp m) =>+    Arr a -> m (Queue a)+initQueueFromBuffer buf = do+    ir <- initRef 0+    let indexQ :: forall m2 . MonadComp m2 => Data Index -> m2 a+        indexQ j = do+          i <- unsafeFreezeRef ir+          getArr buf $ calcIndex i j+        putQ :: forall m2 . MonadComp m2 => a -> m2 ()+        putQ a = do+          i <- unsafeFreezeRef ir+          setArr buf i a+          setRef ir ((i+1) `mod` len)+        withQ :: forall m2 b . (Syntax b, MonadComp m2) => (Pull a -> m2 b) -> m2 b+        withQ f = do+          i   <- unsafeFreezeRef ir+          vec <- unsafeFreezeArr buf+          f $ backPermute (\_ -> calcIndex i) vec+    return Queue {..}+  where+    len = length buf+    calcIndex i j = (len+i-j-1) `mod` len++-- | Create a new cyclic queue initialized by the given vector (which also+-- determines the size)+initQueue :: (Manifestable m vec a, Finite vec, Syntax a, MonadComp m)+    => vec  -- ^ Initial content (also determines the queue size)+    -> m (Queue a)+initQueue init = do+    buf <- newArr $ length init+    manifestStore buf init+    initQueueFromBuffer buf++-- | Create a new cyclic queue of the given length without initialization+newQueue :: (Syntax a, MonadComp m) => Data Length -> m (Queue a)+newQueue l = newArr l >>= initQueueFromBuffer++initQueueFromBuffer2 :: forall m a . (Syntax a, MonadComp m)+    => Data Length  -- ^ Queue size, must be <= half the buffer size+    -> Arr a        -- ^ Buffer+    -> m (Queue a)+initQueueFromBuffer2 len buf = do+    ir <- initRef 0+    let indexQ :: forall m2 . MonadComp m2 => Data Index -> m2 a+        indexQ j = do+          i <- unsafeFreezeRef ir+          getArr buf (len+i-j-1)+        putQ :: forall m2 . MonadComp m2 => a -> m2 ()+        putQ a = do+          i <- unsafeFreezeRef ir+          setArr buf i a+          setArr buf (i+len) a+          setRef ir ((i+1) `mod` len)+        withQ :: forall m2 b . (Syntax b, MonadComp m2) => (Pull a -> m2 b) -> m2 b+        withQ f = do+          i <- unsafeFreezeRef ir+          vec <- unsafeFreezeArr buf+          f $ reverse $ take len $ drop i vec+    return Queue {..}++-- | Create a new cyclic queue. This implementation uses a buffer twice as long+-- as the queue size to avoid modulus operations when accessing the elements.+initQueue2 :: (Pushy m vec a, Finite vec, Syntax a, MonadComp m)+    => vec  -- ^ Initial content (also determines the queue size)+    -> m (Queue a)+initQueue2 init = do+    buf <- newArr (2*len)+    manifestStore buf (init++init)+    initQueueFromBuffer2 len buf+  where+    len = length init++-- | Create a new cyclic queue. This implementation uses a buffer twice as long+-- as the queue size to avoid modulus operations when accessing the elements.+newQueue2 :: (Syntax a, MonadComp m)+    => Data Length  -- ^ Queue size+    -> m (Queue a)+newQueue2 l = do+    buf <- newArr (2*l)+    initQueueFromBuffer2 l buf+

@@ -276,7 +276,7 @@       MarshalFeld (Pull a)   where     type HaskellRep (Pull a) = HaskellRep (Manifest a)-    fwrite hdl = fwrite hdl . toPush+    fwrite hdl = fwrite hdl . toSeq     fread hdl  = (toPull :: Manifest a -> _) <$> fread hdl  data VecChanSizeSpec lenSpec = VecChanSizeSpec (Data Length) lenSpec@@ -299,13 +299,55 @@         untypedReadChanBuf (Proxy :: Proxy a) c 0 len arr         toPull <$> unsafeFreezeArr arr     untypedWriteChan c v = do-        arr <- newArr len+        -- TODO: can we avoid an array copy here if v is already manifest?+        arr <- unsafeThawArr =<< manifestFresh v         untypedWriteChan c len         untypedWriteChanBuf (Proxy :: Proxy a) c 0 len arr       where         len = length v  -- TODO Make instances for other vector types +instance ( Syntax a, BulkTransferable a+         , ContainerType a ~ Arr a+         ) => Transferable (IArr a)+  where+    type SizeSpec (IArr a) = VecChanSizeSpec (SizeSpec a)+    calcChanSize _ (VecChanSizeSpec n m) =+        let hsz = n `Imp.timesSizeOf` (Proxy :: Proxy Length)+            bsz = calcChanSize (Proxy :: Proxy a) m+        in  hsz `Imp.plusSize` (n `Imp.timesSize` bsz)+    untypedReadChan c = do+        len :: Data Length <- untypedReadChan c+        arr <- newArr len+        untypedReadChanBuf (Proxy :: Proxy a) c 0 len arr+        unsafeFreezeArr arr+    untypedWriteChan c v = do+        arr <- unsafeThawArr v+        untypedWriteChan c len+        untypedWriteChanBuf (Proxy :: Proxy a) c 0 len arr+      where+        len = length v++instance ( Syntax a, BulkTransferable a+         , ContainerType a ~ Arr a+         ) => Transferable (Arr a)+  where+    type SizeSpec (Arr a) = VecChanSizeSpec (SizeSpec a)+    calcChanSize _ (VecChanSizeSpec n m) =+        let hsz = n `Imp.timesSizeOf` (Proxy :: Proxy Length)+            bsz = calcChanSize (Proxy :: Proxy a) m+        in  hsz `Imp.plusSize` (n `Imp.timesSize` bsz)+    untypedReadChan c = do+        len :: Data Length <- untypedReadChan c+        arr <- newArr len+        untypedReadChanBuf (Proxy :: Proxy a) c 0 len arr+        pure arr+    untypedWriteChan c v = do+        untypedWriteChan c len+        untypedWriteChanBuf (Proxy :: Proxy a) c 0 len v+      where+        len = length v+ -- | Data structures that are 'Pull'-like (i.e. support '!' and 'length') class    (Indexed vec, Finite vec, IndexedElem vec ~ a) => Pully vec a instance (Indexed vec, Finite vec, IndexedElem vec ~ a) => Pully vec a@@ -655,6 +697,14 @@       where         len = length vec +instance (MonadComp m1, m1 ~ m2) => Pushy m1 (Seq m2 a) a+  where+    toPush (Seq len init) = Push len $ \write -> do+      next <- init+      for (0,1,Excl len) $ \i -> do+        a <- next i+        write i a+ -- | Dump the contents of a 'Push' vector dumpPush     :: Push m a                   -- ^ Vector to dump@@ -735,8 +785,31 @@     v   = toPush vec     len = length v +-- | Convert a vector to a push vector that computes @n@ elements in each step.+-- This can be used to achieve loop unrolling.+--+-- The length of the vector must be divisible by the number of unrolling steps.+unroll :: (Pully vec a, MonadComp m)+    => Length  -- ^ Number of steps to unroll+    -> vec+    -> Push m a+unroll 0 _   = Prelude.error "unroll: cannot unroll 0 steps"+unroll n vec = Push len $ \write -> do+    assert+      ((len `mod` value n) == 0)+      ("unroll: length not divisible by " Prelude.++ show n)+    for (0,n',Excl len) $ \i ->+      Prelude.sequence_+        [ do k <- shareM (i + value j)+             write k (vec!k)+          | j <- [0..n-1]+        ]+  where+    n'  = Prelude.fromIntegral n+    len = length vec  + -------------------------------------------------------------------------------- -- * 2-dimensional push vectors --------------------------------------------------------------------------------@@ -880,6 +953,133 @@   --------------------------------------------------------------------------------+-- * Sequential vectors+--------------------------------------------------------------------------------++-- | Finite sequential vector+--+-- Users interested in infinite streams are referred to the library:+-- <https://github.com/emilaxelsson/feldspar-synch>+data Seq m a+  where+    Seq :: Data Length -> m (Data Index -> m a) -> Seq m a++-- | 'Seq' vector specialized to 'Data' elements+type DSeq m a = Seq m (Data a)++instance Monad m => Functor (Seq m)+  where+    fmap f (Seq len init) = Seq len $ do+      next <- init+      return $ fmap f . next++instance Finite (Seq m a)+  where+    length (Seq len _) = len++instance+    ( Syntax a+    , MarshalHaskell (Internal a)+    , MarshalFeld a+    , m ~ Run+    ) =>+      MarshalFeld (Seq m a)+  where+    type HaskellRep (Seq m a) = HaskellRep (Manifest a)++    fwrite hdl (Seq len init) = do+      next <- init+      fwrite hdl len >> printf " "+      for (0,1,Excl len) $ \i -> next i >>= fwrite hdl >> printf " "++    fread hdl = toSeq . (id :: Manifest _ -> _) <$> fread hdl+      -- Need to go through a temporary array to avoid embedding side-effects in+      -- the resulting vector. E.g. we don't want to duplicate the reads if the+      -- vector is duplicated.++-- | Vectors that can be converted to 'Seq'+class Seqy m vec a | vec -> a+  where+    -- | Convert a vector to 'Seq'+    toSeq :: vec -> Seq m a++-- | A version of 'toSeq' that constrains the @m@ argument of 'Seq' to that of+-- the monad in which the result is returned. This can be a convenient way to+-- avoid unresolved overloading.+toSeqM :: (Seqy m vec a, Monad m) => vec -> m (Seq m a)+toSeqM = return . toSeq++instance (Syntax a, MonadComp m) => Seqy m (Manifest a) a where toSeq = toSeq . toPull+instance (m1 ~ m2)               => Seqy m1 (Seq m2 a) a  where toSeq = id++instance MonadComp m => Seqy m (Pull a) a+  where+    toSeq vec = Seq (length vec) $ return $ \i -> return $ vec!i++zipWithSeq :: (Seqy m vec1 a, Seqy m vec2 b, Monad m) =>+    (a -> b -> c) -> vec1 -> vec2 -> Seq m c+zipWithSeq f vec1 vec2 = Seq (min l1 l2) $ do+    next1 <- init1+    next2 <- init2+    return $ \i -> f <$> next1 i <*> next2 i+  where+    Seq l1 init1 = toSeq vec1+    Seq l2 init2 = toSeq vec2++unfold :: (Syntax b, MonadComp m) => Data Length -> (b -> (b,a)) -> b -> Seq m a+unfold len step init = Seq len $ do+    r <- initRef init+    return $ \_ -> do+      acc <- getRef r+      let (acc',a) = step acc+      setRef r acc'+      return a++mapAccum' :: (Seqy m vec a, Syntax acc, MonadComp m) =>+    (acc -> a -> (acc,b)) -> acc -> vec -> Seq m (acc,b)+mapAccum' step acc0 vec = Seq len $ do+    next <- init+    r    <- initRef acc0+    return $ \i -> do+      a   <- next i+      acc <- getRef r+      let (acc',b) = step acc a+      setRef r acc'+      acc'' <- getRef r+        -- Read from the reference to avoid duplicating `acc'`+      return (acc'',b)+  where+    Seq len init = toSeq vec++mapAccum :: (Seqy m vec a, Syntax acc, MonadComp m) =>+    (acc -> a -> (acc,b)) -> acc -> vec -> Seq m b+mapAccum step acc0 vec = Seq len $ do+    next <- init+    r    <- initRef acc0+    return $ \i -> do+      a   <- next i+      acc <- getRef r+      let (acc',b) = step acc a+      setRef r acc'+      return b+  where+    Seq len init = toSeq vec++-- | This function behaves slightly differently from the standard @scanl@ for+-- lists:+--+-- * The initial element is not included in the output+-- * Thus, the output has the same length as the input+scan :: (Seqy m vec b, Syntax a, MonadComp m) =>+    (a -> b -> a) -> a -> vec -> Seq m a+scan step acc0 = fmap fst . mapAccum' (\acc a -> (step acc a, ())) acc0+  -- The reason for the discrepancy towards `Data.List.scanl` is that it's+  -- generally not possible for a `Seq` of length `l+1` to read elements from a+  -- `Seq` of length `l` without a conditional in the body.++++-------------------------------------------------------------------------------- -- * Writing to memory -------------------------------------------------------------------------------- @@ -897,6 +1097,7 @@ instance ViewManifest (Manifest a) a where viewManifest = Just instance ViewManifest (Pull a) a instance ViewManifest (Push m a) a+instance ViewManifest (Seq m a) a  class ViewManifest vec a => Manifestable m vec a | vec -> a   where@@ -918,12 +1119,11 @@     -- | A version of 'manifest' that allocates a fresh array for the result     manifestFresh :: Syntax a => vec -> m (Manifest a) -    default manifestFresh :: (Pushy m vec a, Syntax a, MonadComp m) =>+    default manifestFresh :: (Finite vec, Syntax a, MonadComp m) =>         vec -> m (Manifest a)     manifestFresh vec = do-        v   <- toPushM vec-        loc <- newArr $ length v-        manifest loc v+        loc <- newArr $ length vec+        manifest loc vec      -- | A version of 'manifest' that only stores the vector to the given array     -- ('manifest' is not guaranteed to use the array)@@ -943,7 +1143,9 @@  instance MonadComp m             => Manifestable m (Pull a) a instance (MonadComp m1, m1 ~ m2) => Manifestable m1 (Push m2 a) a+instance (MonadComp m1, m1 ~ m2) => Manifestable m1 (Seq m2 a) a + class ViewManifest2 vec a | vec -> a   where     -- | Try to cast a vector to 'Manifest2' directly@@ -975,11 +1177,10 @@     -- | A version of 'manifest2' that allocates a fresh array for the result     manifestFresh2 :: Syntax a => vec -> m (Manifest2 a) -    default manifestFresh2 :: (Pushy2 m vec a, Syntax a, MonadComp m) =>+    default manifestFresh2 :: (Finite2 vec, Syntax a, MonadComp m) =>         vec -> m (Manifest2 a)     manifestFresh2 vec = do-        v   <- toPushM2 vec-        loc <- newArr (numRows v * numCols v)+        loc <- newArr (numRows vec * numCols vec)         manifest2 loc vec      -- | A version of 'manifest2' that only stores the vector to the given array

@@ -187,6 +187,8 @@ simplifyUp (SymP _ Cond :$ LitP _ False :$ _ :$ f) = f simplifyUp (SymP _ Cond :$ c :$ t :$ f) | equal t f = t +simplifyUp (SymP _ BitCompl :$ (SymP _ BitCompl :$ a)) = a+ -- simplifyUp (SymP _ ForLoop :$ LitP _ 0 :$ init :$ _) = init   -- This triggers the bug: <https://ghc.haskell.org/trac/ghc/ticket/11336>. The   -- line below is a workaround:@@ -351,7 +353,7 @@       -- Any unary `Tuple` symbol must be a selector (because there are no       -- 1-tuples).     sharable (SymP _ I2N :$ _) _ = False-    sharable (SymP _ (ArrIx _) :$ _) _ = False+--     sharable (SymP _ (ArrIx _) :$ _) _ = False     sharable _ _ = True  -- | Optimize a Feldspar expression

@@ -0,0 +1,109 @@+module Feldspar.Processing.Filters where++++import Prelude ()++import Feldspar+import Feldspar.Data.Vector+import Feldspar.Data.Queue++++recurrenceI+    :: ( Pushy m fvec a+       , Finite fvec+       , Seqy m ivec a+       , Syntax a, Syntax b+       , MonadComp m+       )+    => fvec           -- ^ Initial input buffer+    -> ivec           -- ^ Input vector+    -> (Pull a -> b)  -- ^ Step function, producing one output from previous inputs+    -> Seq m b        -- ^ Output vector+recurrenceI ii vec body = Seq len $ do+    next <- init+    buf  <- initQueue2 ii+    return $ \i -> do+      a <- next i+      iff (length ii /= 0)+        (putQ buf a)+        (return ())+      withQ buf (return . body)+  where+    Seq len init = toSeq vec++recurrenceIO+    :: ( Pushy m fvec a+       , Finite fvec+       , Seqy m ivec a+       , Pushy m bvec b+       , Finite bvec+       , Syntax a, Syntax b+       , MonadComp m+       )+    => fvec                     -- ^ Initial input buffer+    -> ivec                     -- ^ Input vector+    -> bvec                     -- ^ Initial output buffer+    -> (Pull a -> Pull b -> b)  -- ^ Step function, producing one output from+                                --   previous inputs and outputs+    -> Seq m b                  -- ^ Output vector+recurrenceIO ii vec io body = Seq len $ do+    next <- init+    ibuf <- initQueue2 ii+    obuf <- initQueue2 io+    return $ \i -> do+      a <- next i+      iff (length ii /= 0)+        (putQ ibuf a)+        (return ())+      b <- withQ ibuf $ \ib ->+             withQ obuf $ \ob ->+               shareM $ body ib ob+                 -- Sharing important since `b` is shared+      iff (length io /= 0)+        (putQ obuf b)+        (return ())+      return b+  where+    Seq len init = toSeq vec++-- | FIR filter+fir :: (Pully fvec a, Seqy m ivec a, Syntax a, Num a, MonadComp m)+    => fvec     -- ^ Filter coefficients+    -> ivec     -- ^ Input vector+    -> Seq m a  -- ^ Output vector+fir bs inp = recurrenceI (replicate (length bs) 0) inp $ \i ->+    scProd bs i++-- | IIR filter+iir :: ( Pully bvec a, Pully fvec a+       , Seqy m ivec a, Syntax a+       , Fractional a+       , MonadComp m+       )+    => a        -- ^ First feedback coefficient+    -> bvec     -- ^ Remaining feedback coefficients+    -> fvec     -- ^ Feedforward coefficients+    -> ivec     -- ^ Input vector+    -> Seq m a  -- ^ Output vector+iir a0 as bs inp = recurrenceIO+    (replicate (length bs) 0)+    inp+    (replicate (length as) 0)+    (\i o -> 1 / a0 * (scProd bs i - scProd as o))++-- | FIR filter for 'Pull' vectors+--+-- This version avoids creating a queue for previous inputs (since they are all+-- available anyway).+--+-- Note that each input element is referred many times, so the input should+-- normally be a 'Manifest'. In particular, this means that it is usually not a+-- good idea to compose 'firPull' without writing to memory in between.+firPull :: (Pully vec1 a, Pully vec2 a, Syntax a, Num a)+    => vec1    -- ^ Filter coefficients+    -> vec2    -- ^ Input vector+    -> Pull a  -- ^ Output vector+firPull bs = map (scProd bs . reverse) . tail . inits+

@@ -227,6 +227,21 @@ connectStdIO :: (MarshalFeld a, MarshalFeld b) => (a -> Run b) -> Run () connectStdIO f = (readStd >>= f) >>= writeStd +-- | Connect a Feldspar function between serializable types to @stdin@/@stdout@.+-- The input/output will be in the form of a list as recognized by 'toHaskell' /+-- 'fromHaskell' (i.e. the length followed by the elements in sequence).+--+-- The function will be mapped over the input list in a lazy manner.+streamStdIO :: (MarshalFeld a, MarshalFeld b) => (a -> Run b) -> Run ()+streamStdIO f = do+    n :: Data Length <- readStd+    writeStd n+    for (0,1,Excl n) $ \_ ->+      connectStdIO $ \a -> printf " " >> f a+  -- TODO Better to continue until EOF, but I wasn't able to make this work.+  -- Presumably, one needs to check for EOF inside each `fread` and signal that+  -- somehow.+ -- | A version of 'marshalled' that takes 'ExternalCompilerOpts' as additional -- argument marshalled' :: (MarshalFeld a, MarshalFeld b)@@ -261,4 +276,29 @@          -- ^ Function that has access to the compiled executable as a function     -> IO c marshalled = marshalled' def def++-- | A version of 'marshalledStream' that takes 'ExternalCompilerOpts' as+-- additional argument+marshalledStream' :: (MarshalFeld a, MarshalFeld b)+    => CompilerOpts+    -> ExternalCompilerOpts+    -> (a -> Run b)  -- ^ Function to compile+    -> (([HaskellRep a] -> IO [HaskellRep b]) -> IO c)+         -- ^ Function that has access to the compiled executable as a function+    -> IO c+marshalledStream' opts eopts f body =+    withCompiled' opts eopts (streamStdIO f) $ \g ->+      body $ \is -> do+        parse toHaskell <$> g (fromHaskell is)++-- | Compile a function and make it available as an 'IO' function. The compiled+-- function will be applied repeatedly over the list of inputs producing a list+-- of outputs. Note that compilation only happens once, even if the function is+-- used many times in the body.+marshalledStream :: (MarshalFeld a, MarshalFeld b)+    => (a -> Run b)  -- ^ Function to compile+    -> (([HaskellRep a] -> IO [HaskellRep b]) -> IO c)+         -- ^ Function that has access to the compiled executable as a function+    -> IO c+marshalledStream = marshalledStream' def def 

@@ -26,6 +26,7 @@ import qualified Tut5_Matrices              as Tut5 import qualified Tut6_Testing               as Tut6 import qualified Tut7_ImperativeProgramming as Tut7+import qualified Tut8_SequentialVectors     as Tut8 import qualified Concurrent import DFT import FFT@@ -46,8 +47,8 @@     st <- newStore $ length v     f st v -fftS  = wrapStore fft  :: DManifest (Complex Double) -> _-ifftS = wrapStore ifft :: DManifest (Complex Double) -> _+fftS u = wrapStore (flip fft u)  :: DManifest (Complex Double) -> _+ifftS  = wrapStore (flip ifft 1) :: DManifest (Complex Double) -> _  prop_fft_dft dft' fft' = QC.monadicIO $ do     n   :: Int              <- QC.pick $ QC.choose (2,5)@@ -63,12 +64,22 @@     out2 <- QC.run $ fi out1     QC.assert (inp ~= out2) +prop_fib fb1 fb2 = QC.monadicIO $ do+    n   <- QC.pick $ QC.choose (0,40)+    fs1 <- QC.run $ fb1 n+    fs2 <- QC.run $ fb2 n+    QC.assert (fs1 Prelude.== fs2)+ main =     marshalledM (return . dft)  $ \dft'  ->     marshalledM (return . idft) $ \idft' ->-    marshalledM fftS            $ \fft'  ->+    marshalledM (fftS 1)        $ \fft1  ->+    marshalledM (fftS 2)        $ \fft2  ->     marshalledM ifftS           $ \ifft' -> +    marshalled (return . Tut8.fibSeq) $ \fb1 ->+    marshalled (\n -> return $ Pull n Tut2.fib) $ \fb2 ->+       defaultMain $ testGroup "tests"         [ testCase "Tut1"       Tut1.testAll         , testCase "Tut2"       Tut2.testAll@@ -77,10 +88,13 @@         , testCase "Tut5"       Tut5.testAll         , testCase "Tut6"       Tut6.testAll         , testCase "Tut7"       Tut7.testAll+        , testCase "Tut8"       Tut8.testAll         , testCase "Concurrent" Concurrent.testAll-        , testProperty "fft_dft"  $ prop_fft_dft dft' fft'+        , testProperty "fft1_dft" $ prop_fft_dft dft' fft1+        , testProperty "fft2_dft" $ prop_fft_dft dft' fft2         , testProperty "dft_idft" $ prop_inverse dft' idft'-        , testProperty "fft_ifft" $ prop_inverse fft' ifft'+        , testProperty "fft_ifft" $ prop_inverse fft1 ifft'+        , testProperty "fib"      $ prop_fib fb1 fb2         ]    where

@@ -3,9 +3,10 @@  {-# OPTIONS_GHC -fno-warn-partial-type-signatures #-} -import qualified Prelude+import qualified Prelude as P -import Data.List (genericIndex, genericLength)+import Data.List (genericIndex, genericLength, mapAccumL)+import qualified Data.List as List  import qualified Test.QuickCheck.Monadic as QC @@ -15,11 +16,18 @@  import Feldspar.Run import Feldspar.Data.Array+import Feldspar.Data.Vector+import Feldspar.Data.Queue+import Feldspar.Processing.Filters   -generic_prop_marshalHaskell a = parse toHaskell (fromHaskell a) Prelude.== a+--------------------------------------------------------------------------------+-- * Marshalling+-------------------------------------------------------------------------------- +generic_prop_marshalHaskell a = parse toHaskell (fromHaskell a) P.== a+ prop_marshalHaskell_Int8       a = generic_prop_marshalHaskell (a :: Int8) prop_marshalHaskell_Int16      a = generic_prop_marshalHaskell (a :: Int16) prop_marshalHaskell_Int32      a = generic_prop_marshalHaskell (a :: Int32)@@ -39,8 +47,14 @@  property_marshalFeld f a = QC.monadicIO $ do     b <- QC.run $ f a-    QC.assert (a Prelude.== b)+    QC.assert (a P.== b) +++--------------------------------------------------------------------------------+-- * Nested indexing+--------------------------------------------------------------------------------+ -- | Indexing in a 3-dimensional structure nestIndex     :: ( DIArr Int32@@ -56,10 +70,10 @@     forAll genLenIx $ \(l1,i1) ->     forAll genLenIx $ \(l2,i2) ->     forAll genLenIx $ \(l3,i3) ->-      let l = Prelude.fromIntegral (l1*l2*l3)+      let l = P.fromIntegral (l1*l2*l3)       in  forAll (vector l) $ \(as :: [Int32]) -> QC.monadicIO $ do             a <- QC.run $ f (as,(l1,l2,l3),(i1,i2,i3))-            QC.assert (a Prelude.== genericIndex as (i1*l2*l3 + i2*l3 + i3))+            QC.assert (a P.== genericIndex as (i1*l2*l3 + i2*l3 + i3))   where     genLenIx = do         l <- choose (1,5)@@ -81,16 +95,155 @@     forAll genLenIx $ \(l1,i1) ->     forAll genLenIx $ \(l2,i2) ->     forAll (choose (1,5)) $ \l3 ->-      let l = Prelude.fromIntegral (l1*l2*l3)+      let l = P.fromIntegral (l1*l2*l3)       in  forAll (vector l) $ \(as :: [Int32]) -> QC.monadicIO $ do             bs <- QC.run $ f (as,(l1,l2,l3),(i1,i2))-            QC.assert (genericLength bs Prelude.== l3)+            QC.assert (genericLength bs P.== l3)   where     genLenIx = do         l <- choose (1,5)         i <- choose (0,l-1)         return (l,i) +++--------------------------------------------------------------------------------+-- * Queues+--------------------------------------------------------------------------------++-- | Interpreter for queue operations+--+-- Instructions are read as a list of triples:+--+-- * @(0,i,_)@ means get the element at index @i@ and print it to @stdout@+-- * @(1,_,p)@ put @p@ into the queue+interpQueue :: (Syntax a, MarshalFeld a) => Queue a -> Run ()+interpQueue q = streamStdIO $ \(cmd :: Data Word8, ip, pp) -> do+    r <- newRef+    iff (cmd==0)+      (indexQ q ip >>= setRef r)+      (return ())+    iff (cmd==1)+      (putQ q pp >> indexQ q 0 >>= setRef r) $+      (return ())+    getRef r++-- | Make a queue interpreter using 'newQueue'+interpQ1+    :: Data Length  -- ^ Queue size+    -> Run ()+interpQ1 s = do+    q :: Queue (Data Index) <- newQueue s+    interpQueue q++-- | Make a queue interpreter using 'newQueue2'+interpQ2+    :: Data Length  -- ^ Queue size+    -> Run ()+interpQ2 s = do+    q :: Queue (Data Index) <- newQueue2 s+    interpQueue q++myMarshalled :: (MarshalHaskell a, MarshalHaskell b, MonadRun m)+    => m ()+    -> ((a -> IO b) -> IO c)+    -> IO c+myMarshalled prog body =+    withCompiled prog $ \f ->+      body (fmap (parse toHaskell) . f . fromHaskell)++-- | @(command, index parameter, put parameter)@+type QCMD a = (Word8, Index, a)++genQCMD :: Arbitrary a+    => Length  -- ^ Queue size+    -> Gen [QCMD a]+genQCMD len = do+    NonNegative (n :: Int) <- arbitrary+    P.reverse <$> go n 0 []+  where+    go 0 s prog = return prog+    go n 0 prog = do+      a <- arbitrary+      go (n-1) 1 ((1,0,a):prog)+    go n s prog = do+      c <- choose (0,1)+      i <- choose (0,s-1)+      a <- arbitrary+      let c' = P.fromIntegral c+      go (n-1) (P.min len (s+c')) ((c,i,a):prog)++-- | Reference implementation of a queue interpreter+interpQ_ref :: [QCMD a] -> [a]+interpQ_ref = snd . mapAccumL interp []+  where+    interp q (0,i,_) = (q, genericIndex q i)+    interp q (1,_,a) = (a:q, a)++-- | Test that two queue implementations give the same result+property_queue len qprog1 qprog2 = QC.monadicIO $ do+    cmd <- QC.pick (genQCMD len)+    os1 <- QC.run $ qprog1 cmd+    os2 <- QC.run $ qprog2 cmd+    QC.assert (os1 P.== os2)++++--------------------------------------------------------------------------------+-- * Filters+--------------------------------------------------------------------------------++scProd_ref :: Num a => [a] -> [a] -> a+scProd_ref a b = P.sum $ P.zipWith (*) a b++-- | Reference implementation of FIR filter+fir_ref :: Num a => [a] -> [a] -> [a]+fir_ref cs inp =+    [scProd_ref cs is | is <- P.map P.reverse $ P.tail $ List.inits inp]++fir' :: DPull Int32 -> DPull Int32 -> Run (DSeq Run Int32)+fir' cs = return . fir cs++firPull' :: DPull Int32 -> DPull Int32 -> Run (DPull Int32)+firPull' cs = return . firPull cs++property_fir_ref fr cs as = QC.monadicIO $ do+    bs <- QC.run $ fr (cs,as)+    QC.assert (bs P.== fir_ref cs as)++property_fir_firPull fr frPull cs as = QC.monadicIO $ do+    bs1 <- QC.run $ fr (cs,as)+    bs2 <- QC.run $ frPull (cs,as)+    QC.assert (bs1 P.== bs2)++-- | Reference implementation of IIR filter+iir_ref :: Num a => [a] -> [a] -> [a] -> [a]+iir_ref as bs inp = outp+  where+    inps  = P.map P.reverse $ P.tail $ List.inits inp+    outps = P.map P.reverse $ P.tail $ List.inits (0:outp)++    outp = [scProd_ref bs is - scProd_ref as os | (is,os) <- P.zip inps outps]++-- | Same as 'iir' in "Feldspar.Processing.Filters", but without the+-- 'Fractional' constraint, to avoid rounding errors when testing+iirInt :: DPull Int32 -> DPull Int32 -> DPull Int32 -> DSeq Run Int32+iirInt as bs inp = recurrenceIO+    (replicate (length bs) 0)+    inp+    (replicate (length as) 0)+    (\i o -> scProd bs i - scProd as o)++property_iir_ref ir as bs is = QC.monadicIO $ do+    os <- QC.run $ ir (as,bs,is)+    QC.assert (os P.== iir_ref as bs is)++++--------------------------------------------------------------------------------+-- * Main+--------------------------------------------------------------------------------+ main =     marshalled (return :: _ -> Run (Data Int8))               $ \f_Int8 ->     marshalled (return :: _ -> Run (Data Int16))              $ \f_Int16 ->@@ -112,9 +265,16 @@     marshalled nestIndex       $ \nestIndex_c ->     marshalled nestIndexLength $ \nestIndexLength_c -> +    myMarshalled (interpQ1 5) $ \(iq1 :: [QCMD Index] -> IO [Index]) ->+    myMarshalled (interpQ2 5) $ \iq2 ->++    marshalled (return . uncurry (fir     :: DPull Int32 -> DPull Int32 -> DSeq Run Int32)) $ \fr ->+    marshalled (return . uncurry (firPull :: DPull Int32 -> DPull Int32 -> DPull Int32)) $ \frPull ->+    marshalled (return . \(as,bs,is) -> iirInt as bs is) $ \ir ->+       defaultMain $ testGroup "tests"         [ $testGroupGenerator-        , testGroup "marshalFeld"+        , testGroup "marshalling"             [ testProperty "marshalFeld Int8"           $ property_marshalFeld f_Int8             , testProperty "marshalFeld Int16"          $ property_marshalFeld f_Int16             , testProperty "marshalFeld Int32"          $ property_marshalFeld f_Int32@@ -132,9 +292,18 @@             , testProperty "marshalFeld Pair"           $ property_marshalFeld f_Pair             , testProperty "marshalFeld Nested"         $ property_marshalFeld f_Nested             ]-        , testGroup "misc"+        , testGroup "nested indexing"             [ testProperty "nestIndex"       $ property_nestIndex nestIndex_c             , testProperty "nestIndexLength" $ property_nestIndexLength nestIndexLength_c+            ]+        , testGroup "queues"+            [ testProperty "queue vs. ref" $ property_queue 5 (return . interpQ_ref) iq2+            , testProperty "queue 1 vs. 2" $ property_queue 5 iq1 iq2+            ]+        , testGroup "filters"+            [ testProperty "fir vs. ref"     $ property_fir_ref fr+            , testProperty "fir vs. firPull" $ property_fir_firPull fr frPull+            , testProperty "iir vs. ref"     $ property_iir_ref ir             ]         ]