Hi,

this is my first post in your forum and to begin I wish to compliment you for your excellent work in Numerical Recipes 3rd edition.

I have two issues/questions concerning your code. I am trying to obtain the correlation function \xi\left(r\right) out of a model power spectrum (i.e. an array P\left(k\right) of values between kmin and kmax).
The relevant (cosmological) equation is

\xi (r)=\frac{1}{2\pi ^2}\int_{0}^{\infty}k^2 P(k)\cdot \frac{sin(kr)}{kr}\cdot dk.

I tried the dftint approach to compute the integral between kmin and kmax with some success, but I could not find any reference for the error estimation. How is it related to MPOL (the degree of the interpolationg polynomials for the DFT), M (the number of subintervals in which the integration range is divided) and NDFT?

Secondly, I tried to use the routine qromb in chapter 4 to integrate a function defined by a functor but it did not work. My functor involved calling a method from an interp_1d object, here it is the code (by the way, if you have a more efficient way to accomplish this, please let me know):


struct Interpolation_Wrapper {

Spline_interp* my_func;
double w;
Interpolation_Wrapper(VecDoub xx, VecDoub yy) {
static Spline_interp my_funcc(xx,yy);
my_func=&my_funcc;
}
double operator()(const double x) { return my_func->interp(x) * sin(w*x); }
};

Int N=10;
VecDoub xx(N), yy(N);
for(int i=0; i<N; i++) {
xx[i]=i;
yy[i]=i*i;
}

Interpolation_Wrapper my_wrapper(xx,yy);

cout << qromb(my_wrapper, 0, 1) << endl;


All I got was the following errors at compilation time:

prova_functor.c++: In function ‘int main()’:
prova_functor.c++:44: error: ‘main()::Interpolation_Wrapper’ uses local type ‘main()::Interpolation_Wrapper’
prova_functor.c++:44: error: trying to instantiate ‘template<class T> Doub qromb(T&, Doub, Doub, Doub)’
prova_functor.c++:44: error: no matching function for call to ‘qromb(main()::Interpolation_Wrapper&, Doub&, Doub&, Doub&)’

My first thought was that I was getting something wrong in the functor's definition. Hence, I tried with something simpler:


struct Square {
Doub operator()(const Doub x) {
return x*x;
}
};

Square g;

cout << qromb(g, 0, 1) << endl;


I nevertheless got the error:
prova_functor.c++: In function ‘int main()’:
prova_functor.c++:44: error: ‘main()::Square’ uses local type ‘main()::Square’
prova_functor.c++:44: error: trying to instantiate ‘template<class T> Doub qromb(T&, Doub, Doub, Doub)’
prova_functor.c++:44: error: no matching function for call to ‘qromb(main()::Square&, Doub&, Doub&, Doub&)’

In what am I wrong?

Thank you for your attention,

Guido

prova_functor.c++: In function ‘int main()’:
prova_functor.c++:44: error: ‘main()::Square’ uses local type ‘main()::Square’
In what am I wrong?

It looks like you defined the struct inside your main() function.

Try this:

//
//Whatever includes you need. I use the following
//
#include "../code/nr3.h"
#include "../code/interp_1d.h"
#include "../code/quadrature.h"
#include "../code/romberg.h"

struct Square {
Doub operator()(const Doub x) {
return x*x;
}
};

int main()
{
Square g;
cout << "qromb(g, 0, 1) = " << qromb(g, 0, 1) << endl;
return 0;
}


My output (with GNU g++ version 4.1.2)


qromb(g, 0, 1) = 0.333333


Regards,

Dave

Thank you very much Dave, that was indeed the problem. Is there any specific reason for which gcc won't compile if I declare the struct inside the main function?

By the way, I manage to run the integration using the Square struct, but the program hangs out when evaluating qromb(my_wrapper, 0, 1). I know the wrapper functor is working properly, since a


my_wrapper.w = 1;
cout << my_wrapper(3) << endl;


gives 3^2 * sin (3) = 1.2701 as an output. Did I do something wrong in the declaration of the struct? This is what I have used:


struct Interpolation_Wrapper {

Spline_interp* my_func;
double w;
Interpolation_Wrapper(VecDoub xx, VecDoub yy) {
static Spline_interp my_funcc(xx,yy);
my_func=&my_funcc;
}
double operator()(const double x) { return my_func->interp(x) * sin(w*x); }
};


Thank you again for your help.

Regards,

Guido

...Did I do something wrong in the declaration of the struct?...

Well, before feeding stuff to an integration routine, I would make sure that the interpolation was working by looking at more than one point.

Are you creating the spline approximation based on data values as x goes from 1 to 10 and then using that in some integration routine as x goes from zero to 1? That doesn't make much sense to me, but that's what it sounds like.

Anyhow, if there is a question about the functor (or anything else), I would test just that part.


//
// Use a spline interpolation of points from x^2
// to approximate a function x^2 * sin(omega*x)
//
// davekw7x
//

#include "../code/nr3.h"
#include "../code/interp_1d.h"

using namespace std;

struct Spline_wrapper {

Spline_interp* my_func;
Doub w;
Spline_wrapper(VecDoub xx, VecDoub yy) {
static Spline_interp my_funcc(xx,yy);
my_func = &my_funcc;
}
double operator()(const double x) {
return my_func->interp(x) * sin(w*x);
}
};

int main()
{
Int i;
Int nsamps = 11;
Doub omega = 1.0;
VecDoub xx(nsamps);
VecDoub yy(nsamps);

//
// Generate array for interpolation of y = x*x
//
// Values of x go from xmin through xmax
//
Doub xmin = 1.0;
Doub xmax = 10.0;
Doub deltax = (xmax - xmin)/(nsamps - 1.0);
Doub x;

for (i = 0, x = xmin ; i < nsamps; i++, x+=deltax) {
xx[i] = x;
yy[i] = x*x;
}

Spline_wrapper spline(xx,yy);
spline.w = omega;
cout << "Compare interpolated with actual at data points:"
<< endl<< endl;

cout << setw(6) << "x"
<< setw(23) << "interpolated"
<< setw(14) << "actual" << endl;

cout << scientific;
for (Int i=0; i < nsamps; i++) {
cout << setw(12) << xx[i]
<< setw(17) << spline(xx[i])
<< setw(17) << yy[i]*sin(omega * xx[i])
<< endl;
}
cout << endl << endl;

// Now, compare interpolated values for other values of
// x.
// If values of x are outside the range that was
// used to create the interpolation polynomial,
// well...
//
cout << "Compare interpolated values between data points:"
<< endl << endl;

cout << setw(6) << "x"
<< setw(23) << "interpolated"
<< setw(14) << "actual"
<< setw(16) << "error"
<< endl;

Doub yactual;
Doub yinterp;
Doub yerror;

Int npoints = 11;
xmin = 0.0;
xmax = 1.0;
deltax = (xmax-xmin) / (npoints - 1.0);
for (i = 0, x = xmin; i < npoints; i++, x+=deltax) {
yactual = x*x*sin(omega*x);
yinterp = spline(x);
yerror = yinterp - yactual;
cout << setw(12) << x
<< setw(17) << yactual
<< setw(17) << yinterp
<< setw(17) << yerror
<< endl;
}

return 0;
}


Output:


Compare interpolated with actual at data points:

x interpolated actual
1.000000e+00 1.110245e+00 8.414710e-01
1.900000e+00 3.416143e+00 3.416143e+00
2.800000e+00 2.626307e+00 2.626307e+00
3.700000e+00 -7.253457e+00 -7.253457e+00
4.600000e+00 -2.102650e+01 -2.102650e+01
5.500000e+00 -2.134259e+01 -2.134259e+01
6.400000e+00 4.773855e+00 4.773855e+00
7.300000e+00 4.531977e+01 4.531977e+01
8.200000e+00 6.325472e+01 6.325472e+01
9.100000e+00 2.642454e+01 2.642454e+01
1.000000e+01 -5.440211e+01 -5.440211e+01


Compare interpolated values between data points:

x interpolated actual error
0.000000e+00 0.000000e+00 -0.000000e+00 -0.000000e+00
1.000000e-01 9.983342e-04 1.097355e-02 9.975213e-03
2.000000e-01 7.946773e-03 4.394006e-02 3.599329e-02
3.000000e-01 2.659682e-02 9.902720e-02 7.243038e-02
4.000000e-01 6.230693e-02 1.764148e-01 1.141078e-01
5.000000e-01 1.198564e-01 2.762866e-01 1.564302e-01
6.000000e-01 2.032713e-01 3.987650e-01 1.954937e-01
7.000000e-01 3.156667e-01 5.438310e-01 2.281643e-01
8.000000e-01 4.591079e-01 7.112318e-01 2.521239e-01
9.000000e-01 6.344948e-01 9.003795e-01 2.658847e-01
1.000000e+00 8.414710e-01 1.110245e+00 2.687740e-01

Of course, if that is not what you are doing, then change xmin and xmax to whatever your range of integration is going to be. Change npoints to a larger number and see what happens to the error terms between the data points. Etc., etc.


Regards,

Dave

Footnote: Before feeding interpolation results to an integration routine (or anything else, for that matter), I would probably plot the interpolated values for a large number of points over the desired range. Of course, I would only have some data points, not the actual underlying function to compare, but I would, at the very least, "eyeball" the results to see whether the interpolation seems to be reasonable.

Then, and only then, could I justify further investigation to decide whether this gives reasonable results in the integration routine (or whatever).

Hi Dave,

thank you for your answer.

My (now solved) issue was probably due to a gcc or nr3 bug.
Have you tried to input in your code
Int nsamps = 10
instead of Int nsamps = 11?

I did, and the output I obtained with gcc (version 4.0.1, 4.2, 4.4.0) is a series of -inf:

Compare interpolated with actual at data points:

x interpolated actual
1.000000e+00 5.918249e-01 8.414710e-01
2.000000e+00 -inf 3.637190e+00
3.000000e+00 -inf 1.270080e+00
4.000000e+00 inf -1.210884e+01
5.000000e+00 inf -2.397311e+01
6.000000e+00 inf -1.005896e+01
7.000000e+00 -inf 3.219234e+01
8.000000e+00 -inf 6.331893e+01
9.000000e+00 -inf 3.338160e+01
1.000000e+01 inf -5.440211e+01


Compare interpolated values between data points:

x interpolated actual error
0.000000e+00 0.000000e+00 nan nan
1.000000e-01 9.983342e-04 -inf -inf
2.000000e-01 7.946773e-03 -inf -inf
3.000000e-01 2.659682e-02 -inf -inf
4.000000e-01 6.230693e-02 -inf -inf
5.000000e-01 1.198564e-01 -inf -inf
6.000000e-01 2.032713e-01 -inf -inf
7.000000e-01 3.156667e-01 -inf -inf
8.000000e-01 4.591079e-01 -inf -inf
9.000000e-01 6.344948e-01 -inf -inf
1.000000e+00 8.414710e-01 -inf -inf


The same with nsamps=20, 30 but curiously nsamps = 40, 50, 60... works. Using the Intel C++ compiler the program works perfecty with any nsamps value I tried.
The line qromb(my_wrapper, a, b) was freezing my computer independently of "a" and "b" because of the -inf's.

However, my main concern regards the functor: do you think I adopted the correct way to implement the Spline_interp::interp function in the functor? Using the static declaration seemed to me a somewhat dirty move :)

Thank you very much for being so helpful.

Regards,

Guido

P.S. Does somebody kindly have any suggestion to my first issue concerning the error estimate of the dftint routine?

...
Have you tried to input in your code
Int nsamps = 10
instead of Int nsamps = 11?

Output with g++ version 4.1.2 20071124 (Red Hat 4.1.2-42) on my Centos 5.2 workstation:

Compare interpolated with actual at data points:

x interpolated actual
1.000000e+00 1.110245e+00 8.414710e-01
1.900000e+00 3.416143e+00 3.416143e+00
2.800000e+00 2.626307e+00 2.626307e+00
3.700000e+00 -7.253457e+00 -7.253457e+00
4.600000e+00 -2.102650e+01 -2.102650e+01
5.500000e+00 -2.134259e+01 -2.134259e+01
6.400000e+00 4.773855e+00 4.773855e+00
7.300000e+00 4.531977e+01 4.531977e+01
8.200000e+00 6.325472e+01 6.325472e+01
9.100000e+00 2.642454e+01 2.642454e+01
1.000000e+01 -5.440211e+01 -5.440211e+01


Compare interpolated values between data points:

x interpolated actual error
0.000000e+00 0.000000e+00 -0.000000e+00 -0.000000e+00
1.111111e-01 1.368921e-03 1.354872e-02 1.217980e-02
2.222222e-01 1.088384e-02 5.426420e-02 4.338036e-02
3.333333e-01 3.635497e-02 1.223357e-01 8.598071e-02
4.444444e-01 8.492965e-02 2.180169e-01 1.330872e-01
5.555556e-01 1.627825e-01 3.415361e-01 1.787536e-01
6.666667e-01 2.748310e-01 4.929757e-01 2.181446e-01
7.777778e-01 4.244839e-01 6.721254e-01 2.476415e-01
8.888889e-01 6.134297e-01 8.783168e-01 2.648871e-01
1.000000e+00 8.414710e-01 1.110245e+00 2.687740e-01


Now that I actually look at things, I see the data under the "interpolated" and "acual" columns in the second table are swapped. Oops.

The corrected printout would be

Compare interpolated with actual at data points:

x interpolated actual
1.000000e+00 1.110245e+00 8.414710e-01
1.900000e+00 3.416143e+00 3.416143e+00
2.800000e+00 2.626307e+00 2.626307e+00
3.700000e+00 -7.253457e+00 -7.253457e+00
4.600000e+00 -2.102650e+01 -2.102650e+01
5.500000e+00 -2.134259e+01 -2.134259e+01
6.400000e+00 4.773855e+00 4.773855e+00
7.300000e+00 4.531977e+01 4.531977e+01
8.200000e+00 6.325472e+01 6.325472e+01
9.100000e+00 2.642454e+01 2.642454e+01
1.000000e+01 -5.440211e+01 -5.440211e+01


Compare interpolated values between data points:

x interpolated actual error
0.000000e+00 -0.000000e+00 0.000000e+00 -0.000000e+00
1.111111e-01 1.354872e-02 1.368921e-03 1.217980e-02
2.222222e-01 5.426420e-02 1.088384e-02 4.338036e-02
3.333333e-01 1.223357e-01 3.635497e-02 8.598071e-02
4.444444e-01 2.180169e-01 8.492965e-02 1.330872e-01
5.555556e-01 3.415361e-01 1.627825e-01 1.787536e-01
6.666667e-01 4.929757e-01 2.748310e-01 2.181446e-01
7.777778e-01 6.721254e-01 4.244839e-01 2.476415e-01
8.888889e-01 8.783168e-01 6.134297e-01 2.648871e-01
1.000000e+00 1.110245e+00 8.414710e-01 2.687740e-01


At least the two tables are consistent. My point is still that using a spline function outside the interval over which the interpolation was performed is a Bad Thing, and I didn't know whether you were actually doing that.

Anyhow, now we see how valuable it is to have someone else compare results running the exact same program.

The results from compiling with g++ version 3.4.6 20060404 (Red Hat 3.4.6-4) were identical.

My interp_1d.h and nr3.h files are directly from the Numerical Recipes Version 3 CD and were not modified by me in any way.


The line qromb(my_wrapper, a, b) was freezing my computer independently of "a" and "b" because of the -inf's.That's why I like to show the actual program that I ran rather than saying something like, "Well, I did something like that one time and it worked." It's also why I like to test lower functions, like the interpolation, before getting to the "good stuff." (Top-down design, bottom-up implementation.)

Using the static declaration seemed to me a somewhat dirty move.Well you can't pass the address of a non-static class member function to Spline_interp() since the non-static class member function has an extra (implicit-but-invisible-to-the programmer) parameter: (this). The result is that the signature of a non-static class member function is incompatible with that of the function required by Spline_interp.

The address of a non-static class function can't be cast to the address of a "regular" function. Period.

An alternative to the use of a static function might be not to have a class at all. You would have a "regular" function with the variable w having file scope (a "global" variable). Then you can initialize the global in the source code or, in main(), you can set value of w before using the function.

Based on emotion rather than hard, cold logic, I will say that I personally consider use of a static function less onerous than using a global variable, but I wouldn't argue with anyone about it. Chacun Ã* son goût

Regards,

Dave

Hi Dave,

thank you for the useful informations.

By the way, do you have any idea on how to estimate the error in the dftint routine?

Cheers,

Guido