// PROGRAM genetic

#include "TrueBASIC.h"

#define MAXLAT 21
#define MAXL2  401
#define MAXPOP 1001

int main();
void initial(char s[][MAXL2], int *L, int *L2, int T[][3], int *npop, 
             int *nrecombine, int *nmutation, int *ngeneration);
void recombine(char s[][MAXL2], int L2, int npop, int *ntot);
void mutate(char s[][MAXL2], int L2, int npop, int *ntot);
void convert(char* a, int L, int ss[][MAXLAT]);
void energy(int L, int ss[][MAXLAT], int T[][3], double *E);
void selection(char s[][MAXL2], int L, int L2, int T[][3], int npop, 
               int ntot, double Eselect[]);
void showoutput(char s[][MAXL2], int npop, double Eselect[]);
void strsubcat(char s1[], char s2[], int p1, int p2);

int main()
{
    int ntot, L, L2, T[MAXL2][3], i, igeneration, iswap, ngeneration, nmutation, 
        npop, nrecombine, Itmp1_;
    double Eselect[MAXPOP];
    char s[MAXPOP][MAXL2];
	for(Itmp1_ = 0; Itmp1_ < MAXPOP; ++Itmp1_) {
		Eselect[Itmp1_] = 0;
			s[Itmp1_][0] = '\0';
	}

    initial(s, &L, &L2, T, &npop, &nrecombine, &nmutation, &ngeneration);
    for(igeneration = 1; igeneration <= ngeneration; ++igeneration)
    {
        ntot = npop;
        for(iswap = 1; iswap <= nrecombine; ++iswap)
            recombine(s, L2, npop, &ntot);
        for(i = 1; i <= nmutation; ++i)
            mutate(s, L2, npop, &ntot);
        selection(s, L, L2, T, npop, ntot, Eselect);
    }
    showoutput(s, npop, Eselect);

    return 0;
}

void initial(char s[][MAXL2], int *L, int *L2, int T[][3], int *npop, 
             int *nrecombine, int *nmutation, int *ngeneration)
{
    int i, ipop, j;
    char Stmp1_[_LBUFF_];

    rnd(-1);
    printf("string (lattice) size = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", L);
    (*L2) = (*L)*(*L);
    printf("population number = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", npop);
    // L は表現型の大きさ
    printf("number of recombinations per generation = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", nrecombine);
    printf("number of mutations per generation = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", nmutation);
    printf("number of generations = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", ngeneration);
    // 遺伝子型の数をランダムに生成する
    for(ipop = 1; ipop <= (*npop); ++ipop)
    {
        strcpy(s[ipop], "");
        for(i = 0; i < (*L2); ++i)
        {
            if(rnd(0) > 0.5)
                s[ipop][i] = '1';
            else
                s[ipop][i] = '0';
        }
		s[ipop][i] = '\0';
    }
    // ランダムな結合網 Tij を生成する
    for(i = 1; i <= (*L2); ++i)
    {
        for(j = 1; j <= 2; ++j)
        {
            if(rnd(0) > 0.5)
                T[i][j] = 1;
            else
                T[i][j] = -1;
        }
    }
}

void recombine(char s[][MAXL2], int L2, int npop, int *ntot)
// 2つの文字列（遺伝子型）を選んで組換えを行う
{
    int i, j, pos, size;
    int pbc;
    char s1[MAXL2], s2[MAXL2];

    i = (int)(npop*rnd(0)) + 1;
    do
    {
        j = (int)(npop*rnd(0)) + 1;
    } while(i == j);
    size = (int)(rnd(0)*(L2/2.0)) + 1;
    pos = (int)(rnd(0)*L2) + 1;
    strcpy(s1, s[i]);
    strcpy(s2, s[j]);
    if(pos + size <= L2)
    {
        strcpy(s[(*ntot) + 1], "");
        strsubcat(s[(*ntot) + 1], s1, 1, pos-1);
        strsubcat(s[(*ntot) + 1], s2, pos, pos+size);
        strsubcat(s[(*ntot) + 1], s1, pos+size+1, L2);
        strcpy(s[(*ntot) + 2], "");
        strsubcat(s[(*ntot) + 2], s2, 1, pos-1);
        strsubcat(s[(*ntot) + 2], s1, pos, pos+size);
        strsubcat(s[(*ntot) + 2], s2, pos+size+1, L2);
    }
    else                                    // 周期的境界条件を適用する
    {
        pbc = pos + size - L2;
        strcpy(s[(*ntot) + 1], "");
        strsubcat(s[(*ntot) + 1], s2, 1, pbc);
        strsubcat(s[(*ntot) + 1], s1, pbc+1, pos-1);
        strsubcat(s[(*ntot) + 1], s2, pos, L2);
        strcpy(s[(*ntot) + 2], "");
        strsubcat(s[(*ntot) + 2], s1, 1, pbc);
        strsubcat(s[(*ntot) + 2], s2, pbc+1, pos-1);
        strsubcat(s[(*ntot) + 2], s1, pos, L2);
    }
    (*ntot) += 2;
}

void mutate(char s[][MAXL2], int L2, int npop, int *ntot)
{
    int i, pos;
    char c;

    i = (int)(rnd(0)*npop) + 1;
    pos = (int)(rnd(0)*L2) + 1;
    c = s[i][pos-1];
    if(c == '1')
        c = '0';
    else
        c = '1';
    strcpy(s[(*ntot) + 1], "");
    strsubcat(s[(*ntot) + 1], s[i], 1, pos - 1);
    s[(*ntot) + 1][pos - 1] = c;
    s[(*ntot) + 1][pos] = '\0';
    strsubcat(s[(*ntot) + 1], s[i], pos + 1, L2);
    ++(*ntot);
}

void convert(char* a, int L, int ss[][MAXLAT])
{
    int i, j, n;

    // 0 と 1 の文字列を2次元のスピン配列に変換する
    // つまり遺伝子型を表現型に変換する
    for(i = 1; i <= L; ++i)
    {
        for(j = 1; j <= L; ++j)
        {
            n = (j-1)*L + i;
            if(a[n-1] == '1')
                ss[i][j] = 1;
            else
                ss[i][j] = -1;
        }
    }
}

void energy(int L, int ss[][MAXLAT], int T[][3], double *E)
{
    int i, ip, j, jp, n;

    (*E) = 0;
    for(i = 1; i <= L; ++i)
    {
        ip = i + 1;
        if(ip > L)
            ip = 1;
        for(j = 1; j <= L; ++j)
        {
            jp = j + 1;
            if(jp > L)
                jp = 1;
            n = (j-1)*L + i;
            (*E) -= T[n][1]*ss[i][j]*ss[ip][j] - T[n][2]*ss[i][j]*ss[i][jp];
        }
    }
}

void selection(char s[][MAXL2], int L, int L2, int T[][3], int npop, int ntot, 
               double Eselect[])
{
    int i, ipop, ss[MAXLAT][MAXLAT], Itmp1_;
    double E, Elist[MAXPOP], Esum;
    char save[MAXPOP][MAXL2];

    Esum = 0;
    Elist[0] = 0;
    for(i = 1; i <= ntot; ++i)
    {
        convert(s[i], L, ss);
        energy(L, ss, T, &E);
        Esum += - E + 2*L2;                   // Esum > 0 への寄与
        Elist[i] = Esum;
    }
    for(Itmp1_ = 0; Itmp1_ < MAXPOP; ++Itmp1_)
        strcpy(save[Itmp1_], s[Itmp1_]);
    // 新しい種族を選択する
    for(ipop = 1; ipop <= npop; ++ipop)
    {
        E = Esum*rnd(0);
        i = 0;
        do
        {
            ++i;
        } while(!(E < Elist[i]));           // エネルギーに応じて選択する
        strcpy(s[ipop], save[i]);
        Eselect[ipop] = Elist[i - 1] - Elist[i] + 2*L2;
    }
}

void showoutput(char s[][MAXL2], int npop, double Eselect[])
{
    int ipop;

    for(ipop = 1; ipop <= npop; ++ipop)
        printf("%12d %12s %12f\n", ipop, s[ipop], Eselect[ipop]);
    printf("\n");
}

void strsubcat(char s1[], char s2[], int p1, int p2)
{
    char Stmp1_[MAXL2];
    int i;
    if(p1 > p2) return;
    for(i = 0; i <= p2 - p1; ++i)
        Stmp1_[i] = s2[p1-1+i];
    Stmp1_[i] = '\0';
    strcat(s1, Stmp1_);
}