A: number of elements: 560
B: number of elements: 560
B: number of elements: 0
B: number of elements: 0
B: number of elements: 0
C: number of elements: 0

cout << "A: number of elements: " << combinations.size() << endl;
float sum_vect[120][2];
for (int i = 0; i < combinations.size(); ++i) {
    for (int j = 0; j < 2; ++j) {
        sum_vect[i][j] = 0;
    }
}
cout << "B: number of elements: " << combinations.size() << endl;
cout << "B: number of elements: " << combinations.size() << endl;
cout << "B: number of elements: " << combinations.size() << endl;
cout << "B: number of elements: " << combinations.size() << endl;

for (int i = 0; i < combinations.size(); ++i) {
    combination = combinations.at(i);
    for (int j = 0; j < order; ++j) {
        sum_vect[i][0] += (float)virtual_pos[combination[j]][0];
        sum_vect[i][1] += (float)virtual_pos[combination[j]][1];
    }
}

vector<int> optimal_ind;
cout << "C: number of elements: " << combinations.size() << endl;

void nchoosek_helper(int offset, int n, int k, vector<vector<int>> &combinations, vector<int> combination) {
        if (k == 0) {
            combinations.push_back(combination);
            return;
        }
        for (int i = offset; i <= n - k; ++i) {
            combination.push_back(i);
            nchoosek_helper(i+1, n, k-1, combinations, combination);
                combination.pop_back();
        }
    }


double euclidean_norm(double dist1, double dist2){
    return sqrt(pow(dist1,2) + pow(dist2,2));
}

//weirdest: look at A B C
vector<vector<char> > step37::CUDADriver::get_access_pattern(int order){
    vector<vector<char> > result;
    vector<vector<int> > combinations;
    vector<int> combination;

    nchoosek_helper(0, 16 ,order, combinations, combination);
    cout << "number of combinations: " << combination.size() << endl;


    //                    //mapping lexical index 1-16 to 2D array
    int virtual_pos [16][2];
    for (int i = 0; i < 16; ++i) {
        virtual_pos[i][0] = i%4 * order; //write x
        virtual_pos[i][1] = (int)ceil(i/4) * order; //write y
        cout << "mapping " << i << "to (" << i%4 * order << "'" << (int)ceil(i/4) * order<< ")" << endl;
    }

    cout << "A: number of elements: " << combinations.size() << endl;
    float sum_vect[120][2];
    for (int i = 0; i < combinations.size(); ++i) {
        for (int j = 0; j < 2; ++j) {
            sum_vect[i][j] = 0;
        }
    }
    cout << "B: number of elements: " << combinations.size() << endl;
    cout << "B: number of elements: " << combinations.size() << endl;
    cout << "B: number of elements: " << combinations.size() << endl;
    cout << "B: number of elements: " << combinations.size() << endl;

    for (int i = 0; i < combinations.size(); ++i) {
        combination = combinations.at(i);
        for (int j = 0; j < order; ++j) {
            sum_vect[i][0] += (float)virtual_pos[combination[j]][0];
            sum_vect[i][1] += (float)virtual_pos[combination[j]][1];
        }
    }

    vector<int> optimal_ind;
    cout << "C: number of elements: " << combinations.size() << endl;

    cout << "main loop"<< endl;
    for (int i = order; i < order*2; ++i) {
        for (int j = order; j < order*2; ++j) {
            int pos [2];
            //        pos[0] = i;
            //        pos[1] = j;
            pos[0] = j;
            pos[1] = i;

            cout << "current position: (" << j << "," << i << ")" << endl;

            float min_len = std::numeric_limits<float>::infinity(); //minimum length of combined vector
            float min_sum_ind = std::numeric_limits<float>::infinity(); //minimum sum of individual vectors
            int min_idx = -1;

            for (int k = 0; k < combinations.size(); ++k) {
                int curr_vect [2];
                curr_vect[0] = sum_vect[k][0] - pos[0] * order;
                curr_vect[1] = sum_vect[k][1] - pos[1] * order;

                float curr_len = euclidean_norm(curr_vect[0], curr_vect[1]);

                float min_sum_tmp = 0;
                combination = combinations[k];

                for (int l = 0; l < order; ++l) {
                    min_sum_tmp += euclidean_norm(virtual_pos[combination.at(l)][0] - pos[0],
                            virtual_pos[combination.at(l)][1]- pos[1]);
                }

                if (i==4&&j==4){
                    cout << "  ind sum: " << min_sum_tmp << "   len: " << curr_len <<  "    sv: (" << sum_vect[k][0] << "," << sum_vect[k][1] <<
                            ")  cv: (" << curr_vect[0] << ","  << curr_vect[1] << ")" <<endl;

                }

                if (min_len > curr_len ||
                        min_len == curr_len && min_sum_tmp < min_sum_ind){
                    min_len = curr_len;
                    min_idx = k;
                    min_sum_ind = min_sum_tmp;
                }
            }

            //        cout <<

            cout << "pushing minimal idx " << min_idx << endl;
            optimal_ind.push_back(min_idx);
        }
    }

    cout << "main loop done"<< endl;

    //unpack optimal combinations into relative movements
    vector<char> optimal_x((int)pow(order,3));
    vector<char> optimal_y((int)pow(order,3));

    cout << "number of elements: " << combinations.size() << endl;
    for (int i = 0; i <(int)pow(order,2); ++i) {
        cout << "optimal idx: " << optimal_ind.at(i) << endl;
        combination = combinations.at(optimal_ind.at(i));
        for (int j = 0; j < order; ++j) {
            int lex_idx = combination.at(j); //some index between 0 and 15 from 4x4 grid

            //mvt range in grid relative to thread position: -1 to +
            int relative_x = -1 + lex_idx % 4;
            int relative_y = -1 + (int) floor(lex_idx / 4);

            optimal_x[i * order + j] = relative_x;
            optimal_y[i * order + j] = relative_y;
        }
    }

    //DEBUG print
    for (int i = 0; i < (int)pow(order,2); ++i) {
        combination = combinations.at(optimal_ind.at(i));
        cout << "combination: " << i << "   ";
        for (int j = 0; j < order; ++j) {
            cout << combination.at(j) << " ";
        }
        cout <<  endl;
    }

    result.push_back(optimal_x);
    result.push_back(optimal_y);

    for (int i = 0; i < optimal_x.size(); ++i) {
        cout << (int)optimal_x.at(i) << " " << endl;
    }

    cout << "optimal sizes: " << optimal_x.size() << endl;
    cout << "optimal sizes: " << optimal_y.size() << endl;


     cout << "result size: " << result.size() << endl;

    return result;
}

A: number of elements: 560

cout << "A: number of elements: " << combinations.size() << endl;
float sum_vect[120][2];
for (int i = 0; i < combinations.size(); ++i) {
   for (int j = 0; j < 2; ++j) {
      sum_vect[i][j] = 0;
    }
}