How can I use the Bessel function in the residual function?

159 views
Skip to first unread message

Alessandro Gentilini

unread,
Sep 21, 2015, 3:22:01 PM9/21/15
to Ceres Solver
Hello,
I am following this sample (with Ceres Solver version 1.10.0 on Linux):


but instead of fitting the exp function I would like to fit the Bessel functions of the first kind with integer order equal to 1.

My first tentative was to just replace 'exp' with 'j1' but I got the following compiler error:

/home/ag/projects/ceres-hello-world/curve_fitting.cc: In instantiation of ‘bool MyResidual::operator()(const T*, const T*, T*) const [with T = ceres::Jet<double, 2>]’:
/usr/local/include/ceres/internal/variadic_evaluate.h:167:26:   required from ‘static bool ceres::internal::VariadicEvaluate<Functor, T, N0, N1, 0, 0, 0, 0, 0, 0, 0, 0>::Call(const Functor&, const T* const*, T*) [with Functor = MyResidual; T = ceres::Jet<double, 2>; int N0 = 1; int N1 = 1]’
/usr/local/include/ceres/internal/autodiff.h:283:45:   required from ‘static bool ceres::internal::AutoDiff<Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>::Differentiate(const Functor&, const T* const*, int, T*, T**) [with Functor = MyResidual; T = double; int N0 = 1; int N1 = 1; int N2 = 0; int N3 = 0; int N4 = 0; int N5 = 0; int N6 = 0; int N7 = 0; int N8 = 0; int N9 = 0]’
/usr/local/include/ceres/autodiff_cost_function.h:218:25:   required from ‘bool ceres::AutoDiffCostFunction<CostFunctor, kNumResiduals, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>::Evaluate(const double* const*, double*, double**) const [with CostFunctor = MyResidual; int kNumResiduals = 1; int N0 = 1; int N1 = 1; int N2 = 0; int N3 = 0; int N4 = 0; int N5 = 0; int N6 = 0; int N7 = 0; int N8 = 0; int N9 = 0]’
/home/ag/projects/ceres-hello-world/curve_fitting.cc:146:1:   required from here
/home/ag/projects/ceres-hello-world/curve_fitting.cc:111:49: error: cannot convert ‘ceres::Jet<double, 2>’ to ‘double’ for argument ‘1’ to ‘double j1(double)’
     residual[0] = T(y_) - j1(m[0] * T(x_) + c[0]);
                                                 ^
make[2]: *** [CMakeFiles/curve_fitting.dir/curve_fitting.cc.o] Error 1
make[1]: *** [CMakeFiles/curve_fitting.dir/all] Error 2
make: *** [all] Error 2



I then debug the original sample (with 'exp' and I realized that 'exp' is not the C++ library function but it is defined inside the namespace ceres in this file: 

/home/ag/projects/ceres-solver-1.10.0/include/ceres/jet.h

And, as you know, 'exp' is defined as:

// exp(a + h) ~= exp(a) + exp(a) h
template <typename T, int N> inline
Jet<T, N> exp(const Jet<T, N>& f) {
  const T tmp = exp(f.a);
  return Jet<T, N>(tmp, tmp * f.v);
}

So I tried to add  in my code a similar definition for j1:

namespace ceres {
// j0 is the Bessel functions of the first kind with integer order equal to 0
// j1 is the Bessel functions of the first kind with integer order equal to 1
// jn is the Bessel functions of the first kind with integer order equal to n

// In general, f(a + h) ~= f(a) + f'(a) h, via the chain rule.
// j1(a + h) ~= j1(a) + 0.5 ( j0(a) - j2(a) ) h
template <typename T, int N> inline
Jet<T, N> j1(const Jet<T, N>& f) {
  return Jet<T, N>(j1(f.a), T(0.5) * ( j0(f.a)-jn(2,f.a) ) * f.v);
}
}

and now I have the following compiler error (see attached code):

/home/ag/projects/ceres-hello-world/curve_fitting.cc: In instantiation of ‘ceres::Jet<T, N> ceres::j1(const ceres::Jet<T, N>&) [with T = double; int N = 2]’:
/home/ag/projects/ceres-hello-world/curve_fitting.cc:111:49:   required from ‘bool MyResidual::operator()(const T*, const T*, T*) const [with T = ceres::Jet<double, 2>]’
/usr/local/include/ceres/internal/variadic_evaluate.h:167:26:   required from ‘static bool ceres::internal::VariadicEvaluate<Functor, T, N0, N1, 0, 0, 0, 0, 0, 0, 0, 0>::Call(const Functor&, const T* const*, T*) [with Functor = MyResidual; T = ceres::Jet<double, 2>; int N0 = 1; int N1 = 1]’
/usr/local/include/ceres/internal/autodiff.h:283:45:   required from ‘static bool ceres::internal::AutoDiff<Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>::Differentiate(const Functor&, const T* const*, int, T*, T**) [with Functor = MyResidual; T = double; int N0 = 1; int N1 = 1; int N2 = 0; int N3 = 0; int N4 = 0; int N5 = 0; int N6 = 0; int N7 = 0; int N8 = 0; int N9 = 0]’
/usr/local/include/ceres/autodiff_cost_function.h:218:25:   required from ‘bool ceres::AutoDiffCostFunction<CostFunctor, kNumResiduals, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>::Evaluate(const double* const*, double*, double**) const [with CostFunctor = MyResidual; int kNumResiduals = 1; int N0 = 1; int N1 = 1; int N2 = 0; int N3 = 0; int N4 = 0; int N5 = 0; int N6 = 0; int N7 = 0; int N8 = 0; int N9 = 0]’
/home/ag/projects/ceres-hello-world/curve_fitting.cc:146:1:   required from here
/home/ag/projects/ceres-hello-world/curve_fitting.cc:101:26: error: no matching function for call to ‘j1(const double&)’
   return Jet<T, N>(j1(f.a), T(0.5) * ( j0(f.a)-jn(2,f.a) ) * f.v);
                          ^
/home/ag/projects/ceres-hello-world/curve_fitting.cc:101:26: note: candidate is:
/home/ag/projects/ceres-hello-world/curve_fitting.cc:100:11: note: template<class T, int N> ceres::Jet<T, N> ceres::j1(const ceres::Jet<T, N>&)
 Jet<T, N> j1(const Jet<T, N>& f) {
           ^
/home/ag/projects/ceres-hello-world/curve_fitting.cc:100:11: note:   template argument deduction/substitution failed:
/home/ag/projects/ceres-hello-world/curve_fitting.cc:101:26: note:   mismatched types ‘const ceres::Jet<T, N>’ and ‘const double’
   return Jet<T, N>(j1(f.a), T(0.5) * ( j0(f.a)-jn(2,f.a) ) * f.v);
                          ^
make[2]: *** [CMakeFiles/curve_fitting.dir/curve_fitting.cc.o] Error 1
make[1]: *** [CMakeFiles/curve_fitting.dir/all] Error 2
make: *** [all] Error 2

What am I missing?

Thank you very much.
Best regards,
Alessandro Gentilini
curve_fitting.cc

Sameer Agarwal

unread,
Sep 21, 2015, 3:34:07 PM9/21/15
to Ceres Solver
Alessandro,
The problem as you mentioned is that we do not have overloads for Bessel functions that work on Jet<T, N>. 

The overload that you are defining is fine,
template <typename T, int N> inline
Jet<T, N> j1(const Jet<T, N>& f) {
  return Jet<T, N>(j1(f.a), T(0.5) * ( j0(f.a)-jn(2,f.a) ) * f.v);
}

except that the compiler seems unable to find the definition of j1, j2 and jn. Where do you expect it to pull them from? 

Sameer



--
You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/86c5da67-b6d2-4c8e-9d30-d9bc0d9cafc1%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Alessandro Gentilini

unread,
Sep 21, 2015, 3:50:19 PM9/21/15
to Ceres Solver


On Monday, September 21, 2015 at 9:34:07 PM UTC+2, Sameer Agarwal wrote:
Alessandro,
The problem as you mentioned is that we do not have overloads for Bessel functions that work on Jet<T, N>. 

The overload that you are defining is fine,
template <typename T, int N> inline
Jet<T, N> j1(const Jet<T, N>& f) {
  return Jet<T, N>(j1(f.a), T(0.5) * ( j0(f.a)-jn(2,f.a) ) * f.v);
}

except that the compiler seems unable to find the definition of j1, j2 and jn. Where do you expect it to pull them from? 

It seems to me that including <cmath> is enough... This snippet works fine (see also https://ideone.com/xqRlTf ):

#include <cmath>
#include <iostream>

int main(int,char**)
{
std::cout << j0(0.) << " " << j1(0.) << " " << jn(2,0.) << "\n";
return 0;
}

-- 
Alessandro

William Rucklidge

unread,
Sep 21, 2015, 4:11:05 PM9/21/15
to ceres-...@googlegroups.com
Maybe you need to provide forwarding definitions of j0, j1, jn like the ones in jet.h:

inline double exp     (double x) { return std::exp(x);      }

-wjr


To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/839365f8-4672-4de7-8e97-e237d3e602db%40googlegroups.com.

Alessandro Gentilini

unread,
Sep 21, 2015, 4:20:59 PM9/21/15
to Ceres Solver


On Monday, September 21, 2015 at 10:11:05 PM UTC+2, wjr wrote:
Maybe you need to provide forwarding definitions of j0, j1, jn like the ones in jet.h:

inline double exp     (double x) { return std::exp(x);      }

Thank you!
With the following snippet my code compiles with success... now I have a Segmentation fault but that's a different problem ;-)

namespace ceres {

inline double j1     (double x) { return j1(x);      }//‘j1’ is not a member of ‘std’

Alessandro Gentilini

unread,
Sep 21, 2015, 4:25:22 PM9/21/15
to Ceres Solver


On Monday, September 21, 2015 at 10:20:59 PM UTC+2, Alessandro Gentilini wrote:


On Monday, September 21, 2015 at 10:11:05 PM UTC+2, wjr wrote:
Maybe you need to provide forwarding definitions of j0, j1, jn like the ones in jet.h:

inline double exp     (double x) { return std::exp(x);      }

Thank you!
With the following snippet my code compiles with success... now I have a Segmentation fault but that's a different problem ;-)

I am afraid that my "solution" is an infinite recursion...

Sameer Agarwal

unread,
Sep 21, 2015, 4:28:23 PM9/21/15
to Ceres Solver
Alessandro, 
Could you start by naming the overload for jets 

Bessel0, Bessel1 and BesselN and then calling j0, j1 and jn inside it. That way you should not need the recursive definition.

Sameer


To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/96e8d6ab-b2d9-423f-834c-a1ce9366a002%40googlegroups.com.

Alessandro Gentilini

unread,
Sep 21, 2015, 4:35:50 PM9/21/15
to Ceres Solver


On Monday, September 21, 2015 at 10:28:23 PM UTC+2, Sameer Agarwal wrote:
Alessandro, 
Could you start by naming the overload for jets 

Bessel0, Bessel1 and BesselN and then calling j0, j1 and jn inside it. That way you should not need the recursive definition.

Sameer

Thank you Sameer, I will have a look tomorrow.

Keir Mierle

unread,
Sep 21, 2015, 8:26:48 PM9/21/15
to ceres-...@googlegroups.com
On Mon, Sep 21, 2015 at 1:20 PM, Alessandro Gentilini <agent...@gmail.com> wrote:


On Monday, September 21, 2015 at 10:11:05 PM UTC+2, wjr wrote:
Maybe you need to provide forwarding definitions of j0, j1, jn like the ones in jet.h:

inline double exp     (double x) { return std::exp(x);      }

Thank you!
With the following snippet my code compiles with success... now I have a Segmentation fault but that's a different problem ;-)

namespace ceres {

inline double j1     (double x) { return j1(x);      }//‘j1’ is not a member of ‘std’

I think you want an explicit std::j1 here

inline double j1(double x) {return std::j1(x);}

You can see the other Ceres functions in jet.h add an explicit std:: to deal with this as well.
 

// j0 is the Bessel functions of the first kind with integer order equal to 0
// j1 is the Bessel functions of the first kind with integer order equal to 1
// jn is the Bessel functions of the first kind with integer order equal to n

// In general, f(a + h) ~= f(a) + f'(a) h, via the chain rule.
// j1(a + h) ~= j1(a) + 0.5 ( j0(a) - j2(a) ) h
template <typename T, int N> inline
Jet<T, N> j1(const Jet<T, N>& f) {
  return Jet<T, N>(j1(f.a), T(0.5) * ( j0(f.a)-jn(2,f.a) ) * f.v);
}
}

You will need to define the other ones as well (j0, jn) in addition to j1.

If you finish this and get it working, I would consider taking such a patch. However, I have to say, the name "j1" is terrible; I've certainly named loop counters with j1 j2 etc in the past.
 

Alessandro Gentilini

unread,
Sep 23, 2015, 4:25:12 PM9/23/15
to Ceres Solver


On Tuesday, September 22, 2015 at 2:26:48 AM UTC+2, Keir Mierle wrote:


On Mon, Sep 21, 2015 at 1:20 PM, Alessandro Gentilini <agent...@gmail.com> wrote:


On Monday, September 21, 2015 at 10:11:05 PM UTC+2, wjr wrote:
Maybe you need to provide forwarding definitions of j0, j1, jn like the ones in jet.h:

inline double exp     (double x) { return std::exp(x);      }

Thank you!
With the following snippet my code compiles with success... now I have a Segmentation fault but that's a different problem ;-)

namespace ceres {

inline double j1     (double x) { return j1(x);      }//‘j1’ is not a member of ‘std’

I think you want an explicit std::j1 here

inline double j1(double x) {return std::j1(x);}

You can see the other Ceres functions in jet.h add an explicit std:: to deal with this as well.

I think j1 is not a member of 'std', at least on my machine where I am using gcc version 4.9.3 (Debian 4.9.3-1).
j0, j1 and jn are not members of 'std' in Windows too.
 
 

// j0 is the Bessel functions of the first kind with integer order equal to 0
// j1 is the Bessel functions of the first kind with integer order equal to 1
// jn is the Bessel functions of the first kind with integer order equal to n

// In general, f(a + h) ~= f(a) + f'(a) h, via the chain rule.
// j1(a + h) ~= j1(a) + 0.5 ( j0(a) - j2(a) ) h
template <typename T, int N> inline
Jet<T, N> j1(const Jet<T, N>& f) {
  return Jet<T, N>(j1(f.a), T(0.5) * ( j0(f.a)-jn(2,f.a) ) * f.v);
}
}

You will need to define the other ones as well (j0, jn) in addition to j1.

If you finish this and get it working, I would consider taking such a patch. However, I have to say, the name "j1" is terrible; I've certainly named loop counters with j1 j2 etc in the past.

I get it working on Debian and on Windows 7 and I will be happy to contribute to Ceres, I submitted the patch https://ceres-solver-review.googlesource.com/7060

On Debian I built it with gcc version 4.9.3 (Debian 4.9.3-1) and on Windows with Visual Studio 2012.

I agree: j0, j1 and jn are not so good as variable names, anyway they are used in gcc and Visual Studio...


Thank you for your help.

-- 
Alessandro Gentilini

Sameer Agarwal

unread,
Sep 23, 2015, 4:46:21 PM9/23/15
to Ceres Solver
Thanks for the patch Alessandro, a review is in your mailbox.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/5a016505-b9a4-49a7-b62a-06f45c8cda3c%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages