I have an Nx3 tensor C. N is not known at compile time (the N 3-coordinates are being read from a user provided data file). I'd like to create an NxNx3 tensor (call it D) with all the vectorial differences of the coordinates. So, an element of D would look like D(i, j) = C(i) - C(j). While these differences are signed, I only need them going in one direction (it doesn't matter which, just that it be consistent). I'm looking for a way to only evaluate the unique elements without just writing a double for-loop. The documentation I've found for the Eigen TensorSymmetry module does not make it totally apparent how to do this, or whether it is even possible. It is worth saying that I was able to write a numpy-esque difference of differences for one dimension using:
Eigen::array<int, 2> flip({1,0});
Eigen::array<int, 2> bc({1,N});
d = c.broadcast(bc) - c.broadcast(bc).shuffle(flip);
Here c is a N length tensor of rank 1, and d is an NxN size tensor of rank 2.
I've tried various permutations of this for an NxNx3 tensor permuting or expanding the Nx3 input data but can't seem to find anything that looks right. What I'm working with right now is an expression like:
Eigen::array<int, 3> expand({N, N, 3});
Eigen::array<int, 3> trans({1,0,2});
Tensor<double, 2> C(N, 3);
Tensor<double, 3> D(N, N, 3);
// fill up C from file-object like thing foo not relevant to this problem
for (auto i = 0; i < N; i++)
for (auto j = 0; j < 3; j++)
C(i, j) = foo.getCoords()[j];
D = C.broadcast(expand) - C.broadcast(expand).shuffle(trans);
The above will compile, but produces an obviously wrong result (doesn't even have dimensions expected, which is peculiar...). I've tried various permutations of options for expand, for example expand({N,3,1}) and expand({1,N,3}), and also for trans, but there's clearly just something I'm missing here; nothing has worked yet.
As I mentioned above I also don't see how or whether I could exploit the antisymmetry of D to avoid some redundant subtractions.
Any advice would be appreciated.
My somewhat unsatisfactory solution to my own problem has been to write a for-loop. This doesn't take advantage of any symmetries, or allow for lazy evaluation, but it does appear to give the expected results.
Eigen::array<int, 2> bcChip({N, 3});
for (auto i = 0; i < N; i++){
D.chip(i, 0) = C - C.chip(i, 0).eval().broadcast(bcChip);
}
Note that if someone were able to provide a slicker/better implementation of this that did take advantage of Eigen Tensor's features I'd accept their answer instead.