> Right. That's why I'm saying we need to go back and rethink what we've 
> already done! What we really need to do to nail this down is trace 
> through the uses implied by the functionals in agrad/autodiff.hpp I've 
> convinced myself (and I think Daniel) that we'll never have a case 
> where fvar<T1> and fvar<T2> go into an operator with different T1 and 
> T2 types. The proof is by induction. The function call at the top 
> level has all the same agrad types (var, fvar<double>, 
> fvar<fvar<double>>, fvar<var>, or whatever). Then no operation mixes 
> anything other than scalars in, which have obvious promotions because 
> they have zero derivatives. If so, I'd like to simplify all of the 
> operations in fvar to reflect this. We can't just do the above because 
> there's an ambiguity in promotion of fvars. Consider this: foo(const 
> fvar<fvar<T>>& a, const fvar<T>& b) If we write them in dual notation, 
> then a is going to look like this: a = ((a.val, a.dx1), (a.dx2, 
> a.dx1.dx2)) Depending on where it came from, b could be either
> (b.val, b.dx1) --> ((b.val, b.dx1), (?, ?)) or
> (b.val, b.dx2) --> ((b.val, ?), (b.dx2, ?)) But we don't have enough 
> information to fill in the ? from the values on the left-hand-side. 
> This is different for constants. If we have a double x, then it's 
> trivially: (x,0) as fvar<double> ((x,0), (0, 0)) as fvar<fvar<double>> 
> etc. The base case we need to handle explicitly by implementing a 
> version for double and for int. I think we're going to need an 
> is_arithmetic metaprogram to make sure the type getting promoted isn't 
> itself an auto-diff type. Sorry this is getting so complex.
I was able to get the idea of the problem here, although not the 
solution. It's a little harder to test and understand fwd diff as it's 
not yet used in Stan. Basically you something like this:

     template <typename T>
     inline
     fvar<typename stan::return_type<T>::type>
     log_diff_exp(const fvar<T>& x1, const fvar<T>& x2) {
       using stan::math::log_diff_exp;
       using stan::math::NOT_A_NUMBER;
       using std::exp;
       if (x1.val_ <= x2.val_)
         return fvar<T>(NOT_A_NUMBER, NOT_A_NUMBER);
       return fvar<T>(log_diff_exp(x1.val_, x2.val_),
                      x1.d_ / (1 - exp(x2.val_ - x1.val_) ) + x2.d_ / (1 
- exp(x1.val_ - x2.val_) ) );
     }

???
  (not tested)


> Great. Could you coordinate with Peter Li before jumping in on this? I 
> just suggested some fixes for Peter Li (on the dev list) for the 
> matrix ops, which really need to be improved for efficiency and 
> numerical stability. - Bob
I was going to open a discussion at stan-dev before proceeding, I just 
already done these changes so I wouldn't forget them later.