Re: [Haskell-cafe] Purely funcional LU decomposition
---BeginMessage--- Rafael Gustavo da Cunha Pereira Pinto wrote: interesting to look at real matrix code that people have written and think about what would be needed in a library to make it easier to write. -- Dan What I miss most is a data structure with O(1) (amortized) direct access. STArray and STUArray? ---End Message--- ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
---BeginMessage--- Paulo Tanimoto wrote: Pretty cool, thanks for releasing this into the wild. I remember looking into this about a year ago. By the way, have you seen Matt's DSP library? http://haskelldsp.sourceforge.net/ He's got LU and others in there, if my memory serves me. The last release seems to be 2003, so it might be worth emailing him to see what happened and if he has plans for the future. I have cleaned it up and uploaded it to Hackage some time ago: http://hackage.haskell.org/cgi-bin/hackage-scripts/package/dsp You may also be interested in http://hackage.haskell.org/cgi-bin/hackage-scripts/package/numeric-quest http://ofb.net/~frederik/vectro/ http://haskell.org/haskellwiki/Linear_algebra ---End Message--- ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
2009/2/5 Henning Thielemann lemm...@henning-thielemann.de -- Forwarded message -- From: Henning Thielemann thunderb...@henning-thielemann.de To: Rafael Gustavo da Cunha Pereira Pinto rafaelgcpp.li...@gmail.com Date: Thu, 05 Feb 2009 14:43:13 +0100 Subject: Re: [Haskell-cafe] Purely funcional LU decomposition Rafael Gustavo da Cunha Pereira Pinto wrote: interesting to look at real matrix code that people have written and think about what would be needed in a library to make it easier to write. -- Dan What I miss most is a data structure with O(1) (amortized) direct access. STArray and STUArray? Are they? :-) -- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
Rafael Gustavo da Cunha Pereira Pinto wrote: What I miss most is a data structure with O(1) (amortized) direct access. Finger trees get close, O(log(min(i,n-i))): http://hackage.haskell.org/packages/archive/containers/latest/doc/html/Data-Sequence.html (And Theta(1) amortized for all dequeue operations to boot.) One of the changes I thought today was to remove the ++ operator and create a list of lists that I would concat in the last call. Don't do it manually, let an abstract type do the optimization for you. This idea[1] is often known as difference lists, which the estimable Don Stewart has already packaged up for us: http://hackage.haskell.org/cgi-bin/hackage-scripts/package/dlist :) [1] Actually, not quite that idea. The idea instead is to represent a list as a function [a]-[a] and then compose these functions (rather than concatenating the strings). At the very end once you're done, pass in [] and get back the whole list in O(n) time, rather than the O(n^2) that commonly arises from using (++). -- Live well, ~wren ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
On Wed, Feb 4, 2009 at 04:42, Thomas Davie tom.da...@gmail.com wrote: Shinyness indeed – a quick note though, as ghc doesn't support dynamic linking of Haskell code, the above is equivalent to the GPL. I always use LGPL. Anyway, I will keep it that way, as I still have hopes on dynamic linking Haskell becoming a reality Would be lovely if you packaged this up and stuck it on Hackage :) I may do that when I am finished with all the other stuff. I am building a (very primitive) circuit simulator using Haskell. It is a toy I am doing just to learn Haskell (until now I learned how to use Arrays and Parsec)... Bob -- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
Matt's code is pretty comprehensive. His LU implementation is much cleaner, essentially because he used the ijk format, while I used the kij. I'll take a look and e-mail him eventually. Thanks! On Tue, Feb 3, 2009 at 23:26, Paulo Tanimoto tanim...@arizona.edu wrote: Hi Rafael, 2009/2/3 Rafael Gustavo da Cunha Pereira Pinto rafaelgcpp.li...@gmail.com : Hello folks After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell... At first I thought it would be necessary to use a mutable or monadic version of Array, but then I figured out it is a purely interactive process. I am releasing this code fragment as LGPL. Pretty cool, thanks for releasing this into the wild. I remember looking into this about a year ago. By the way, have you seen Matt's DSP library? http://haskelldsp.sourceforge.net/ He's got LU and others in there, if my memory serves me. The last release seems to be 2003, so it might be worth emailing him to see what happened and if he has plans for the future. Regards, Paulo ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe -- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
2009/2/3 Rafael Gustavo da Cunha Pereira Pinto rafaelgcpp.li...@gmail.com: After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell... import Data.Array.IArray ... e_an i j=a!(i,j)-(lik i)*a!(k,j) There are three different representations of arrays in this code: arrays, lists and a functional one where f represents an array with elements (f i j). Looks like a candidate for a stream fusion type thing so the user only writes code using one representation and some RULES switch back and forth between representations behind the scenes. -- Dan ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
On Wed, Feb 4, 2009 at 17:09, Dan Piponi dpip...@gmail.com wrote: 2009/2/3 Rafael Gustavo da Cunha Pereira Pinto rafaelgcpp.li...@gmail.com : After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell... import Data.Array.IArray ... e_an i j=a!(i,j)-(lik i)*a!(k,j) There are three different representations of arrays in this code: arrays, lists and a functional one where f represents an array with elements (f i j). Looks like a candidate for a stream fusion type thing so the user only writes code using one representation and some RULES switch back and forth between representations behind the scenes. -- Dan Those different representations are derived from my (very) low Haskell handicap! :-D 1) I used lists for the intermediate L and U matrices so I wouldn't have the overhead of calling accumArray and assocs everywhere 2) I used the function that returns an array element just to shorten the line length, and to make sure I was typing it right! When I started this, I thought I would end with a couple of folds, but then I turned to the recursive version, so I could implement pivoting on matrix a' more easily. List comprehension also derived from this. I started thinking of a fold for L matrix and then it hit me that it was a list comprehension... Matt's DSP library has one of the most elegant implementations I saw, but it is harder to implement pivoting, which is really important for my application. lu :: Array (Int,Int) Double -- ^ A - Array (Int,Int) Double -- ^ LU(A) lu a = a' where a' = array bnds [ ((i,j), luij i j) | (i,j) - range bnds ] luij i j | ij = (a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k - [1 ..(j-1)] ]) / a'!(j,j) | i=j = a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k - [1 ..(i-1)] ] bnds = bounds a -- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
On Wed, Feb 4, 2009 at 12:57 PM, Rafael Gustavo da Cunha Pereira Pinto rafaelgcpp.li...@gmail.com wrote: Those different representations are derived from my (very) low Haskell handicap! :-D BTW That wasn't intended in any way as a criticism. I think it's interesting to look at real matrix code that people have written and think about what would be needed in a library to make it easier to write. -- Dan ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
On Wed, Feb 4, 2009 at 19:14, Dan Piponi dpip...@gmail.com wrote: On Wed, Feb 4, 2009 at 12:57 PM, Rafael Gustavo da Cunha Pereira Pinto rafaelgcpp.li...@gmail.com wrote: Those different representations are derived from my (very) low Haskell handicap! :-D BTW That wasn't intended in any way as a criticism. I think it's Not taken! interesting to look at real matrix code that people have written and think about what would be needed in a library to make it easier to write. -- Dan What I miss most is a data structure with O(1) (amortized) direct access. One of the changes I thought today was to remove the ++ operator and create a list of lists that I would concat in the last call. -- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Purely funcional LU decomposition
Hello folks
After a discussion on whether is possible to compile hmatrix in
Windows, I decided to go crazy and do a LU decomposition entirely in
Haskell...
At first I thought it would be necessary to use a mutable or
monadic version of Array, but then I figured out it is a purely interactive
process.
I am releasing this code fragment as LGPL.
{-
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published
by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see http://www.gnu.org/licenses/.
-}
import Data.Array.IArray
type Dim=(Int,Int)
lu::Array Dim Double - (Array Dim Double,Array Dim Double)
lu a = (aa l,aa u)
where
(l,u)=lu' a [] []
aa = accumArray (+) 0 (bounds a)
lu'::(Floating e) = Array Dim e - [(Dim,e)] - [(Dim,e)] -
([(Dim,e)],[(Dim,e)])
lu' a l u=if (ui==li) then (
((ui,uj),1.0):l,((ui,uj),a!(ui,uj)):u) else (lu' an (l++ln) (u++un))
where
k=li
((li,lj),(ui,uj))=bounds a
lik i=(a!(i,k)/a!(k,k))
un=[((k,j),a!(k,j)) | j-[lj..uj]]
ln=((lj,lj),1.0):[((i,k),lik i) | i-[li+1..ui]]
an=array ((li+1,lj+1),(ui,uj)) [((i,j),e_an i j) |
i-[li+1..ui],j-[lj+1..uj]]
e_an i j=a!(i,j)-(lik i)*a!(k,j)
Still needed:
1) Partial pivoting
2) Profiling... Lots!
3) Parallelize
--
Rafael Gustavo da Cunha Pereira Pinto
Electronic Engineer, MSc.
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
Hi Rafael, 2009/2/3 Rafael Gustavo da Cunha Pereira Pinto rafaelgcpp.li...@gmail.com: Hello folks After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell... At first I thought it would be necessary to use a mutable or monadic version of Array, but then I figured out it is a purely interactive process. I am releasing this code fragment as LGPL. Pretty cool, thanks for releasing this into the wild. I remember looking into this about a year ago. By the way, have you seen Matt's DSP library? http://haskelldsp.sourceforge.net/ He's got LU and others in there, if my memory serves me. The last release seems to be 2003, so it might be worth emailing him to see what happened and if he has plans for the future. Regards, Paulo ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Purely funcional LU decomposition
On 3 Feb 2009, at 22:37, Rafael Gustavo da Cunha Pereira Pinto wrote: Hello folks After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell... At first I thought it would be necessary to use a mutable or monadic version of Array, but then I figured out it is a purely interactive process. I am releasing this code fragment as LGPL. Shinyness indeed – a quick note though, as ghc doesn't support dynamic linking of Haskell code, the above is equivalent to the GPL. Would be lovely if you packaged this up and stuck it on Hackage :) Bob___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
