2D array of vectors

[PhysicsIsFun]

Hi guys,

can someone give me an easy to understand way on how to do this?
I have N² 2d vectors and I would like to store them in a N x N-array.
So, mathematically I want a simple N x N matrix with each element being a pair of doubles.

I have found this for 2d arrays, can I modify this somehow?

 int** ary = new int*[rowCount];
for(int i = 0; i < rowCount; ++i)
    ary[i] = new int[colCount];

Regards!

[Thomas1965]

In C++ you should not use dynamic arrays.
Better to use an vector<vector<double>> or even better create a class Array2D - sth. like this perhaps:

#include <vector>
#include <algorithm>
#include <cassert>

/**
 * class Array2D - Implementation of an 2D Array
 **/

template<typename T>
class Array2D
{
public:
	Array2D(size_t rows, size_t cols): m_RowCount(rows), m_ColCount(cols)
	{
		assert(rows > 0 && cols > 0);
		m_Items.resize(rows * cols);
		std::fill(m_Items.begin(), m_Items.end(), T());
	}
	Array2D(size_t rows, size_t cols, const T& val) : m_RowCount(rows), m_ColCount(cols)
	{
		assert(rows > 0 && cols > 0);
		m_Items = new T[rows * cols];
		std::fill(m_Items.begin(), m_Items.end(), val);
	}

	size_t rows() const 
	{ 
		return m_RowCount; 
	}
	
	size_t cols() const 
	{ 
		return m_ColCount; 
	}

	T& operator()(const size_t row, const size_t col)
	{
		assert(row > 0 && col > 0);
		size_t index = row * m_ColCount + col;
		
		return m_Items[index];
	}
	
	const T& operator()(const size_t row, const size_t col) const
	{
		assert(row > 0 && col > 0);
		size_t index = row * m_ColCount + col;

		return m_Items[index];
	}

	T* operator[](size_t row)
	{
		assert(row < m_RowCount);
		return &m_Items[row * m_ColCount];
	}

	const T* operator[](size_t row) const
	{
		assert(row < m_RowCount);

		return &m_Items[row * m_ColCount];
	}

	
private:
	std::vector<T> m_Items;
	size_t  m_RowCount, m_ColCount;
};

CAUTION: This code is just an example of how things can be done.
It is neither the only way nor the best way.
It is not properly tested and might not be correct.

[MikeyBoy]

Sure. Instead of making your elements int, make them std::pair<double, double>

http://en.cppreference.com/w/cpp/utility/pair

Or, if you prefer, create a struct whose members are the two values you're storing, and make that the element type.

EDIT: Ninja'd. And, as Thomas suggests, using std::vector is better than dynamically allocated C-style arrays. To avoid confusion, note that the std::vector type in C++ has nothing to do with the mathematical meaning of "vector"; it's essentially an object like an array, that manages its own memory so that you don't have to.

[PhysicsIsFun]

I think I don't know enough yet to understand that code. I am not even there at classes and objects and so on and neither do I know what a 'dynamic' array is.

I got my problem solved with the example I found, just changing the data type. It should be fine for the simulation stuff I am doing, but I will keep learning more about C++...

Thanks anyway!

[MikeyBoy]

"dynamic array" means that you are dynamically allocating the memory for the arrays, which is what you're doing in the code you posted. This is an old-fashioned way of doing things, and is error-prone, as you have to manage the memory yourself, and make sure it's properly cleaned up afterwards.

Seriously, if you're new to C++, and you don't understand the implications of managing your dynamically-allocated memory, you really should do yourself a favour and learn how to use std::vector as soon as possible. You'll save yourself a whole lot of headaches.

[jonnin]

Ive seen dynamic array to mean "can change in size", and/or "size unknown at compile time" and / or 'uses dynamic memory somehow'. I prefer the first (size can change).

struct and class at this level can just be thought of as lumping multiple pieces of data into a single user defined type. there is more to it, but that is a start.

[PhysicsIsFun]

Hi MikeyBoy and jonnin,

thanks for your explanation.
Unfortunately I am doomed to learning by doing, I am writing simulations for computational physics and aside from a few basics, I don't know anything about C/C++.
I try to read through a C++ book, but I am not that far there yet^^"

I finished my programm and during runtime, there seems to be an "out of range" error for a vector. After commenting-out other parts, it seems the error comes from the part which I wanted to do in this thread.

I created a 2d-Array whose elements are vectors.
Then there is lots of code.

pos is a vector, image_pos is a copy of pos. Anything else like L or rx are just double variables.

Can someone see easily, whether there is an 'out of range' error?
The error message is like:
"throwing instance 'out of range', what() vector::M_range_check__n (which is 0)>=this->size() (which is 0)

So some vector has 0 size?

  //matrix storing the 2d forces between all particles
  vector <double> **F_ij=new vector <double> *[pos.size()];
  for(int i=0;i<pos.size();i++){
    F_ij[i]=new vector <double>[pos.size()];

  }

  //filling the matrix
  for(int i=0;i<pos.size();i++){
    F_ij[i][i]={0,0};  //no force of particle on itself
    for(int j=0; j<i;j++){
        F_ij[i][j]={-F_ij[j][i].at(0),-F_ij[j][i].at(1)}; //Newton's 3rd law
    }
    for(int j=i+1; j<pos.size();j++){
        rx=image_pos.at(i).at(0)-image_pos.at(j).at(1); //distances between particles
        ry=image_pos.at(i).at(1)-image_pos.at(j).at(1);
        if(rx>=L/2){
                rx-=L;
                image_pos.at(j).at(0)+=L;   //imaging particle to neighbouring box
        }
        else if(rx<=-L/2){
                rx+=L;
                image_pos.at(j).at(0)-=L;
        }
        if(ry>=L/2){
                ry-=L;
                image_pos.at(j).at(1)+=L;
        }
        else if(ry<=-L/2){
                ry+=L;
                image_pos.at(j).at(1)-=L;
        }
        r=sqrt(pow(rx,2)+pow(ry,2));
        if(r<=R){                           //calculate forces
                f_abs=48*pow(r,-14)-24*pow(r,-8);
                fx=f_abs*(image_pos.at(i).at(0)-image_pos.at(j).at(0));
                fy=f_abs*(image_pos.at(i).at(1)-image_pos.at(j).at(1));
                F_ij[i][j]={fx, fy};

                Epot+=4*(pow(r,-12)-pow(r,-6))+1;
        }
    }

Btw. it deals with N interacting particles in a box under periodic boundary conditions, just in case one is interested :)

[lastchance]

Line 15
rx=image_pos.at(i).at(0)-image_pos.at(j).at(1); //distances between particles
should be
rx=image_pos.at(i).at(0)-image_pos.at(j).at(0); //distances between particles

Are your forces on lines 36 and 37 consistent with what you have just done to rx and ry? Also, if f = -dV/dr (as it usually is) then your potential has the wrong sign.

Neither of these will cause your range-check error. Maybe you need to cout what pos.size() and image_pos.size() are before this code snippet: they might not be what you hoped.

[PhysicsIsFun]

Ah, thanks, of course!
I will change it.

pos is a vector or N=144 double pairs (the x and y coordinates of N particles).
I declared it with vector <vector <double>> pos=readin_pos();

using a read-in function for a binary file with 144 lines of the form: double x double y

vector<vector <double>> readin_pos(){
    vector <vector <double>> pos;
    ifstream inFile;
    inFile.open("positions.txt");
    if(inFile.fail()){                  // Error test
        cerr<< "Error opening file"<<endl;
        exit(1);
    }
    vector <double> pos_xy;
    double pos_x;
    double pos_y;
    while(!inFile.eof()){
            inFile>> pos_x;
            inFile>> pos_y;
            pos_xy={pos_x, pos_y};
            pos.push_back(pos_xy);
    }
    inFile.close();
    pos.pop_back();
    return pos;
}

This function seems to work.

[PhysicsIsFun]

Now, with cout commands after the end of every loop, it seems the error occurs when i=1 in the first j-loop-iteration.

//matrix storing the 2d forces between all particles
  vector <double> **F_ij=new vector <double> *[pos.size()];
  for(int i=0;i<pos.size();i++){
    F_ij[i]=new vector <double>[pos.size()];

  }

  //filling the matrix
  for(int i=0;i<pos.size();i++){
    F_ij[i][i]={0,0};  //no force of particle on itself
    cout<<"diagonal element for i="<<i<<" works"<<endl;
    for(int j=0; j<i;j++){
        F_ij[i][j]={-F_ij[j][i].at(0),-F_ij[j][i].at(1)};  //Newton's 3rd law
        cout <<"Newton's 3rd law for j="<<j<<" works"<<endl;
    }
    for(int j=i+1; j<pos.size();j++){
        rx=image_pos.at(i).at(0)-image_pos.at(j).at(0); //distances between particles
        ry=image_pos.at(i).at(1)-image_pos.at(j).at(1);
        if(rx>=L/2){
                rx-=L;
                image_pos.at(j).at(0)+=L;   //imaging particle to neighbouring box
        }
        else if(rx<=-L/2){
                rx+=L;
                image_pos.at(j).at(0)-=L;
        }
        if(ry>=L/2){
                ry-=L;
                image_pos.at(j).at(1)+=L;
        }
        else if(ry<=-L/2){
                ry+=L;
                image_pos.at(j).at(1)-=L;
        }
        r=sqrt(pow(rx,2)+pow(ry,2));
        if(r<=R){                           //calculate forces
                f_abs=48*pow(r,-14)-24*pow(r,-8);
                fx=f_abs*(image_pos.at(i).at(0)-image_pos.at(j).at(0));
                fy=f_abs*(image_pos.at(i).at(1)-image_pos.at(j).at(1));
                F_ij[i][j]={fx, fy};

                Epot+=4*(pow(r,-12)-pow(r,-6))+1;
        }
        cout<<"j: "<<j<<endl;
    }
    cout <<"i: "<< i<< "  "<<pos.size()<<endl;
  }

[PhysicsIsFun]

Ahh, I have the mistake.

The if() in line 36 decides wether F_ij gets filled. If it doesn't happen, there will be an error in line 13, since the matrix has no elements yet.

I thought the vectors in the array would all be autoset to {0,0}, but it doesn't look like it.
Can I somehow add the initialization in line 4 so that all vectors are {0,0}?

[jonnin]

debugging 101 :)
you found it to be having trouble in the j loop, so now figure out where in there it goes bad. keep adding print statements.

this might help you.

void debugger(int dummy)
{
cout <<"got here " << dummy << endl;
cout.flush();
}

use:
debugger(0);
code
debugger(1); //2,3,4.... etc
...

the dummy parameter will match the print statement to find the issue.

these types of bugs are hard to spot, you are usually better off finding it with print statements than trying to eyeball spot it. There are lots of good pairs of eyes in here but even so 5-10 min of printing will answer before most of the crew here even sees the post...

[jonnin]

Yes, you can initialize the vector.

vector<double> example(size, 0.0); //sets all to 0.

do that to the innermost data vector (the external one is a vector of vectors... no value that you can tamper with).

[PhysicsIsFun]

Thanks for the tipp @Jonnin,

I see debugging is a skill in itself.
If the code does not run afterwards, I will try this method.

Can you answer the last question concerning the initialisation of the vectors in the array? Or do I have to write another loop.

[PhysicsIsFun]

Like this?

//matrix storing the 2d forces between all particles
  vector <double> **F_ij=new vector <double> *[pos.size()];
  for(int i=0;i<pos.size();i++){
    F_ij[i]=new vector <double>(2,0.0) [pos.size()];

[jonnin]

whats with the pointers?

this is a 2d vector 5x5 initialized to zero:

vector< vector< int > > x (5, vector<int> ( 5, 0 ) );

if you want to do something to it with pointers, we can do that, but at that point you have something... weird.

BTW this is one of those places where c++ syntax breaks down. you might expect it to be

vector< vector< int >(5,0) > x (5);
but that is not correct.

[John87Connor]

there are 2 ways to implement NxN matrices:

1. you use a container that manages the memory for you (as described by Thomas1965 in the 2nd post here)

2. you dont let the memory be managed dynamically by containers (like std::vector or so)

consider this:

// your element of 2 components (you could also use std::pair<...>)
template<typename T>
struct Element { T First, Second; };

// a vector (or column in the matrix)
template<typename T, unsigned int dimensions>
struct Vector {
private:
    T data[dimensions];
public:
    T& operator[](unsigned int idx) { return data[idx]; }
    const T& operator[](unsigned int idx) const { return data[idx]; }
};

// the final 2D matrix
const unsigned int numberofrows= 3;
const unsigned int numberofcolumns= 4;
typedef Vector<Vector<Element<double>, numberofcolumns>, numberofrows> Matrix;

with that, you can create an array of matrices that are consecutive layed down in memory, otherwise with the 1. aproach you'd have all the matrices in N different memory locations:

unsigned int N = 50;
std::vector<Array2D> firstversion(N); // uses N different memory locations
std::vector<Matrix> secondversion(N); // uses 1 memory location

the disadvantage of te 2. version is that you cant create too large NxM matrices on the stack (stack overflow). have a look at this library:
https://glm.g-truc.net/0.9.2/api/a00001.html