Santos-Pacheco-Lenaerts model

/* 繋ぎ変えメカニズムをもつジレンマモデル
   Santos-Pacheco-Lenaerts model */
/* 描画ライブラリとして plplotが必要 */

#include "plConfig.h"
#include "plplot.h"
#include "plevent.h"
#include <unistd.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define AGENT_NUM 100
#define LINK_NUM 500

#define XMAX 4.0  /* glplotのウィンドウサイズ */
#define YMAX 4.0
#define BGCOLOR "FFFFFF"  /* glplotの背景色 */
#define TAU_A 1.0
#define TAU_E 1.0

typedef struct agent_type{
  int attri;
  double payoff;
  int linkto[AGENT_NUM];
  int c_max;

  double c;
} AGENT;
double x[AGENT_NUM][1],y[AGENT_NUM][1];


int main(){
  double te,ta,y,w,z,c,M,K,O,T,S,P,R,Z,Q,pl,p,r; /* 変数定義としては悪い見本！ まねしないように*/
  int i=0,n,m,k,j,o,l,q,f,v,B,L,h,mk,b,x,lm,d,e,V,C,F,D,ok,ck,ak,cd;
  int myIDX_in_partner_linkto, partnerIDX_in_my_linkto,IDX_new_in_partner_linkto;
  AGENT a[AGENT_NUM];
  int iteration,checkFlag, me, partner, both_unsatisfied=0;
  double param_w, prob_choose_strategy;  /* W */

  srand((unsigned) time(NULL));
  /* T: Temptation, P: punishiment, S: sucker's, R: Reward for mutual coop. */

  param_w = TAU_E/TAU_A;
  printf("W= %f\n", param_w);

  prob_choose_strategy = 1.0 / (1+ param_w);
  R=1;P=0; /* ペアの利得定数 */

  T=0.0;
  while(T < 2.0){
    T+=0.2;
    printf("\n利得G(T):%.2f\n",T);

    S=-1.0;
    while(S <1.0){
      S+=0.2;
      printf("利得S:%.2f  ",S);


      for(i=0; i<AGENT_NUM;i++){
        a[i].c_max=-1;
      }
      for (i=0;i<AGENT_NUM;i++){
        a[i].attri =(int)(0.5 + (double)random()/RAND_MAX); // CかDか

        do {  // 各ノードについて最初のリンク１本をつくる リンクが重複しないように
          checkFlag=0;
          n = (int) random()%AGENT_NUM;
          for(j=0;j<=a[i].c_max;j++) if(a[i].linkto[j]==n) checkFlag++;
        } while(checkFlag!=0);
        a[i].c_max++;
        a[i].linkto[a[i].c_max] = n;
        a[n].c_max++;
        a[n].linkto[a[n].c_max] = i;
        if (1.0<=a[i].attri){e++;}

      }

      for(i=0;i<LINK_NUM-AGENT_NUM;i++){
        n = random()%AGENT_NUM;
        m = n;

        do {
          checkFlag = 0;
          m = random()%AGENT_NUM;
          for(j=0;j<=a[n].c_max;j++) if(a[n].linkto[j]==m) checkFlag++;
          for(j=0;j<=a[m].c_max;j++) if(a[m].linkto[j]==n) checkFlag++;
        } while(checkFlag!=0);

        a[n].c_max++;
        a[n].linkto[a[n].c_max]=m;
        a[m].c_max++;
        a[m].linkto[a[m].c_max]=n;
      }
      /* 以上、初期のランダムグラフ作成 */
      /*for(ck=0;ck<=AGENT_NUM;ck++){printf("つながりの数:%d;%d\n",ck,a[ck].c_max);}
       */
      //      printf("SGAP %f,%f,%f,%f\n",S,G,A,P);

      /*　初期状態の描画 */
      //                        drawInit();
      //                        drawAll(a);
      /*******************/

      /* 時間発展の繰り返し */
      for(iteration=0;iteration<1000000;iteration++) {/*printf("%d回目\n",i+1);*/
        /* ひとつのノードをランダムに選び出す*/
        k = rand()%AGENT_NUM;

        if(a[k].c_max==0) continue; /* リンクが１本になったらゲームはしない */

        /* それとつながった相手をランダムに選ぶ */
        partnerIDX_in_my_linkto  = rand() % a[k].c_max;
        f = a[k].linkto[partnerIDX_in_my_linkto];

        /*  printf("Aの戦略:%d  ", a[k].attri );*/

        /* attri=1: cooperator,  attri=0: defector */
        a[k].payoff=0;
        for(j=0;j<=a[k].c_max;j++){
          if(a[k].attri == 1 &&  a[a[k].linkto[j]].attri == 0)
            { a[k].payoff+=S;}
          else if(a[k].attri == 0 && a[a[k].linkto[j]].attri==1)
            { a[k].payoff+=T;}
          else if(a[k].attri==0 &&  a[a[k].linkto[j]].attri==0)
            { a[k].payoff+=P;}
          else if(a[k].attri==1 &&  a[a[k].linkto[j]].attri==1)
            { a[k].payoff+=R;
            } else {
            printf("error %d %d %d %d %d\n",k,j,a[j].c_max,a[k].attri,a[j].attri);
          }
        }
        r=a[k].payoff;
        /* printf("Aの総利得:%.001f\n",a[k].payoff);*/

        /* printf("Bの戦略:%d  ", a[f].attri );*/

        a[f].payoff=0;
        for(v=0;v<=a[f].c_max;v++){
          if(a[f].attri == 1 &&  a[a[f].linkto[v]].attri ==0)
            { a[f].payoff+=S;}
          else if(a[f].attri ==0 && a[a[f].linkto[v]].attri==1)
            { a[f].payoff+=T;}
          else if(a[f].attri==0 &&  a[a[f].linkto[v]].attri==0)
            { a[f].payoff+=P;}
          else if(a[f].attri==1 &&  a[a[f].linkto[v]].attri==1)
            { a[f].payoff+=R;
            } else {
            printf("error %d %d %d %d %d\n",f,v,a[f].c_max,a[f].attri,a[v].attri);
          }
        }
        p=a[f].payoff;
        /* printf("Bの総利得%.001f\n",a[f].payoff);*/

        y=-0.005*(p-r);
        if(abs(y)>100.0)        y=100.0;
        w=1.0/(1.0+exp(y));  /* w = 1.0/(1 + exp(-\beta * (Pi(B) - Pi(A)) */



        if((double)random()/(double)RAND_MAX < prob_choose_strategy) {
          /****** 戦略の変更を選択  *******/
          if((double)random()/(double)RAND_MAX < w) {
            a[k].attri = a[f].attri;  /* attri=0なら1に，1なら0に */
          }

        } else if(a[f].c_max!=0) {
          /***** つなぎ替えを選択 *******/
          both_unsatisfied =0;
          if(a[k].attri ==1 ) { /*自分がCooperatior */
            if( a[f].attri == 1) continue; /* 相手がCooperatorなら何もしない．両方満足 */
            else {
              /* 相手がDefectorなら，自分のつなぎ替えを試みる */
              me = k; partner = f;
            }
          } else { /* 自分がDefectorの場合 */
            if(a[f].attri == 1) continue; /* 相手がCooperatorなら何もしない */
            else {
              /*　どちらもDefector(不満)の場合は，どちらかが確率的につなぎ替え */
              both_unsatisfied = 1;
              if((double)random()/(double)RAND_MAX < w) {
                me = k; partner = f;
              } else {
                me = f; partner = k; /* partnerの方を入れ替えるため，meとpartnerを交代 */
                for(b=0; b <= a[me].c_max; b++) {
                  if(a[me].linkto[b]==partner) {
                    partnerIDX_in_my_linkto = b;
                    break;
                  }
                }
                if(b>a[me].c_max) {printf("見つからない!!%d, %d",b,a[me].c_max);
                  exit(0);}
              }
            }
          }
          /* 以降，つなぎ替えを試みる方をme，切られる相手をpartnerと記す */
          v = 0; cd=0;
          /*　新しい友達探し */
          while(v < a[partner].c_max*2) {  /* 無限ループにおちいらいないために　c_maxの２倍繰り返したらやめる　*/
            IDX_new_in_partner_linkto = rand() % a[partner].c_max;
            B = a[partner].linkto[IDX_new_in_partner_linkto];
            if(B!=me) { /* 再び自分が選ばれることを避ける */
              for (b=0; b<=a[me].c_max; b++){ /* すでに自分の友達かどうかチェック．リンクの重複を避ける */
                if(B == a[me].linkto[b]) cd++;
              }
              if(cd==0) break; /* 避けられた -> 新しい友達候補がみつかったらループから抜ける */
            }
            v++;
          }

          if(v < a[partner].c_max*2) { /* 新たな友達候補が見つかったときだけ */

            if(both_unsatisfied==1 || (double)random()/(double)RAND_MAX < w){ /* 確率的につなぎ替え */

              for(b=0; b<=a[partner].c_max;b++) {
                if(a[partner].linkto[b] == me) {
                  myIDX_in_partner_linkto = b;
                  break;
                }
              }
              /* test 出力 */
              /* printf("\nme=%d, f=%d,B=%d\n", me, partner,B);
                 for(j=0;j<=a[me].c_max;j++) printf("%d, ",a[me].linkto[j]);
                 printf("\n");
                 for(j=0;j<=a[partner].c_max;j++) printf("%d, ",a[partner].linkto[j]);
              */

              if(b>a[partner].c_max) {
                printf("相手のlinktoに自分の登録が見つかりませんでした．me=%d\n",me);
                //                                                              for(j=0;j<=a[partner].c_max;j++) printf("%d, ",a[partner].linkto[j]);
                exit(0);
              }
              /*相手のlinktoから自分の登録を抹消 */
              L = myIDX_in_partner_linkto;
              while(L<a[partner].c_max) {
                a[partner].linkto[L]=a[partner].linkto[L+1];
                L++;
              }
              a[partner].c_max--;

              /* 自分のlinktoに新しい友達を登録 (上書きにより既存の相手を抹消) */
              a[me].linkto[partnerIDX_in_my_linkto] = B;

              /* 新しい友達のlinktoへ自分を登録 */
              a[B].linkto[a[B].c_max+1]=me;
              a[B].c_max++;
            }
            //                                                  f = B; /** ???? これは何？　*/
          }
          /*BににつながったCをランダムに選ぶ*/
        } /* つなぎ替えの終わり */
        if(i%100==0) { /* 描画 */
          //                                    printf("i=%d\n",i);
          //                                    drawAll(a);
        }
      } /* １回の試行(r <= p) の終わり */

      /*printf("%d\n",e);*/ /*for(ak=0;ak<=AGENT_NUM;ak++){printf("数%d;%d\n",ak,a[ak].c_max);}*/
      for (b=0;b<AGENT_NUM;b++){if( a[b].attri==1){Z++;}}

      printf("%.5f \n",(double)Z/AGENT_NUM);
      Z=0;
    }
  }

  plend();
  exit(0);
}

int drawInit(AGENT a[]){
  int i,j;
  plsdev("xwin");
  plsetopt("geometry","700x700");
  plsetopt("bg",BGCOLOR);
  //  plinit();
  plstar(1,1);

  for(i=0;i<AGENT_NUM; i++){
    x[i][0] = (double) random()/RAND_MAX*XMAX;  /* AGENT_NUM個のノード位置 */
    y[i][0] = (double) random()/RAND_MAX*YMAX;
  }

}

int drawAll(AGENT a[], int new_page) {

  int i,j;
  int col[2] = {1,9}; /* マークの色 */

  plenv(0,XMAX,0,YMAX,1,-2);

  plcol0(1);
  plssym(0.,2.0);

  for(i=0;i<AGENT_NUM; i++){
    plcol0(col[a[i].attri]);
    plpoin(1,x[i],y[i],a[i].attri+16);
  }

  plcol0(5);
  for(i=0;i<AGENT_NUM;i++){
    j=0;
    while(j<=a[i].c_max){
      pljoin(x[i][0],y[i][0],x[a[i].linkto[j]][0],y[a[i].linkto[j]][0]);
      j++;
    }
  }
  /*  plend(); */
  return(0);
}

/* 変化するノードを目立つように大きくし，前の注目ノードを元に戻す */
int drawChangeNode(int n, int attri, int prev_n, int prev_attri) {
  int col[2] = {1,9}; /* マークの色 */

  plssym(0,2.0); /* 一回り大きく */

  plcol0(col[attri]);
  plpoin(1,x[n],y[n],attri+16);

  plcol0(15);  /* 15は白 */
  plpoin(1,x[prev_n],y[prev_n],prev_attri+16); /* 前のを普通の大きさに戻す */

  plcol0(col[prev_attri]);
  plpoin(1,x[prev_n],y[prev_n],prev_attri+16);

}

/* リンク線を消す */
int plotDelLink(int m, int del_n) {/*２つのAGENT番号を受け取って消す */
  int col[2] = {1,9}; /* マークの色 */
  plcol0(15);
  pljoin(x[m][0],y[m][0],x[del_n][0],y[del_n][0]);
}

/* 新たにリンク線を張る */
int plotNewLink(int m, int new_n) { /*  */
  int col[2] = {1,9}; /* マークの色 */
  plcol0(5);
  pljoin(x[m][0],y[m][0],x[new_n][0],y[new_n][0]);
}