repchar(outfp," ",lchrom-5);
        fprintf(outfp,"適應度\n");
        repchar(outfp,"-",80);
        skip(outfp,1);
        writepop(outfp);
        repchar(outfp,"-",80);
        skip(outfp,1);
     }
    fprintf(outfp," 第 %d 代統計: \n",gen);
    fprintf(outfp," 總交叉操作次數 = %d, 總變異操作數 = %d\n",ncross,nmutation);
    fprintf(outfp," 子群中最小適應度：%f 最大適應度：%f  平均適應度 %f\n", min,max,A[subpop]);
    fprintf(outfp," 子群中發現最佳個體 =>  所在代數: %d  ", bestfit.generation);
    fprintf(outfp," 適應度：%f  染色體: ", bestfit.fitness);
    writechrom((&bestfit)->chrom);
    fprintf(outfp," 對應的變量值: %f", bestfit.varible);
	skip(outfp,1);
	fprintf(outfp," 迄今發現最佳個體   =>  所在分層: %d 所在子群: %d  所在代數: %d ", runA,bestfit.subpopulation,genA);
	fprintf(outfp," 適應度：%f  染色體: ", maxA);
    writechrom((&bestfit)->chromA);
	fprintf(outfp," 對應的變量值: %f", varibleA);
	skip(outfp,1);
	fprintf(outfp," 迄今發現最小適應度 => %f ",minA);
    skip(outfp,1);
    repchar(outfp,"-",80);
    skip(outfp,1);  
}

void writepop()
{
    struct individual *pind;
    int j;
    for(j=0; j<popsize; j++)
    {
        fprintf(outfp,"%3d)  ",j+1);
        /* 當前代個體 */
        pind = &(oldpop[subpop][j]);
        writechrom(pind->chrom);
        fprintf(outfp,"  %8f | ", pind->fitness);
        /* 新一代個體 */
        pind = &(newpop[subpop][j]);
        fprintf(outfp,"(%2d,%2d)   %2d   ",
        pind->parent[0], pind->parent[1], pind->xsite);
        writechrom(pind->chrom);
        fprintf(outfp,"  %8f\n", pind->fitness);
    }
}

void writechrom(chrom)           /* 輸出染色體編碼 */
unsigned *chrom;
{
    int j, k, stop;
    unsigned mask = 1, tmp;
    for(k = 0; k < chromsize; k++)
    {
        tmp = chrom[k];
        if(k == (chromsize-1))
            stop = lchrom - (k*(8*sizeof(unsigned)));
        else
            stop =8*sizeof(unsigned);
        for(j = 0; j < stop; j++)
        {
            if(tmp&mask)
                fprintf(outfp,"1");
            else
                fprintf(outfp,"0");
            tmp = tmp>>1;
        }
    }
}

void preselect()
{
    int j;
    sumfitness = 0;
    for(j = 0; j < popsize; j++) sumfitness += oldpop[subpop][j].fitness;
}

void preselectA()
{
	int j;
	sumA = 0.0;
	for(j = 0; j< subpopnum; j++) sumA += A[j];
}

int select()                    /* 輪盤賭選擇*/
{
    extern float randomperc();
    float sum, pick;
    int i;
    pick = randomperc();
    sum = 0;
    if(sumfitness != 0)
    {
        for(i = 0; (sum < pick) && (i < popsize); i++)
            sum += oldpop[subpop][i].fitness/sumfitness;
    }
    else
        i = rnd(1,popsize);
    return(i-1);
}

int selectA()
{
	extern float randomperc();
    double sum;
	float pick;
    int i;
    pick = randomperc();
    sum = 0.0;
    if(sumA != 0)
    {
        for(i = 0; (sum < pick) && (i < subpopnum); i++)
            sum += A[i]/sumA;
    }
    else
        i = rnd(1,subpopnum);
    return(i-1);
}

void statistics(pop)  /* 計算種群統計數據 */
struct individual *pop;
{
    int i, j;
    sumfitness = 0.0;
    min = pop[0].fitness;
    max = pop[0].fitness;
    /* 計算最大、最小和累計適應度 */
    for(j = 0; j < popsize; j++)
    {
        sumfitness = sumfitness + pop[j].fitness;            
        if(pop[j].fitness > max) max = pop[j].fitness;        
        if(pop[j].fitness < min) min = pop[j].fitness;
        /* new global best-fit individual */
        if(pop[j].fitness > bestfit.fitness)
	  {
	   for(i = 0; i < chromsize; i++)
	    bestfit.chrom[i]      = pop[j].chrom[i];
            bestfit.fitness    = pop[j].fitness;
            bestfit.varible   = pop[j].varible;  
            bestfit.generation = gen;
	  }
      }
    /* 計算平均適應度 */
	if(bestfit.fitness>=maxA)
	{
		maxA=bestfit.fitness;
		for(i = 0; i < chromsize; i++)
			bestfit.chromA[i] = bestfit.chrom[i];
		varibleA = bestfit.varible;
		bestfit.subpopulation = subpop+1;
		runA=run;
		genA=bestfit.generation;
	}
	if(min<=minA) minA = min;
    avg = sumfitness/popsize;
	if(subpop!=(subpopnum-1))
	A[subpop]=avg;
}

void title()
{
  printf("分層遺傳算法\n");
}

void repchar (outfp,ch,repcount)
FILE *outfp;
char *ch;
int repcount;
{
    int j;
    for (j = 1; j <= repcount; j++) fprintf(outfp,"%s", ch);
}

void skip(outfp,skipcount)
FILE *outfp;
int skipcount;
{
    int j;
    for (j = 1; j <= skipcount; j++) fprintf(outfp,"\n");
}

void objfunc(critter)            /* 計算適應度函數值 */
struct individual *critter;
{
    unsigned mask=1;
    unsigned bitpos;
    unsigned tp;
    double pow(), bitpow ;
    int j, k, stop;
    critter->varible = 0.0;
    for(k = 0; k < chromsize; k++)
    {
        if(k == (chromsize-1))
            stop = lchrom-(k*(8*sizeof(unsigned)));
        else
            stop =8*sizeof(unsigned);
        tp = critter->chrom[k];
        for(j = 0; j < stop; j++)
        {
            bitpos = j + (8*sizeof(unsigned))*k;
            if((tp&mask) == 1)
            {
                bitpow = pow(2.0,(double) bitpos);
                critter->varible = critter->varible + bitpow;
            }
            tp = tp>>1;
        }
    }
    critter->varible =-1+critter->varible*3/(pow(2.0,(double)lchrom)-1);
    critter->fitness =critter->varible*sin(critter->varible*10*atan(1)*4)+2.0;
}

void  mutation(unsigned *child)   /*變異操作*/
{
    int j, k, stop;
    unsigned mask, temp = 1;
    for(k = 0; k < chromsize; k++)
    {
        mask = 0;
        if(k == (chromsize-1))
            stop = lchrom - (k*(8*sizeof(unsigned)));
        else
            stop = 8*sizeof(unsigned);
        for(j = 0; j < stop; j++)
        {
            if(flip(pmutation))
            {
                mask = mask|(temp<<j);
                nmutation++;
            }
        }
        child[k] = child[k]^mask;
    }
}


int crossover (unsigned *parent1, unsigned *parent2, unsigned *child1, unsigned *child2)
/* 由兩個父個體交叉產生兩個子個體 */
{
    int j, jcross, k;
    unsigned mask, temp;
    if(flip(pcross))
    {
        jcross = rnd(1 ,(lchrom - 1));/* Cross between 1 and l-1 */
        ncross++;
        for(k = 1; k <= chromsize; k++)
        {
            if(jcross >= (k*(8*sizeof(unsigned))))
            {
                child1[k-1] = parent1[k-1];
                child2[k-1] = parent2[k-1];
            }
            else if((jcross < (k*(8*sizeof(unsigned)))) && (jcross > ((k-1)*(8*sizeof(unsigned)))))
            {
                mask = 1;
                for(j = 1; j <= (jcross-1-((k-1)*(8*sizeof(unsigned)))); j++)
                {
                    temp = 1;
                    mask = mask<<1;
                    mask = mask|temp;
                }
                child1[k-1] = (parent1[k-1]&mask)|(parent2[k-1]&(~mask));
                child2[k-1] = (parent1[k-1]&(~mask))|(parent2[k-1]&mask);
            }
            else
            {
                child1[k-1] = parent2[k-1];
                child2[k-1] = parent1[k-1];
            }
        }
    }
    else
    {
        for(k = 0; k < chromsize; k++)
        {
            child1[k] = parent1[k];
            child2[k] = parent2[k];
        }
        jcross = 0;
    }
    return(jcross);
}

void advance_random()          /* 產生55個隨機數 */
{
    int j1;
    double new_random;
    for(j1 = 0; j1 < 24; j1++)
    {
        new_random = oldrand[j1] - oldrand[j1+31];
        if(new_random < 0.0) new_random = new_random + 1.0;
        oldrand[j1] = new_random;
    }
    for(j1 = 24; j1 < 55; j1++)
    {
        new_random = oldrand [j1] - oldrand [j1-24];
        if(new_random < 0.0) new_random = new_random + 1.0;
        oldrand[j1] = new_random;
    }
}

int flip(float prob)          /* 以一定概率產生0或1 */
{
    float randomperc();
    if(randomperc() <= prob)
        return(1);
    else
        return(0);
}

void randomize()            /* 設定隨機數種子并初始化隨機數發生器 */
{
    int j1;
    for(j1=0; j1<=54; j1++)
      oldrand[j1] = 0.0;
    jrand=0;
    warmup_random(randomseed);
}

double randomnormaldeviate()         /* 產生隨機標準差 */
{
    double sqrt(), log(), sin(), cos();
    float randomperc();
    double t, rndx1;
    if(rndcalcflag)
    {   rndx1 = sqrt(- 2.0*log((double) randomperc()));
        t = 6.2831853072 * (double) randomperc();
        rndx2 = rndx1 * sin(t);
        rndcalcflag = 0;
        return(rndx1 * cos(t));
    }
    else
    {
        rndcalcflag = 1;
        return(rndx2);
    }
}

float randomperc()            /*與庫函數random()作用相同, 產生[0,1]之間一個隨機數 */
{
    jrand++;
    if(jrand >= 55)
    {
        jrand = 1;
        advance_random();
    }
    return((float) oldrand[jrand]);
}

int rnd(low, high)           /*在整數low和high之間產生一個隨機整數*/
int low,high;
{
    int i;
    float randomperc();
    if(low >= high)
        i = low;
    else
    {
        i = (randomperc() * (high - low + 1)) + low;
        if(i > high) i = high;
    }
    return(i);
}


void warmup_random(float random_seed)       /* 初始化隨機數發生器*/
{
    int j1, ii;
    double new_random, prev_random;

    oldrand[54] = random_seed;
    new_random = 0.000000001;
    prev_random = random_seed;
    for(j1 = 1 ; j1 <= 54; j1++)
    {
        ii = (21*j1)%54;
        oldrand[ii] = new_random;
        new_random = prev_random-new_random;
        if(new_random<0.0) new_random = new_random + 1.0;
        prev_random = oldrand[ii];
    }
    advance_random();
    advance_random();
    advance_random();
    jrand = 0;
}


main()           /*  主程序  */
{
    struct individual **temp;
    FILE   *fopen();
    void   title();
    char   *malloc();
        if((outfp = fopen("output4.txt","w")) == NULL)
        {
           fprintf(stderr,"Cannot open output file %s\n","output4.txt");
            exit(-1);
        }
     title();
     initialize();
     maxruns=2;
     /*printf("輸入遺傳算法執行次數(1-5):");
     scanf("%d",&maxruns);*/
   for(run=1; run<=maxruns; run++)
   {
	 
	 //subpopnum=20;
	 if(run==1)
	 {
		//initialize();
		bestfit.subpopulation = 0;
		varibleA = 0;
	 }
	 randomseed = 0;
	 ncross=0;
	 nmutation=0;
	 for(subpop=0;subpop<subpopnum;subpop++)
	 {
        /* 初始化隨機數發生器 */
		randomseed += (float)(1)/(float)(subpopnum+run);
		randomize();
    /* 初始化全局計數變量和一些數值*/
		nmutation = 0;
		ncross = 0;
		bestfit.fitness = 0.0;
		bestfit.generation = 0;
		
    /* 初始化種群，并統計計算結果 */
		if(run==1)
			initpop();
		statistics(oldpop[subpop]);
			for(gen=0; gen<maxgen; gen++)
			{
				generation();
            /* 計算新一代種群的適應度統計數據 */
				statistics(newpop[subpop]);
				temp = oldpop;
				oldpop = newpop;
				newpop = temp;
			}
			report();
			if(subpop==(subpopnum-1))
			{
				ncrossA=0;
				nmutationA=0;
				fprintf(outfp,"\n第 %d / %d 次分層優化\n", run,maxruns);
			/* 選擇交叉變異 */
				CMTransA();
			/* 計算新一代種群的適應度統計數據 */
			
				for(shu = 0; shu<subpopnum; shu++)
					statistics(newpop[shu]);
				reportA();

				temp = oldpop;
				oldpop = newpop;
				newpop = temp;
			}
	 }
   }
	 freeall();
     printf("程序已經運行完畢，請查看結果！");
	 fclose(outfp);
	 getch();
}