Re: [Haskell-cafe] Acquiring a random set of a specific size (w/o dups) from a range of Ints

[Patrick Perry]

Correction: the version in monte-carlo takes time O(n), but it only consumes k 
random numbers.  The streaming algorithm consumes n random numbers.

On Jun 14, 2011, at 12:04 PM, Patrick Perry wrote:

> There are time/space tradeoffs in sampling without replacement.  The version 
> in monte-carlo takes space O(n) and time O(k), for all k.  I chose this 
> algorithm instead of a streaming algorithm, with takes space O(k) and time 
> O(n).
> 
> If k is much less than n, you can improve a bit.  Here is a version taking 
> space O(k) and time O(k^2), *provided k << n*;  when k is large relative to 
> n, the time increases to roughly O(n log n).  
> 
> import Control.Monad( liftM )
> import Control.Monad.MC
> import Data.List( nub )
> 
> sampleSmallIntSubset :: MonadMC m => Int -> Int -> m [Int]
> sampleSmallIntSubset n k = (take k . nub) `liftM` repeatMC (sampleInt n)
> 
> sampleSmallSubset :: MonadMC m => [a] -> Int -> m [a]
> sampleSmallSubset xs k = (map (xs !!)) `liftM` sampleSmallIntSubset (length 
> xs) k
> 
> -- strict version of the above
> sampleSmallSubset' :: MonadMC m => [a] -> Int -> m [a]
> sampleSmallSubset' xs k = do
>ys <- sampleSmallSubset xs k
>(length ys) `seq` return ys -- by forcing the length, we make sure we are 
> done with the RNG when we return
> 
> 
> If you're wondering about the strict version, the following example should be 
> instructive:
> 
> Let's note the first three numbers generated by the RNG with seed 0:
> 
>ghci> (replicateM 3 $ uniform 0 1) `evalMC` mt19937 0
>[0.999741748906672,0.16290987539105117,0.28261780529282987]
> 
> Now, let's look at the first number we get after a call to sampleSmallSubset, 
> if we don't evaluate the subset:
> 
>ghci> sampleSmallSubset [1..1] 2 >> uniform 0 1) `evalMC` mt19937 0
>0.999741748906672
> 
> Apparently, the call to sampleSmallSubset doesn't consume any random numbers. 
>  Now, let's try the strict version:
> 
>ghci> (sampleSmallSubset' [1..1] 2 >> uniform 0 1) `evalMC` mt19937 0
>0.28261780529282987
> 
> This time, the call to sampleSmallSubset' consumes two random numbers, even 
> if we throw away the result.
> 
> Since you're just starting out, I recommend using the strict versions of the 
> functions.  They are less space-efficient than the lazy versions, but they 
> are guaranteed to be safe.
> 
> 
> 
> 
>> From: Jonas Almström Duregård 
> 
>> Subject: Re: [Haskell-cafe] Acquiring a random set of a specific size (w/o 
>> dups) from a range of Ints
>> To: "michael rice" <
>> nowgate at yahoo.com
>>> 
>> Cc: "Felipe Almeida Lessa" <
>> felipe.lessa at gmail.com>, haskell-cafe at haskell.org
>> 
>> Date: Tuesday, June 14, 2011, 5:17 AM
>> 
>>> 
>> Shuffle [1..20], then take 5?
>> 
>>> 
>> Yes, so simple, I'm embarrassed I didn't think of it.
>> 
>> 
>> That works well for small numbers, but I'm guessing it will evaluate the 
>> entire list so it should not be used for large inputs. If you have a large 
>> interval and use a relatively small part of it, the following function 
>> should be significantly faster (it builds a random permutation lazily):


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] Acquiring a random set of a specific size (w/o dups) from a range of Ints

[Patrick Perry]

There are time/space tradeoffs in sampling without replacement.  The version in 
monte-carlo takes space O(n) and time O(k), for all k.  I chose this algorithm 
instead of a streaming algorithm, with takes space O(k) and time O(n).

If k is much less than n, you can improve a bit.  Here is a version taking 
space O(k) and time O(k^2), *provided k << n*;  when k is large relative to n, 
the time increases to roughly O(n log n).  

import Control.Monad( liftM )
import Control.Monad.MC
import Data.List( nub )

sampleSmallIntSubset :: MonadMC m => Int -> Int -> m [Int]
sampleSmallIntSubset n k = (take k . nub) `liftM` repeatMC (sampleInt n)

sampleSmallSubset :: MonadMC m => [a] -> Int -> m [a]
sampleSmallSubset xs k = (map (xs !!)) `liftM` sampleSmallIntSubset (length xs) 
k

-- strict version of the above
sampleSmallSubset' :: MonadMC m => [a] -> Int -> m [a]
sampleSmallSubset' xs k = do
ys <- sampleSmallSubset xs k
(length ys) `seq` return ys -- by forcing the length, we make sure we are 
done with the RNG when we return


If you're wondering about the strict version, the following example should be 
instructive:

Let's note the first three numbers generated by the RNG with seed 0:

ghci> (replicateM 3 $ uniform 0 1) `evalMC` mt19937 0
[0.999741748906672,0.16290987539105117,0.28261780529282987]

Now, let's look at the first number we get after a call to sampleSmallSubset, 
if we don't evaluate the subset:

ghci> sampleSmallSubset [1..1] 2 >> uniform 0 1) `evalMC` mt19937 0
0.999741748906672

Apparently, the call to sampleSmallSubset doesn't consume any random numbers.  
Now, let's try the strict version:

ghci> (sampleSmallSubset' [1..1] 2 >> uniform 0 1) `evalMC` mt19937 0
0.28261780529282987

This time, the call to sampleSmallSubset' consumes two random numbers, even if 
we throw away the result.

Since you're just starting out, I recommend using the strict versions of the 
functions.  They are less space-efficient than the lazy versions, but they are 
guaranteed to be safe.




> From: Jonas Almström Duregård 

> Subject: Re: [Haskell-cafe] Acquiring a random set of a specific size (w/o 
> dups) from a range of Ints
> To: "michael rice" <
> nowgate at yahoo.com
> >
> Cc: "Felipe Almeida Lessa" <
> felipe.lessa at gmail.com>, haskell-cafe at haskell.org
> 
> Date: Tuesday, June 14, 2011, 5:17 AM
> 
> >
>  Shuffle [1..20], then take 5?
> 
> >
>  Yes, so simple, I'm embarrassed I didn't think of it.
> 
> 
> That works well for small numbers, but I'm guessing it will evaluate the 
> entire list so it should not be used for large inputs. If you have a large 
> interval and use a relatively small part of it, the following function should 
> be significantly faster (it builds a random permutation lazily):

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: ieee754 version 0.7.1

[Patrick Perry]

Upon popular request, I've renamed the "ieee" package to "ieee754".  The new 
hackage page is http://hackage.haskell.org/package/ieee754 .



On Sep 19, 2010, at 11:16 PM, Conrad Parker wrote:

> On 20 September 2010 11:18, Patrick Perry  wrote:
>>> Given that IEEE is actually a standards body and they have many
>>> standards, wouldn't it be more appropriate to call this library
>>> ieee754?
>> 
>> If it seems important to people, I'd be happy to change the name.  I'm
>> not religious about these things.  Will it clutter up hackage, though?
> 
> I reckon it's worth making it obvious that this library does 754 and
> not, say, 1394 or 802.11 ;-) On the other hand if you intend on
> expanding the package to implement every IEEE standard ... (j/k)
> 
> Anyway, good work. Does this have any overlap with
> data-binary-ieee754? There was some recent discussion here about the
> encoding speed in that package.
> 
> Conrad.

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] ANN: ieee version 0.7

[Patrick Perry]

> I needed "real" IEEE754 binary support for round-trip parsing, where
> (for example) maintaining the particular bit pattern of a NaN is
> important. For 99% of people, the "unsafe" method will work fine.

How does a C-style cast not preserve the bit pattern of a NaN?  Again,
sorry if this is a stupid question.
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] ANN: ieee version 0.7

[Patrick Perry]

I'm looking at data-binary-754 right now, and it seems pretty
complicated.  Why don't they just cast the values bitwise to integers
and then serialize those?  Forgive my naivete-- I don't know much
about binary encoding issues.

I went the route of implementing everything in C and then using the
FFI.  There's a lot of lot of bit-twiddling involved when you work
with the guts of IEEE754, which is a lot easier to do in C.  See
http://github.com/patperry/hs-ieee/blob/master/cbits/feqrel_source.c
for a truly ugly example.  I ported this code from the Tango math
library for D (
http://www.dsource.org/projects/tango/browser/trunk/tango/math/IEEE.d
).  It's original author, Don Clugston, claims that the function is
about as fast as a ">" comparison.  Haskell's great, but I don't think
it could get nearly as fast for a function like this.


Patrick

On Sun, Sep 19, 2010 at 11:16 PM, Conrad Parker  wrote:
> On 20 September 2010 11:18, Patrick Perry  wrote:
>>> Given that IEEE is actually a standards body and they have many
>>> standards, wouldn't it be more appropriate to call this library
>>> ieee754?
>>
>> If it seems important to people, I'd be happy to change the name.  I'm
>> not religious about these things.  Will it clutter up hackage, though?
>
> I reckon it's worth making it obvious that this library does 754 and
> not, say, 1394 or 802.11 ;-) On the other hand if you intend on
> expanding the package to implement every IEEE standard ... (j/k)
>
> Anyway, good work. Does this have any overlap with
> data-binary-ieee754? There was some recent discussion here about the
> encoding speed in that package.
>
> Conrad.
>
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] ANN: ieee version 0.7

[Patrick Perry]

> Given that IEEE is actually a standards body and they have many
> standards, wouldn't it be more appropriate to call this library
> ieee754?

If it seems important to people, I'd be happy to change the name.  I'm
not religious about these things.  Will it clutter up hackage, though?


Patrick
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: ieee version 0.7

[Patrick Perry]

ieee is a Haskell library for dealing with IEEE floating numbers.  It
was originally written to make testing with floating point values less
painful.  The library provides an "approximate equality" type class,
AEq, with approximate equality operator "~==".  One property of "~=="
is that nearby floating point numbers are deemed equivalent, so that,
for example, "1 ~== 1.0001" evaluates to True.

Documentation is on hackage: http://hackage.haskell.org/package/ieee

Changes since the last release:

* Add IEEE type class with instances for Double, Float, CDouble, and CFloat

* Add functions for getting/setting NaN payloads

* Add succIEEE/predIEEE for advancing up and down the IEEE number line
  (ported from Tango's nextUp and nextDown)

* Add bisectIEEE for midpoints of two numbers on the IEEE number line
  (ported from Tango's ieeeMean)

* Add identicalIEEE for exact (bitwise) equality of IEEE numbers

* Add copySign for setting the sign bit of an IEEE number

* Add sameSignificandBits for seeing how many significand bits of two
  IEEE numbers agree, ported from Tango's feqrel

* Add nan, infinity, maxFinite, minNormal constants for IEEE numbers

* Add maxNum and minNum

* Rename maxF and minF to maxNaN and minNaN

* Switch to a simpler "~==" comparison for complex numbers

* Make "~==" comparison use sameSignificandBits for IEEE types

* Make "===" comparison use bitwise equality for IEEE types

* Remove old "eqRel" comparisons.

* Remove old epsilon' and delta constants

* Remove (RealFloat a) => AEq (Complex a) instance in favor of explicit
  instances for Complex {Double,Float,CDouble,CFloat}
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: monte-carlo version 0.3

[Patrick Perry]

monte-carlo is a haskell library providing a monad and associated monad 
transformer for computing with quasi-random numbers.  It provides a high-level 
interface to the distributions supported by the GNU Scientific Library 
(currently just uniform, normal, exponential, Poisson, Levy, but others are 
easy to add).  You can see examples of the library in action at 
http://github.com/patperry/hs-monte-carlo/tree/master/examples/ .   Hackage 
documentation is at http://hackage.haskell.org/package/monte-carlo .

A few blog posts describe an older version of the library, but beware the API 
has changed slightly, and these examples require some modification:

http://ptrckprry.com/blog/2008/08/26/a-monte-carlo-monad-for-haskell/ (see 
examples/Pi.lhs)
http://ptrckprry.com/blog/2008/12/31/monte-carlo-poker-odds/ (see 
examples/Poker.hs)

The changes since the last release are listed below.  I welcome feedback, bug 
reports, and contributions.


Patrick


Changes since version 0.2.

* Add strict versions of sampleSubset, sampleIntSubset, and shuffleInt.

* Port to vector-0.6.0.

* Add Exponential and Levy alpha-Stable distributions.

* Add Summary.Bool for indicators.

* Move Summary to Data.Summary

* Introduce `repeatMC`, which produces an infinite (lazy) stream of values, and
  `replicateMC`, which produces a lazy list of specified length.

* Remove `repeatMC/repeatMCWith`.

* Build fix for 6.8.2 from Robert Gunst.

* The function `sample`, `sampleWithWeights`, `sampleSubset`, and
  `shuffle` no longer require that you explicitly pass in the length.

* The pure RNG is now a newtype, so you can't use the functions from
  GLS.Random.Gen on it anymore.
  
* The internals of the monad have been cleaned up.  IO is used internally
  instead of `seq` calls and `unsafePerformIO` everywhere.  This results in
  a modest performance boost.

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] QuickCheck2 question

[Patrick Perry]

I'm having some trouble with QuickCheck2 and Control.Applicative.  
Specifically, I have the following functions (slightly simplified):

copyVector :: IOVector -> Vector -> IO ()
freezeVector :: IOVector -> IO (Vector)

I have a test, part of which looks like

monadicIO $ do
run $ copyVector <$> pure dst <*> freezeVector src
-- more stuff...

When I try running the test (using test-framework-quickcheck2), the 
`copyVector` function never gets called.

If I re-write the test as

monadicIO $ do
src' <- run $ freezeVector
run $ copyVector dst src'
-- etc.

then everything works fine (i.e., the call to `copyVector` gets executed).  

Does anyone have any clues as to what is going on?  I suspect the problem is 
related to Test.QuickCheck.Monadic using `unsafePerformIO` without a `{-# 
NOINLINE -#}` pragma [1], but I'm not completely sure.  Adding such a pragma to 
`monadicIO` and `run` in Monadic.hs does not fix the problem.

I'm using Quickcheck-2.1.0.2 and GHC 6.10.1.  

Thanks in advance for any help,


Patrick

[1] 
http://hackage.haskell.org/packages/archive/QuickCheck/2.1.0.2/doc/html/src/Test-QuickCheck-Monadic.html
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] fast Eucl. dist. - Haskell vs C

[Patrick Perry]

Hi Kenneth,

I wrote a benchmark similar to yours using the haskell blas library I  
wrote (latest version is on github at http://github.com/patperry/blas/tree/master 
, an older version is on hackage).


The pure Haskell implementation is not very good, and the problem  
seems to be repeated boxing/unboxing of doubles.  It may be possible  
to avoid this by writing a custom "fold" function for Vectors, but I  
haven't tried it.  Using a vectorized subtraction function along with  
the BLAS1 "dnrm2" function, I was able to get performance about  
2.8-5.4x worse than C.  This isn't great, but it's not terrible,  
either.  I didn't have the performance problems you ran in to using  
the FFI.  My benchmark includes:


distVectors1: a native implementation using foldl' and zipWith
distVectors2: a vectorize implementation, defined as @distVectors2 x y  
= norm2 (x-y)@

distVectors3: a custom loop in haskell
distVectors4: an FFI call to C

All 4 functions use the "Vector" type, provided by the "blas" package.

I compiled with "-O -fexcess-preceision" and ran timings for vectors  
of size 10, 100, and 1000.  Here are the results:


Macintosh-62:~/Code/blas/examples (euclidean)$ ./Euclidean 10
n: 10
* distVectors1: .
  1.366ns per iteration / 732029.06 per second.
* distVectors2: .
  1.286ns per iteration / 34.02 per second.
* distVectors3: ..
  0.950ns per iteration / 1052847.93 per second.
* distVectors4: ...
  0.452ns per iteration / 2213752.56 per second.
Macintosh-62:~/Code/blas/examples (euclidean)$ ./Euclidean 100
n: 100
* distVectors1: ..
  12.496ns per iteration / 80025.79 per second.
* distVectors2: 
  2.149ns per iteration / 465311.61 per second.
* distVectors3: ..
  8.728ns per iteration / 114572.53 per second.
* distVectors4: ..
  0.597ns per iteration / 1675765.65 per second.
Macintosh-62:~/Code/blas/examples (euclidean)$ ./Euclidean 1000
n: 1000
* distVectors1: ..
  125.920ns per iteration / 7941.52 per second.
* distVectors2: ..
  10.677ns per iteration / 93661.99 per second.
* distVectors3: ...
  85.533ns per iteration / 11691.42 per second.
* distVectors4: 
  1.962ns per iteration / 509668.17 per second.

The "distVectors2" function is computing the norm in a more  
numerically-stable way.  It is also allocating a temporary array to  
store x-y.


I haven't looked over your benchmarking code too thoroughly, but parts  
of it seem suspect to me.  Namely, you are using "mapM" instead of  
"mapM_".  Also, it seems like you are including GC times in your  
benchmarks.


Attached is my code.  You will need the latest "blas" from github to  
compile it.  You also might want to look into hmatrix.  I have very  
little time for haskell code right now (I'm trying to finish my Ph.D.  
in the next month).  This summer I plan to work more on this stuff,  
and I'll probably make a release in July.



Patrick



Euclidean.hs
Description: Binary data





dist.c
Description: Binary data


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] Re: Matrices

[Patrick Perry]

Hi Cetin,

This is probably the easiest way:

> (m - es)^2/es
listMatrix (2,2)  
[0.6163270413689804,0.600918865334756,0.33217626255600896,0.3238718559921087 
]


This will create 2 temporary arrays.  Alternatively,

> let res = listMatrix (shape m) [ (o-e)^2 / e | o <- colElems m | e  
<- colElems es ] where colElems = concatMap elems . cols

> res
listMatrix (2,2)  
[0.6163270413689804,0.600918865334756,0.33217626255600896,0.3238718559921087 
]


In a later version of the library, "colElems" will probably be built- 
in.  This version has the disadvantage that it doesn't use BLAS calls.


If you *really* want to get tricky and eliminate temporary arrays:

> runSTMatrix $ do
res <- newCopyMatrix m
subMatrix res es
modifyWith (^2) res
divMatrix res es
return res


Now, to get the sum:

> sumAbs $ vectorFromMatrix res
1.8732940252518542

or

> sum $ elems res
1.8732940252518542

The first version will usually be more efficient, since it calls the  
BLAS1 function "dasum".



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] hmatrix, Windows and GCC

[Patrick Perry]

Well, I guess I am not the only one!

This blog show exactly what I am looking for!

http://quantile95.com/2008/10/31/ann-blas-bindings-for-haskell-version-06/


If you're looking to implement other linear algebra algorithms in  
Haskell, there's enough at


http://github.com/patperry/lapack

to implement a QR decomposition with column pivoting with fairly  
minimal effort.  Householder reflections and permutation matrices are  
already supported. You will also need the latest BLAS bindings, also  
available at github.  Note that these libraries are not compatible  
with hmatrix.


What would *really* be awesome would be an implementation of an SVD  
algorithm: take a bidiagonal matrix and produce a (lazy) list of  
Givens rotations that diagonalizes the matrix.  That would probably be  
a 2-4 week project, but would be really useful.  Take a look at the  
code and references at http://netlib.org/lapack/explore-html/zbdsqr.f.html 
 if you're interested.


The LAPACK version, zbdsqr, takes "B" and factors it as "B = Q S  
P^H".  Optionally, it sets


U  := U Q
VT := P^H VT
C  := Q^H C

Annoyingly, there is no way to specify an input matrix "D" and set  
"D := D P".


It'd be really cool to have a Haskell function of type

svdBidiagonal :: (WriteBanded a m) => a (n,p) e -> m  
[(Givens,Givens)]


The type signature means that "a" is a mutable banded matrix of shape  
"(n,p)" with elements of type "e" that can be modified in monad "m".   
I'm envisioning that the result is a lazy list, "gs", and that no  
computation gets done until the "gs" list is traversed.  So, if  
someone just wants singular values, then they would do "liftM length .  
svdBidiagonal".  If they want to update "U := U Q", then they would do  
something like:


   mapM_ (applyGivens u) $ liftM (fst . unzip) $ svdDidiagonal a

The user has the option to update as many matrices they want in  
whatever way they want.  This is only possible with laziness.  It is a  
perfect example of the "glue" John Hughes talks about in "Why  
Functional Programming Matters".



Patrick



On Wed, Jan 28, 2009 at 08:21, Rafael Gustavo da Cunha Pereira Pinto <
RafaelGCPP.Linux at gmail.com> wrote:

> I was planning to recompile everything (ATLAS, LAPACK and GHC  
included)
> this weekend, so I can have a similar environment on Windows and  
Linux...

> Having to "borrow" libraries
>
> Since I am married, this means it will actually happen on some  
weekend till

> 2010.
>
>
> What I really would like to try is a (purely?) functional approach  
to
> create a (P)LU decomposition of a matrix. I am not too much  
worried (at
> first) with performance or memory constraints, since I only want  
to see how
> beautiful it gets (or not!).  (This one might happen somewhere in  
this

> century...)
>
>
> Thanks anyway
>
>
>
> On Wed, Jan 28, 2009 at 07:57, allan  wrote:
>
>> Hi
>>
>> The INSTALL file in the hmatrix repository has some very clear
>> instructions for installation on Windows.
>> http://perception.inf.um.es/~aruiz/darcs/hmatrix/INSTALL

>>
>> However note this section at the bottom:
>> "Unfortunately the lapack dll supplied by the R system does not  
include

>> zgels_, zgelss_, and zgees_, so the functions depending on them
>> (linearSolveLS, linearSolveSVD, and schur for complex data)
>> will produce a "non supported in this OS" runtime error."
>>
>> Of course linearSolve is exactly what you will be wanting so this  
won't

>> work for you.
>> I ran into exactly this problem myself. I actually didn't get as  
far as a

>> run-time error as I got a linker error.
>>
>> I don't have any solution for you though, sorry.
>>
>> regards
>> allan
>>
>>
>>
>>
>> Rafael Gustavo da Cunha Pereira Pinto wrote:
>>
>>>
>>>
>>>   Hi all,
>>>
>>> I am writing a program that uses hmatrix for solving some linear  
systems.
>>> The hmatrix package depends on BLAS, which, in turn, depend on  
GCC 4.2 to be

>>> built (at least ATLAS does).
>>>
>>> GHC 6.10 for Windows is pre-packaged with GCC 3.4.5, and it  
leaves me

>>> with the impression that I would have incompatible ABIs.
>>>
>>> My questions:
>>>
>>> 1) Why GHC 6.10 still uses GCC 3.4.5 in Windows? I know mingw  
considers

>>> GCC 4.2 to be alpha, but, lets face it, 4.2 is almost obsolete!
>>> 2) Is it possible to rebuild GHC 6.10, using Windows and GCC  
4.2? Is

>>> there any guide for doing this?
>>> 3) Has any of you tried hmatrix on Windows? How did you do it?
>>>
>>> Thanks,
>>>
>>> Rafael
>>>
>>> --
>>> Rafael Gustavo da Cunha Pereira Pinto
>>> Electronic Engineer, MSc.
>>>
>>>
>>>  


>>>
>>> ___
>>> Haskell-Cafe mailing list
>>> Haskell-Cafe at haskell.org
>>> http://www.haskell.org/mailman/listinfo/haskell-cafe
>>>
>>
>>
>>
>> --
>> The University of Edinburgh is a charitable body, registered in
>> Scotland, with

[Haskell-cafe] Question about allocaArray and friends

[Patrick Perry]

In C, if you try to "alloca" too much memory, then the stack gets  
overwritten and bad things happen.  Does GHC exhibit the same behavior  
with allocaArray and the like?  Is there a way to find out how much is  
safe to allocate?


Thanks in advance for the help,


Patrick


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: monad-interleave 0.1

[Patrick Perry]

My two favorite functions in Haskell are "unsafeInterleaveIO" and  
"unsafeInterleaveST".  I'm surprised there isn't a type class for  
them.  I just fixed this by adding the "monad-interleave" package to  
hackage:


http://hackage.haskell.org/cgi-bin/hackage-scripts/package/monad-interleave

The package adds a "Control.Monad.Interleave" module.  Here is the  
entirety of the module:


\begin{code}
-- | Monads that have an operation like 'unsafeInterleaveIO'.
class Monad m => MonadInterleave m where
-- | Get the baton from the monad without doing any computation.
unsafeInterleave :: m a -> m a

instance MonadInterleave IO where
unsafeInterleave = unsafeInterleaveIO
{-# INLINE unsafeInterleave #-}

instance MonadInterleave (ST s) where
unsafeInterleave = unsafeInterleaveST
{-# INLINE unsafeInterleave #-}

instance MonadInterleave (Lazy.ST s) where
unsafeInterleave = Lazy.unsafeInterleaveST
{-# INLINE unsafeInterleave #-}
\end{code}

Use it if you need it.


Patrick


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] Question about touchForeignPtr

[Patrick Perry]

Thanks for your help, Duncan.

On Jan 12, 2009, at 6:10 AM, Duncan Coutts wrote:


Do the (ForeignPtr e) and the (Ptr e) point to the same thing? They

appear to be related because you dereference p but touch f.

It used to be the ForeignPtr was slower to dereference than a Ptr  
and so

caching it like this could make it faster (but then requires lots of
careful thought about touchForeignPtr). These days ForeignPtr is  
just as

fast and so you can use withForeignPtr and never have to worry about
touchForeignPtr. You'll notice ByteString works this way.


Often they point to the same thing, but not always.  If they pointed to
the same thing, I wouldn't store p separately, but often p has a nonzero
offset from f.  You either have to store the offset or cache the  
pointer.

I found that it was both faster and simpler to cache p rather than the
offset.


Patrick
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] Question about touchForeignPtr

[Patrick Perry]

I have the following code:

IOVector n e = IOVector !ConjEnum !Int (ForeignPtr e)! (Ptr e)! Int!
newtype Vector n e = IOVector n e

unsafeAtVector :: Vector n e -> Int -> e
unsafeAtVector (Vector (IOVector c _ f p inc)) i =
let g = if c == Conj then conjugate else id
in inlinePerformIO $ do
   e  <- peekElemOff p (i*inc)
   io <- touchForeignPtr f
   let e' = g e
   e' `seq` io `seq` return e'
{-# INLINE unsafeAtVector #-}


The Ptr, 'p' is derived from the ForeignPtr, 'f'. For some offset,  
'o', it

is defined as:

p = unsafeForeignPtrToPtr f `advancePtr` o

The "touchForeignPtr" is there to keep the garbage collector from  
deallocating
the memory before we have a chance to read 'e'.  My question is the  
following:
Is the `seq` on `io` necessary (from a safety standpoint)?  Or am I  
just being

paranoid?

Thanks in advance for any help,


Patrick




___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] ANN: Haskell BLAS bindings version 0.7

[Patrick Perry]

Thanks for the suggestion!

Done and uploaded to hackage.


Patrick

On Jan 10, 2009, at 9:07 AM, Eugene Kirpichov wrote:


Why don't you put them into a separate non-hidden Unsafe module and
provide documentation for it? Users who don't care simply won't look
there, whereas users who do care (for whom you are providing access)
will still have a possibility to do so.

2009/1/10 Patrick Perry :
The reason I want to expose modules but hide the documentation is  
because
there are a lot of "unsafeX" functions I want to provide access to,  
but
which 99% of users don't care about.  The array library does the  
same thing.



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe



___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] ANN: Haskell BLAS bindings version 0.7

[Patrick Perry]

The reason I want to expose modules but hide the documentation is  
because there are a lot of "unsafeX" functions I want to provide  
access to, but which 99% of users don't care about.  The array library  
does the same thing.



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] ANN: Haskell BLAS bindings version 0.7

[Patrick Perry]

Hi David,

The problem is with Hackage, not with haddock.


Patrick

On Jan 10, 2009, at 3:38 AM, David Waern wrote:


2009/1/10 Patrick Perry :

Here's the haddock documentation; I'm not sure if Hackage honors {-#
OPTIONS_HADDOCK hide #-} when it displays the exposed modules:

http://quantile95.com/blas/


It should, so If it doesn't then please tell us about it. We have a  
trac at:


 http://trac.haskell.org/haddock/

David


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: Haskell BLAS bindings version 0.7

[Patrick Perry]

Here's the haddock documentation; I'm not sure if Hackage honors {-#  
OPTIONS_HADDOCK hide #-} when it displays the exposed modules:


http://quantile95.com/blas/


Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: Haskell BLAS bindings version 0.7

[Patrick Perry]

New version!

The blas package is a set of bindings to the BLAS (Basic Linear Algebra
Subprograms) library.

On Hackage:

http://hackage.haskell.org/cgi-bin/hackage-scripts/package/blas-0.7

What's new:

* Get rid of most functional dependencies in favor of type families.
  There is one remaining functional dependency that cannot be gotten
  rid of until GHC implements equality constraints in superclass
  contexts.
* Put the right superclass constraints in  ReadMatrix/ReadBanded.
* Fix freeze/thaw functions.  They used to be awkward to use.
* Fix a bug in getting a row view of a banded matrix.
* Make sure no NaNs/Infinities/Denormals are used when testing.
* Documentation.
* Export functions for writing QuickCheck tests.
* Remove Perm and Diag.  These will be in the LAPACK bindings.
* Get rid of "UnsafeIOToM" type class.
* Lots of INLINE everywhere.  This will improve performance (?)
* Switch to autoconf for build system.

Requisite blog post:

http://quantile95.com/2009/01/09/new-haskell-blas-bindings/

Not much example code.  Here's some old stuff:

 http://github.com/patperry/blas/tree/master/examples/LU.hs

Thanks for listening.  Any feedback is always welcome.


Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: monte-carlo-0.2, gsl-random-0.2.3

[Patrick Perry]

I've released a new version of the monte-carlo packages for haskell.   
Here are the highlights for monte-carlo:


Changes in 0.2:

* More general type class, MonadMC, which allows all the functions to  
work

  in both MC and MCT monads.

* Functions to sample from discrete distributions.

* Functions to sample subsets

For a quick tutorial, see my blog post at 
http://quantile95.com/2008/12/31/monte-carlo-poker-odds/

Happy New Year, everyone!


Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] How to think about this? (profiling)

[Patrick Perry]

I agree with everyone else who has said that the better solution is  
Lemmih's.  It is simple, fast, and does not use much memory.


On the other hand, here is a more faithful implementation of what you  
were trying to write.  To use mutable arrays, you need to work in  
either the IO or the ST monad.  The advantage of ST is the function  
"runST", which is roughly equivalent to "unsafePerformIO", but is much  
safer.


I hope this helps,


Patrick

\begin{code}

import Control.Monad
import Control.Monad.ST
import Data.Array.ST
import System.Environment

-- | Create an array for storing the results of foo 0 through foo n.
-- Initialize foo 0, foo 1, and foo 2, but set all other values to
-- 'missing'.
newMemo :: Int -> ST s (STUArray s Int Int)
newMemo n = do
arr <- newArray (0,n) missing :: ST s (STUArray s Int Int)
writeArray arr 0 0
writeArray arr 1 1
writeArray arr 2 2
return arr

-- | Code to indicate a missing value in the memo array.
missing :: Int
missing = -1

-- | Compute the value of the function at @i@, using the memo-ed results
-- in the array.
fooWithMemo :: STUArray s Int Int -> Int -> ST s Int
fooWithMemo arr i = do
  x <- readArray arr i
  if x /= missing
  then return x
  else do
  r <- liftM3 (\ a b c -> a + b + c)
   (fooWithMemo arr $ i - 1)
   (fooWithMemo arr $ i - 2)
   (fooWithMemo arr $ i - 3)
  writeArray arr i r
  return r

foo :: Int -> Int
foo n = runST $ do
arr <- newMemo n
fooWithMemo arr n

main = do
  (n:_)  <- liftM (map read) getArgs
  print $ foo n

\end{code}

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] unsafeThawByteArray#

[Patrick Perry]

I've noticed that there is no longer an "unsafeThawByteArray#"  
function in GHC.Prim.

Looking through the darcs history, I found a patch with description:

[project @ 2000-03-13 12:11:43 by simonmar]
simonmar**2313121144
 Remove unsafeThawByteArray# primop (which was a no-op), and use
 unsafeCoerce# instead.
]

Is the following code safe?

data IntArray = IntArray   !Int (ByteArray#)
data STIntArray s = STIntArray !Int (MutableByteArray# s)

unsafeThaw :: IntArray -> ST s (STIntArray s)
unsafeThaw (IntArray n arr#) =
ST $ \s1# ->
let coerceArray :: State# s -> MutableByteArray# s
coerceArray _ = unsafeCoerce# arr#
marr# = coerceArray s1#
marr  = STIntArray n marr# in
(# s1#, marr #)

Thanks,


Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: permutation-0.2

[Patrick Perry]

Hi Everyone,

I've uploaded a new version of the permutation library to hackage.   
Here is the package description:


This library includes data types for storing permutations.  It
implements pure and impure types, the latter which can be modified
in-place.  The main utility of the library is converting between
the linear representation of a permutation to a sequence of swaps.
This allows, for instance, applying a permutation or its inverse
to an array with O(1) memory use.
.
Much of the interface for the library is based on the permutation
functions in the GNU Scientific Library (GSL).

Here's the url:
http://hackage.haskell.org/cgi-bin/hackage-scripts/package/permutation


Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] Re: Statically dimension-checked hmatrix

[Patrick Perry]

What is the situation regarding statically dimension-checked linear
algebra libraries? It seems that Frederik Eaton was working on one in
2006 (see the paper "Statically typed linear algebra in Haskell"), and
he produced the Vectro library from this, based on GSLHaskell.

Are there any more recent efforts into this, particularly using the
new TFs? If not, I might have a go at it, as a thin wrapper for
hmatrix.


The BLAS bindings I wrote use phantom types to make sure the  
dimensions are consistent.  See, for example the functions to get row  
and column views of a matrix:


http://hackage.haskell.org/packages/archive/blas/0.6/doc/html/Data-Matrix-Dense-Class.html#v%3ArowViews

Also, the multiplication routines:

http://hackage.haskell.org/packages/archive/blas/0.6/doc/html/BLAS-Matrix-Immutable.html

I've found that phantom types are a good trade off: they provide  
enough type safety to catch common mistakes without being too much of  
a hassle to use.  Another benefit (which I didn't anticipate) is that  
they help a *lot* with documentation.


I'm skeptical of value that comes from stronger typing.  At some point  
the types become too much of a hassle to be worth using.  Even phantom  
types get in the way sometimes and require use of either GADTs,  
unsafeCoerce, or both, which is annoying.


I'd be interested to hear about any developments you make in this area.


Patrick

p.s. If you're interested in trying it out, the version of blas on  
hackage doesn't compile with ghc-6.10.1, but the version in the darcs  
repository at http://stat.stanford.edu/~patperry/code/blas does.


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: blas version 0.6

[Patrick Perry]

New version of BLAS bindings out.  Now with support for the ST monad!   
This breaks backwards compatibility, unfortunately.


More info (and some sample code) here: 
http://quantile95.com/2008/10/31/ann-blas-bindings-for-haskell-version-06/

Also, if you want to help, please let me know.  There is plenty of  
work to do.



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] ANN: gsl-random 0.1 and monte-carlo-0.1

[Patrick Perry]

Thanks for the heads up, Don.  I fixed the file in version 0.1.1.

I also changed the default CBLAS to the one that comes with the GSL.   
Anyone who cares at all about performance will want to configure the  
package to use the best CBLAS on their system.  If this is ATLAS, make  
sure you pass the "-fatlas" argument when you configure the package.



Patrick

On Aug 28, 2008, at 9:56 AM, Don Stewart wrote:


patperry:

Hi everyone,

I've started on bindings for the random number generators and random
distributions provided by the gsl.  The package is available here:
http://hackage.haskell.org/cgi-bin/hackage-scripts/package/gsl-random

I've also written a monad and transformer for doing monte carlo
computations that uses gsl-random internally.  Here is that package:
http://hackage.haskell.org/cgi-bin/hackage-scripts/package/monte- 
carlo


For a quick tutorial in the latter package, see my blog:
http://quantile95.com/2008/08/27/a-monte-carlo-monad-for-haskell/

There is also a more complicated example in the "examples" directory.
Currently, only normal, uniform, and poisson random variables are
supported.  I have no plans to add anything else unless I need it,  
but
I will happily accept patches if someone else is willing to do the  
work.


One thing you may need to watch out for is that gsl-random needs to
link with cblas.  If your cblas is not called "cblas", you may have  
to

edit "gsl-random.cabal" to get things to work correctly.  To use the
cblas that comes with the gsl, change "cblas" to "gslcblas".  To use
ATLAS, change "cblas" to "cblas atlas".


Arch packages available here,

   http://aur.archlinux.org/packages.php?ID=19416
   http://aur.archlinux.org/packages.php?ID=19417

Note I had to patch gsl-random's .cabal file to not explicitly look  
for

extra libs in /opt/local/lib. The exact path (if any) should be cabal
(or ./configure's) job, I think.

-- Don




___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: gsl-random 0.1 and monte-carlo-0.1

[Patrick Perry]

Hi everyone,

I've started on bindings for the random number generators and random  
distributions provided by the gsl.  The package is available here:

http://hackage.haskell.org/cgi-bin/hackage-scripts/package/gsl-random

I've also written a monad and transformer for doing monte carlo  
computations that uses gsl-random internally.  Here is that package:

http://hackage.haskell.org/cgi-bin/hackage-scripts/package/monte-carlo

For a quick tutorial in the latter package, see my blog:
http://quantile95.com/2008/08/27/a-monte-carlo-monad-for-haskell/

There is also a more complicated example in the "examples" directory.   
Currently, only normal, uniform, and poisson random variables are  
supported.  I have no plans to add anything else unless I need it, but  
I will happily accept patches if someone else is willing to do the work.


One thing you may need to watch out for is that gsl-random needs to  
link with cblas.  If your cblas is not called "cblas", you may have to  
edit "gsl-random.cabal" to get things to work correctly.  To use the  
cblas that comes with the gsl, change "cblas" to "gslcblas".  To use  
ATLAS, change "cblas" to "cblas atlas".



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] Bug with QuickCheck 1.1 and GHC 6.8.2

[Patrick Perry]

'split' has type

split :: StdGen -> (StdGen, StdGen)

It takes a generator and produces two independent generators.  The old  
version of variant calls 'split' n times to get a new generator, and  
it only calls "split" on the second result.  The new version calls  
split (log n) times.  It gets a unique generator for each value of n  
by using both results from split.


For example,

rands 2 r0 = let (r1, r2) = split r0
(r3, r4) = split r1
   in r4

rands 3 r0 = let (r1,r2) = split r0
(r3',r4') = split r2
   in r4'


Patrick

On Aug 14, 2008, at 4:57 AM, Johan Tibell wrote:

On Thu, Aug 14, 2008 at 12:15 PM, Patrick Perry  
<[EMAIL PROTECTED]> wrote:

Actually, a much better solution is:

variant :: Int -> Gen a -> Gen a
variant v (Gen m) = Gen (\n r -> m n (rands r v))
where
rands r0 0 = r0
rands r0 n = let (r1,r2) = split r0
 (n',s)  = n `quotRem` 2
 in case s of
  0 -> rands r1 n'
  _ -> rands r2 n'


It's not really clear to me how this works. Could you please explain?

Cheers,

Johan


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] Bug with QuickCheck 1.1 and GHC 6.8.2

[Patrick Perry]

Actually, a much better solution is:

variant :: Int -> Gen a -> Gen a
variant v (Gen m) = Gen (\n r -> m n (rands r v))
 where
  rands r0 0 = r0
  rands r0 n = let (r1,r2) = split r0
   (n',s)  = n `quotRem` 2
   in case s of
0 -> rands r1 n'
_ -> rands r2 n'


Patrick


On Aug 14, 2008, at 2:17 AM, Johan Tibell wrote:

On Thu, Aug 14, 2008 at 1:58 AM, Patrick Perry  
<[EMAIL PROTECTED]> wrote:

variant :: Int -> Gen a -> Gen a
variant v (Gen m) = Gen (\n r -> m n (rands r !! v'))
where
v' = abs (v+1) `mod` 1
rands r0 = r1 : rands r2 where (r1, r2) = split r0


Note that if you have a uniformly distributed random value in [0, n)
and take its value mod k you don't end up with another random
uniformly distributed value unless n `mod` k == 0. Consider n = 3 and
k = 2. What you can do is to throw away all random numbers larger than
n - (n `mod` k) and just generate a new number. This will terminate
with a high probability if n >> k.

Cheers,

Johan


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] Bug with QuickCheck 1.1 and GHC 6.8.2

[Patrick Perry]

The bug is in the "variant" function in QuickCheck.  I replaced

variant :: Int -> Gen a -> Gen a
variant v (Gen m) = Gen (\n r -> m n (rands r !! (v+1))
 where
  rands r0 = r1 : rands r2 where (r1, r2) = split r0

with

variant :: Int -> Gen a -> Gen a
variant v (Gen m) = Gen (\n r -> m n (rands r !! v'))
 where
  v' = abs (v+1) `mod` 1
  rands r0 = r1 : rands r2 where (r1, r2) = split r0

and now everything works fine.  "1" seems like a reasonable value  
here, but one could hard-code a lower or higher value as well.


Can someone make sure this gets fixed in the next version of  
QuickCheck?  I'm not sure who the maintainer is.



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] Bug with QuickCheck 1.1 and GHC 6.8.2

[Patrick Perry]

I'm running into problems with generating an arbitrary function of  
type Double -> Double.  Specifically, the following code:


{-# LANGUAGE PatternSignatures #-}
import Test.QuickCheck
import Text.Show.Functions

prop_ok (f :: Double -> Double) =
f (-5.5) `seq` True

prop_bug (f :: Double -> Double) =
f (-5.6) `seq` True

main = do
putStr "prop_ok:\t"  >> quickCheck prop_ok
putStr "prop_bug:\t" >> quickCheck prop_bug


On an Intel Core 2 Duo running Mac OS 10.5.4 with GHC 6.8.2 the output  
I get is


prop_ok:OK, passed 100 tests.
prop_bug:   Test: Prelude.(!!): negative index

On Intel Xeon running RedHat with GHC 6.8.2 both tests hang.

Has anyone seen this before?  Is it fixed in QuickCheck2?

Thanks,


Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: BLAS bindings for haskell, version 0.5

[Patrick Perry]

Hey everyone,

I've put together a new release of the Haskell BLAS bindings, now  
available on hackage.


Here are the new features:

* Add Banded matrix data type, as well as Tri Banded and Herm Banded.
* Add support for trapezoidal dense matrices (Tri Matrix (m,n) e, where
  m is not the same as n).  Note that trapezoidal banded matrices are
  *NOT* supported.
* Add Diag matrix data type for diagonal matrices.
* Add Perm matrix data type, for permutation matrices.
* Enhance the RMatrix and RSolve type classes with an API that allows
  specifying where to store the result of a computation.
* Enhance the IMatrix, RMatrix, ISolve, and RSolve type classes to add
  "scale and multiply" operations.
* Remove the scale parameter for Tri and Herm matrix data types.
* Flatten the data types for DVector and DMatrix.
* Some inlining and unpacking performance improvements.

As far as what to expect in version 0.6, I plan to add support for  
operations
in the ST monad.  This is going to require a pretty big code  
reorganization and
will cause quite a few API breakages, but I think it's worth the  
pain.  The

next release will also come with a tutorial and examples.  After the big
code reorganization, I'll get started on LAPACK bindings.

Please let me know if you are using the library.  I'm really  
interested in
what people like and don't like.  If you think that some functionality  
is
missing, let me know.  If you think the API is awkward in certain  
places,

let me know that, too.


Patrick
Cross-posted to:
http://quantile95.com/2008/08/11/ann-blas-bindings-for-haskell-version-05/


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] Re: blas bindings, why are they so much slower the C?

[Patrick Perry]

Yeah, I think that's where most of the performance gains came from.  I  
also added a re-write rule for unsafeGet dot (since it doesn't matter  
if the arguments are conjugated or not if the vectors are real)  that  
shaved off about a tenth of a second.



Patrick

On Jul 24, 2008, at 4:26 PM, Don Stewart wrote:


patperry:

Last month Anatoly Yakovenko published some disturbing numbers about
the Haskell BLAS bindings I wrote being significantly slower than
using plain C.  I wanted to let everyone know that I've closed the
performance gap, and now for doing ten million dot products, the
overhead for using Haskell instead of C is about 0.6 seconds on my
machine, regardless of the size of the vectors.  The next version  
will

incorporate the changes.  If you can't wait for a formal release, the
darcs repository is at http://www-stat.stanford.edu/~patperry/code/blas/

Anyone interested in more details can check out my blog:
http://quantile95.com/2008/07/24/addressing-haskell-blas-performance-issues/

Thanks everyone for the input on this (especially Anatoly).  If any
else finds any performance discrepancies, please let me know and I
will do whatever I can to fix them.



Great work, Patrick!

So if I read correctly, the main change was to flatten the
representation (and thus in loops the vector's structure will be
unpacked and kept in registers, which isn't possible for sum types).

-- Don


___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] Re: blas bindings, why are they so much slower the C?

[Patrick Perry]

Last month Anatoly Yakovenko published some disturbing numbers about  
the Haskell BLAS bindings I wrote being significantly slower than   
using plain C.  I wanted to let everyone know that I've closed the  
performance gap, and now for doing ten million dot products, the  
overhead for using Haskell instead of C is about 0.6 seconds on my  
machine, regardless of the size of the vectors.  The next version will  
incorporate the changes.  If you can't wait for a formal release, the  
darcs repository is at http://www-stat.stanford.edu/~patperry/code/blas/


Anyone interested in more details can check out my blog:
http://quantile95.com/2008/07/24/addressing-haskell-blas-performance-issues/

Thanks everyone for the input on this (especially Anatoly).  If any  
else finds any performance discrepancies, please let me know and I  
will do whatever I can to fix them.



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] BLAS Solve Example

[Patrick Perry]

Sorry Darrin, the BLAS library only includes matrix multiplication and  
solving triangular systems.  To solve a general system, you would need  
to use LAPACK, but there aren't any bindings for that library yet.  I  
would suggest you take a look at the hmatrix package, which includes a  
lot more linear algebra than blas does..



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: ieee-0.2

[Patrick Perry]

ieee is a library that provides approximate comparison of floating  
point numbers based, NaN-aware minimum and maximum, and a type class  
for approximate comparisons.


This version fixes a bug in the comparison implementation for Maybe.

http://hackage.haskell.org/cgi-bin/hackage-scripts/package/ieee
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] ANN: blas-0.4.1

[Patrick Perry]

There were a few hiccups in the release of 0.4, so this version is to  
make installation go a little smoother for people.  There are now  
installation instructions (the INSTALL file) and configuration  
settings for some common CBLAS (ATLAS, MKL, vecLib, GSL).  The default  
is now to assume that ATLAS is installed, so hopefully most people  
won't have to fuss too much with this.


Alberto Ruiz reported two small issues with precision in the unit  
tests, and those have now been fixed.



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] Re: Patrick Perry's BLAS package

[Patrick Perry]

a function like foo :: (BLAS1 e) => Matrix (m,n) e ->
Matrix (n,k) e -> Int -> Vector m e foo a b i = let x =
row b i in a <*> x will not type-check. (”<*>” is the
function to multiply a matrix by a vector. Everything is
ok if you replace “row” by “col”.) This feature has caught
a few bugs in my code.


If I understand this correctly, the compiler can catch
dimension mismatches so that using `col' will result in a
compilation error when m and k are different, is it so?


Yes, the compiler infers the type of expression to be

Matrix(m,n) e -> Matrix (l,m) e -> Int -> Vector n e

which is different from the declared type.  The error will be  
something like "Expected type Vector m e, but got type Vector a1 e  
instead".



I haven't look through the code, yet.  But it looks like
there are certain levels of similarities between blas and
hmatrix.  Is it possible for these two libraries to
cooperate well with each other?  (I didn't look at HBlas, so
can't say much about that.)


This is more work than you might think.  The data structures are  
different, many of the function names are different, and the module  
namespaces overlap.  I wouldn't recommend any switch from hmatrix to  
blas right now unless you really need mutable types-- hmatrix has a  
lot more linear algebra functionality.



Apart from some warnings, the library compiles fine in my
system.  But there is a minor issue about the library it
links against when `./Setup test'.  I need to use `-lcblas'
instead of `-lblas' to get it to link to correct libraries.
I don't know other people's system.  But in my system,
Gentoo Linux, I use blas library provided by atlas, and
libblas.so is a fortran library and libcblas.so is for C.


I don't know of a good way to get around this.  It seems like every  
system has a different convention for the location and name of the  
cblas libraries.  So, everyone has to edit the "extra-libraries" and  
the "extra-lib-dirs" field in the blas.cabal file.  Does anyone know  
of a better way of doing this?



Patrick

___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

[Haskell-cafe] Re: Patrick Perry's BLAS package

[Patrick Perry]

Wow, thanks for noticing, Alberto!  For anyone interested, I put up a  
formal announcement describing the bindings a little bit here:


http://quantile95.com/

I just registered the domain yesterday, so it may take a few days to  
resolve the DNS magic.  Here's the text of the announcement:


I’m really happy that people seem to be interested in the library.  
Alberto, in particular, is the primary author of hmatrix, another  
haskell linear algebra library (which I stole a few ideas from), so if  
he endorses it, that means a lot to me.


So, Yet Another Linear Algebra Library? I’ve already mentioned  
hmatrix. There’s also another one called HBlas. Why would anyone want  
a third? Here are my reasons:


	* Support for both immutable and mutable types. Haskell tries to make  
you use immutable types as much as possible, and indeed there is a  
very good reason for this, but sometimes you have a 100MB matrix, and  
it just isn’t very practical to make a copy of it every time you  
modify it. hmatrix only supports immutable types, and HBlas only  
supports mutable ones. I wanted both.
	* Access control via phantom types. When you have immutable and  
mutable types, it’s very annoying to have separate functions for each  
type. Do I want to have to call “numCols” for immutable matrices and  
“getNumCols” for mutable ones, even though both functions are pure,  
and both do exactly the same thing? No. If I want to add an immutable  
matrix to a mutable one, to I want to first call “unsafeThaw” on the  
immutable one to cast it to be mutable? No. With the phantom type  
trick, you can get around this insanity. Jane Street Capital has a  
very good description of how this works.
	* Phantom types for matrix and vector shapes. This is a trick I  
learned from darcs. It means that the compiler can catch many  
dimension-mismatch mistakes. So, for instance, a function like

foo :: (BLAS1 e) =>
Matrix (m,n) e -> Matrix (n,k) e -> Int -> Vector m e
foo a b i = let x = row b i in a <*> x
will not type-check. (”<*>” is the function to multiply a matrix by a  
vector. Everything is ok if you replace “row” by “col”.) This feature  
has caught a few bugs in my code.


	* Taking the conjugate transpose (”herm”) of a matrix is an O(1)  
operation. This is similar to hmatrix, where taking the transpose is  
O(1). As BLAS and LAPACK (mostly) support this, it makes no sense to  
copy a matrix just to work with the conjugate transpose. Why conjugate  
transpose instead of just transpose? Because the former is a far more  
common operation. This is why the “‘” operator in MATLAB is conjugate  
transpose. The drawback for this feature is that BLAS and LAPACK do  
not support it everywhere. In particular, QR decomposition with  
pivoting is going to be a huge pain in the ass to support for herm-ed  
matrices.
	* Support for triangular and hermitian views of matrices. This is a  
feature of BLAS that no one seems to support (not even MATLAB). In  
addition to the “Matrix” type, there are “Tri Matrix” and “Herm  
Matrix” types that only refer to the upper- or lower-triangular part  
of the matrix.
Hopefully the features above are compelling enough to make people want  
to use the library. These bindings have been a lot of work. For me to  
come up with the feature list above, I’ve already gone through a few  
iterations of dramatic re-writes (hence the version number). Of  
course, I always welcome suggestions for how to make it better.


What’s next? In the immediate future, I plan to add banded matrices.  
I’ve already written a good chunk of code for this, but it isn’t very  
well tested, so I decided to leave it out of the release. I’m also  
going to add permutation matrices. I don’t have plans to add support  
for packed triangular matrices, but if someone else wanted to do that,  
I would be happy to include it. The same goes for symmetric complex  
matrices.


LAPACK support is on the horizon, but that may take awhile. Also, I  
probably won’t do more than SVD, QR, and Cholesky, since those are all  
I need. Expect a preliminary announcement by the end of the summer.


This work would not have been possible without looking at the other  
excellent linear algebra libraries out there. In particular the GNU  
Scientific Library was the basis for much of the design. I also drew  
inspiration from hmatrix and the haskell array libraries.


Please let me know if you have any success in using the library, and  
if you have any suggestions for how to make it better.



Patrick

On Jun 6, 2008, at 5:36 AM, Alberto Ruiz wrote:


Hello all,

I have just noticed that yesterday this fantastic package has been  
uploaded to hackage:


http://hackage.haskell.org/cgi-bin/hackage-scripts/package/blas-0.4

We finally have a high quality library for numeric linear algebra.  
This is very good news for the Haskell community.


Patrick, many thanks for your excellent work. Do you have similar  
plans for LAPACK?


Alberto


_