Re: [eigen] calculating on ranges

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

OOps there is a big constness problem, I marked y as const but then I modify it.
However, in order to allow passing a temporary xpr as y, I have to
make it const. So again the solution is a const_cast, yes that's ugly
but the C++ language leaves us no choice if we want to pass a
temporary segment() as y.

template<typename VectorType>
void addGaussianToCurve(
const MatrixBase<VectorType>& x,
const MatrixBase<VectorType>& _y,
typename VectorType::Scalar a,
typename VectorType::Scalar b,
typename VectorType::Scalar c)
{
 VectorType& y = _y.const_cast_derived();
 y += a * ( -(x.cwise() - b).cwise().square() / c ).cwise().exp();
}

Cheers,
Benoit

2009/1/27 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> Hi,
>
> Yes, Eigen's lazy evaluation means that if you implement this with
> Eigen it will also evaluate the function only on the range that's
> needed.
>
> Compared to your C function, I don't know if you would prefer the
> Eigen way or not. Normally, Eigen would bring the advantage of
> vectorization, except that here you use exp and that is not yet
> vectorized (it could be if that were important. we would use Pade
> approximants if there are no SSE instructions for that.). The other
> advantage of using Eigen is that it will automatically unroll the loop
> if the size (here 10) is known at compile-time and the total cost
> isn't too big (you can tune EIGEN_UNROLLING_LIMIT).
>
> Here is a basic Eigen implementation:
>
> #include <Eigen/Array>
>
> template<typename VectorType>
> void addGaussianToCurve(
> const MatrixBase<VectorType>& x,
> const MatrixBase<VectorType>& y,
> typename VectorType::Scalar a,
> typename VectorType::Scalar b,
> typename VectorType::Scalar c)
> {
>  y += a * ( -(x.cwise() - b).cwise().square() / c ).cwise().exp();
> }
>
> Now if you have vectors VectorXd x,y; you can call your function like this:
> addGaussianToCurve(
> x.segment(start,length),
> y.segment(start,length),
> a,b,c);
>
> Now if you know length at compile time you can do:
> addGaussianToCurve(
> x.segment<length>(start),
> y.segment<length>(start),
> a,b,c);
> And it will automatically unroll the loop if appropriate (see above).
>
> Again if vectorization matters it's no big deal, we just need a packet
> version of exp(). You'd then benefit from it automatically.
>
> Cheers,
> Benoit
>