Unrepeatable 2D vector problem

I'm writing a small program to calculate a physics problem but I'm having problems erasing elements from a 2D array. I'm compiling with XCode (so GDB).

The problem always comes in the walker.erase(walker.begin()+loacationInArray) part (in the void diffuseWalkers function), and always after a different number of function calls. So sometimes, the 'thermalisation' loop runs maybe 50 times, other times it runs all the way through. I usually get an EXC_BAD_ACCESS error, and very occasionally a malloc error.

Any help would be much appreciated as I've tried everything and really can't see what I'm doing wrong. I've posted the code below. A header file called 'heliumHeader.h' contains some mathematical functions that aren't really important for my problem.

    #include "heliumHeader.h"


    void copyVectorInformation(vector<double>& walker, vector<double>& walkerTemp);
    void diffuseWalker(double beta, double a, double alpha, int locationInArray, vector<double>& walker, double dt, 
           double& trialEnergy, int& numberOfWalkers, int targetNumberOfWalkers);
    using namespace std;

    int main(){


 srand(time(NULL));

 double s=0.; //inter-molecular distance
 double beta = 0.15;
 double alpha = 2.;
 double dt=0.1; 
 int numberOfWalkers = 1; //number of particles
 int targetNumberOfWalkers = numberOfWalkers;

 //2D-vectors to hold information
 vector<vector<double> > walker(numberOfWalkers,6); 

 //solve for a
 double a = 1.;
 for(int i=0; i<20;i++)
  a = 1./(1.+exp(-s/a));

 //set up sums
 double localEnergy, trialEnergy, localEnergy2, trialEnergy2;
 localEnergy = trialEnergy = localEnergy2 = trialEnergy2 = 0.;

 //put the walkers randomly in space & get the energy of that configuration
 for(int i=0; i<walker.size(); i++){
  for(int j=0; j<6; j++)
   walker[i][j]=randomPositiveNegative();
  localEnergy += calculateEnergy(beta, alpha, a, walker[i]);
 }

 localEnergy /= numberOfWalkers;

 double beforeEnergy = localEnergy;
 cout << "local energy of random " << localEnergy << endl;
 trialEnergy = -2.903;

 //move the walkers
 for(int thermalisationCounter = 1; thermalisationCounter<1000; thermalisationCounter++){
  for(int i=0; i<walker.size(); i++)
   diffuseWalker(beta, a, alpha, i, walker[i], dt, trialEnergy, numberOfWalkers, targetNumberOfWalkers);
  cout << thermalisationCounter << endl;
 }

 //recalculate the local energy
 for(int i=0; i<walker.size(); i++){
  for(int j=0; j<6; j++)
   localEnergy += calculateEnergy(beta, alpha, a, walker[i]);
 }

 localEnergy /= numberOfWalkers;

 for(int numberOfSteps = 1; numberOfSteps<1000; numberOfSteps++){
  for(int i=0; i<walker.size(); i++)
   diffuseWalker(beta, a, alpha, i, walker[i], dt, trialEnergy, numberOfWalkers, targetNumberOfWalkers);
  cout << numberOfSteps << endl;
 }

 //get initial energy of random positiions
 for(int i=0; i<walker.size(); i++)
  localEnergy += calculateEnergy(beta, alpha, a, walker[i]);

 localEnergy /= numberOfWalkers;

 cout << "before energy " << beforeEnergy << " local energy " << localEnergy << 
 " trial energy " << trialEnergy << " number " << numberOfWalkers << endl;

 return 0;
    }


    void copyVectorInformation(vector<double>& walker, vector<double>& walkerTemp){
     for(int i=0; i<6; i++)
      walkerTemp[i] = walker[i];
    }//end of copyVectorInformation

    void diffuseWalker(double beta, double a, double alpha, int locationInArray, vector<double>& walker, double dt, 
           double& trialEnergy, int& numberOfWalkers, int targetNumberOfWalkers){

 vector<double> driftFunction(6);
 vector<double> walkerTemp(6); //temporary store of information

 //copy the walker information over
 copyVectorInformation(walker, walkerTemp);

 //get the drift functions
 calculateDriftFunctions(beta, alpha, a, walker, driftFunction);

 //get the old local energy
 double preMoveLocalEnergy = calculateEnergy(beta, alpha, a, walker);

 //move the walker
 for(int j=0; j<6;j++)
  walker[j] += 0.5*dt*driftFunction[j] + randomGauss()*sqrt(dt);

 //caclulate the local energy of the new position
 double postMoveLocalEnergy = calculateEnergy(beta, alpha, a, walker);

 //calculate the weight, q and branching ration s
 double q = exp(-dt*(0.5*(postMoveLocalEnergy + preMoveLocalEnergy) - trialEnergy));
 double s = q+randomPositive();

 if(int(s)<1){
  //kill the walker 
  walker.erase(walker.begin()+locationInAr开发者_Go百科ray);
  numberOfWalkers--;
 }
 else{
  //reproduce walker int(s) number of times
  for(int k=0; k<int(s); k++){
   walker.push_back(walker[locationInArray]);
   numberOfWalkers++;
  }
 }

 //update the trial energy
 trialEnergy += 0.2*log((double)targetNumberOfWalkers/(double)numberOfWalkers);

    }//end of diffuse walkers

You have:

 if(int(s)<1){
  //kill the walker 
  walker.erase(walker.begin()+locationInArray);
  numberOfWalkers--;
 }

assuming

(int(s) < 1)

is always true (always erasing) then

 //move the walkers
 for(int thermalisationCounter = 1; thermalisationCounter<1000; thermalisationCounter++){
  for(int i=0; i<walker.size(); i++)
   diffuseWalker(beta, a, alpha, i, walker[i], dt, trialEnergy, numberOfWalkers, targetNumberOfWalkers);
  cout << thermalisationCounter << endl;
 }

will remove all diagonal elements. Is it that what you intended? Maybe you mixed row index with column index (or 1D with 2D)?

Anyway, the problem is that locationInArray is too big and you are erasing an item that is no longer in the vector.

I think you are erasing from the wrong vector. In the main function walker is a 2D vector. In diffuseWalker you are erasing from the walker that has been passed into the function which is a 1D vector of doubles. You use locationInArray and decrement numberOfWalkers here which doesn't make sense for the 1D walker - it looks like you want to actually erase from the 2D vector. I'd suggest using a different name for the 2D and 1D vectors to avoid confusion like this.

I think your problem lies in the fact that vectors are dynamic in their size. It looks like you are in a loop when you call .erase. If you erase something in the middle, and continue cycling through the elements, I not surprised that random weirdness ensues.

Yeah - thanks very much. It's true, I was deleting the diagnoal elements of the array, which I discovered by doing a much smaller toy example in another program.

I solved it by passing the function the full 2D vector, so that when I erase walker(walker.begin()+locationInArray) the iterator refers to the first index, 'i' of the walker[i][j] vector. Before, when I just passed it as 1D vector (essentially the 'j' index piece), I was erasing a location that didn't exist if 'locationInArray >=6' as 0 <= j <= 5.

Many thanks for your answers.