Monte Carlo simulation
Main code to run simulations about evolutionary game theory in C.
Packages
1
2
3
4
5
6
7
8
9
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "pointers.h"
#include "regular_lattices.h"
Definitions
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
//------------------------------------------------------
#define fps 0000.01000000005 // fps^-1
//#define densidade_arquivo
#define densidade_terminal
/********************************************************************
*** Global variables ***
********************************************************************/
const int tmax=9999;//5000;
long long int N;
int L;
int G; // =4 (rede hexagonal) =5 (rede quadrada / kagome) =7 (rede triangular / quadrada 3D) =9 (rede moore / quadrada 4D)
double RUIDO;
/********************************************************************
*** GSL ***
********************************************************************/
unsigned long rseed;
const gsl_rng_type * T;
gsl_rng * rand_vec;
//#define DEBUG
void set_gsl_rng(void)
{
#ifdef DEBUG
rseed=1681248046;
#else
rseed=time(NULL);
#endif
gsl_rng_env_setup();
T = gsl_rng_default;
rand_vec = gsl_rng_alloc (T);
gsl_rng_set (rand_vec, rseed);
}
Initial distributions
1
2
3
4
5
6
7
8
9
10
11
12
13
14
/********************************************************************
*** distribuição dos estados ***
********************************************************************/
void calculo_ci_rand(int strategy[N], int **viz)
{
int n;
for(n=0; n < N; n++)
{
double aa = gsl_rng_uniform(rand_vec);
if(aa < 0.5 ) {strategy[n] = 1;}
else {strategy[n] = 0;}
}
}
Payoff calculation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
/********************************************************************
*** Payoff ***
********************************************************************/
void calculo_payoff_pd ( double *payoff, double r, double gama, double delta, int x, int strategy[N], int **viz )
{
// x: sitio central do grupo que o jogo ocorre
// y: sitio que esta jogando
int nd=0;
int nc=0;
int np=0;
int i=0;
for(i=0;i<G;i++) // vizinho 0 e o proprio sitio
{
switch(strategy[viz[x][i]])
{
case 0: ++nd; break;
case 1: ++nc; break;
case 2: ++np; break;
default:
fprintf(stderr,"ERRO calculo vizinhos\n");
fflush(stderr);
}
}
double S = -r;
double P = 0;
double R = 1;
double T = 1+r;
payoff[1] = R*(nc-1) + S*nd;
payoff[0] = T*nc + P*nd;
}
void calculo_payoff_pgg ( double *payoff, double r, double gama, double delta, int x, int strategy[N], int **viz )
{
// x: sitio central do grupo que o jogo ocorre
// y: sitio que esta jogando
int nd, nc, np, i, j, focal_node;
payoff[1] = 0;
payoff[0] = 0;
for(j=0;j<G;j++) // set only first term for FPGG (loop between groups)
{
nd = 0;
nc = 0;
np = 0;
focal_node = viz[x][j];
for(i=0;i<G;i++) // loop inside a group
{
switch(strategy[viz[focal_node][i]])
{
case 0: ++nd; break;
case 1: ++nc; break;
case 2: ++np; break;
default:
fprintf(stderr,"ERRO calculo vizinhos\n");
fflush(stderr);
}
}
payoff[1] += (r/G)*nc - 1;
payoff[0] += (r/G)*nc;
}
}
Update rule
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/********************************************************************
*** Update rule ***
********************************************************************/
void update_rule( int x, int vizinho, int strategy[N], double Px, double Py , int t)
{
double Wxy = 1.0/(1.0 + exp(-(Py-Px)/RUIDO));
double l = gsl_rng_uniform(rand_vec);
//printf("Px=%lf Py=%lf delta=%lf W=%lf \n",Px,Py,Py-Px,Wxy);
//troca de estado
if(Wxy > l)
{
strategy[x] = strategy[vizinho];
}
}
Monte Carlo Step
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
/********************************************************************
*** MCS ***
********************************************************************/
void calculo_mcs(double *payoff, int strategy[N], int **viz, double r, double gama, double delta, int t)
{
int n;
for(n=0; n<N; n++)
{
int x = (int) N*gsl_rng_uniform(rand_vec);
int sitio = x;
int y = 1 + (int) (G-1)*gsl_rng_uniform(rand_vec);
int vizinho = viz[x][y];
calculo_payoff_pgg(payoff,r,gama,delta,sitio,strategy,viz);
double Px = payoff[strategy[sitio]];
calculo_payoff_pgg(payoff,r,gama,delta,vizinho,strategy,viz);
double Py = payoff[strategy[vizinho]];
//printf("%d %f %f\n", t, Px, Py);
update_rule(x,vizinho,strategy,Px,Py,t);
}
}
Densities calculation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
/********************************************************************
*** Densidades ***
********************************************************************/
int calculo_densidades(int strategy[N], double *payoff, int **viz, double r, double gama, double delta, int t, FILE *fp)
{
int k;
int ND=0;
int NC=0;
int NP=0;
for(k=0;k<N;k++)
{
switch (strategy[k])
{
case 1: ++NC; break;
case 2: ++NP; break;
case 0: ++ND; break;
default:
fprintf(stderr,"ERRO - tipo de estrategia\n");
fflush(stderr);
}
}
#ifdef densidade_terminal
printf("%d %lf %lf\n", t, (double)NC/(N), (double)ND/(N));
#endif
#ifdef densidade_arquivo
fprintf(fp,"%d %lf %lf\n", t, (double)NC/(N), (double)ND/(N));
#endif
// Return 1 if any strategy dominates, 0 otherwise
// used to stop simulation in the first scenario
if (NC == N || ND == N || NP == N) {
return 1;
} else {
return 0;
}
}
Main
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
/********************************************************************
*** Main ***
********************************************************************/
/*
- initial conditions
- time dynamics
- calculate payoffs
- update strategy
- compute densities
*/
void main(int argc, char *argv[])
{
//------------------------------------------------------------------------------
// gcc dilema.c -lm -lgsl -lgslcblas -O3 ---- //***
// ./a.out r gama delta L ruido //***
//***
double payoff[2], gama, delta, r; //***
int **viz; //***
//***
r=atof(argv[1]); //***
gama=atof(argv[2]); //***
delta=atof(argv[3]); //***
L=atoi(argv[4]); //***
RUIDO=atof(argv[5]); //***
//***
set_gsl_rng(); // gsl //***
//***
//gera arquivo //***
char filename[200]; //***
FILE *fp; //***
#ifdef densidade_arquivo
sprintf(filename,"dados_r%f_g%f_d%f_seed%ld.txt",r,gama,delta,rseed); //***
#endif //***
fp=fopen(filename,"w"); //abre .dat //***
//***
// -------------------------------------- //***
// specify depending on regular_lattices.h //***
N = L*L; //***
G = 5; //***
viz = create_2d_int_pointer_h(N,G); //***
square_lattice(viz,L); //***
int strategy[N]; // sempre depois de declarar N //***
//-------------------C.I.-------------------------------------------------------
calculo_ci_rand(strategy,viz);
calculo_densidades(strategy,payoff,viz,r,gama,delta,0,fp);
//-------------------MCS--------------------------------------------------------
int t;
for(t=1; t < tmax; t++)
{
calculo_mcs(payoff,strategy,viz,r,gama,delta,t);
calculo_densidades(strategy,payoff,viz,r,gama,delta,t,fp);
// stop the simulation if one strategy dominates
int stop = calculo_densidades(strategy, payoff, viz, r, gama, delta, t, fp);
if (stop == 1)
{
while (t < tmax)
{
calculo_densidades(strategy, payoff, viz, r, gama, delta, t, fp);
t++;
}
break;
}
}//MCS
fclose(fp);
free_2d_int_pointer(viz,N,G);
} //main
Other functions
Gnuplot snapshots
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
void snap_gnuplot(int state[N], int topologia[N], double *investimento, int label[N], int t)
// ./a.out | gnuplot
{
// plota de cima pra baixo, esquerda pra direita
// gnuplot inverte baixo e cima
int i,j;
printf("set title \"MCS = %d\" \n",t);
printf("set autoscale keepfix\n");
printf("set palette model RGB\n");
printf("unset border\n");
//printf("unset colorbox\n");
printf("unset xtics\n");
printf("unset ytics\n");
//printf("set palette defined ( 0 \"dark-red\", 0.5 \"light-red\", 1 \"#0000B3\", 2 \"#000057\")\n");
printf("set palette defined ( 0 \"red\", 1 \"blue\", 2 \"yellow\")\n");
//printf("set palette defined ( 0 \"#a6611a\", 0.5 \"#dfc27d\", 1 \"#80cdc1\", 2 \"#018571\")\n");
printf("set cbrange[0:2]\n");
printf("set xrange[0:%d]\n",L);
printf("set yrange[0:%d]\n",L);
printf("set size square\n");
printf("plot \"-\" matrix with image\n");
for(i=0;i<L;i++)
{
for(j=0;j<L;j++)
{
printf("%d ",state[j+i*L]);
}
/*for(j=0;j<L;j++)
{
printf("%d ",state[j+i*L]);
}*/
printf("\n");
}
printf("\n");
printf("e\n"); printf("pause(%lf)\n",fps);
}
Gnuplot gif
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
void snap_gif(int state[N], int topologia[N], int t)
{
int i,j;
for(i=0;i<L;i++)
{
for(j=0;j<L;j++)
{
printf("%d ",state[j+i*L]);
}
/*for(j=0;j<L;j++)
{
printf("%d ",topologia[j+i*L]);
}*/
printf("\n");
}
printf("\n");
}
Hoshen-Kopelman
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
void calculo_percolacao( int state[N], int **viz, int label[N])
{
int largest_label = 0;
int n,i;
for(n=0;n<N;n++)
{
label[n] = state[n]*(n+1);
}
for(n=0;n<N;n++)
{
if(state[n]==1)
{
if(state[viz[n][1]] == 0 && state[viz[n][4]] != 0)
{
int a = label[n];
label[n] = label[viz[n][4]];
for(i=0;i<N;i++)
{
if(label[i]==a)
{label[i]=label[n];}
}
}
else if(state[viz[n][1]] != 0 && state[viz[n][4]] == 0)
{
int a = label[n];
label[n] = label[viz[n][1]];
for(i=0;i<N;i++)
{
if(label[i]==a)
{label[i]=label[n];}
}
}
else if(state[viz[n][1]] != 0 && state[viz[n][4]] != 0)
{
int a = label[viz[n][1]];
int b = label[n];
label[viz[n][1]] = label[viz[n][4]];
label[n] = label[viz[n][4]];
for(i=0;i<N;i++)
{
if(label[i]==a || label[i]==b)
{label[i]=label[n];}
}
}
}
}
}
Other Initial distributions
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
void calculo_ci_square1(int strategy[N], int **viz)
{
int n;
for(n=0; n < N; n++)
{
int i,j;
j=n%L;
i=n/L;
//if(j<L/2){strategy[j+i*L] = COOPERATOR; investimento[j+i*L] = 0.5;}
//else{strategy[j+i*L] = COOPERATOR; investimento[j+i*L] = 2.0;}
strategy[j+i*L] = 0; investimento[j+i*L] = 0.0;
//if(j==L/3 || j==2*L/3 || i==L/3 || i==2*L/3){strategy[j+i*L] = COOPERATOR;investimento[j+i*L] = 0.;}
if(j>L/3 && j<2*L/3 && i>L/3 && i<2*L/3){strategy[j+i*L] = 1;investimento[j+i*L] = 1.;}
}
}
void calculo_ci_square2(int strategy[N], int **viz)
{
int n;
for(n=0; n < N; n++)
{
int i,j;
j=n%L;
i=n/L;
strategy[n] = 0;
if(j>1*L/5 && j<2*L/5 && i>1*L/5 && i<2*L/5){strategy[j+i*L] = 1;investimento[j+i*L] = 0.55;}
if(j>3*L/5 && j<4*L/5 && i>1*L/5 && i<2*L/5)
{
double temp = gsl_rng_uniform(rand_vec);
if(temp < 0.5) {investimento[j+i*L] = 3.25;strategy[j+i*L] = 1;}
else{investimento[j+i*L] = 3.25;strategy[j+i*L] = 1;}
}
if(j>1*L/5 && j<2*L/5 && i>3*L/5 && i<4*L/5)
{
double temp = gsl_rng_uniform(rand_vec);
if(temp < 1./3) {strategy[j+i*L] = 1;investimento[j+i*L] = 0.;}
else if(temp < 2./3) {strategy[j+i*L] = 1;investimento[j+i*L] = 2.45;}
else{strategy[j+i*L] = 1;investimento[j+i*L] = 3.25;}
}
if(j>3*L/5 && j<4*L/5 && i>3*L/5 && i<4*L/5){strategy[j+i*L] = 1;investimento[j+i*L] = 2.45;}
}
}
void calculo_ci_stripes(int strategy[N], int **viz)
{
int n;
for(n=0; n < N; n++)
{
int i,j;
j=n%L;
i=n/L;
if(i<L/4) {strategy[j+i*L] = 0;investimento[j+i*L] = 0;}
if(i>=L/4 && i<2*L/4) {strategy[j+i*L] = 1;investimento[j+i*L] = 0.4;}
if(i>=2*L/4 && i<3*L/4) {strategy[j+i*L] = 1;investimento[j+i*L] = 0.7;}
if(i>=3*L/4) {strategy[j+i*L] = 1;investimento[j+i*L] = 0.0;}
}
}
This post is licensed under CC BY 4.0 by the author.