Options

The bacterium of nonlinear PDE

I plan to write one post in the series Fluid Flows and Infinite-Dimensional Manifolds about the nonlinear aspects of the Burgers equation. This equation turns up as the geodesic equation of the diffeormorphism group of the circle.

It is also the simplest example of a nonlinear equation, and it has solutions in closed forms. So, when Isaac Held writes about the fruit fly of climate models, the Burgers Equation is kind of the bacterium of nonlinear partial differential equations. If you are a coder without much experience in numerical mathematics and programming, this is a great place to start.

I would like to compare the solution in closed form in one (spatial) dimension to the approximation using the Fourier-Galerkin approximation. If we could get an online animation of the solutions, that would be a big big plus. I will try to get this done during the next two months, but cannot promise anything.

This thread is for anyone who would like to try to program this him/herself. Here are my thoughts about how this should be done so far:

  • Java is a pain when you need complex numbers already, because it does not come with this datatype and it does not have operator overloading and structures. So I won't use it.

  • C# shares one problem with Java, namely that there are no libraries for number crunching available, those are all in Fortran or in C++, so I won't use it either (or any language from the .NET platform),

  • Scala's overall tool chain was really disappointing when last I looked, so I won't use that either,

  • JavaScript is a language to create interactive web pages, I don't think it is appropriate for anything else,

  • closed source software like Matlab or Mathematica should not be used in science,

  • Sage is certainly a possibility, it uses C++ and Python,

  • Python is a toy, it is an interpreted scripting language which lacks a lot of concepts (type safety, for example) that would make it eligible for large software projects like Java or C++,

  • FORTRAN is known in academia only, at least I think so. Most coders out there won't know it, which includes me. And I am not inclined to spend months to learn it.

So I will try to use C++, but am still undecided which libraries I should use. Plus I have to invest some work to get fluent in C++, which I haven't used for a decade.

If you think I'm wrong and it is much easier to get the job done in a different way, we can talk about that here. If you agree with me and would like to help me, we can also talk about that here. Right now I am looking for the best way to include libaries for FFT and for the modified Bessel functions into an existing project, but I haven't checked in anything yet to our repository.

If there are several successful solutions to the posed problem, I will probably write a blog post about all of them, comparing them.

I will start today with stubs for said blog post and a page for C++ (doesn't exist right now), and will notify you here every other week or so, so that you can check my progress :-)

Comments

  • 1.

    Not necessarily disagreeing with your conclusion to use C++ but:

    If you had wanted to use a Matlab-style language, there's the open source Octave and (newer so less mature) Julia that you could use.

    Regarding your point about typesafety, as someone who is very familiar with the area and has even written a type-inference engine for a hindley-milner style type system, I'm not sure that, particularly in pure science, typechecking for "safety" is all that important. Science generally isn't something where you either burn millions of retail CDs or are running "control systems", the two big areas where compile time safety is really important.

    On the other hand, I think a language with inferred types is beneficial for science because the compiler/runtime can specialise the code for the types being used, which increases the execution speed.

    Comment Source:Not necessarily disagreeing with your conclusion to use C++ but: If you had wanted to use a Matlab-style language, there's the open source Octave and (newer so less mature) Julia that you could use. Regarding your point about typesafety, as someone who is very familiar with the area and has even written a type-inference engine for a hindley-milner style type system, I'm not sure that, particularly in pure science, typechecking for "safety" is all that important. Science generally isn't something where you either burn millions of retail CDs or are running "control systems", the two big areas where compile time safety is really important. On the other hand, I think a language with inferred types is beneficial for science because the compiler/runtime can specialise the code for the types being used, which increases the execution speed.
  • 2.

    Not necessarily disagreeing with your conclusion to use C++ but...

    Yeah, I guessed that :-) In our previous discussions, I took the position of promoting Java because you already were on the C++ side. It's a bad habit of mine: I often argue for a position that I don't necessarily fully agree with, but that is underrepresented in the ongoing discussion. This has lead to a lot of misunderstandings in the past, so I try to get rid of it or at least explain it better.

    If you had wanted to use a Matlab-style language, there's the open source Octave and (newer so less mature) Julia that you could use.

    Right, Octave is missing on my list. I think I need full control of the implementation at the abstraction level of a 4th generation programming language, so I won't use Octave, either, but maybe others know better :-)

    Regarding your point about typesafety, as someone who is very familiar with the area and has even written a type-inference engine for a hindley-milner style type system, I'm not sure that, particularly in pure science, typechecking for "safety" is all that important. Science generally isn't something where you either burn millions of retail CDs or are running "control systems", the two big areas where compile time safety is really important.

    I don't have any experience with big software projects in science like climate models. I'm simply talking about my own limited experience which tells me that I would not like to take responsibility for a big project written in Python. I hope that this thread could shed some light on different experience and best practices like this one.

    Comment Source:<blockquote> <p> Not necessarily disagreeing with your conclusion to use C++ but... </p> </blockquote> Yeah, I guessed that :-) In our previous discussions, I took the position of promoting Java because you already were on the C++ side. It's a bad habit of mine: I often argue for a position that I don't necessarily fully agree with, but that is underrepresented in the ongoing discussion. This has lead to a lot of misunderstandings in the past, so I try to get rid of it or at least explain it better. <blockquote> <p> If you had wanted to use a Matlab-style language, there's the open source Octave and (newer so less mature) Julia that you could use. </p> </blockquote> Right, Octave is missing on my list. I think I need full control of the implementation at the abstraction level of a 4th generation programming language, so I won't use Octave, either, but maybe others know better :-) <blockquote> <p> Regarding your point about typesafety, as someone who is very familiar with the area and has even written a type-inference engine for a hindley-milner style type system, I'm not sure that, particularly in pure science, typechecking for "safety" is all that important. Science generally isn't something where you either burn millions of retail CDs or are running "control systems", the two big areas where compile time safety is really important. </p> </blockquote> I don't have any experience with big software projects in science like climate models. I'm simply talking about my own limited experience which tells me that I would not like to take responsibility for a big project written in Python. I hope that this thread could shed some light on different experience and best practices like this one.
  • 3.

    I had some problems to use the GSL, because my compiler could not locate the library at link time - the GSL itself, that is its header files, are included in my Cygwin installation. Someday I will have to find out how to configure Eclipse so that this works. This is one aspect of C++ that is really very difficult for me, because in Java all you have to do is to put your libraries on your Classpath. I was successful with the inclusion of boost for the modified bessel functions of the first kind, so that I could write a class that calculates the exact solution. I have to write tests and all tbat stuff, of course. For the spectral approximation, I need an FFT and an ODE solver. I found the library ALGLIB (http://www.alglib.net/) which has everything I need and was successful with including it in my project, and compile and use it.

    So I think I have all the infrastructure I need in place and can concentrate on programming. There is not much to do, but if you have to look up how to declare constructors and stuff, that still takes a while ;-)

    Java is really a faux amis when you write C++, because the syntax is similar, but semantics can be quite different.

    Comment Source:I had some problems to use the GSL, because my compiler could not locate the library at link time - the GSL itself, that is its header files, are included in my Cygwin installation. Someday I will have to find out how to configure Eclipse so that this works. This is one aspect of C++ that is really very difficult for me, because in Java all you have to do is to put your libraries on your Classpath. I was successful with the inclusion of boost for the modified bessel functions of the first kind, so that I could write a class that calculates the exact solution. I have to write tests and all tbat stuff, of course. For the spectral approximation, I need an FFT and an ODE solver. I found the library ALGLIB (http://www.alglib.net/) which has everything I need and was successful with including it in my project, and compile and use it. So I think I have all the infrastructure I need in place and can concentrate on programming. There is not much to do, but if you have to look up how to declare constructors and stuff, that still takes a while ;-) Java is really a faux amis when you write C++, because the syntax is similar, but semantics can be quite different.
  • 4.
    edited June 2012

    I ran the calculation of the exact solution for some parameters and it looks good. I will clean up the code and someday soon check it in. I also uploaded a picture of a solution on Burgers equation. The best way to do QA on it is the implementation of a different algorithm: The spectral approximation. This will be next.

    Comment Source:I ran the calculation of the exact solution for some parameters and it looks good. I will clean up the code and someday soon check it in. I also uploaded a picture of a solution on [[Burgers equation]]. The best way to do QA on it is the implementation of a different algorithm: The spectral approximation. This will be next.
  • 5.

    Tim wrote

    I had some problems to use the GSL, because my compiler could not locate the library at link time - the GSL itself, that is its header files, are included in my Cygwin installation.

    I presume you know that on "Unixy platforms" (incl. cygwin) at the lowest level it all happens in terms of flags to the compiler invocations

    -I/path/to/a/directory/with/somerelevantheaders

    and including the paths to actual libraries after object files (eg, /here/is/my/library.so) when invoking linking.

    Someday I will have to find out how to configure Eclipse so that this works.

    Yeah, quite how the Eclipse GUI lets you set project structure I something that's scared me off from using Eclipse. (Well, that plus it's memory footprint...)

    Comment Source:Tim wrote > I had some problems to use the GSL, because my compiler could not locate the library at link time - the GSL itself, that is its header files, are included in my Cygwin installation. I presume you know that on "Unixy platforms" (incl. cygwin) at the lowest level it all happens in terms of flags to the compiler invocations -I/path/to/a/directory/with/somerelevantheaders and including the paths to actual libraries after object files (eg, /here/is/my/library.so) when invoking linking. > Someday I will have to find out how to configure Eclipse so that this works. Yeah, quite how the Eclipse GUI lets you set project structure I something that's scared me off from using Eclipse. (Well, that plus it's memory footprint...)
  • 6.

    I think something should be said here about dynamic versus static type checking, and strict versus loose type checking.

    Dynamic versus static type systems

    Languages such as Python, Javascript, PHP, Lua, Matlab and Octave are dynamically typed. This means that the types of variables are not declared to the compiler; variables acquire their types at run-time. This dynamic nature means the compiler can do almost nothing to check types, or check that variables are defined before use. These checks occur at run-time, when variables are used. One problem with a dynamic type system is what I call "lurking typos". Because type errors and undefined variables are not detected at compile time, it is possible to run code that contains errors in parts of the code that are not exercised, e.g. a branch of a conditional or a function definition. This means that, for anything other than trivial scripts and interactive sessions, extra care is needed to properly exercise all parts of the code during testing. Another well-known disadvantage of a dynamic type system is computing inefficiency: when running a loop over many items, type checks are done for each iteration, which can easily take far more compute time than the main computation.

    Languages such as C, C++, Java and Haskell are statically typed. In most statically-typed languages, the types of variables must be declared in the code, and the compiler can then check that types are used correctly, and that variables have been declared before use. In languages such as Haskell, the compiler can infer almost all types from the context where variables are used, so most explicit declarations can be omitted. One big advantage of static typing is that the programmer can convey much useful information to the compiler, so a good compiler can 1) detect potential nonsense, and 2) optimise the compiled target for speed/space. If the language and compiler are very good, then we may approach the ideal of "if it compiles, it works", so avoiding much tedious testing. The main downsides are that static type systems tend to be a good deal more complex than dynamic ones, and that there is a tendency for verbosity where the type of everything has to be declared.

    Strict versus loose type systems

    If a type system is strict, then there are restricted ways in which variables may be employed in expressions, depending on their types. For example, adding a number and a string would not be allowed, because it is nonsense. In a loose type system, such expressions may have a meaning. For example, the number could be implicitly converted to a string, and concatenated with the string variable, yielding a string as the result.

    While it is common for dynamically-typed languages to be more loosely typed than statically-typed languages, the ideas of dynamicism and looseness should not be conflated. In my opinion, Python's type system is quite strict, compared to, say, Perl or PHP. And it would be quite possible for some miscreant to overload operator+() in C++ to allow addition of strings and numbers.

    Comment Source:I think something should be said here about dynamic versus static type checking, and strict versus loose type checking. ###Dynamic versus static type systems### Languages such as Python, Javascript, PHP, Lua, Matlab and Octave are dynamically typed. This means that the types of variables are not declared to the compiler; variables acquire their types at run-time. This dynamic nature means the compiler can do almost nothing to check types, or check that variables are defined before use. These checks occur at run-time, when variables are used. One problem with a dynamic type system is what I call "lurking typos". Because type errors and undefined variables are not detected at compile time, it is possible to run code that contains errors in parts of the code that are not exercised, e.g. a branch of a conditional or a function definition. This means that, for anything other than trivial scripts and interactive sessions, extra care is needed to properly exercise all parts of the code during testing. Another well-known disadvantage of a dynamic type system is computing inefficiency: when running a loop over many items, type checks are done for each iteration, which can easily take far more compute time than the main computation. Languages such as C, C++, Java and Haskell are statically typed. In most statically-typed languages, the types of variables must be declared in the code, and the compiler can then check that types are used correctly, and that variables have been declared before use. In languages such as Haskell, the compiler can infer almost all types from the context where variables are used, so most explicit declarations can be omitted. One big advantage of static typing is that the programmer can convey much useful information to the compiler, so a good compiler can 1) detect potential nonsense, and 2) optimise the compiled target for speed/space. If the language and compiler are _very_ good, then we may approach the ideal of "if it compiles, it works", so avoiding much tedious testing. The main downsides are that static type systems tend to be a good deal more complex than dynamic ones, and that there is a tendency for verbosity where the type of everything has to be declared. ###Strict versus loose type systems### If a type system is _strict_, then there are restricted ways in which variables may be employed in expressions, depending on their types. For example, adding a number and a string would not be allowed, because it is nonsense. In a _loose_ type system, such expressions may have a meaning. For example, the number could be implicitly converted to a string, and concatenated with the string variable, yielding a string as the result. While it is common for dynamically-typed languages to be more loosely typed than statically-typed languages, the ideas of dynamicism and looseness should not be conflated. In my opinion, Python's type system is quite strict, compared to, say, Perl or PHP. And it would be quite possible for some miscreant to overload `operator+()` in C++ to allow addition of strings and numbers.
  • 7.

    I will try to keep that in mind - I think I have confused both concepts occasionally. When I wrote my "final judgement of Python" at the beginning of this thread, I mostly thought about the "dynamic class system" of Python. Good for little scripts, saves you some typing. But in large projects, you don't lose time typing, you loose time chasing bugs and trying to understand legacy code.

    Comment Source:I will try to keep that in mind - I think I have confused both concepts occasionally. When I wrote my "final judgement of Python" at the beginning of this thread, I mostly thought about the "dynamic class system" of Python. Good for little scripts, saves you some typing. But in large projects, you don't lose time typing, you loose time chasing bugs and trying to understand legacy code.
  • 8.

    Glyn wrote:

    The main downsides are that static type systems tend to be a good deal more complex than dynamic ones, and that there is a tendency for verbosity where the type of everything has to be declared.

    Haskell's ghc compiler has the most advanced type checker to come out of the research world. Type declarations are mostly optional.

    Comment Source:Glyn wrote: > The main downsides are that static type systems tend to be a good deal more complex than dynamic ones, and that there is a tendency for verbosity where the type of everything has to be declared. Haskell's ghc compiler has the most advanced type checker to come out of the research world. Type declarations are mostly optional.
  • 9.

    David wrote:

    I presume you know that on "Unixy platforms" (incl. cygwin) at the lowest level it all happens in terms of flags to the compiler invocations -I/path/to/a/directory/with/somerelevantheaders and including the paths to actual libraries after object files (eg, /here/is/my/library.so) when invoking linking.

    Yes I know, but it is better to assume that I don't. In this case I had no patience to try to configure eclipse.

    Comment Source:David wrote: <blockquote> <p> I presume you know that on "Unixy platforms" (incl. cygwin) at the lowest level it all happens in terms of flags to the compiler invocations -I/path/to/a/directory/with/somerelevantheaders and including the paths to actual libraries after object files (eg, /here/is/my/library.so) when invoking linking. </p> </blockquote> Yes I know, but it is better to assume that I don't. In this case I had no patience to try to configure eclipse.
  • 10.
    edited June 2012

    Jim wrote:

    Haskell's ghc compiler has the most advanced type checker to come out of the research world. Type declarations are mostly optional.

    Hic Rhodos, hic salta! Since you are impressed by Haskell, try to solve the Burgers equation finger exercise, and implement a computation of both the solution in closed form and its spectral approximation (of the 1D problem with periodic boundary conditions and the sine function as initial condition).

    Since the last time I said "Hic Rhodos, hic salta", the other guy came back and told me "I have found Rhodos on the map, but where is Salta?", here is the explanation from Wikipedia:

    In ancient times there was a Roman saying: "Hic Rhodus, hic salta!"—"Rhodes is here, here perform your jump", an admonition to prove one's idle boasts by deed rather than talk. It comes from an Aesop's fable called "The Boastful Athlete", and was cited by Hegel and Marx.

    Of course, "idle boasts" is a lot harsher that what I intend to say, but I think you get the point :-) ("salta" is the imperative "jump!" of the verb saltare).

    Comment Source:Jim wrote: <blockquote> <p> Haskell's ghc compiler has the most advanced type checker to come out of the research world. Type declarations are mostly optional. </p> </blockquote> Hic Rhodos, hic salta! Since you are impressed by Haskell, try to solve the [[Burgers equation]] finger exercise, and implement a computation of both the solution in closed form and its spectral approximation (of the 1D problem with periodic boundary conditions and the sine function as initial condition). Since the last time I said "Hic Rhodos, hic salta", the other guy came back and told me "I have found Rhodos on the map, but where is Salta?", here is the explanation from Wikipedia: <blockquote> <p> In ancient times there was a Roman saying: "Hic Rhodus, hic salta!"—"Rhodes is here, here perform your jump", an admonition to prove one's idle boasts by deed rather than talk. It comes from an Aesop's fable called "The Boastful Athlete", and was cited by Hegel and Marx. </p> </blockquote> Of course, "idle boasts" is a lot harsher that what I intend to say, but I think you get the point :-) ("salta" is the imperative "jump!" of the verb saltare).
  • 11.
    edited March 2013

    Tim wrote:

    Of course, "idle boasts" is a lot harsher that what I intend to say, but I think you get the point :-)

    I fail to see how qualifying Glyn's statement that statically-typed languages have a:

    tendency for verbosity where the type of everything has to be declared.

    has anything to do with boasting, so no, I don't get the point. Whatever it is, I can't guess how it might be likely to encourage me to help work on Burgers equation.

    My efforts at trying to program equations from Ch 3 of North have been repeatedly postponed as I have been spending a lot of time helping Glyn with various interface and stochastic issues.

    I am now continuing work on these slab ocean models (although nobody has suggested any interactive models to build) and don't have the mileage for anything else at the moment.

    When I get round to writing a wiki page, I will point people at Dan Piponi's Eleven reasons to use Haskell as a mathematician and Jerzy Karczmarcuk's papers on functional scientific programming. His most recent paper is here. If I do this I hope it will not be taken as some claim to be able to program anything I've read.

    Comment Source:Tim wrote: > Of course, "idle boasts" is a lot harsher that what I intend to say, but I think you get the point :-) I fail to see how qualifying Glyn's statement that statically-typed languages have a: > tendency for verbosity where the type of everything has to be declared. has anything to do with boasting, so no, I don't get the point. Whatever it is, I can't guess how it might be likely to encourage me to help work on Burgers equation. My efforts at trying to program equations from Ch 3 of North have been repeatedly postponed as I have been spending a lot of time helping Glyn with various interface and stochastic issues. I am now continuing work on these slab ocean models (although nobody has suggested any interactive models to build) and don't have the mileage for anything else at the moment. When I get round to writing a wiki page, I will point people at Dan Piponi's [Eleven reasons to use Haskell as a mathematician](http://blog.sigfpe.com/2006/01/eleven-reasons-to-use-haskell-as.html) and Jerzy Karczmarcuk's [papers](https://users.info.unicaen.fr/~karczma/arpap/) on functional scientific programming. His most recent paper is [here](http://arxiv.org/html/1109.0323v1). If I do this I hope it will not be taken as some claim to be able to program anything I've read.
  • 12.
    edited June 2012

    Sorry, I hope I did not offend you in any way. Let me try it another way:

    My efforts at trying to program equations from Ch 3 of North have been repeatedly postponed as I have been spending a lot of time helping Glyn with various interface and stochastic issues. I am now continuing work on these slab ocean models (although nobody has suggested any interactive models to build) and don't have the mileage for anything else at the moment.

    That's great! I'm sure we can make a blog post showing the results from those models, interactive or not. I'll certainly try to help as far as I am able.

    When I get round to writing a wiki page, I will point people at Dan Piponi's ~'11 reasons a mathematician should use Haskell' and Jerzy Karczmarcuk's work on functional scientific programming.

    I could start a stub on Haskell and we see where we go from there...

    If I do this I hope it will not be taken as some claim to be able to program anything I've read.

    No, of course not - my point was that the Burgers equation project is a good litmus test if a programming language / framework is helpful for climate models, because you need to

    • compute special functions (Bessel functions, in this case, spherical harmonics, for sure, for more general applications),

    • use standard algorithms like FFT and ODE solvers,

    • produce input and output that can be processed by some visualization software,

    • maybe do some of the steps in parallel.

    So the question is: Is Haskell useful for this kind of stuff or does "mathematician should use Haskell" refer to some more special context? (NB: the last time I did anything with a functional programming language was playing around a bit with Gofer 10 years ago. Then it did not seem to be useful for anything but torture students.)

    Comment Source:Sorry, I hope I did not offend you in any way. Let me try it another way: <blockquote> <p> My efforts at trying to program equations from Ch 3 of North have been repeatedly postponed as I have been spending a lot of time helping Glyn with various interface and stochastic issues. I am now continuing work on these slab ocean models (although nobody has suggested any interactive models to build) and don't have the mileage for anything else at the moment. </p> </blockquote> That's great! I'm sure we can make a blog post showing the results from those models, interactive or not. I'll certainly try to help as far as I am able. <blockquote> <p> When I get round to writing a wiki page, I will point people at Dan Piponi's ~'11 reasons a mathematician should use Haskell' and Jerzy Karczmarcuk's work on functional scientific programming. </p> </blockquote> I could start a stub on Haskell and we see where we go from there... <blockquote> <p> If I do this I hope it will not be taken as some claim to be able to program anything I've read. </p> </blockquote> No, of course not - my point was that the [[Burgers equation]] project is a good litmus test if a programming language / framework is helpful for climate models, because you need to * compute special functions (Bessel functions, in this case, spherical harmonics, for sure, for more general applications), * use standard algorithms like FFT and ODE solvers, * produce input and output that can be processed by some visualization software, * maybe do some of the steps in parallel. So the question is: Is Haskell useful for this kind of stuff or does "mathematician should use Haskell" refer to some more special context? (NB: the last time I did anything with a functional programming language was playing around a bit with Gofer 10 years ago. Then it did not seem to be useful for anything but torture students.)
  • 13.
    edited June 2012

    I agree with Tim that coding up a realistic test case such as solutions to Burgers' Equation is a good way to compare programming languages. When I started investigating many different programming languages, my test case was a linear regression analysis of some real-life experimental data.

    My experience with Haskell versus, say, C++ and Python is that the core numeric code can be expressed very clearly in Haskell, especially when there are recursive definitions in the maths, which is very often then case when analysing or generating chunks of data. You get the terseness of Python, but with static typing and compilation to native code, as in C++.

    The first problem I encountered with Haskell was with IO. Haskell is not, in my experience, particularly suited to munching on text and writing large amounts of text to files. It happens that my linear regression test problem required reading data from a large log file, in an inconvenient format, that could contain errors. Python lapped this up, as one might expect from a scripting language. Coding the same in Haskell was a fair bit more awkward. One difficulty is the odd type construction and syntax required to get Haskell to deal with a mutable outside world.

    Mutable state is bread-and-butter in most languages in common use, but can at best only be emulated in a pure functional language like Haskell. For many computations, a natural way to proceed is to update the contents of objects such as local variables, arrays and records at each stage in the computation, and to express this as a sequence of operations. This computing model will be familiar to anyone who has programmed in languages like C, C++, Java, and the common scripting languages. In Haskell, one may be able to recast the computation to avoid mutable state, e.g. by using recursion, or one might find or contruct some special data types that emulate mutable state. This involves additional programming effort.

    I suppose it can be argued that this awkwardness in expressing certain types of computation in Haskell arises from preconceptions of how things should be done, derived from previous experience in imperative and object-orientated languages. I would be the first to admit that this would apply to me. So perhaps I would not have encountered such problems if my first language was a functional language.

    Comment Source:I agree with Tim that coding up a realistic test case such as solutions to Burgers' Equation is a good way to compare programming languages. When I started investigating many different programming languages, my test case was a linear regression analysis of some real-life experimental data. My experience with Haskell versus, say, C++ and Python is that the core numeric code can be expressed very clearly in Haskell, especially when there are recursive definitions in the maths, which is very often then case when analysing or generating chunks of data. You get the terseness of Python, but with static typing and compilation to native code, as in C++. The first problem I encountered with Haskell was with IO. Haskell is not, in my experience, particularly suited to munching on text and writing large amounts of text to files. It happens that my linear regression test problem required reading data from a large log file, in an inconvenient format, that could contain errors. Python lapped this up, as one might expect from a scripting language. Coding the same in Haskell was a fair bit more awkward. One difficulty is the odd type construction and syntax required to get Haskell to deal with a mutable outside world. Mutable state is bread-and-butter in most languages in common use, but can at best only be emulated in a pure functional language like Haskell. For many computations, a natural way to proceed is to update the contents of objects such as local variables, arrays and records at each stage in the computation, and to express this as a sequence of operations. This computing model will be familiar to anyone who has programmed in languages like C, C++, Java, and the common scripting languages. In Haskell, one may be able to recast the computation to avoid mutable state, e.g. by using recursion, or one might find or contruct some special data types that emulate mutable state. This involves additional programming effort. I suppose it can be argued that this awkwardness in expressing certain types of computation in Haskell arises from preconceptions of how things should be done, derived from previous experience in imperative and object-orientated languages. I would be the first to admit that this would apply to me. So perhaps I would not have encountered such problems if my first language was a functional language.
  • 14.
    edited June 2012

    Sorry for coming in on this a bit late: I've been away.

    Anyway, I think I ought to say that I think Glyn's [not Jim as I originally wrote] discussion of dynamic vs static typing missed that there's (at least) three levels:

    1. Every variable has exactly one (possibly polymorphic) type that's known at compile time.

    2. Every variable has exactly one (possibly polymorphic) type that may not be known at compile time but only at run-time.

    3. It's valid for a variable, depending on the execution path, to have one of several different types at a given point in the program.

    As an example

    if cond: x=1
    else: x="message"
    print(x)
    

    is a trivially simple example of this that's a valid Python program that the run-time must run correctly.

    2 is what Julia uses (at least away from the very top-level command line) while 3 applies to Python, Javascript, etc. With 2 the Julia guys seem to be saying you get most of the speed advantages of 1 through using a JIT, but without the static assertion of types all being consistent at "compile time". 3 is the case you really don't want with scientific code.

    Regarding the question of whether "mathematicians should only use Haskell in certain contexts": I really like Haskell and use it for a lot of things. However, a big issue with Haskell is that it's a lazy functional language, ie at each evaluation step the runtime has to figure out what's the minimum thing it must do in order to make progress in the computation. On the one hand this is great in that it's possible to write lots of simple, intuitive programs that give the correct result even when there are "problems" in irrelevant parts of the input. (Eg, you can do calculations which in the code apply operations over infinite lists and providing the output is only affected by a finite portion of the lists you'll get a result, unlike in other languages where you'd have to explicitly code a careful evaluation strategy in your user-level code.)

    However, this comes at a cost since quite often things that will eventually need to be evaluated won't "need" to be evaluated at the time they are created. For example, pseudo-code

    let x=[1..100] ; y=fold (+) x
    in printLn (show y)
    

    will (without extra "strictness annotations") create a big expression tree

    (100+(99+(98+(97+....))))
    

    during the first pass over x, and only evaluate it when it's used. In contrast, "eager" functional languages and imperative languages would accumulate the sum so y is never more than a machine integer. This means that lazy functional languages often have much more unpredictable time and memory behaviour, with seemingly inconsequential changes in the code dramatically altering the memory usage, relative to other languages.

    So personally I'd use Haskell for a "complex, subtle, input-dependent" algorithm like a shot, but I wouldn't use it for anything like the core engine of an ODE solver where speed is the dominant concern and the algorithm is relatively direct and independent of the input.

    Comment Source:Sorry for coming in on this a bit late: I've been away. Anyway, I think I ought to say that I think Glyn's [not Jim as I originally wrote] discussion of dynamic vs static typing missed that there's (at least) three levels: 1. Every variable has exactly one (possibly polymorphic) type that's known at compile time. 2. Every variable has exactly one (possibly polymorphic) type that may not be known at compile time but only at run-time. 3. It's valid for a variable, depending on the execution path, to have one of several _different_ types at a given point in the program. As an example if cond: x=1 else: x="message" print(x) is a trivially simple example of this that's a valid Python program that the run-time must run correctly. 2 is what Julia uses (at least away from the very top-level command line) while 3 applies to Python, Javascript, etc. With 2 the Julia guys seem to be saying you get most of the speed advantages of 1 through using a JIT, but without the static assertion of types all being consistent at "compile time". 3 is the case you really don't want with scientific code. Regarding the question of whether "mathematicians should only use Haskell in certain contexts": I really like Haskell and use it for a lot of things. However, a big issue with Haskell is that it's a _lazy_ functional language, ie at each evaluation step the runtime has to figure out what's the _minimum_ thing it _must_ do in order to make progress in the computation. On the one hand this is great in that it's possible to write lots of simple, intuitive programs that give the correct result even when there are "problems" in irrelevant parts of the input. (Eg, you can do calculations which in the code apply operations over infinite lists and providing the output is only affected by a finite portion of the lists you'll get a result, unlike in other languages where you'd have to explicitly code a careful evaluation strategy in your user-level code.) However, this comes at a cost since quite often things _that will eventually need to be evaluated_ won't "need" to be evaluated at the time they are created. For example, pseudo-code let x=[1..100] ; y=fold (+) x in printLn (show y) will (without extra "strictness annotations") create a big expression tree (100+(99+(98+(97+....)))) during the first pass over x, and only evaluate it when it's used. In contrast, "eager" functional languages and imperative languages would accumulate the sum so y is never more than a machine integer. This means that lazy functional languages often have much more unpredictable time and memory behaviour, with seemingly inconsequential changes in the code dramatically altering the memory usage, relative to other languages. So personally I'd use Haskell for a "complex, subtle, input-dependent" algorithm like a shot, but I wouldn't use it for anything like the core engine of an ODE solver where speed is the dominant concern and the algorithm is relatively direct and independent of the input.
  • 15.
    edited June 2012

    David wrote:

    Jim's discussion of dynamic vs static typing

    Not me guv, Glyn wot dun it.

    If I can type the correct words in the correct order. Static versus dynamic is not the same as strong versus weak typing.

    hth

    Comment Source:David wrote: > Jim's discussion of dynamic vs static typing Not me guv, Glyn wot dun it. If I can type the correct words in the correct order. Static versus dynamic is not the same as strong versus weak typing. hth
  • 16.
    edited June 2012

    David wrote:

    a big issue with Haskell is that it's a lazy functional language, ie at each evaluation step the runtime has to figure out what's the minimum thing it must do in order to make progress in the computation.

    In Jerzy's papers with code examples of automatic differentiation and differential forms (linked at #12 above), it's precisely laziness which allows Haskell to deal with infinite lists and construct infinite towers of derivatives and series expansions: an example of the concision you refer to.

    Evaluating an infinite list from the end is not a good example and nobody would do it. foldr is mathematically sound but unusable on infinite types but is asymmetric with foldl which is not a correct mathematical definition but only uses what it needs when it's needed on infinite types. I don't think your criticism is borne out today; the price of functional programming is just profiling and this seems to have worked in all cases I've come across.

    Comment Source:David wrote: > a big issue with Haskell is that it's a lazy functional language, ie at each evaluation step the runtime has to figure out what's the minimum thing it must do in order to make progress in the computation. In Jerzy's papers with code examples of automatic differentiation and differential forms (linked at #12 above), it's precisely laziness which allows Haskell to deal with infinite lists and construct infinite towers of derivatives and series expansions: an example of the concision you refer to. Evaluating an infinite list from the end is not a good example and nobody would do it. foldr is mathematically sound but unusable on infinite types but is asymmetric with foldl which is not a correct mathematical definition but only uses what it needs when it's needed on infinite types. I don't think your criticism is borne out today; the price of functional programming is just profiling and this seems to have worked in all cases I've come across.
  • 17.
    edited June 2012

    Here's the Haskell GSL wrapper package function (a wiki page by a mathematician explaining this and optimal algorithms seems appropriate :).

    'besselJn :: Int -> Double -> Either String (Double, Double)
    
    besselJn n x = unsafePerformIO $
    alloca $ \gslSfPtr -> do
        c_deactivate_gsl_error_handler
    
        status <- c_besselJn (fromIntegral n) (realToFrac x) gslSfPtr
        if status == 0
            then do
                GslSfResult val err <- peek gslSfPtr
                return $ Right (realToFrac val, realToFrac err)
           else do
               error <- c_error_string status
               error_message <- peekCString error
               return $ Left ("GSL error: "++error_message)'
    

    or if this is sub-optimal:

    -- | yn() computes the Bessel function of the second kind for the
    -- integer Bessel0 n for the positive integer value x (expressed as a
    -- double).
    --
    yn :: Int -> Double -> Double
    yn x y = realToFrac (c_yn (fromIntegral x) (realToFrac y))
    

    is in cmath wrapper package which wraps what it says.

    You don't have to have beliefs to reap the benefits, just test :)

    Comment Source:Here's the [Haskell GSL wrapper package](http://en.wikibooks.org/wiki/Haskell/FFI#Implementing_the_Bessel_Function) function (a wiki page by a mathematician explaining this and optimal algorithms seems appropriate :). 'besselJn :: Int -> Double -> Either String (Double, Double) besselJn n x = unsafePerformIO $ alloca $ \gslSfPtr -> do c_deactivate_gsl_error_handler status <- c_besselJn (fromIntegral n) (realToFrac x) gslSfPtr if status == 0 then do GslSfResult val err <- peek gslSfPtr return $ Right (realToFrac val, realToFrac err) else do error <- c_error_string status error_message <- peekCString error return $ Left ("GSL error: "++error_message)' or if this is sub-optimal: -- | yn() computes the Bessel function of the second kind for the -- integer Bessel0 n for the positive integer value x (expressed as a -- double). -- yn :: Int -> Double -> Double yn x y = realToFrac (c_yn (fromIntegral x) (realToFrac y)) is in [cmath wrapper package](http://hackage.haskell.org/package/cmath-0.3) which wraps what it says. You don't have to have beliefs to reap the benefits, just test :)
  • 18.

    Apologies about mixing up the names: I'm a bit jetlagged. You're right about strong vs weak typing, the point I was more making was that "strong, dynamic" may be as almost as efficient as "strong, static" (but without compile-time assurances), it's "weak" that really screws up efficiency.

    I think laziness is absolutely wonderful except in the one case that strict evaluation would work correctly and you really, really care about speed. This blog article, which I see is 4 years old now..., contains the line

    Which is practically identical to the naive list version, we simply replaced foldl’ with foldlU, and [n .. m] with an explict enumeration (as list enumerator syntax isn’t available).
    

    I just about understand what he's doing after spending an MSc thesis using functional programming, but I just wonder whether "if you're not at the SPJ/Lennart/Don" level of understanding if this stuff is as obvious.

    (If I could stomach the verbose syntax I might get around to trying eager language O'Caml.)

    Comment Source:Apologies about mixing up the names: I'm a bit jetlagged. You're right about strong vs weak typing, the point I was more making was that "strong, dynamic" may be as almost as efficient as "strong, static" (but without compile-time assurances), it's "weak" that really screws up efficiency. I think laziness is absolutely wonderful except in the _one_ case that strict evaluation would work correctly and you really, really care about speed. This [blog article, which I see is 4 years old now...,](http://donsbot.wordpress.com/2008/06/04/haskell-as-fast-as-c-working-at-a-high-altitude-for-low-level-performance/) contains the line Which is practically identical to the naive list version, we simply replaced foldl’ with foldlU, and [n .. m] with an explict enumeration (as list enumerator syntax isn’t available). I just about understand what he's doing after spending an MSc thesis using functional programming, but I just wonder whether "if you're not at the SPJ/Lennart/Don" level of understanding if this stuff is as obvious. (If I could stomach the verbose syntax I might get around to trying eager language O'Caml.)
  • 19.
    edited June 2012

    David Tweed wrote:

    Anyway, I think I ought to say that I think Jim's [Glyn's] discussion of dynamic vs static typing missed that there's (at least) three levels: ...

    Thanks David, that is a much clearer explanation than mine.

    This means that lazy functional languages often have much more unpredictable time and memory behaviour ...

    This fact pains me a good deal, actually. I believe in the value of code analysis, which for me means figuring out what the code actually does, versus intentions and expectations. I am used to using code analysis for debugging in imperative languages such as C, where it is usually possible to track down where and when the code fails to meet expectations, by examining the results of intermediate computations. In Haskell, there is no "where and when".

    (If I could stomach the verbose syntax I might get around to trying eager language O'Caml.)

    I have used OCaml a bit. It implements most of the expected functional programming toolbox, e.g. recursion, closures, generic functions, partial function application, Hindley-Milner type inference. etc. It also has fairly direct syntax and semantics for mutable data and object-orientated programming. I generally find that OCaml's runtime behaviour in terms of speed and memory usage is very good. I would not say that OCaml is especially verbose, except maybe for the fact that different operators have to be used for different numeric types, e.g. (+) for integers and (.+) for floating point.

    And by the way, OCaml is primarily French, not Irish. There is no apostrophe after the 'O' ;-)

    Comment Source:David Tweed wrote: > Anyway, I think I ought to say that I think Jim's [Glyn's] discussion of dynamic vs static typing missed that there's (at least) three levels: ... Thanks David, that is a much clearer explanation than mine. > This means that lazy functional languages often have much more unpredictable time and memory behaviour ... This fact pains me a good deal, actually. I believe in the value of code analysis, which for me means figuring out what the code actually does, versus intentions and expectations. I am used to using code analysis for debugging in imperative languages such as C, where it is usually possible to track down where and when the code fails to meet expectations, by examining the results of intermediate computations. In Haskell, there is no "where and when". > (If I could stomach the verbose syntax I might get around to trying eager language O'Caml.) I have used OCaml a bit. It implements most of the expected functional programming toolbox, e.g. recursion, closures, generic functions, partial function application, Hindley-Milner type inference. etc. It also has fairly direct syntax and semantics for mutable data and object-orientated programming. I generally find that OCaml's runtime behaviour in terms of speed and memory usage is very good. I would not say that OCaml is especially verbose, except maybe for the fact that different operators have to be used for different numeric types, e.g. (+) for integers and (.+) for floating point. And by the way, OCaml is primarily French, not Irish. There is no apostrophe after the 'O' ;-)
Sign In or Register to comment.