0

I'm writing up a bunch (~ a dozen) algorithms that iteratively process a vector in the following pattern:

ArrayXd prev;
ArrayXd curr;

prev = some_initial_vector();
for (i = 0; i < N; ++i) {
    // do lots of stuff, many lines to compute curr
    // use lots of algorithm specific variables initialised above
    ...
    prev = curr;
}
return curr;

I would like to have a way of returning the entire history of the values of curr as well, as rows in an ArrayXXd.

I have tried solving this by writing two classes that exhibit a curr handle, one as an ArrayXd & the other as a Block<ArrayXXd, 1, -1>, but that failed as it's not possible to reassign Blocks.

What is a good way of solving this problem? Maybe I could store the Blocks or the ArrayXds themselves in a std::vector and then convert that to an ArrayXXd at the end.

Edit: Added sample input, output

struct NotAccumulator {
    typedef ArrayXd return_type;
    ArrayXd curr;
    ArrayXd prev;
    void record () {}
    return_type result() {
        return prev;
    }
};

struct RowAccumulator {
    typedef ArrayXXd return_type;
    ArrayXd curr;
    ArrayXd prev;
    RowAccumulator(const uint N) {
        history.reserve(N);
    }

    void record () {
        history.push_back(curr);
    }
    return_type result () {
        uint rows = history.size();
        uint cols = history[0].size();
        ArrayXXd result_matrix (rows, cols);
        for(uint i = 0; i < rows; ++i) {
            result_matrix.row(i) = Map<ArrayXd> (history[i].data(), cols);
        }
        return result_matrix;
    }

private:
    std::vector<ArrayXd> history;
};

template <typename Accumulator>
typename Accumulator::return_type add_one(const ArrayXd & start, const uint how_many_times, Accumulator & u) {
    u.prev = start;
    for (uint i = 0; i < how_many_times; ++i) {
        u.curr = 1 + u.prev;
        u.record();
        u.prev = u.curr;
    }
    return u.result();
}

ArrayXd start (3);
start << 1, 0, -1;
NotAccumulator notAccumulator;
RowAccumulator rowAccumulator (5);
cout << add_one(start, 5, notAccumulator) << endl;
// outputs 6 5 4
cout << add_one(start, 5, rowAccumulator) << endl;
// outputs 2 1 0\n 3 2 1\n ... 6 5 4
Community
  • 1
  • 1
user357269
  • 1,835
  • 14
  • 40

1 Answers1

0

If you want to avoid the data copy as much as possible, you could allocate the N-row ArrayXXd in advance and avoid using return.

Buffers like prev and curr are unnecessary too. The following code also assumes the number of iterations N is known in advance.

void GenerateHistory(int N, ArrayXXd* result) {
    result.row(0) = gen_row0_without_prev();
    for(int curr=1; curr<N; curr++) {
        result.row(curr) = gen_row_with_prev_row(result.row(curr-1));
    }
}

You could call this function with

int N = 100;
int m = 20;
ArrayXXd history(N, m);
GenerateHistory(N, &history);

EDIT

Let's put the copying issue aside. Here's an updated code that should fulfill the requirement according to your code sample.

ArrayXXd add_one(const ArrayXd& start, int num_iterations, bool acc) {
    if (acc) {
        ArrayXXd result = start;
        for(int i=0; i<num_iterations; i++) {
            result = result + 1;
        }
        return result;
    } else {
        ArrayXXd result(num_iterations, start.cols());
        result.row(0)=start+1;
        for(int i=1; i<num_iterations; i++) {
            result.row(i)=result.row(i-1);
        }
        return result;
    }
}

ArrayXd start (3);
start << 1, 0, -1;
bool acc = true;
cout << add_one(start, 5, !acc) << endl;
// outputs 6 5 4
cout << add_one(start, 5, acc) << endl;
// outputs 2 1 0\n 3 2 1\n ... 6 5 4
kangshiyin
  • 9,681
  • 1
  • 17
  • 29
  • Thank you for your answer. What you wrote doesn't seem to solve my problem as I would like to be able to return _either_ the whole history (as you do in your answer) _or_ just the last value. I have updated my question, hopefully it's clearer now – user357269 Aug 01 '16 at 21:23
  • Correct me if I'm wrong, but putting `result` as an input doesn't lead to less copying. A similar transformation is done by the compiler under the name _Return Value Optimisation_ automatically – user357269 Aug 01 '16 at 21:26
  • You could use `ArrayXXd` with only 1 row for the NotAcc case. You don't need a separate return type. There will be a copy if you assign the return to a variable like `r=add_one(…)` – kangshiyin Aug 01 '16 at 21:40
  • Instead of the `ArrayXXd*` you could use a `Eigen::Ref`. No copy will occur. – Avi Ginsburg Aug 02 '16 at 03:49
  • @AviGinsburg both OK. Using pointers indicates the parameter is a output as suggested by some coding style guide. – kangshiyin Aug 02 '16 at 06:11