#19 Dependency-free implementation of PDL::Fit::Polynomial

closed
None
5
2010-07-29
2009-11-03
David Mertens
No

The package PDL::Fit::Polynomial depends on Slatec, but that's silly because Slatec has its own polynomial fitter called polyfit. This means that I can't fit my data to a line unless I use PDL::Fit::Linfit, which does linear fitting for arbitrary functions.

This is particularly frustrating since a polynomial fitter is not too hard to implement by hand. In the very least, the documentation should point the user to PDL::Fit::Linfit as a Slatec-free option.

Discussion

  • David Mertens
    David Mertens
    2009-11-03

    Arrgh, even PDL::Fit::Linfit depends on Slatec! This means I can't do something as simple as fitting data to a line without a FORTRAN compiler! Unless GSL has some fitting routines, but once again, I'll have to use an external dependency. Grrrr.

     
  • Chris Marshall
    Chris Marshall
    2009-11-03

    We could use an f2c version of the slatec fortran code in the PDL tree.
    As for GSL, it is open source, modular, and C based, so there would be
    no problems including an internal subset of the library in PDL (or even all).

    --Chris

     
  • Chris Marshall
    Chris Marshall
    2009-11-03

    A clean up of the dependencies and functionality in PDL by module would be great.
    It would be nice to indicate what is provided where and when there are multiple alternatives.

    One example is all the FFT stuff in PDL: inconsistent calling sequences, undocumented
    implementation information, no not inplace FFT in the base PDL, no support for complex data with adjacent Re and Im parts,...

     
  • David Mertens
    David Mertens
    2010-02-08

    Good news: I have implemented a PDL::PP version of the Householder diagonalization scheme. I've started using it for my own research. The module still needs a lot of work, including documentation, testing, and a nicer wrapper, but at the moment I have a working PDL::PP function that will do linear fitting (for arbitrary functions).

     
  • David Mertens
    David Mertens
    2010-02-08

    My first very rough cut. It uses Inline::Pdlpp. I have looked into using Module::Build, but I can't figure it out, so my next step is to rewrite this using Extutils::MakeMaker (using one of the pp util scripts) and go from there.

     
    Attachments
  • Chris Marshall
    Chris Marshall
    2010-02-20

    The current PDL source tree has an f2c conversion of the SLATEC
    fortran sources. All that is needed is to get that to compile instead
    of the fortran version and a "fortran free" SLATEC would be available.

     
  • David Mertens
    David Mertens
    2010-02-20

    I posted on CPAN but forgot to post here that I've created the module PDL::Fit::Householder. It does general linear fitting using Householder transformations under the hood. I really like the interface (I actually wrote it from scratch and never compared it with PDL::Fit::Polynomial) since it handles general functions, polynomials, and hand-crafted piddles all in the same function call. I'd like develop this package into a general PDL::Fit package that interfaces a whole bunch of different methods under-the-hood, based on some heuristics or based on user specification. Anyway, if anybody wants to see my proposed work, it's on CPAN.

     
  • David Mertens
    David Mertens
    2010-02-20

    I noticed the f2c of the Slatec source and played around with it a little. Though I wasn't able to get it to work, I would like to see this happen. I think it would be nice to have PDL::Fit::Polynomial (or a more general PDL::Fit::Linear) interface with any number of diagonalization schemes, be it Slatec, GSL, or hand-written PDL. The advantage of hand-written PDL stuff is that we can handle bad values in ways that may be more difficult using Slatec or GSL.

     
  • Chris Marshall
    Chris Marshall
    2010-02-20

    The f2c conversion was vanilla only and there may be some runtime IO issues and so
    forth to be handled. Thanks for taking a look.

    I think the idea of a general PDL::Xxxx namespace that can support different
    implementations (e.g., Term::Readline::Gnu, Term::Readline::Perl or our own
    PDL::Image2d...) is a good one. One thing I would like to see is a PDL::FFT that
    supports FFT, FFTW2 and FFTW3 under the hood.

    As for hand-rolled PDL bad value stuff, it seems like it should be possible to
    implement a general bad value strategy for elementwise functions where
    the only requirement is an existing element operation (provided by a library
    such as GSL, Slatec,...) and a choice of bad value support method. How
    many common ones are there for elementwise operations? Binary ops?

     
  • I wrote PDL::Fit::Polynomial, iI wanted it to be as pure PDL as possible to avoid dependencies. There are no for() loops, everything is achieved with slice constructs and it threads nicely over the external arguments. However it DOES call the matrix inversion part of SLATEC. Maybe someone could write a pure PDL matrix inversion routine? This would be of general use.

     
  • David Mertens
    David Mertens
    2010-02-21

    PDL::Fit::Householder implements a pure PDL::PP matrix inversion routine. :)

     
  • David Mertens
    David Mertens
    2010-02-21

    I'm sorry, that wasn't a very useful comment, nor entirely correct. Strictly speaking, PDL::Fit::Householder doesn't actually invert a matrix. Rather, in $A x $coefs = $y, it converts $A (a not necessarily square matrix) to an upper-triangular matrix and properly modifies $y so that you can use back-substitution on $A to solve for $coefs.

     
  • Chris Marshall
    Chris Marshall
    2010-07-25

    As far as I can tell the inv() routine in PDL::MatrixOps
    is an all PDL routine for matrix inversion. Is it sufficient
    for this processing?

     
  • Chris Marshall
    Chris Marshall
    2010-07-27

    This has been fixed in PDL git by replacing the use of
    the PDL::Slatec matinv() routine by the PDL::MatrixOps
    routine inv() for matrix inversion. This removes the
    dependency on Slatec for this routine. The ticket is
    marked Pending and will be closed in 2weeks if no
    further updates are made.

     
  • Chris Marshall
    Chris Marshall
    2010-07-27

    • assigned_to: nobody --> marshallch
    • status: open --> pending
     
  • David Mertens
    David Mertens
    2010-07-27

    This won't work because inv() only works with square matrices, but for fitting you need to diagonalize (actually upper-triangularize) rectangular matrices. Unfortunately I have a routine that does this, but I don't think it threads properly and I'm pretty sure it's inefficient. It's in my Householder module on CPAN. It needs a lot of cleaning up before it would be a candidate for inclusion in PDL. Basically, I want to make my Householder routine to return something that an LU back-substitution solver can use, and then I'll propose it for inclusion.

     
  • David Mertens
    David Mertens
    2010-07-27

    • status: pending --> open
     
  • Chris Marshall
    Chris Marshall
    2010-07-27

    Please add a test to poly.t that shows this.
    Matrix inverse is not defined for rectangular
    matrices (at least in the normal sense).

    Just because inv() doesn't work, maybe the
    lu_decomp and lu_backsub could be used
    directly. Hard to tell without a test case.

    --Chris

     
  • Chris Marshall
    Chris Marshall
    2010-07-27

    I cooked up a few more test cases in the pdl2 shell.
    No errors with using inv() instead of matinv(). A
    cursory look at the code doesn't show any way
    that the matrix being inverted would not be square.

     
  • Chris Marshall
    Chris Marshall
    2010-07-29

    Marking this closed again. If there is a problem with the new
    implementation, please open a bug tracker ticket for that.
    This implementation does resolve the feature request. Thanks.

     
  • Chris Marshall
    Chris Marshall
    2010-07-29

    • status: open --> closed