User:Dllu/BML

From Wikipedia, the free encyclopedia
Sample output for a large lattice of 2000 by 2000, with a density of 31%.

Here's some C++ code for dealing with the Biham-Middleton-Levine traffic model. Note: this code was only used to generate the pictures in the said article, and not the videos. I used MATLAB for the videos instead. That said, this code is far faster (about 10x faster) and more elegant.

Just like how Rule 184 consists of a one-dimensional array of cells, each containing a binary value (0 or 1), the BML Traffic Model consists of a two-dimensional array of cells, each containing a trinary value (0, 1, or 2). Hence the "secret sauce" in the code, which really isn't that secret. It just looks cool.

Also, system("PAUSE") only works in Windows. Just delete that line if you want to compile it in something else.

The code shown here should output some bitmap images. Enjoy! Please suggest any improvements in the talk page.

A more interactive real-time Processing version is available at Open Processing.

Code[edit]

//Purpy Pupple's code for dealing with the Biham-Middleton-Levine Traffic Model!
//http://en.wikipedia.org/wiki/user:Purpy_Pupple/BML
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <vector>
#define one true //ones move vertically
#define two false//twos move horizontally
#define OUTFILE "bml"
#define MAX(x,y) (x>y?x:y)
#define MIN(x,y) (x<y?x:y)
int PERCENT=32;
const int H_SIZE=512;//height
const int W_SIZE=512;//width
const int statesone[27]={0,0,0,0,1,1,2,2,2,1,1,1,0,1,1,2,2,2,0,0,0,0,1,1,2,2,2,};//secret sauce
const int statestwo[27]={0,0,0,1,1,1,0,2,2,0,0,0,1,1,1,0,2,2,2,2,2,1,1,1,0,2,2,};//more secret sauce
std::vector<int>mobilityone; //mobility
std::vector<int>mobilitytwo; //mobility
int X[H_SIZE][W_SIZE];//to store the current state before it is one's turn to move
int Y[H_SIZE][W_SIZE];//to store the current state before it is two's turn to move

void printgraph();
void printimg(bool oneortwo, int iterations);

void trafficroll(bool oneortwo){
	//"rolls" traffic depending on whose turn it is.
	int a,b,c;//the value of the spot before, at, and after the position you're currently at.
	int mobility=0;//mobility
	if(oneortwo==one){
		for(int i=0;i<H_SIZE;i++){
			for(int j=0;j<W_SIZE;j++){
				a=X[(i-1+H_SIZE)%H_SIZE][j]; b=X[i][j]; c=X[(i+1)%H_SIZE][j];
				Y[i][j]=statesone[c+3*b+9*a];//awesome base three arithmetic!
				if(Y[i][j]>X[i][j]) mobility++;
			}
		}
		mobilityone.push_back(mobility);
	}else{
		for(int i=0;i<H_SIZE;i++){
			for(int j=0;j<W_SIZE;j++){
				a=Y[i][(j-1+W_SIZE)%W_SIZE]; b=Y[i][j]; c=Y[i][(j+1)%W_SIZE];
				X[i][j]=statestwo[c+3*b+9*a];//awesome base three arithmetic!
				if(X[i][j]>Y[i][j]) mobility++;
			}
		}
		mobilitytwo.push_back(mobility);
	}
	return;
}

void trafficset(){
	//for initializing the X. Always start with X, not Y!
	int n, xn, yn;
	for(int j=0;j<H_SIZE;j++){
		for(int k=0;k<W_SIZE;k++){
			X[j][k]=0;//initialize everything to zero!
		}
	}
	for(int i=0;i<PERCENT*H_SIZE*W_SIZE/200;i++){
		if(H_SIZE*W_SIZE<32768){//I feel like optimizing this by calling rand() less if the size is small enough. But it doesn't do much.
			do{
				n=rand()%(H_SIZE*W_SIZE);
			}while(X[n/W_SIZE][n%W_SIZE]!=0);//find a random place that's empty
			X[n/W_SIZE][n%W_SIZE]=1;         //set that as a one
			do{
				n=rand()%(H_SIZE*W_SIZE);
			}while(X[n/W_SIZE][n%W_SIZE]!=0);//find a random place that's empty
			X[n/W_SIZE][n%W_SIZE]=2;         //set that as a two
		}else{
			do{
				xn=rand()%(H_SIZE);
				yn=rand()%(W_SIZE);
			}while(X[xn][yn]!=0);//find a random place that's empty
			X[xn][yn]=1;         //set that as a one
			do{
				xn=rand()%(H_SIZE);
				yn=rand()%(W_SIZE);
			}while(X[xn][yn]!=0);//find a random place that's empty
			X[xn][yn]=2;         //set that as a two
		}
	}
	return;
}

void print(bool oneortwo){
	//for printing out the ascii diagram in command line!
	for(int j=0;j<H_SIZE;j++){
		for(int k=0;k<W_SIZE;k++){
			oneortwo?
				printf("%c",Y[j][k]==0?' ':(Y[j][k]==1?'|':'-'))://two
				printf("%c",X[j][k]==0?' ':(X[j][k]==1?'|':'-'));//one
		}
		printf("\n");
	}
	printf("\n\n");
}

int main(){
	//for(PERCENT=26;PERCENT<40;PERCENT+=1){ //uncomment this if you want to batch produce images for many different percents
		//std::vector<long> hashes;
		srand ( time(NULL) );//seed random number generator
		trafficset();//initialize the thing
		//print(two);
		printimg(two,0);
		int start=clock();
		int i=0,j=0;
		for(i=0;i<32000;i++){
			trafficroll(one);
			trafficroll(two);
			//if(i%200==0) printimg(two,i); //something like this is useful for making movies.
		}
		int end=clock();
		//print(two);
		printimg(two,i);
		printgraph();
		while(!mobilityone.empty()){
			mobilityone.pop_back();
		}
		while(!mobilitytwo.empty()){
			mobilitytwo.pop_back();
		}
		printf("\nElapsed time is %f seconds. %d loops iterated. \n",(end-(float)start)/CLOCKS_PER_SEC, i);
		//system("PAUSE");
	//}
	return 0;
}

//From now on you will find that the majority of the code deals with printing bitmaps.

void printimg(bool oneortwo, int iterations){
	//print out the lattice to a bitmap file.
	//the code to output bitmaps was taken from http://en.wikipedia.org/wiki/User:Evercat/Buddhabrot.c
	unsigned int headers[13];
	FILE * outfile;
	int extrabytes;
	int paddedsize;
	int x; int y; int n;

	extrabytes = 4 - ((W_SIZE * 3) % 4); 

	char filename[200];
	
	sprintf(filename, "%s x %d y %d p %d iterated %d.bmp", OUTFILE, H_SIZE, W_SIZE, PERCENT, iterations);

	if (extrabytes == 4) extrabytes = 0;

	paddedsize = ((W_SIZE * 3) + extrabytes) * H_SIZE;

	// Headers...
	// Note that the "BM" identifier in bytes 0 and 1 is NOT included in these "headers".
                     
	headers[0]  = paddedsize + 54;      // bfSize (whole file size)
	headers[1]  = 0;                    // bfReserved (both)
	headers[2]  = 54;                   // bfOffbits
	headers[3]  = 40;                   // biSize
	headers[4]  = W_SIZE;  // biWidth
	headers[5]  = H_SIZE; // biHeight

	// Would have biPlanes and biBitCount in position 6, but they're shorts.
	// It's easier to write them out separately (see below) than pretend
	// they're a single int, especially with endian issues...

	headers[7]  = 0;                    // biCompression
	headers[8]  = paddedsize;           // biSizeImage
	headers[9]  = 0;                    // biXPelsPerMeter
	headers[10] = 0;                    // biYPelsPerMeter
	headers[11] = 0;                    // biClrUsed
	headers[12] = 0;                    // biClrImportant

	outfile = fopen(filename, "wb");

	//
	// Headers begin...
	// When printing ints and shorts, we write out 1 character at a time to avoid endian issues.
	//

	fprintf(outfile, "BM");

	for (n = 0; n <= 5; n++){
	   fprintf(outfile, "%c", headers[n] & 0x000000FF);
	   fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
	   fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
	   fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
	}

	// These next 4 characters are for the biPlanes and biBitCount fields.

	fprintf(outfile, "%c", 1);
	fprintf(outfile, "%c", 0);
	fprintf(outfile, "%c", 24);
	fprintf(outfile, "%c", 0);

	for (n = 7; n <= 12; n++){
	   fprintf(outfile, "%c", headers[n] & 0x000000FF);
	   fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
	   fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
	   fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
	}

	//
	// Headers done, now write the data...
	//
	for(y = H_SIZE - 1; y >= 0; y--){
		for(x = 0; x <= W_SIZE - 1; x++){
			if(oneortwo){
				if(Y[x][y]==0){
					fprintf(outfile, "%c", 255); fprintf(outfile, "%c", 255); fprintf(outfile, "%c", 255);
				}else if(Y[x][y]==1){
					fprintf(outfile, "%c", 0); fprintf(outfile, "%c", 0);   fprintf(outfile, "%c", 255);
				}else if(Y[x][y]==2){
					fprintf(outfile, "%c", 255);   fprintf(outfile, "%c", 0);   fprintf(outfile, "%c", 0);
				}
			}else{
				if(X[x][y]==0){
					fprintf(outfile, "%c", 255); fprintf(outfile, "%c", 255); fprintf(outfile, "%c", 255);
				}else if(X[x][y]==1){
					fprintf(outfile, "%c", 0); fprintf(outfile, "%c", 0);   fprintf(outfile, "%c", 255);
				}else if(X[x][y]==2){
					fprintf(outfile, "%c", 255);   fprintf(outfile, "%c", 0);   fprintf(outfile, "%c", 0);
				}
			}
		}
		if (extrabytes){     // See above - BMP lines must be of lengths divisible by 4.
			for (n = 1; n <= extrabytes; n++){
				fprintf(outfile, "%c", 0);
			}
		}
	}
	printf("file printed: %s\n", filename); 
	fclose(outfile);
	return;
}

void printgraph(){
	//output a graph that's approximately 5000 by 2000 (but varies) that shows the mobility vs time
	//the code to output bitmaps was taken from http://en.wikipedia.org/wiki/User:Evercat/Buddhabrot.c
	unsigned int headers[13];
	FILE * outfile;
	int extrabytes;
	int paddedsize;
	int x; int y; int n;

	int scale_height=(H_SIZE*W_SIZE*PERCENT/200)/2000;
	int scale_width=MIN(mobilityone.size(),mobilitytwo.size())/5000;

	int graph_height=H_SIZE*W_SIZE*PERCENT/200/scale_height;
	int graph_width=MIN(mobilityone.size(),mobilitytwo.size())/scale_width;

	extrabytes = 4 - ((graph_width * 3) % 4); 

	char filename[200];
	
	sprintf(filename, "%s x %d y %d p %d MOBILITY.bmp", OUTFILE, H_SIZE, W_SIZE, PERCENT);

	if (extrabytes == 4) extrabytes = 0;

	paddedsize = ((graph_width * 3) + extrabytes) * graph_height;

	// Headers...
	// Note that the "BM" identifier in bytes 0 and 1 is NOT included in these "headers".
                     
	headers[0]  = paddedsize + 54;      // bfSize (whole file size)
	headers[1]  = 0;                    // bfReserved (both)
	headers[2]  = 54;                   // bfOffbits
	headers[3]  = 40;                   // biSize
	headers[4]  = graph_width;  // biWidth
	headers[5]  = graph_height; // biHeight

	// Would have biPlanes and biBitCount in position 6, but they're shorts.
	// It's easier to write them out separately (see below) than pretend
	// they're a single int, especially with endian issues...

	headers[7]  = 0;                    // biCompression
	headers[8]  = paddedsize;           // biSizeImage
	headers[9]  = 0;                    // biXPelsPerMeter
	headers[10] = 0;                    // biYPelsPerMeter
	headers[11] = 0;                    // biClrUsed
	headers[12] = 0;                    // biClrImportant

	outfile = fopen(filename, "wb");

	//
	// Headers begin...
	// When printing ints and shorts, we write out 1 character at a time to avoid endian issues.
	//

	fprintf(outfile, "BM");

	for (n = 0; n <= 5; n++){
	   fprintf(outfile, "%c", headers[n] & 0x000000FF);
	   fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
	   fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
	   fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
	}

	// These next 4 characters are for the biPlanes and biBitCount fields.

	fprintf(outfile, "%c", 1);
	fprintf(outfile, "%c", 0);
	fprintf(outfile, "%c", 24);
	fprintf(outfile, "%c", 0);

	for (n = 7; n <= 12; n++){
	   fprintf(outfile, "%c", headers[n] & 0x000000FF);
	   fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
	   fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
	   fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
	}

	//
	// Headers done, now write the data...
	//
	for(y = 0; y<= graph_height - 1;y++){
		for(x = 0; x <= graph_width - 1; x++){
			int valuer=255;
			int valueg=255;
			int valueb=255;
			for(int i=0;i<scale_height;i++){
				for(int j=0;j<scale_width;j++){
					if(mobilityone[x*scale_width+j]==y*scale_height+i){
						valueg=0;valueb=0;
					}
					if(mobilitytwo[x*scale_width+j]==y*scale_height+i){
						valueg=0;valuer=0;
					}
				}
			}
			fprintf(outfile, "%c", valueb);   fprintf(outfile, "%c", valueg);   fprintf(outfile, "%c", valuer);
		}
		if (extrabytes){     // See above - BMP lines must be of lengths divisible by 4.
			for (n = 1; n <= extrabytes; n++){
				fprintf(outfile, "%c", 0);
			}
		}
	}
	printf("file printed: %s\n", filename); 
	fclose(outfile);
	return;
}

Gallery of images[edit]

Known bugs or issues[edit]

  • If the mobility is 1, sometimes it disappears in the mobility graph.
  • The mobility graph is a scatter plot - consequently the dots are not joined and may look faint.

Real-time Processing version[edit]

A version has been ported to Processing. Click here.

Future improvements[edit]

  • Terminate program immediately when it reaches global jam or free flow (I did this with my MATLAB script) so that it doesn't waste time doing nothing.