Re: Haskell in Scientific Computing?
Dave Tweed wrote: But there's a lot of problems, probably more in the hazy region between science engineering, where `numerically intensive' algorithms are developed which don't look anything like existing classical techniques. Here the issue is to generate CORRECT results REASONABLY QUICKLY, ie, the time has to be within a factor of 3-4 times of a C implementation but this slowdown is acceptable if you are more confident your infant algorithm is correctly implemented in your infant code. Let me give an example here that demonstrates a need for speed and one reason why some people still use FORTRAN. One of my pervious jobs was developing new receiver techniques for digital communication systems. Given a description of the communication channel we were working with (point-to-point mincrowave, cellular, etc) we could generate our theoretical error rate curves bases on system parameters. Our job was see how close we could get to the curve; ie. there isn't a concept of correct results, just degrees of success. We would develop a solution and then run a simulation of this solution against a set of system parameters (generally signal to noise ratio). For a 99% confidence in our results, we would run our simulation until we got 100 bit-errors. For theoretical errors rates to 1e-5, this would require about 1000 input points (I think I did my math right, it's been a while since I did this...). Needless to say, these would take a while. Even a 1.5 times slowdown would mean that we could run less simulations per day. For numerically intensive applications like this where running time is measure in hours, speed really does matter. A friend of mine uses FORTRAN just for this reason. Looking back, I would have really liked to have used Haskell. Infinite lists map very nicely into DSP and digital comm methods; I was able to boil down a failrly ugly C funtion to 3 lines of Haskell using infinite lists. -- Matt Donadio ([EMAIL PROTECTED]) | 43 Leopard Rd, Suite 102 Sr. Software Engineer | Paoli, PA 19301-1552 Image Signal Processing, Inc. | Phone: +1 610 407 4391 http://www.isptechinc.com | FAX: +1 610 407 4405
Re: Functional DSP
Philip Wadler wrote: By the way, FFTW although written in C, was generated by a functional program written in Caml. You can find more details on the web page I asked the authors (Steven Johnson and Matteo Frigo from MIT) a few months ago why they chose Objective Caml over SML, Haskell, and Scheme. The answer was bascially because they knew it already and it ran well on one of their older laptops. I think what it more interesting is that fact that the code uses explicit recursion (as oppose to "traditional" loop based algorithms) and all of the codelets/butterflys are executed via byte code. On top of that, it managed to beat most other FFT packages available today. --Matt Donadio ([EMAIL PROTECTED])
Re: Haskell 98 draft report
Hans Aberg wrote: Thinking of it, "round" should probably be viewed as a method to convert a float to another float of less precision (and not a conversion to an integer) To be picky, rounding a fixed point value to less bits is a very common procedure (at least it is in the DSP world) to acount for bit-growth in an algorithm. --Matt Donadio ([EMAIL PROTECTED])
Debugging Arrays
Hi all, I am having a tough time debugging an array problem, and I was wondering if anyone had some pointers. I am working on a complicated function that relies on a few mutually recursive arrays (AR parameter estimation using Burg's Method). I am getting a runtime error, index: Index out of range, when I test out the function. For the life of me, I cannot figure out where the error is. One complication is that not all of the arrays elements are defined for the index bounds. I have tried several simplifications to isolate the arrays, but when I simplify too far, I run into type problems. Is there an easy way to trace execution? Does anyone have any advice on how to go about debugging this? I am using Hugs 98, November 2002 if it matters. Thanks. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Arrays vs Lists and Specialization
Hi all, I'm sorry if this topic has been rehashed a lot, but I poked around in the mailing list archive, and didn't find what I was looking for. I currently have some free time on my hands, and have been implementing some digital signal processing and spectral/frequency estimation algorithms along with the needed matrix routines in Haskell. For those unfamiliar with this field, most algorthms are defined in textbooks in terms of indexing through discrete sequences. For example the implementation of cross-correlation rxy x y k | k = 0 = sum [ x!(i+k) * (conjugate (y!i)) | i - [0..(n-1-k)] ] | k 0 = conjugate (rxy y x (-k)) where n = snd (bounds x) + 1 looks very similar to the textbook definition if one uses arrays. A definition using lists would probably use drop, map, and zipWith, and look nothing like the definitions found in the standard texts. Many spectral estimation routines are defined in terms of special matrices (ie, Toeplitz, etc). Arrays defined recursively by list comprehensions make it easy to implement algorithms like Levinson-Durbin recursion, and they look very similar to the mathematical definitions: levinson r p = (array (1,p) [ (k, a!(p,k)) | k - [1..p] ], realPart (rho!p)) where a = array ((1,1),(p,p)) [ ((k,i), ak k i) | k - [1..p], i - [1..k] ] rho = array (1,p) [ (k, rhok k) | k - [1..p] ] ak 1 1 = -r!1 / r!0 ak k i | k==i = -(r!k + sum [ a!(k-1,l) * r!(k-l) | l - [1..(k-1)] ]) / rho!(k-1) | otherwise = a!(k-1,i) + a!(k,k) * (conjugate (a!(k-1,k-i))) rhok 1 = (1 - (abs (a!(1,1)))^2) * r!0 rhok k = (1 - (abs (a!(k,k)))^2) * rho!(k-1) This could be rewritten for lists, but would probably need to be defined in terms of an aux. recursive function, which destroys the simplicity of the above definition. OK, my question then has to do with the efficiency of lists versus arrays. Do the latest compilers handle handle arrays efficiently, or are lists really the way to go? If there is a performace difference, is it generally big enough to warrant rewriting algorithms? A related question is how is specilization handled in arrays with lazy evaluation. In the definition of levinson above, the a array is defined in terms of the ak function. By doing this, you save some horizontal space, but it also unburdens the programmer from tracking the recursive dependencies. a!(k,k) is needed before a!(i,j) can be calculated, but lazy evaluation takes care of this. If the above function is specialized for r::Array Int (Complex Double) and p::Int, would I be correct to say that the final value of the function would be unboxed, but all intermediate values wouldn't? Now, in some cases, a user may need all of the model orders from 1..p. This is handled easilly enough by just changing the first line. to levinson r p = (a, fmap realPart rho) Would the a matrix in the tuple be unboxed with specilization? If anyone is interesting in what I have put together, I will be making everything public sometime next week. I have a lot of algorithms implemented, but I need to clean up the documentation a bit (well, a lot). Thanks. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Class Question
Hi all, I have yet another question. I think I have too much time on my hands. I have two functions: rf :: (RealFloat a, Integral b) = [a] - b - Complex a rf x k = ... cf :: (RealFloat a, Integral b) = [Complex a] - b - Complex a cf x k = ... I would like to add these to the class system so I can create an overloaded version, f. From my understanding of the Report, the following should work, but I get an error, which I really don't understand. class (Fractional a) = Floating a where f :: (Integral b) = [a] - b - a instance (RealFloat a) = Floating a where f x k = rf x k instance (RealFloat a) = Floating (Complex a) where f x k = cf x k hugs98 spits out ERROR f.lhs:59 - Syntax error in instance head (constructor expected) where line 59 is the first instance declaration. Can anyone give me some guidance, or hints, here? Thanks a million. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Bizarre Haskell Problem
Hi all, I am currently a having a bizarre Haskell problem, and was wondering if anyone had any suggestions. This is a snippet from a source file: foo n = [ a ^* i | i - [0..(n-2)] ] where i ^* j = (i ^ j) `mod` n a = generator n rader :: Array Int (Complex Double) - Int - Array Int (Complex Double) rader f n = foo a n -- DEBUG: should be f' where h = listArray (0,n-2) [ f!(a ^* (n-(1+n'))) | n' - [0..(n-2)] ] g = listArray (0,n-2) [ w (a ^* n') | n' - [0..(n-2)] ] f' = array (0,n-1) ((0, sum [ f!i | i - [0..(n-1)] ]) : [ (a ^* i, f!0 + sum [ h!j * g!((i-j)`mod`(n-1)) | j - [0..(n-2)] ]) | i - [0..(n-2)] ]) w i = cis (-2 * pi * fromIntegral i / fromIntegral n) i ^* j = (i ^ j) `mod` n a = generator n Under hugs and ghc, calling 'foo' and 'rader' with the proper arguments will give me different results under certain circumtances. In rader, n is the number of elements in the array. If n is 23, then everything is OK. If n = 23, then rader returns the wrong result. The bizarre thing is that if I comment out the definitions of h and f' in rader, then it returns the correct results. Hugs also give different reduction counts depending on whether h and f' are commented out or not. If I add any more definitions that reference f, then rader misbehaves. What is even more bizarre is that if I copy everything into a separate file, then it work for all n. Any suggestions or hints, other than just use a separate file? Separate files isn't really an option, because the above is simplified a bit, and the final version would have mutually recursive modules, which hugs can't handle. Thanks. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
ANNOUNCE: DSP Libraries for Haskell
Hi all, I would like to annouce the availability of Haskell libraries for digital signal processing, spectral estiation, and frequency estimation. This project is a reincarnation of a project started several years ago, but went dormant for a while. The URL is http://users.snip.net/~donadio/haskell/ The home pages lists everything that is available. Of particular interest is: -- Module for LU decomposition of matrices (LU solver, matrix inversion, etc) -- Cholesky solver -- FFT module that is not limited to N=2^v. It includes Radix-2 decimation in time FFT for N=2^v Radix-2 decimation in frequency FFT for N=2^v Cooley-Tukey FFT for composite N Prime factor algorithm for composite N Rader's algorithm using direct convolution for prime N Rader's algorithm using FFT convolution fpr prime N Hard-coded FFT for a few small values of N Inverse FFT defined in terms of the FFT 2N-point real FFT computed with N-point FFT More algorithms will be appearing over time, as well as some demo applications. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Re: time since the epoch
Hi all, I'm jumping into this a bit late, but I have some good info about time. Two sites that I know about that have good tutorials and white papers are http://www.datum.com/res_technical.html http://www.truetime.com/DOCSn/TTreferencematerial.html In particlar, the paper Timing and Time Code Reference Book is very interesting. It explains several different time standards that are currently used, and tells you why UTC needs leap seconds. Pretty much the whole world runs on UTC. All of the common time distribution systems use UTC. Technically, GPS doesn't, but the GPS signal includes the correction to UTC. I understand the argument for using TAI. Maybe internally the libray should use TAI, but default to giving the user UTC? -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Re: Learning Haskell moved
I would add some significant papers, such as John Hughes' Why functional programming matters, etc. I've added that paper to the introduction. Other suggestions are welcome. Conception, Evolution, and Application of Functional Programming Languages, Paul Hudak, ACM Computing Surveys, Volume 21, Number 3, pp.359-411, 1989 This would probably be good, too. --Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Re: compiling
Mike T. Machenry wrote: I am having a problem compiling my code. Usually I run it with ghci -package data -fglashow-exts Main.hs Main declares a main function and imports all my other files. when I try to ghc it to compile I get that it can't find an interface file for each file in my project. How do I compile something? Also it can't find an interface for Data. How do I make an interface file for that? Read the ghc docs, but try ghc -o Main --make -package data -fglashow-exts Main.hs -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
FFI Tutorial / Examples
Hi all, I may be being a bit dense about this, but I am having some trouble understanding how to use FFI, especially with respect to interfacing Haskell lists/arrays to C arrays. For example, say I have the C functions void foo (double *input, double *output, int N); double bar (double *input, int N); and I want to create an FFI interface and have the resulting type signatures be foo :: Array Int Double - Array Int Double bar :: Array Int Double - Double where the bounds of the arrays are (0,N-1), and both foo and bar are pure. I read through the FFI docs, but I am still confused about how to do this. Can anyone point me to an FFI tutorial, or some examples? I have a feeling that once I see some examples using lists and arrays that things will fall into place. Thanks. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Looking for Libraries
Hi all, This is probably a long shot, but I am looking for a few libraries and don't want to put effort into something that has already been done. Is there a high level Haskell graphics library that would give functionality similar to gnuplot? I know I could build one myself, but I hate graphics programming. Does anyone have an FFI interface to LAPACK? Does anyone have an FFI interface to a free linear program solver? Or, a Haskell implmentation of the two-step simplex method? (I sold my copy of Chvatal's book in college for beer money). Thanks. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Re: Looking for Libraries
Ferenc Wagner wrote: Is there a high level Haskell graphics library that would give functionality similar to gnuplot? Why not simply USE gnuplot? Or plotutils? They have simple textual interfaces, do good work, and are fairly standard tools (on a Unix system, at least). I would like to be able to have plotting capabilities directly from a Haskell program rather than using a spawned process and temporary files, especially for interactive programs. This approach would probably work well enough, though. Does anyone have an FFI interface to LAPACK? Have a look at http://www.haskell.org/pipermail/haskell/2002-June/009833.html Duh. I thought I looked through the archives for this, but I guess I didn't. Thanks. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Re: main modules in GHC, deriving differences between GHC and Hugs
Graham Klyne wrote: GHC seems to require a 'main' module in a file in order to generate an exe file. This makes it awkward to create unit test programs, because I generally create one (to run stand-alone) for each major module. Now I want to create a master test module that runs all the individual module tests. But if the module tests are all main modules it seems I cannot bring them all together into a larger program. Have I overlooked any way to create an executable program from any module containing a main function of the appropriate type? The easiest way to handle this is to run all the source through the C preprocessor, and put #ifdef's around main and the module name. Something like #ifdef UNIT_TEST module Main where #else module Foo where #endif foo x y = x + y tests = [ foo 2 3 == 5, foo 3 4 /= 6 ] #ifdef UNIT_TEST main = print $ and tests #else foo_test = and tests #endif I haven't done this with Haskell, but I have done it with a lot of my C libraries. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Re: User-Defined Operators
Wolfgang Jeltsch wrote: A related problem is that I cannot see a way to define a new log-like function (as Lamport names them), i.e., a function with a name consisting of several letters which have to be set in upright font with no spaces between them. Examples are log, min, max, sin, cos and tan. Check out the AMS-LaTeX package. I think it has a macro to solve this. It also includes a zillion new symbols/operators. http://www.ams.org/tex/amslatex.html If you have TeTeX installed as your TeX system, then it should be included. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell
Bizarre Haskell Problem
Hello, I am having a bizarre Haskell problem that I am having difficulty debugging. I am not positive this is a compiler problem, but my results are not making any sense. I have attached a few source files which compiled with ghc-5.04.2 running under Win95. The files were compiled as: ghc -c DFT.lhs ghc -c FFT1.lhs ghc -c FFT2.lhs ghc -c Main1.hs ghc -c Mail2.hs ghc -o test1 Main1.o FFT1.o DFT.o ghc -o test2 Main2.o FFT2.o DFT.o Running test1 gives the following results (the last line is wrong): foo 19:[1,2,4,8,16,13,7,14,9,18,17,15,11,3,6,12,5,10] rader1 19: [1,2,4,8,16,13,7,14,9,18,17,15,11,3,6,12,5,10] foo 23:[1,5,2,10,4,20,8,17,16,11,9,22,18,21,13,19,3,15,6,7,12,14] rader1 23: [1,5,2,10,4,20,8,17,16,11,9,22,18,21,1,4,8,18,22,6,19,2] Running test2 gives the following results (these are the results I expect): foo 19:[1,2,4,8,16,13,7,14,9,18,17,15,11,3,6,12,5,10] rader1 19: [1,2,4,8,16,13,7,14,9,18,17,15,11,3,6,12,5,10] foo 23:[1,5,2,10,4,20,8,17,16,11,9,22,18,21,13,19,3,15,6,7,12,14] rader1 23: [1,5,2,10,4,20,8,17,16,11,9,22,18,21,13,19,3,15,6,7,12,14] The only difference bewteen the sources is that in FFT1.lhs, lines 215 and 217 are present, while in FFT2.lhs, they are commented out. Note that these two lines reference the parameter, f. Also note that rader1 calls foo, so I am confused as to how that can produce different results, as test1 shows. FFT1.rader1 works for all n = 19, but fails for n = 23 (n has to be prime, however). Also, if I copy the offending code from FFT1.lhs to a separate file, then I get the results I expect, but this is a less than ideal solution. ghc-5.04.2 was installed with the Windows Installer from the website. If you play with the code, foo n should produce a permutation of the sequence [1..(n-1)] for all prime n (ie, foo n produces the permutation for a generator of the Galois field n). I appologize in advance if this is a bug on my part, but based on what I am seeing, I am getting results that should not happen. Thanks. -- Matthew Donadio ([EMAIL PROTECTED]) FFT2.lhs Description: haskellprogram Main2.hs Description: haskellprogram DFT.lhs Description: haskellprogram FFT1.lhs Description: haskellprogram Main1.hs Description: haskellprogram
FFI Help
Hi all, I am just starting to experiment with FFI, and am running into a problem. I was to create an FFI to the lgamma(3) found in many of the newer libm implementations. My code follows the sig. The lgamma function works. The gamma function core dumps (I am using ghc 5.04.3) on me. gdb reports a SIGSEGV in signgam(), but I'm not sure why. I believe that I need to use the monad because signgam is only valid after lgamma returns. Does anyone have an idea what I am doing wrong? Thanks. -- Matthew Donadio [EMAIL PROTECTED] module Gamma (gamma, lgamma) where import System.IO.Unsafe foreign import ccall math.h lgamma lgammaC :: Double - IO Double foreign import ccall math.h signgam signgamC :: IO Int lgamma :: Double - Double lgamma x = unsafePerformIO $ lgammaC x gamma :: Double - Double gamma x = unsafePerformIO $ gammaIO x gammaIO :: Double - IO Double gammaIO x = do lg - lgammaC x s - signgamC return $ fromIntegral s * exp lg ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: Green Card and Exceptions
Alastair Reid wrote:
The %fail statements (described in the last few paragraphs of
http://www.haskell.org/greencard/downloads/greencard-latest/type-sig.html)
consist of two C expressions. For example:
%fail {f == NULL} {errstring(errno)}
The first is a test for failure.
The second is an expression which returns a C string.
If the test expression fails, the string expression is evaluated and used to
generate a UserError.
OK, I think I misread the manual. Sect 7.6 talks about functions with
side effects, so I assumed that the function had to have type (IO a) to
use %fail.
It would also be nice to be able to generate a different error instead of
UserError. We'd need to specify the type and the exception constructor so a
plausible syntax would be:
%fail {f == NULL} (UserError (string {errstring(errno)}))
[Detail: should it be a Haskell98 IOError constructor or a non-standard but
widely implemented exception constructor? Should it be a function or a
constuctor?]
In the case of what I am doing, I'm not sure if IOError really make
sense philosophically. The failures I need are underflow, overflow,
loss of precision, etc. Since IOError is a type synonym for
IOException, then perhaps accepting an Exception constructor is
appropriate.
To keep compatibility with old libraries it may be wise to keep %fail as
is, and have a new directive %throw that accepts an Exception
constructor, and uses either throw or throwIO.
On the other hand, now that I know that I can use %fail with pure
functions, I can make that work.
As GreenCard maintainer, I've got to ask:
- How many users of GreenCard are still out there?
New GreenCard user. In my case, I need access to C land for typedefs
and macros. I could write my own stubs, but GreenCard saves me this
step.
- Are you developing new libraries or just maintaining the ones you've got?
New library.
- Is there a demand for new features?
A more generic %fail mechanism?
Thanks for the response.
--
Matthew Donadio ([EMAIL PROTECTED])
___
Haskell-Cafe mailing list
[EMAIL PROTECTED]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Testing was Re: main modules in GHC, deriving differences between GHC and Hugs
(redirected to haskell-cafe) Hal Daume III wrote: Yes, but there's a problem with this solution. Namely, if Foo.hs takes a long time to compile, then you can't leverage having already created Foo.o and Foo.hi when making the Main. Yeah, it's not perfect, but I think we just have different methodologies for testing. Typically, I do most of my development with ghci or hugs. Each module will have a single variable that represents the tests tests :: [Bool] tests = [ test1, test2, test3 ... ] and then I define a variable called test :: Bool test = and tests so I can just load a module and either evaluate tests or test to check things out. This generalizes to importing and testing several modules at once (as long as I take care of name conflists). This only works for simple modules, though. For more complicated ones, I have a pretty tester. For example, my Haskell FFT library is collection of mutually recursive modules. I have a specialzed test function for this that tests a range of lengths. It looks someone like main = testfft n1 n2 testfft n1 n2 = sequence $ map test1fft [n1..n2] test1fft n = do putStr $ show n ++ :\t putStr $ if ok then OK\n else ERROR\n where ok = and [ test1 n, test2 n, test3 n ] This way I can compile the code, and run the executable as ffttest 2 2048 | grep ERROR and I am confident that I get full coverage of the algorithm. I am always to hear about other methods of automated testing. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: Poll: How to respond to homework questions
Tom Pledger wrote: I'm curious about what the people on this list consider appropriate, as responses to homework questions. Even if there isn't a consensus, it may be interesting to see how opinion is divided. Please consider the following. (A) Give a perfect answer. (B) Give a subtly flawed answer. (C) Give an obfuscated answer. (D) Give a critique of what the questioner has tried so far. (E) Give relevant general advice without answering the specific question. I mean them to include giving references, e.g. a Wiki URL to relevant general advice would count as (E). I think this depends on the situation. There is a big difference between I am having some trouble with this homework problem. This is what I did. Could someone give me some tips? Thanks and How do I write a map function in Haskell? For the first case, I would vote for D and/or E as appropriate. For the second case, I vote for (F) Ignore. -- Matthew Donadio ([EMAIL PROTECTED]) ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
