Segmentation fault, while running a Monte-Carlo-random-walk simulation in c++

Question

I am trying to simulate a random walk of 2000 particles, while the one boundary has the ability to make particles bound on that and merely perform a biased step.

There are of course probabilities for binding unbinding etc... Below I have the whole code. However I get segfault error. I put some print statements in the code to see where the issue lies. But nothing. What I found strange though, is that although seed is fixed, the length of the output statement determined the loop, where code crushed.

I am totally inexperienced in these issues, so if you have any idea on what I could do, would be appreciated.

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <fstream>
#include <string>
using namespace std;


const int pi=6;   
const int epsilon=10;


const int X=3000; 
const int Y=30; //length
const int time_steps=100000;
const int N=2000; // number of molecules

int kinesins[N][3]={0};//[X,Y,bound or not]
int grid[X][Y][2]={0};

void place_kinesins(){
    for (int i=0; i<N;i++){
        int x= rand()%X;
        int y= (rand()%(Y-2))+2;
        if (grid[x][y][0]==0){
            kinesins[i][0]=x;
            kinesins[i][1]=y;
            kinesins[i][2]=0;
            grid[x][y][0]=1;
        }else{i--;}
    }   
}
void create_boundaries(){
    for(int i=0;i<Y;i++){
        grid[0][i][1]=-1;
        grid[X-1][i][1]=-3;
    
    }
    for (int i=0; i<X; i++){
        grid[i][Y-1][1]=-2;
        }
}
void create_filament(){ //in order to create binding affinity.
    for(int i=0; i<X;i++){
        grid[i][1][1]=pi;   
    }
}

void step(int kinesin, int x_step, int y_step){
    int x=kinesins[kinesin][0];
    int y=kinesins[kinesin][1];
    int x_end=x+x_step;
    int y_end=y+y_step;
    if (grid[x_end][y_end][0]==0){      
        grid[x][y][0]=0;
        kinesins[kinesin][0]=x_end;
        kinesins[kinesin][1]=y_end;
        grid[x_end][y_end][0]=1;
    }

}                                   

void bound(int kinesin){
    int roll=rand()%10000 ;
    if (roll<epsilon){
        kinesins[kinesin][2]=0;
        step(kinesin,0,1);
    }else{
        if (roll%63==0){                        //controls the binding rate speed
            step(kinesin, 1,0);
        };
    }
}

void unbound(int kinesin){
    cout<<"1";
    int x= kinesins[kinesin][0];
    int y= kinesins[kinesin][1];
    
    int type= grid[x][y][1];
    
    switch(type){
        case 0:{
        cout<<"2";
            int roll=rand()%4;
            
            switch(roll){
                
                case 0:
                    step(kinesin,-1,0);
                    break;
                case 1:
                    step(kinesin,1,0);
                    break;
                case 2:
                    step(kinesin,0,1);
                    break;
                case 3:             
                    step(kinesin,0,-1);
                    break;
            }
            break;
            }
        case -1:
            step(kinesin,1,0);
            break;
        case -2:
            step(kinesin,0,-1);
            break;
        case -3:
            step(kinesin,-1,0);
            break;
        default:
            int roll=rand()%10000;
            if(roll<grid[x][y][1]){kinesins[kinesin][2]=1;}
            else{ if(roll%2==0){step(kinesin,0,1);}}
    }
}



void kinesin_move(int kinesin){
    cout<<" "<<kinesins[kinesin][0]<<kinesins[kinesin][1];
    if (kinesins[kinesin][2]==0){ 
        unbound(kinesin);
    }else{
    cout<<"3";
        bound(kinesin);
    }

}



        
        
void simulation(){
    
    
    
    for(int j=7000; j<time_steps;j++){
    cout<<endl<< j<<"  "<<endl;
        
        for (int kin=0; kin<N; kin++){
            cout<<kin;
            kinesin_move(kin);
            cout<<"E  " ;
            }   
                
    }
    
}




void programm(){
    srand(1);
    create_boundaries();
    create_filament();
    cout<<"Filament done"<<endl;
    place_kinesins();
    cout<<"Kines placed"<<endl;
    simulation();

}


int main(){
    programm();
    return 0;
}

Answer 1

Problem:

In the function step you're accessing the array grid out of its bonds, which produces Undefined Behaviour and the segmentation fault.

This can be proven adding an assert before if (grid[x_end][y_end][0]==0):

   if(!(x_end < X && x_end >= 0))
        std::cerr << x_end << std::endl;
    assert(x_end < X && x_end >= 0);
    if (grid[x_end][y_end][0]==0){
        grid[x][y][0]=0;
        kinesins[kinesin][0]=x_end;
        kinesins[kinesin][1]=y_end;
        grid[x_end][y_end][0]=1;
    }

Output:

3000
Assertion failed: x_end < X && x_end >= 0 main.cpp line 57

This application has requested the Runtime to terminate it in an unusual way.

Solution:

You will have to check the arguments for step won't make it go out of bonds before each call.

Additional information:

1. using namespace std; is considered a bad practice (More info here).
2. rand is not uniformly distributed. You may want use the <random> library instead (More info here).
3. Global variables are bad (More info here).
4. You may want to use std::array instead of raw, C-style arrays.