- Haskell-Cafe - Haskell.org

A Performance Puzzle
by Gregory Wright 04 Aug '18

04 Aug '18

Hi, Something Haskell has lacked for a long time is a good medium-duty linear system solver based on the LU decomposition. There are bindings to the usual C/Fortran libraries, but not one in pure Haskell. (An example "LU factorization" routine that does not do partial pivoting has been around for years, but lacking pivoting it can fail unexpectedly on well-conditioned inputs. Another Haskell LU decomposition using partial pivoting is around, but it uses an inefficient representation of the pivot matrix, so it's not suited to solving systems of more than 100 x 100, say.) By medium duty I mean a linear system solver that can handle systems of (1000s) x (1000s) and uses Crout's efficient in-place algorithm. In short, a program does everything short of exploiting SIMD vector instructions for solving small subproblems. Instead of complaining about this, I have written a little library that implements this. It contains an LU factorization function and an LU system solver. The LU factorization also returns the parity of the pivots ( = (-1)^(number of row swaps) ) so it can be used to calculate determinants. I used Gustavson's recursive (imperative) version of Crout's method. The implementation is quite simple and deserves to be better known by people using functional languages to do numeric work. The library can be downloaded from GitHub: https://github.com/gwright83/luSolve The performance scales as expected (as n^3, a linear system 10 times larger in each dimension takes a 1000 times longer to solve): Benchmark luSolve-bench: RUNNING... benchmarking LUSolve/luFactor 100 x 100 matrix time 1.944 ms (1.920 ms .. 1.980 ms) 0.996 R² (0.994 R² .. 0.998 R²) mean 1.981 ms (1.958 ms .. 2.009 ms) std dev 85.64 μs (70.21 μs .. 107.7 μs) variance introduced by outliers: 30% (moderately inflated) benchmarking LUSolve/luFactor 500 x 500 matrix time 204.3 ms (198.1 ms .. 208.2 ms) 1.000 R² (0.999 R² .. 1.000 R²) mean 203.3 ms (201.2 ms .. 206.2 ms) std dev 3.619 ms (2.307 ms .. 6.231 ms) variance introduced by outliers: 14% (moderately inflated) benchmarking LUSolve/luFactor 1000 x 1000 matrix time 1.940 s (1.685 s .. 2.139 s) 0.998 R² (0.993 R² .. 1.000 R²) mean 1.826 s (1.696 s .. 1.880 s) std dev 93.63 ms (5.802 ms .. 117.8 ms) variance introduced by outliers: 19% (moderately inflated) Benchmark luSolve-bench: FINISH The puzzle is why the overall performance is so poor. When I solve random 1000 x 1000 systems using the linsys.c example file from the Recursive LAPACK (ReLAPACK) library -- which implements the same algorithm -- the average time is only 26 ms. (I have tweaked ReLAPACK's dgetrf.c so that it doesn't use optimized routines for small matrices. As near as I can make it, the C and haskell versions should be doing the same thing.) The haskell version runs 75 times slower. This is the puzzle. My haskell version uses a mutable, matrix of unboxed doubles (from Kai Zhang's matrices <https://hackage.haskell.org/package/matrices> library). Matrix reads and writes are unsafe, so there is no overhead from bounds checking. Let's look at the result of profiling: Tue Jul 31 21:07 2018 Time and Allocation Profiling Report (Final) luSolve-hspec +RTS -N -p -RTS total time = 7665.31 secs (7665309 ticks @ 1000 us, 1 processor) total alloc = 10,343,030,811,040 bytes (excludes profiling overheads) COST CENTRE MODULE SRC %time %alloc unsafeWrite Data.Matrix.Dense.Generic.Mutable src/Data/Matrix/Dense/Generic/Mutable.hs:(38,5)-(39,38) 17.7 29.4 basicUnsafeWrite Data.Vector.Primitive.Mutable Data/Vector/Primitive/Mutable.hs:115:3-69 14.7 13.0 unsafeRead Data.Matrix.Dense.Generic.Mutable src/Data/Matrix/Dense/Generic/Mutable.hs:(34,5)-(35,38) 14.2 20.7 matrixMultiply.\.\.\ Numeric.LinearAlgebra.LUSolve src/Numeric/LinearAlgebra/LUSolve.hs:(245,54)-(249,86) 13.4 13.5 readByteArray# Data.Primitive.Types Data/Primitive/Types.hs:184:30-132 9.0 15.5 basicUnsafeRead Data.Vector.Primitive.Mutable Data/Vector/Primitive/Mutable.hs:112:3-63 8.8 0.1 triangularSolve.\.\.\ Numeric.LinearAlgebra.LUSolve src/Numeric/LinearAlgebra/LUSolve.hs:(382,45)-(386,58) 5.2 4.5 matrixMultiply.\.\ Numeric.LinearAlgebra.LUSolve src/Numeric/LinearAlgebra/LUSolve.hs:(244,54)-(249,86) 4.1 0.3 primitive Control.Monad.Primitive Control/Monad/Primitive.hs:152:3-16 3.8 0.0 basicUnsafeRead Data.Vector.Unboxed.Base Data/Vector/Unboxed/Base.hs:278:813-868 3.3 0.0 basicUnsafeWrite Data.Vector.Unboxed.Base Data/Vector/Unboxed/Base.hs:278:872-933 1.5 0.0 triangularSolve.\.\ Numeric.LinearAlgebra.LUSolve src/Numeric/LinearAlgebra/LUSolve.hs:(376,33)-(386,58) 1.3 0.1 <snip> A large amount of time is spent on the invocations of unsafeRead and unsafeWrite. This is a bit suspicious -- it looks as if these call may not be inlined. In the Data.Vector.Unboxed.Mutable library, which provides the underlying linear vector of storage locations, the unsafeRead and unsafeWrite functions are declared INLINE. Could this be a failure of the 'matrices' library to mark its unsafeRead/Write functions as INLINE or SPECIALIZABLE as well? On the other hand, looking at the core (.dump-simpl) of the library doesn't show any dictionary passing, and the access to matrix seem to be through GHC.Prim.writeDoubleArray# and GHC.Prim.readDoubleArray#. If this program took three to five times longer, I would not be concerned, but a factor of seventy five indicates that I've missed something important. Can anyone tell me what it is? Best Wishes, Greg

4 6

Re: [Haskell-cafe] A Performance Puzzle
by dominic＠steinitz.org 04 Aug '18

04 Aug '18

> Something Haskell has lacked for a long time is a good medium-duty > linear system solver based on the LU decomposition. There are bindings > to the usual C/Fortran libraries, but not one in pure Haskell. (An > example "LU factorization" routine that does not do partial pivoting has > been around for years, but lacking pivoting it can fail unexpectedly on > well-conditioned inputs. Another Haskell LU decomposition using partial > pivoting is around, but it uses an inefficient representation of the > pivot matrix, so it's not suited to solving systems of more than 100 x > 100, say.) > > By medium duty I mean a linear system solver that can handle systems of > (1000s) x (1000s) and uses Crout's efficient in-place algorithm. In > short, a program does everything short of exploiting SIMD vector > instructions for solving small subproblems. > > Instead of complaining about this, I have written a little library that > implements this. It contains an LU factorization function and an LU > system solver. The LU factorization also returns the parity of the > pivots ( = (-1)^(number of row swaps) ) so it can be used to calculate > determinants. I used Gustavson's recursive (imperative) version of > Crout's method. The implementation is quite simple and deserves to be > better known by people using functional languages to do numeric work. > The library can be downloaded from GitHub: > https://github.com/gwright83/luSolve <https://github.com/gwright83/luSolve> Great news :) Dominic Steinitz dominic(a)steinitz.org http://idontgetoutmuch.wordpress.com Twitter: @idontgetoutmuch

1 1

Need some help figuring out how to use the repl with reflex using only stack
by Michael Litchard 04 Aug '18

04 Aug '18

I've created a minimal example to show you what I've tried. Below is also a truncated error log. Could someone help me figure out how to get ghci working with this example? https://github.com/mlitchard/reflex-todo stack repl --stack-yaml=stack-ghc.yaml Truncated error log Error: -- While building custom Setup.hs for package jsaddle-dom-0.9.2.0 using: /home/mlitchard/.stack/setup-exe-cache/x86_64-linux/Cabal-simple_mPHDZzAJ_1.24.2.0_ghc-8.0.2 --builddir=.stack-work/dist/x86_64-linux/Cabal-1.24.2.0 build lib:jsaddle-dom --ghc-options " -ddump-hi -ddump-to-file" Process exited with code: ExitFailure 1 Logs have been written to: /home/mlitchard/projects/git/reflex-todo/.stack-work/logs/jsaddle-dom-0.9.2.0.log Preprocessing library jsaddle-dom-0.9.2.0... [587 of 599] Compiling JSDOM.Custom.XMLHttpRequest ( src/JSDOM/Custom/XMLHttpRequest.hs, .stack-work/dist/x86_64-linux/Cabal-1.24.2.0/build/JSDOM/Custom/XMLHttpRequest.o ) /home/mlitchard/projects/git/reflex-todo/.stack-work/downloaded/Vheiln5kqwE0/src/JSDOM/Custom/XMLHttpRequest.hs:39:46: error: • Could not deduce (Control.Monad.Catch.MonadThrow Language.Javascript.JSaddle.Types.JSM) arising from a use of ‘throwM’ from the context: MonadDOM m bound by the type signature for: throwXHRError :: MonadDOM m => Maybe XHRError -> m () at src/JSDOM/Custom/XMLHttpRequest.hs:38:1-53 • In the second argument of ‘(.)’, namely ‘throwM’ In the second argument of ‘maybe’, namely ‘(liftDOM . throwM)’ In the expression: maybe (return ()) (liftDOM . throwM) /home/mlitchard/projects/git/reflex-todo/.stack-work/downloaded/Vheiln5kqwE0/src/JSDOM/Custom/XMLHttpRequest.hs:42:20: error: • No instance for (Control.Monad.Catch.MonadMask Language.Javascript.JSaddle.Types.JSM) arising from a use of ‘bracket’ • In the first argument of ‘(.)’, namely ‘bracket aquire id’ In the expression: bracket aquire id . const In an equation for ‘withEvent’: withEvent aquire = bracket aquire id . const /home/mlitchard/projects/git/reflex-todo/.stack-work/downloaded/Vheiln5kqwE0/src/JSDOM/Custom/XMLHttpRequest.hs:45:31: error: • Could not deduce (Control.Monad.Catch.MonadCatch Language.Javascript.JSaddle.Types.JSM) arising from an operator section from the context: MonadDOM m bound by the type signature for: send' :: MonadDOM m => XMLHttpRequest -> Maybe JSVal -> m () at src/JSDOM/Custom/XMLHttpRequest.hs:44:1-62 • In the expression: (`onException` abort self) In the second argument of ‘($)’, namely ‘(`onException` abort self) $ do { result <- liftIO newEmptyMVar; r <- withEvent (onAsync self Generated.error . liftIO $ putMVar result (Just XHRError)) $ withEvent (onAsync self abortEvent . liftIO $ putMVar result (Just XHRAborted)) $ withEvent (onAsync self load . liftIO $ putMVar result Nothing) $ do { void $ case mbVal of { Nothing -> ... Just val -> ... }; liftIO $ takeMVar result }; throwXHRError r }’ In the expression: liftDOM $ (`onException` abort self) $ do { result <- liftIO newEmptyMVar; r <- withEvent (onAsync self Generated.error . liftIO $ putMVar result (Just XHRError)) $ withEvent (onAsync self abortEvent . liftIO $ putMVar result (Just XHRAborted)) $ withEvent (onAsync self load . liftIO $ putMVar result Nothing) $ do { void $ case mbVal of { Nothing -> ... Just val -> ... }; liftIO $ takeMVar result }; throwXHRError r }

4 4

How to turn Integer into Proxy (n :: KnownNat) with n equal to the integer?
by David Banas 04 Aug '18

04 Aug '18

I’d like to use Data.Vector.Sized.slice’, but it wants the starting position and length parameters as Proxys, not Integers. How can I convert an Integer into a Proxy (n :: KnownNat) with n equal to the integer?

2 1

[random-fu] sampleRVar with StdGen in RWS?
by Sergiu Ivanov 03 Aug '18

03 Aug '18

Dear Café, I would like to use random-fu to do some pseudo-random simulations for a given StdGen (so that I can run the same simulation multiple times, if needed). The following works: testState :: StdGen -> Int testState = evalState (sampleRVar $ uniform 1 10) The following doesn't: testRWS :: StdGen -> Int testRWS = fst . evalRWS (sampleRVar $ uniform 1 10) () I get <interactive>:2:26: error: • No instance for (MonadRandom (RWST () b0 StdGen Identity)) arising from a use of ‘sampleRVar’ • In the first argument of ‘evalRWS’, namely ‘(sampleRVar $ uniform 1 10)’ In the second argument of ‘(.)’, namely ‘evalRWS (sampleRVar $ uniform 1 10) ()’ In the expression: fst . evalRWS (sampleRVar $ uniform 1 10) () Indeed, I do see a MonadRandom instance for StateT, but none for RWST [0]. Is there a reason to not have a MonadRandom instance for RWST? Am I looking in the wrong place? - Sergiu [0] https://hackage.haskell.org/package/random-source-0.3.0.6/docs/Data-Random-…

2 4

hp2pretty-0.9
by Claude Heiland-Allen 03 Aug '18

03 Aug '18

Hi all, hp2pretty[0] graphs heap profiles as SVG, much like hp2ps graphs as PostScript. I release hp2pretty-0.9[1], with support for detached keys in a separate text file with individual label images, which lets you construct a single key for multiple graphs, and also workaround the issue of truncated names (running off the edge of the image) when using the +RTS -L option. I applied the same principle to the titles, they can now be omitted from the images and added as captions in a surrounding document if desired. The examples/ folder shows how to embed the output in PDF[2] using the LaTeX 'svg' package (which calls out to 'inkscape'). There's also an HTML5[3] example. Because one can now create PDF output from SVG, the vapourware plans for other output formats have been dropped. Get it from Hackage[4], use `cabal unpack hp2pretty` to access examples. Or browse the git repository[5]. Thanks, Claude [0] https://mathr.co.uk/blog/hp2pretty.html [1] https://mathr.co.uk/blog/2018-08-01_hp2pretty-0.9_released.html [2] https://mathr.co.uk/hp2pretty/0.9/report.pdf [3] https://mathr.co.uk/hp2pretty/0.9/report.html [4] https://hackage.haskell.org/package/hp2pretty [5] https://code.mathr.co.uk/hp2pretty -- https://mathr.co.uk

1 0

why different type variables enforce different types?
by Rui Azevedo 02 Aug '18

02 Aug '18

after using haskell for a while i was surprised by this little thing f::a->b f a=a and this is not making sense to me till now, so please can someone enlighten me, why does haskell assume that b can not be the same type of a? after all having two distinct math variables does not mean they can not have the same value. I know i could use f::a->a, but its not the same thing, on this last case I'm requesting them to be the same type, but on the first I'm not forbidding it, is there a math reason behind it? would type inference be rendered impossible if considering the case of them being equal? thanks Rui Azevedo (neu-rah)

3 2

JOB | Permanent Java Developer (Montreal, Canada)
by James Tobin 01 Aug '18

01 Aug '18

Hello, I'm working with an employer that is looking to hire (for their Montreal office) a permanent Java developer that has extensive hands-on Java coding experience. Consequently I had hoped that some members of this mailing list may like to discuss with me further; off-list. I can be reached using "JamesBTobin (at) Gmail (dot) Com". Kind regards, James

4 3

JOB | Permanent Java Development Manager (Montreal, Canada)
by James Tobin 01 Aug '18

01 Aug '18

Hello, I'm working with an employer that is looking to hire (for their Montreal office) a permanent development manager that has extensive hands-on Java coding experience. Consequently I had hoped that some members of this mailing list may like to discuss with me further; off-list. I can be reached using "JamesBTobin (at) Gmail (dot) Com". Kind regards, James

1 0

Missing common function for Maybe types
by Marc Busqué 31 Jul '18

31 Jul '18

Hi! I have two functions: ``` foo :: a -> Maybe b bar :: a -> Maybe c ``` >From which I want to build a higher order function: ``` foobar :: a -> (a -> Maybe b) -> (a -> Maybe c) -> Either b c ``` The implementation I need is: ``` foobar x f g = case (f x) of Nothing -> g x Just y -> Just y ``` I'm a bit surprised that looking at hoogle I don't find a built-in solution for this quite common need for `Maybe` types (or perhaps for any monad). Am I looking in the wrong way? Does it exist a similar abstraction but with a different shape? Thanks in advance, Marc Busqué http://waiting-for-dev.github.io/about/

7 16