Replacing spline from SBProfile Table.cpp

[rmandelb]

Table.cpp in SBProfile uses a modified version of the NR spline. Per Gary's description, "The spline implementation is completely hidden inside the Table<> class so it can be rewritten with no effect on code outside this class." Best option might be to rewrite it ourselves using the rather simple definition of a spline (downside: significant testing required). Otherwise, is there some other open source code that we could use to replace this?

[PaulPrice]

Many of these "pull out NR" issues might be easily resolved by using gsl.

[rmjarvis]

Ah, but then we would need some "pull out gsl" issues... :)

[dkirkby]

Isn't the main issue with NR its restrictive license (which I believe forbids including any NR code in an open source project)? When we discussed wether something like GSL would be useful, it was only in the context of special functions, but now we are also talking about algorithms for splines, 1d integrators, and random numbers. I suppose it is a matter of taste, but we seem to be going to a lot of trouble to avoid using a mainstream open source project that is well documented and tested.

[rmjarvis]

I was mostly joking here. (Hence the smiley.) You're right that the bigger problem with NR is the license. But in general, I think a number of us have a pretty low opinion of both NR and GSL in terms of their algorithm choices and implementations. I don't actually know specifically about the spline routines (for this issue), but the Boost random number generator is much better than the NR one, my integrator is much faster than the simpson's rule one that I replaced, and don't even get me started on either of their linear algebra packages. Basically, in every case where I used a numerical recipes or gsl algorithm, and then looked into the literature, I found much better algorithms were possible. So now I assume that the same is probably true for their other algorithms that I haven't looked into as closely.

[dkirkby]

I guess I am not up on the latest research in 1d integrators, but I have
not had any problems in my limited use of GSL algorithms in code that
wasn't performance critical and it saves me a lot of time to let someone
else take care of the documentation and testing. Is there some part of the
great3 simulation code where integration or splines are likely to be a
bottleneck?

David

On Sat, Feb 25, 2012 at 12:38 PM, rmjarvis <
reply@reply.github.com

wrote:

I was mostly joking here. (Hence the smiley.) You're right that the
bigger problem with NR is the license. But in general, I think a number of
us have a pretty low opinion of both NR and GSL in terms of their algorithm
choices and implementations. I don't actually know specifically about the
spline routines (for this issue), but the Boost random number generator is
much better than the NR one, my integrator is much faster than the
simpson's rule one that I replaced, and don't even get me started on either
of their linear algebra packages. Basically, in every case where I used a
numerical recipes or gsl algorithm, and then looked into the literature, I
found much better algorithms were possible. So now I assume that the same
is probably true for their other algorithms that I haven't looked into as
closely.

---

Reply to this email directly or view it on GitHub:
#4 (comment)

[TallJimbo]

Unless somebody volunteers to do a custom spline implementation, I think adding GSL with the understanding we shouldn't use it unless necessary is probably the right call for now. We really shouldn't have anything related to NR any longer than we have to.

A custom implementation also has the potential to be a lot more user-friendly in C++; we could wrap the GSL C interface in a nice C++ layer too, of course, but for something as simple as splines that may be similar to the amount of work it would take to do a custom C++ implementation.

So even if we do need GSL for now, I'd also eventually like to see us phase it out if we don't run into something harder than splines that we want from it. But getting the NR out should be a higher priority.

[barnabytprowe]

Hi all, I've assigned myself to this issue (for the same reasons as I assigned myself to #3)... Will keep updated with progress.

[gbernstein]

The files include/galsim/Poisson.h and src/Poisson.cpp are not needed for anything in SBProfile. However when I removed them (and the #include "galsim/Poisson.h" in GalSim.h), an attempt to build produces this message:

scons: *** [src/.obj/Poisson.os] Source src/Poisson.cpp' not found, needed by targetsrc/.obj/Poisson.os'.
scons: building terminated because of errors.

I am guessing this is due to the python wrapping of the Poisson methods but don't know where to remove the dependence. Jim, maybe you do?

[TallJimbo]

I think you just need to remove it from src/files.txt

[barnabytprowe]

Mike / @rmjarvis , I was wondering if I could pick your brains? I've pushed a first attempt at a homegrown natural cubic spline algorithm, which uses TMV for the symmetric tridiagonal matrix solution required to find the second derivatives of the interpolant.

I haven't been able to test the output of this yet, as I'm getting these errors from ld:

  "tmv::BandMatrixView<double, 0>::setZero()", referenced from:
      tmv::GenSymBandMatrix<double>::assignToB(tmv::BandMatrixView<double, 0>) constin Table.os
      tmv::GenSymBandMatrix<double>::assignTosB(tmv::SymBandMatrixView<double, 0>) constin Table.os
      tmv::GenBandMatrix<double>::assignToM(tmv::MatrixView<double, 0>) constin Table.os
      void tmv::Copy<double, double>(tmv::GenBandMatrix<double> const&, tmv::BandMatrixView<double, 0>)in Table.os
  "tmv::BandMatrixView<std::complex<double>, 0>::setZero()", referenced from:
      tmv::GenSymBandMatrix<double>::assignToB(tmv::BandMatrixView<std::complex<double>, 0>) constin Table.os
      tmv::GenSymBandMatrix<double>::assignTosB(tmv::SymBandMatrixView<std::complex<double>, 0>) constin Table.os
      tmv::GenBandMatrix<double>::assignToM(tmv::MatrixView<std::complex<double>, 0>) constin Table.os
      tmv::GenBandMatrix<std::complex<double> >::assignToM(tmv::MatrixView<std::complex<double>, 0>) constin Table.os
      tmv::GenSymBandMatrix<std::complex<double> >::assignToB(tmv::BandMatrixView<std::complex<double>, 0>) constin Table.os
      tmv::GenSymBandMatrix<std::complex<double> >::assignTosB(tmv::SymBandMatrixView<std::complex<double>, 0>) constin Table.os
      void tmv::Copy<double, std::complex<double> >(tmv::GenBandMatrix<double> const&, tmv::BandMatrixView<std::complex<double>, 0>)in Table.os
      ...
  "tmv::GenBandMatrix<double>::norm1() const", referenced from:
      tmv::GenBandMatrix<double>::normInf() constin Table.os
      vtable for tmv::GenBandMatrix<double>in Table.os
      vtable for tmv::ConstBandMatrixView<double, 0>in Table.os
      construction vtable for tmv::GenBandMatrix<double>-in-tmv::ConstBandMatrixView<double, 0>in Table.os
      vtable for tmv::BandMatrixView<double, 0>in Table.os
      construction vtable for tmv::GenBandMatrix<double>-in-tmv::BandMatrixView<double, 0>in Table.os
      vtable for tmv::BandMatrix<double, 2>in Table.os
      ...
  "tmv::BandMatrixView<std::complex<double>, 0>::canLinearize() const", referenced from:
      void tmv::DoCopy<double, std::complex<double> >(tmv::GenBandMatrix<double> const&, tmv::BandMatrixView<std::complex<double>, 0>)in Table.os
      void tmv::DoCopy<std::complex<double>, std::complex<double> >(tmv::GenBandMatrix<std::complex<double> > const&, tmv::BandMatrixView<std::complex<double>, 0>)in Table.os
      vtable for tmv::BandMatrixView<std::complex<double>, 0>in Table.os
  "tmv::GenBandMatrix<std::complex<double> >::norm1() const", referenced from:
      tmv::GenBandMatrix<std::complex<double> >::normInf() constin Table.os
      vtable for tmv::GenBandMatrix<std::complex<double> >in Table.os
      vtable for tmv::ConstBandMatrixView<std::complex<double>, 0>in Table.os
      construction vtable for tmv::GenBandMatrix<std::complex<double> >-in-tmv::ConstBandMatrixView<std::complex<double>, 0>in Table.os
      vtable for tmv::BandMatrixView<std::complex<double>, 0>in Table.os
      construction vtable for tmv::GenBandMatrix<std::complex<double> >-in-tmv::BandMatrixView<std::complex<double>, 0>in Table.os
      vtable for tmv::BandMatrix<std::complex<double>, 2>in Table.os
      ...
  "tmv::BandStorageLength(tmv::StorageType, long, long, long, long)", referenced from:
      tmv::BandMatrixView<double, 0>::BandMatrixView(double*, long, long, long, long, long, long, long, tmv::ConjType, long)in Table.os
      tmv::ConstBandMatrixView<double, 0>::ConstBandMatrixView(double const*, long, long, long, long, long, long, long, tmv::ConjType, long)in Table.os
      tmv::SymBandMatrix<double, 4>::SymBandMatrix(tmv::GenSymBandMatrix<double> const&)in Table.os
      tmv::BandMatrix<double, 2>::BandMatrix(tmv::GenBandMatrix<double> const&)in Table.os
      tmv::BandMatrix<double, 0>::BandMatrix(tmv::GenBandMatrix<double> const&)in Table.os
      tmv::BandMatrix<double, 4>::BandMatrix(tmv::GenBandMatrix<double> const&)in Table.os
      tmv::BandMatrixView<std::complex<double>, 0>::BandMatrixView(std::complex<double>*, long, long, long, long, long, long, long, tmv::ConjType, long)in Table.os
      ...

plus many more lines, until...

  "tmv::GenSymMatrix<std::complex<double> >::setDiv() const", referenced from:
      vtable for tmv::GenSymMatrix<std::complex<double> >in Table.os
      vtable for tmv::SymMatrixView<std::complex<double>, 0>in Table.os
      construction vtable for tmv::GenSymMatrix<std::complex<double> >-in-tmv::SymMatrixView<std::complex<double>, 0>in Table.os
  "tmv::GenSymMatrix<std::complex<double> >::cref(long, long) const", referenced from:
      vtable for tmv::GenSymMatrix<std::complex<double> >in Table.os
      vtable for tmv::SymMatrixView<std::complex<double>, 0>in Table.os
      construction vtable for tmv::GenSymMatrix<std::complex<double> >-in-tmv::SymMatrixView<std::complex<double>, 0>in Table.os
  "non-virtual thunk to tmv::GenSymMatrix<std::complex<double> >::setDiv() const", referenced from:
      vtable for tmv::GenSymMatrix<std::complex<double> >in Table.os
      vtable for tmv::SymMatrixView<std::complex<double>, 0>in Table.os
      construction vtable for tmv::GenSymMatrix<std::complex<double> >-in-tmv::SymMatrixView<std::complex<double>, 0>in Table.os
ld: symbol(s) not found for architecture x86_64
collect2: ld returned 1 exit status
scons: *** [lib/libgalsim.0.dylib] Error 1
scons: building terminated because of errors.

Although rather a short few lines of code, this is something of a first foray into proper coding of my own C++ up until now, and so I'm quite prepared for having done something silly. Is there anything in the code above that sticks out at you?

For info, I'm running TMV 0.70 (but installed with SHARED=True). The above doesn't look like any sort of bug to me, just some kind of install/linker problem, but I can certainly try replacing 0.70 with the newer 0.71 if you think it might be of use. I'd be grateful if you have any ideas.

Also, please feel free to comment on the overall code! (That goes for everyone... I need to check it still produces matching output to the old code, but compiling and linking were my first hurdle).

[barnabytprowe]

Very sorry everybody for the multiple messages. Github was telling me there was an error with the message... In the past I'd remedied such problems by altering the amount of marked-down text, and there was quite a bit in the original message. Very sorry if I clogged your inboxes!

[rmjarvis]

We hadn't been linking with -ltmv_symband, since we hadn't been using any symmetric or banded matrix stuff yet. I just switched that in the SConstruct file, so you should be good to go now.

[barnabytprowe]

Ah, thanks Mike, works a charm...

[barnabytprowe]

OK guys, I've completed the switchout and added regression tests. Am going to pull request...

[barnabytprowe]

OK, merged in #155 the branch that fixes this. Closing.