EGT
Main code to run simulations about evolutionary game theory.
Packages
1
2
3
4
5
6
7
8
9
10
11
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "pointers.h"
#include "mc.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
/********************************************************************
*** C.I. ***
********************************************************************/
#define dens_aleat
#define dens_inic_c 0.5
#define dens_inic_p 0.
#define dens_inic_pc 0.
#define dens_inic_dc 0.
//#define dens_quad
//#define dens_quad2
//#define dens_listras
/********************************************************************
*** Plot ***
********************************************************************/
#define fps 0000.01000000005 // fps^-1
//#define densidade_arquivo
#define densidade_terminal
/********************************************************************
*** Jogo ***
********************************************************************/
//#define PRISONERS_DILEMMA
#define PGG_FOCAL
//#define PGG
Global variables
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
//------------------------------------------------------------------------------------
enum tipo_rede {UNIDIMENSIONAL,QUADRADA, CUBICA, QUADRIDIMENSIONAL, HEXAGONAL, KAGOME,
TRIANGULAR, MOORE};
enum tipo_rede REDE_ATUAL;
enum tipo_estrategia {DEFECTOR, COOPERATOR, PUNISHER};
double RUIDO;
double INV_RUIDO;
#define prob_mobil 0.
#define EPSILON (1e-8)
const int tmax=99999;//5000;
double measure_time = 1;
long long int N;
int L;
long int L2;//= L*L;
long int L3;// = L*L*L;
long long int L4;// = L*L*L*L;
int G; // =4 (rede hexagonal) =5 (rede quadrada / kagome) =7 (rede triangular / quadrada 3D) =9 (rede moore / quadrada 4D)
int ND=0;
int NC=0;
int NP=0;
int NPC=0;
int NDC=0;
int nd=0;
int nc=0;
int np=0;
int ndc=0;
int npc=0;
double investimento_total;
double payoff_total_C;
double payoff_total_D;
double payoff_total_resto;
//------------------------------------------------------------------------------------
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);
return;
}
Initial distributions/conditions
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
/********************************************************************
*** distribuição dos estados ***
********************************************************************/
void calculo_ci_estado(int state[N], int **viz, double *investimento)
{
int n;
for(n=0; n < N; n++)
{
#ifdef dens_aleat
double temp = gsl_rng_uniform(rand_vec);
if(temp < dens_inic_c ) state[n] = 1;
else
{
if(temp < dens_inic_c + dens_inic_p) state[n] = 2;
else
{
if(temp < dens_inic_c + dens_inic_p + dens_inic_pc ) state[n] = 3;
else
{
if(temp < dens_inic_c + dens_inic_p + dens_inic_pc + dens_inic_dc) state[n] = 4;
else state[n] = 0;
}
}
}
#endif
#ifdef dens_quad
int i,j;
j=n%L;
i=n/L;
state[n] = 0;
if(j>1*L/5 && j<2*L/5 && i>1*L/5 && i<2*L/5){state[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;state[j+i*L] = 1;}
else{investimento[j+i*L] = 3.25;state[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) {state[j+i*L] = 1;investimento[j+i*L] = 0.;}
else if(temp < 2./3) {state[j+i*L] = 1;investimento[j+i*L] = 2.45;}
else{state[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){state[j+i*L] = 1;investimento[j+i*L] = 2.45;}
#endif
#ifdef dens_quad2
int i,j;
j=n%L;
i=n/L;
//if(j<L/2){state[j+i*L] = COOPERATOR; investimento[j+i*L] = 0.5;}
//else{state[j+i*L] = COOPERATOR; investimento[j+i*L] = 2.0;}
state[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){state[j+i*L] = COOPERATOR;investimento[j+i*L] = 0.;}
if(j>L/3 && j<2*L/3 && i>L/3 && i<2*L/3){state[j+i*L] = 1;investimento[j+i*L] = 1.;}
#endif
#ifdef dens_listras
int i,j;
j=n%L;
i=n/L;
if(i<L/4) {state[j+i*L] = 0;investimento[j+i*L] = 0;}
if(i>=L/4 && i<2*L/4) {state[j+i*L] = 1;investimento[j+i*L] = 0.4;}
if(i>=2*L/4 && i<3*L/4) {state[j+i*L] = 1;investimento[j+i*L] = 0.7;}
if(i>=3*L/4) {state[j+i*L] = 1;investimento[j+i*L] = 0.0;}
#endif
//state[n] = 0;investimento[n] = 0;
}
return;
}
/********************************************************************
*** distribuição da topologia ***
********************************************************************/
void calculo_ci_topologia(int topologia[N])
{
int i;
for(i=0;i<N;i++)
{
topologia[i] = G;
//topologia aleatoria
/*double temp = gsl_rng_uniform(rand_vec);
if( temp < 1/3. ) topologia[n] = QUADRADA; //ou KAGOME
else
{
if( temp < 2/3.) topologia[n] = TRIANGULAR; // ou CUBICA
else
topologia[n] = MOORE;
}*/
//printf("%d\n",topologia[i]);
}
}
//#define M_PI
/********************************************************************
*** distribuição da contribuiçao ***
********************************************************************/
void calculo_ci_investimento(double *investimento, int state[N], double gama)
{
int i;
for(i=0;i<N;i++)
{
double ga1 = gsl_rng_uniform(rand_vec);
double ga2 = gsl_rng_uniform(rand_vec);
double gaussian = 4 + 1*sqrt(-2*log(ga1))*cos(2*M_PI*ga2); // media 3 desvio 1
int k = 1 + (int) 3*gsl_rng_uniform(rand_vec); // X*gsl_rng_uniform(rand_vec) vai de 0 até X-1
double kk = 5*gsl_rng_uniform(rand_vec);
if(k==1){kk=0.5;}
if(k==2){kk=1.0;}
if(k==3){kk=2.0;}
if(state[i] == 0){investimento[i] = 0.0;}
if(state[i] == 1){investimento[i] = 1.0;}
if(state[i] == 2){investimento[i] = 1.0;}
/*int jj,j;
j=i%L;
jj=i/L;
if(jj<L/3){investimento[j+jj*L] = 1;}
if(jj>=L/3 && jj<2*L/3){investimento[j+jj*L] = k;}
if(jj>=2*L/3){investimento[j+jj*L] = 0;}*/
//printf(" %d %d %lf\n",i,state[i],investimento[i]);
}
}
Functions
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
72
73
74
75
76
77
78
79
80
81
82
83
84
/********************************************************************
*** Payoff ***
********************************************************************/
void calculo_payoff ( double *payoff, double r, double gama, double delta, int x, int topologia[N], int state[N], double *investimento, int **viz )
{
// x: sitio central do grupo que o jogo ocorre
// y: sitio que esta jogando
// ***********************************************************************************
// calcula numero de cada estrategia
nd=0;
nc=0;
np=0;
ndc=0;
npc=0;
investimento_total = 0.;
int i=0;
for(i=0;i<topologia[x];i++) // vizinho 0 e o proprio sitio
{
switch(state[viz[x][i]])
{
case 0: ++nd; break;
case 1: ++nc; break;
case 2: ++np; break;
case 3: ++npc; break;
case 4: ++ndc; break;
default:
fprintf(stderr,"ERRO calculo vizinhos\n");
fflush(stderr);
}
investimento_total += investimento[viz[x][i]];
}
// ***********************************************************************************
// calcula payoff
#if defined(PGG) || defined(PGG_FOCAL)
double pool = (r/topologia[x])*investimento_total;
payoff[1] = pool; //- investimento[x];
payoff[0] = pool; //- investimento[x];
#endif
#ifdef PRISONERS_DILEMMA
// lembrar de mudar o nc pra ter ou n autointeraçao
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;
//pure cordination game: P = R > S = T
// 1 = C = left
// 0 = D = right
//Snowdrift / Hawk and Dove / chicken game: T > R > S > P
// 1 = Dove, deviate
// 0 = Hawk, straight
//pure cordination game: P = R > S = T
// 1 = C = left
// 0 = D = right
//stag hunt game: R > T >= S*nd;
// 1 = C = stag
// 0 = D = rabbit
#endif
return;
}
Total payoff
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
/********************************************************************
*** payoff total ***
********************************************************************/
void calculo_payoff_total(double *payoff, int state[N], int **viz, double r, double gama, double delta, int topologia[N], double *investimento)
{
int i;
payoff_total_C = 0.;
payoff_total_D = 0.;
payoff_total_resto = 0.;
for(i=0;i<N;i++)
{
calculo_payoff(payoff,r,gama,delta,i,topologia,state,investimento,viz);
if(state[i]==1){payoff_total_C += payoff[state[i]];}
if(state[i]==0){payoff_total_D += payoff[state[i]];}
else{payoff_total_resto += payoff[state[i]];}
//printf("%d %d %lf %lf %lf %lf\n", state[i], topologia[i], investimento_total, payoff[state[i]], payoff_total_C, payoff_total_D);
}
}
Update rule
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
/********************************************************************
*** Update rule ***
********************************************************************/
void update_rule( int x, int vizinho, int state[N], double Px, double Py ,int topologia[N], double *investimento, 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)
{
state[x] = state[vizinho];
topologia[x] = topologia[vizinho];
investimento[x] = investimento[vizinho]; // tirar pra fixa investimento com local da rede
//state[x] = fabs(state[x]-1); //Ising
}
}
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
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
/********************************************************************
*** MCS ***
********************************************************************/
void calculo_mcs(double *payoff, int state[N], int **viz, double r, double gama, double delta, int topologia[N], double *investimento, int t)
{
int n;
int x,y,sitio,vizinho;
double Px, Py;
int sitio2, vizinho2;
for(n=0; n<N; n++)
{
Px = 0;
Py = 0;
x = (int) N*gsl_rng_uniform(rand_vec);
sitio = x;
y = 1 + (int) (topologia[x]-1)*gsl_rng_uniform(rand_vec);
vizinho = viz[x][y];
if (investimento[vizinho] != investimento[x] || state[vizinho] != state[x])
{
#if defined(PGG_FOCAL)
calculo_payoff(payoff,r,gama,delta,sitio,topologia,state,investimento,viz);
Px = payoff[state[sitio]] - investimento[sitio];
calculo_payoff(payoff,r,gama,delta,vizinho,topologia,state,investimento,viz);
Py = payoff[state[vizinho]] - investimento[vizinho];
#endif
#ifdef PGG
for(y=0;y<topologia[sitio];y++)
{
sitio2 = viz[sitio][y];
calculo_payoff(payoff,r,gama,delta,sitio2,topologia,state,investimento,viz);
Px += payoff[state[sitio]] - investimento[sitio];
//printf("%d %d %d %d %lf \n", sitio, state[sitio],sitio2, state[sitio2],Px);
}
for(y=0;y<topologia[vizinho];y++)
{
vizinho2 = viz[vizinho][y];
calculo_payoff(payoff,r,gama,delta,vizinho2,topologia,state,investimento,viz);
Py += payoff[state[vizinho]] - investimento[vizinho];
//printf("%d %d %d %d %lf \n", vizinho, state[vizinho],vizinho2, state[vizinho2],Py);
}
#endif
#if defined(PRISONERS_DILEMMA)
calculo_payoff(payoff,r,gama,delta,sitio,topologia,state,investimento,viz);
Px = payoff[state[sitio]];
calculo_payoff(payoff,r,gama,delta,vizinho,topologia,state,investimento,viz);
Py = payoff[state[vizinho]];
#endif
update_rule(x,vizinho,state,Px,Py,topologia,investimento,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
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
73
/********************************************************************
*** Densidades ***
********************************************************************/
void calculo_densidades(int state[N], double *investimento, double *payoff, int **viz, int histograma[N], int topologia[N], double r, double gama, double delta, int t, FILE *fp)
{
int k;
//long double mediaC=0.;
//long double mediaC2=0.;
//long double mediaD=0.;
//long double mediaD2=0.;
ND=0;
NC=0;
NP=0;
NPC=0;
NDC=0;
for(k=0;k<N;k++)
{
switch (state[k])
{
case 1: ++NC; break;
case 2: ++NP; break;
case 0: ++ND; break;
case 3: ++NPC; break;
case 4: ++NDC; break;
default:
fprintf(stderr,"ERRO - tipo de estrategia\n");
fflush(stderr);
}
/*if(state[k] == 1)
{
mediaC += investimento[k];
mediaC2 += investimento[k]*investimento[k];
}
if(state[k] == 0)
{
mediaD += investimento[k];
mediaD2 += investimento[k]*investimento[k];
}*/
}
/*double media_C = (double)mediaC/(NC);
double desvio_C = (double)sqrt((0.0000000001+mediaC2-(mediaC*mediaC)/(NC))/NC);
double media_D = (double)mediaD/(ND);
double desvio_D = (double)sqrt((mediaD2-(mediaD*mediaD)/(ND))/ND);
if(NC == 0)
{
media_C = 0.;
desvio_C = 0.;
}
if(ND == 0)
{
media_D = 0.;
desvio_D = 0.;
}
calculo_payoff_total(payoff,state,viz,r,gama,delta,topologia,investimento);
double payoff_total = payoff_total_C + payoff_total_D + payoff_total_resto;
*/
#ifdef densidade_terminal
//printf("%d %lf %lf %lf %lf %lf %lf %lf\n", t, (double)NC/(N), (double)ND/(N), (double)NP/(N), media_C ,desvio_C, payoff_total_C, payoff_total_D);
printf("%d %lf %lf\n", t, (double)NC/(N), (double)ND/(N));
#endif
#ifdef densidade_arquivo
fprintf(fp,"%d %lf %lf %lf %lf %lf\n", t, (double)NC/(N), (double)ND/(N), (double)NP/(N), media_C ,desvio_C);
#endif
}
Snapshots (gnuplot)
just run the code with | gnuplot
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
/********************************************************************
*** Snapshots ***
********************************************************************/
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);
}
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
/********************************************************************
*** Main ***
********************************************************************/
int main(int argc, char *argv[])
{
// gcc dilema.c -lm -lgsl -lgslcblas -O3 ----
// ./a.out r gama delta L ruido
int **viz;
double payoff[5], gama, delta, r;
r=atof(argv[1]);
gama=atof(argv[2]);
delta=atof(argv[3]);
L=atoi(argv[4]);
RUIDO=atof(argv[5]);
set_gsl_rng(); // algo da 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
// --------------------------------------
// especificar dependendo da rede
N = 2*L*L;
G = 4;
viz = create_2d_int_pointer_h(N,G);
honeycomb_lattice(viz,L);
// --------------------------------------
int state[N], topologia[N], histograma[N], label[N];
double investimento[N];
//-------------------C.I.--------------------
calculo_ci_estado(state,viz,investimento);
calculo_ci_topologia(topologia);
calculo_ci_investimento(investimento,state,gama);
calculo_densidades(state,investimento,payoff,viz,histograma,topologia,r,gama,delta,0,fp);
//-------------------MCS----------------------
int t=0;
for(t=1; t < tmax; t++)
{
//snap_gnuplot(state, topologia,investimento,label,t);
/*
if(t==1 || t==40 || t==90 || t==120 || t==300){
printf("set terminal postscript eps enhanced color\n");
printf("set output 'snapshotC1_MCS%d.eps'\n",t);
snap_gnuplot(state, topologia,investimento,label,t);
printf("unset terminal\n");
fprintf(fp, "set terminal postscript eps enhanced color\n");
fprintf(fp, "set output 'snapshotC1_MCS%d.eps'\n",t);
snap_arquivo(state, topologia,investimento,label,t,fp);
}
*/
//--------------------------------------------------------------
calculo_mcs(payoff,state,viz,r,gama,delta,topologia,investimento,t);
calculo_densidades(state,investimento,payoff,viz,histograma,topologia,r,gama,delta,t,fp);
/*
if(t>70000)
{
calculo_percolacao(state,viz,label);
calculo_tamanho_cluster(state,viz,label,investimento,t);
}
*/
//calculo_mobilidade(state,viz);
//----------------------------------------------------------------
measure_time = rint(1.03*t); //printf("%d %lf\n",t, measure_time);
if((NC/(N)==1) || (NP/(N)==1) || (NPC/(N)==1) || (ND/(N)==1) )
{
do
{
if ( t >= measure_time)
{
calculo_densidades(state,investimento,payoff,viz,histograma,topologia,r,gama,delta,t,fp);
measure_time = rint(1.03*measure_time);
}
t=t+1;
}
while(t<tmax);
break;
}
}//MCS
fclose(fp);
return 0;
free_2d_int_pointer(viz,N,G);
} //main
This post is licensed under CC BY 4.0 by the author.