http://arstechnica.com/science/2014/05/scientific-computings...
I have written real time 3D reconstruction and SLAM algorithms in Haskell, so I have had to face this head on.
mmMult :: Monad m
=> Array U DIM2 Double
-> Array U DIM2 Double
-> m (Array U DIM2 Double)
mmMult a b = sumP (Repa.zipWith (*) aRepl bRepl)
where
t = transpose2D b
aRepl = extend (Z :.All :.colsB :.All) a
bRepl = extend (Z :.rowsA :.All :.All) t
(Z :.colsA :.rowsA) = extent a
(Z :.colsB :.rowsB) = extent b
http://www.haskell.org/haskellwiki/Numeric_Haskell:_A_Repa_T... readIOArray m (i,j)
But that is misleading because there are two kinds of arrays in Haskell, and the kind (namely, mutable) I used in my literal answer is the less common kind. In fact, it is so uncommon and so contrary to the spirit of Haskell that the language implementations do not take the trouble to make it fast. Last time I played with the Glasgow Haskell Compiler, for example, a related operation (readIORef) was over 100 times slower than its counterpart in C.Parenthetically, because you need to use "a monad" to perform any sort of mutable operation in Haskell, the translation of
m[i,j] := m[j,i]
into Haskell is not writeIOArray m (i,j) (readIOArray m (j,i))
like someone unfamiliar with monadic style might think, but rather readIOArray m (j,i) >>= writeIOArray m (i,j)
or equivalently do
element <- readIOArray m (j,i)
writeIOArray m (i,j) element
But that is a side issue. The central fact is that even though it is technically possible, it is considered bad style in Haskell to access or update individual array elements.> with Haskell the code would look more like the actual equations you're working from
so at least in matrix algorithms, and finite difference methods for partial differential equations, the textbooks usually explain the algorithms using a notation that refers to matrix elements. So in these cases, an imperative implementation would be close to the math textbook representation, and the Haskell code would not.
And another reason why Haskell will not fair well. A lot of numerical mathematics is linear algebra which involves, almost exclusively, updating array elements.
I'm sure Haskell is good at some things but for many of the types of problems that numerical analysts face Haskell's paradigm just doesn't make sense. Especially when we have languages like Fortran/Matlab which were explicitly designed for computational mathematics.
Look more like the actual equations? Can you give me a concrete example of Haskell code that solves a simple PDE?
As languages go, Haskell out paces (not just in elegance but pragmatism too) most of the mainstream imperative languages.
The tooling for scientific computing does need to catch up though.
Remember that scientists don't care about the quality of their code. In many cases, once the code is complete, it will only run once and then be thrown away. There is no maintenance. Things like technical debt don't exist.
And pragmatism? Haskell doesn't even have loops. Loops with mutable variables inside are a very intuitive way of thinking about things and often a better fit to a problem than recursion. Sure, you can use tail recursion (which is just ugly looping) instead, but you're adding unwanted cognitive overhead.
What I would like is a functional language with enclosed, sealed off loops. Syntactic sugar for tail recursion, basically, but Haskell isn't built with programmers in mind; it's built for CS researchers. Computer language researchers, more specifically.
In any dependently typed language, absolutely yes.
It's very simple to build a function which doesn't have to understand the possible dimension conflicts and lift it to work on this new type, returning an either (or a maybe, if there's only one failure mode) in place of a definite value.
It's also very simple to propagate such errors forward, so they'll short circuit a computation when you have non-matching matrices used in a calculation that's multiple steps.
In Haskell, I don't have to remember to write special functions which guard against this: I write functions that operate on the matrices and add the guarding at the very end. I can ensure that all my calls use the guarding functions, because they have a different type signature.
Trying to do this same thing in Python require that I remember to always use the guarded calls, and doesn't have as clean of an interface to create the guarded functions from standard functions.
Yes they can, it's known as 'type providers'. http://blogs.msdn.com/b/dsyme/archive/2013/01/30/twelve-type...
I can do it in C# for small matrices, by including the dimensions as type bindings (http://bling.codeplex.com), but that is not practical for larger matrices, nor would it work for matrices loaded from IO.
To the downvoters: are you claiming Haskell is dependently typed? If so, why isn't it used in http://hackage.haskell.org/package/matrix-0.2.1/docs/Data-Ma...? Or just one matrix data type that checks dimension lengths via the type system?
Repa for instance does provide type errors based on dimensionality of matrices. This is a package on Hackage.
See http://www.haskell.org/haskellwiki/Numeric_Haskell:_A_Repa_T...
Looking at the linked page, extent isn't a part of a matrix's type signature, so it would be checked dynamically, correct?
This is exactly how Repa works, it uses a Peano encoding of the extent of dimensions to make invalid array operations inexpressible.
Yes. In fact this is one of the earliest examples [0] of what you can do with the ever more numerous extensions to Haskell's type system.
There are BLAS bindings [1] and several native linear-algebra libraries in Hackage [2] that enlist the type-checker to enforce shape constraints at compile time.
0. See McBride (2001), "Faking It": http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.22.2...
1. http://hackage.haskell.org/package/blas
2. Repa, for instance: http://hackage.haskell.org/package/repa, http://hackage.haskell.org/package/repa-algorithms
It seems so
http://stackoverflow.com/questions/8332392/type-safe-matrix-...
The code is usually used once, but it’s modified and “improved” for the next paper. So you add a few new routines here and a few files to read or write there, ... And one day you have a 95 pages fortran 77 program (it’s does’t even follow the fortran 90 standard).
If you need any more input from me to convince you to learn something new, then it's not worth my time or yours.
Like I said the tooling isn't as mature but the language enables me to produce results faster with fewer bugs that is much easier to maintain down the road.
The only numerical computing I've ever really programmed for was manipulation of time series data and some standard statistical methods.
I did it at first in Python with Pandas and numpy. Then moved to Haskell. The language was better but the tooling I had to shore up in places.
Take that if it's useful to you :)
Not when you're doing mathematical research. You can write x86 machine code and the time spent programming is still completely dwarfed by the time spent doing math. Mathematics is full of renowned professors who write F77 by hunting and pecking with one finger. There's no reason for them to use anything "more productive" because they can already code far, far faster than they can think of what to code.
There are extreme outliers for whom this doesn't hold, but they are few and far between.
Productive in what sense and at what cost (in time, especially, but generally)? What benefit is there to learning Haskell when the existing "tooling" in FORTRAN77, C, or C++ is quite adequate for the purpose of research? I mean for the specific purpose of applying it in a research context; not the general sense (i.e. this isn't meant to question the worth learning something new in general).
Some people also use Fortran 95/2003/2008.
That's the funny thing about researchers, some seem extremely obstinate (I suppose in many ways that's a required attribute otherwise they might give up). I wonder if this is why we still have Fortran code around...
The state of the art in numerical computing really isn't to the point we can work from equations because things like sparsity and which solver you use (Algebraic Multi Grid? Conjugate Gradient Descent? ...) matter a lot and the established algorithms and implementations are really fast and well tested.