/*************************************************/
/* システム工学の基礎 5章の問題4の解答プログラム */
/* コンパイル方法 $gcc q5_4.c -lm                */
/*************************************************/

#include <stdio.h>
#include <math.h>

const double myu = 0.01;
const double kai = 0.7; 

/* |i| */
const int babs[4] = { 0, 1, 1, 2};

/* Z行列 Z(j,k) = Z[k][j] */
const int Z[10][4] = {
  0,0,0,2,
  0,0,1,1,
  0,0,2,0,
  0,1,0,1,
  0,1,1,0,
  0,2,0,0,
  1,0,0,1,
  1,0,1,0,
  1,1,0,0,
  2,0,0,0,
};

double Q(int k, int v);
double r(int j, int k);
double (*p)(int, int);
double R(int i, int h);

/* ランダム選択の選択確率 */
double select_rand(int j, int k){
  return Z[k][j]/2.0;
}

/* 1MAX問題におけるルーレット選択確率 */
double select_1max(int j, int k){
  int i;
  double m = 0.0;

  for(i=0; i<4; i++){
    m += Z[k][i] * babs[i];
  }
  if(m == 0.0){
    return select_rand(j ,k);
  }

  return Z[k][j]*babs[j]/m;
}

int main(int argc, char *argv[]){
  int k, v;

  p = select_rand;
  printf("random\n");
  for(k=0; k<10; k++){
    for(v=0; v<10; v++){
      printf("%1.2e%s", Q(k,v), v == 9 ? "\n" : ", ");
    }
  }

  p = select_1max;
  printf("1max\n");
  for(k=0; k<10; k++){
    for(v=0; v<10; v++){
      printf("%1.2e%s", Q(k,v), v == 9 ? "\n" : ", ");
    }
  }

  return 0;
}

double Q(int k, int v){
  int j, M = 2; /* 集団数 M! = 2! = 2 */
  double q = M;

  for(j=0; j<4; j++){
    if(Z[v][j] > 0){
      q *= pow(r(j,k),Z[v][j])/Z[v][j];
    }
  }

  return q;
}

double r(int j, int k){
  int i, h;
  double r = 0.0;

  for(i=0; i<4; i++){
    for(h=0; h<4; h++){
      r += p(j^i,k)*p(j^h,k)*R(i,h);
    }
  }

  return r;
}

double R(int i, int h){
  int L = 2; /* 遺伝子型の長さ L = 2 */
  double M,M1,M2,C,C1,C2;

  M1 = pow(1.0-myu,L-babs[i])*pow(myu,babs[i]);
  M2 = pow(1.0-myu,L-babs[h])*pow(myu,babs[h]);
  M = (1.0-kai)*(M1+M2)/2.0;
  C1 = pow(myu,babs[1|h]+babs[i]-babs[1|i])*pow(1.0-myu,L-babs[i]+babs[1|i]-babs[1|h]);
  C2 = pow(myu,babs[1|i]+babs[h]-babs[1|h])*pow(1.0-myu,L-babs[h]+babs[1|h]-babs[1|i]);
  C = kai*(C1+C2)/2.0;

  return M+C;
}