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.
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.
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.
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
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,...
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).
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.
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.
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.
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.
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.
PDL::Fit::Householder implements a pure PDL::PP matrix inversion routine. :)
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.
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?
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.
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.
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
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.
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.
Log in to post a comment.