davisad
Registered Member

I am using Eigen in my code along with the automatic differentiation (AD) software Sacado (http://trilinos.sandia.gov/packages/doc ... index.html). I actually had some trouble getting the two pieces of software to work together. Sacado basically overloads scalar types in such a way that when you call operations like (*,+,/,,ect...) the corresponding derivative is also computed. The idea being if you write a function that returns a scalar and call it with a Sacado type the function also computes the derivative with no additional implementation. I am over simplifying this explanation but hopefully you get the idea.
I had some trouble with the "?" operator. In order to get Sacado working through an LU solver I had to change line 188 of the file SelfCwiseBinaryOp.h: // tmp = PlainObject::Constant(rows(),cols(), NumTraits<Scalar>::IsInteger ? other : Scalar(1)/other); tmp = PlainObject::Constant(rows(),cols(), NumTraits<Scalar>::IsInteger ? other : Scalar(1/other)); The reason is that in the case when Scalar is a Sacado type the value returned by expressions are often different from the AD types themselves (I'd have to dig into the Sacado code to figure out why). So you have to change the casting to: isInt(x) ? x : ScalarT(1/x); This fix is pretty simple. It has at least solved my problem. I'm not sure if it comes up in other parts of Eigen though. I figured I'd post this and see if anyone else has had this problem. I know there are a couple other AD tools that are using Eigen. I haven't seen it with Sacado though. 
ggael
Moderator

I think Sacado has some expression template mechanisms, so 1/x is probably a complex type. I'll apply such a fix.
Anyway, note that if your computation actually include a matrix inversion (or solving), then it is much better to use highlevel differentiation formulas than tracking the derivatives at the scalar level through the matrix factorization: dA^1 / dx =  A^1 * dA/dx * A^1 and you can still use AD to compute dA/dx, and LU to perform the A^1 products without explicitly forming A^1 
pjsa
Registered Member

Dear Gael,
I also have experienced the same issue with line 188 of the file SelfCwiseBinaryOp.h. If you can add the fix in your repository soon that would be very helpful. 
ggael
Moderator

Here you go:
https://bitbucket.org/eigen/eigen/commits/227098a74040/ changeset: 227098a74040 user: ggael date: 20130223 20:13:21 summary: Avoid problematic ternary operator (viewtopic.php?f=74&t=109486) 
Registered users: Baidu [Spider], Bing [Bot], davidemme, Exabot [Bot], Google [Bot], La Ninje, Majestic12 [Bot], nezumi, nharmony, paviluf, YaCy [Bot], Yahoo [Bot]