/*--------------------------------------------------------------------------

		Genetic Algorithm Base program for MaTX
				9-30.1994,  Shinya Kisui

		Genetic Algorithm for Analog Space

				Last Modified	4-26.1995

	This program is Genetic-Algorithm-Optimizeing-Routine for MaTX.
	Using this program, you can solve 最適値(準最適値)探索 problem.

	Format is as follows,
	
	---------------(example)------------------

		Matrix	GA_opt( List );

		Matrix	Answer;
		List	init_param;

		Answer = GA_opt( init_param );
	
	------------------------------------------

	init_param has 8 parameters,

		init_param = { N_I, N_E, N_P, L_P, R_C, R_M, EXC, SED }
	
		Integer	N_I;	: Number of Individuals
		Integer	N_E;	: Number of Elites
		Integer	N_P;	: Number of Pheno
		Integer	L_P;	: Length of Pheno
		Integer	EXC;	: Number of Exchange-Step
		Real	R_C;	: Rate   of Crossing
		Real	R_M;	: Rate   of Mutation
		Integer	SED;	: Seed   of Random

	returned value is ( N_E  x  N_P ) matrix, which means

		Elite-1	: pheno-1(Real)   pheno-2(Real) ... Pheno-(N_P)(Real)
		Elite-2	: pheno-1(Real)   pheno-2(Real) ... Pheno-(N_P)(Real)
		Elite-3	: pheno-1(Real)   pheno-2(Real) ... Pheno-(N_P)(Real)
						.
		Elite-(N_E).................................................


	You have to make fitnessfunction() in fitnessfunction.mm .
	Please look the sample of fitnessfunction.mm .


	If you have any question at this program, please mail to

				kisui@ctrl.titech.ac.jp

------------------------------------------------------------------------*/

/*------------------------------------------------------*/
/*		GA - Main Routine			*/
/*							*/
/*	Initialize Parameters, Excute Optimize,		*/
/*	Display-Print the result............		*/
/*	And Return the Result for Main-Program		*/
/*							*/
/*------------------------------------------------------*/

Func Matrix GA_opt( init_param )
	List	init_param;
{
	/*------------------------------------------------------*/
	/*		Conditions of the Space			*/
	/*------------------------------------------------------*/

	Integer	Num_Ind;		/*	Number of Individuals	*/
	Integer	Num_Eleat;		/*	Number of Eites		*/
	Integer	Exch;			/*	Number of Exchange	*/
	Integer	step;			/*	Generation		*/
	Real	R_Cros;			/*	Rate of Crossing	*/
	Real	R_Mut;			/*	Rate of Mutation	*/
	Integer	G_Len;			/*	Length of Gene		*/
	Integer	G_Num;			/*	Number of Gene		*/

	/*------------------------------------------------------*/
	/*		Parameters of the Individuals		*/
	/*------------------------------------------------------*/

	List	Indiv;			/*	Individual		*/
	List	Group;			/*	Group of Individuals	*/
	List	Group_Eleat;		/*	Group of Elites		*/
	Array	Gene;			/*	Gene			*/
	Matrix	Pheno;			/*	Phenotype		*/
	Real	Fit;			/*	Fitness			*/

	/*------------------------------------------------------*/
	/*		Parameters for DATA-Logging		*/
	/*------------------------------------------------------*/

	Array	max_log;		/*	Log of Max-Fitness	*/
	Array	avg_log;		/*	Log of Avg-Fitness	*/
	Array	exc_log;		/*	Log of Exchange		*/

	/*------------------------------------------------------*/
	/*		Another Parameters			*/
	/*------------------------------------------------------*/

	Integer	i;
	Matrix	Ans;			/*	Matrix for Return Answer */

	/*------------------------------------------------------*/
	/*		Prototype of Functions 			*/
	/*------------------------------------------------------*/

	List	setup_val();		/*	Parameters Setup	*/
	List	setup();		/*	Initial Condition setup	*/
	List	make_Indiv();		/*	Make the Individuals	*/
	List	mutate();               /*      Excute the Mutation    	*/
	List	crossover();            /*      Excute the Crossing    	*/
        List    calc_phenotype();       /*      Calculate the Phenotype	*/
	
        List    calc_fitness();		/*      Calculate the Fitness  	*/
        List    reproduction();         /*      Reproduction		*/
        List    sorting();              /*      Sorting, save the Elites*/
        Array    log_max();		/*      Save the Max-Fitness	*/
        Array    log_avg();		/*      Save the Avg-Fitness	*/
        Array    log_exc();		/*      Save the Exchange-Number*/

	/*------------------------------------------------------*/
	/*		Program	Start				*/
	/*------------------------------------------------------*/

	{ Num_Ind, Num_Eleat, Exch, R_Cros, R_Mut, G_Len, G_Num } = 
			setup_val( init_param );

	{ Group, Group_Eleat, Indiv, max_log, avg_log, exc_log } = 
			setup( Num_Ind, Num_Eleat, G_Len, G_Num, Exch );

	  Group = 	make_Indiv( Group, Num_Ind );

	for( step = 1 ; step <= Exch ; step++ ){
		Group = calc_phenotype( Group, Num_Ind );
		Group = calc_fitness( Group, Num_Ind );
		{ Group, Group_Eleat }=
			sorting( Group, Group_Eleat, Num_Ind, Num_Eleat );

		max_log = log_max( step, Group_Eleat, max_log );
		avg_log = log_avg( step, Group, Num_Ind, avg_log );
		exc_log = log_exc( step, exc_log );

		Group = reproduction( Group, Num_Ind );
		Group = crossover( Group, Num_Ind, R_Cros );
		Group = mutate( Group, Num_Ind, R_Mut );
		print	step;
	}

        print [ [exc_log][max_log][avg_log] ] >> "max_avg.log";

        Ans = Z( Num_Eleat, G_Num );  
        for( i=1 ; i<= Num_Eleat ; i++ ){
                Indiv   = Group_Eleat( i, List );
                Pheno   = Indiv( 2, Matrix );
                Ans( i, : ) = [Pheno];
        }
        return Ans;
}

/*------------------------------------------------------*/
/*		初期値設定ルーチン			*/
/*							*/
/*	○ 生物数・エリート数				*/
/*	○ 世代交代回数					*/
/*	○ 突然変異率・遺伝子交差率			*/
/*	○ 染色体本数・染色体長さ			*/
/*							*/
/*	の各パラメータを設定・変更するルーチン		*/
/*------------------------------------------------------*/

Func List setup_val( param_list )
	List	param_list;
{
        Integer Num_Ind;
        Integer Num_Eleat;
        Integer Exch;
        Real    R_Cros;
        Real    R_Mut;
        Integer G_Len;
        Integer G_Num;

	Integer	select;
	Integer	no_op, i;

	no_op = param_list( 8, Integer );

	for( i = 1 ; i <= no_op ; i++ ){
		rand();
	}

	Num_Ind		= param_list( 1, Integer );
	Num_Eleat	= param_list( 2, Integer );
	G_Num		= param_list( 3, Integer );
	G_Len		= param_list( 4, Integer );
	R_Cros		= param_list( 5, Real    );
	R_Mut		= param_list( 6, Real    );
	Exch		= param_list( 7, Integer );

	return	{ Num_Ind, Num_Eleat, Exch, R_Cros, R_Mut, G_Len, G_Num };
}

/*------------------------------------------------------*/
/*		仮想生物設定ルーチン			*/
/*							*/
/*	○ 仮想生物群・エリート群			*/
/*							*/
/*	の初期設定を行うルーチン			*/
/*------------------------------------------------------*/

Func List setup( Num_Ind, Num_Eleat, G_Len, G_Num, Exch )
	Integer	Num_Ind;
	Integer	Num_Eleat;
	Integer	G_Len;
	Integer	G_Num;
	Integer	Exch;
{
        List    Indiv;
        List    Group;
        List    Group_Eleat;
	Array	max_log;
	Array	avg_log;
	Array	exc_log;

        Matrix  gene_size;
        Matrix  pheno_size;
        Real    Fit;

	Integer	i;

	gene_size 	= Z( G_Num, G_Len );
	pheno_size  	= Z( 1, G_Num );
	Fit		= 0.0;

	Indiv	   	= { gene_size, pheno_size, Fit };
	Group	   	= makelist( Num_Ind );
	Group_Eleat	= makelist( Num_Eleat );

	for( i=1 ; i<=Num_Ind ; i++ ){
		Group( i ) = Indiv;
	}
	for( i=1 ; i<=Num_Eleat ; i++ ){
		Group_Eleat( i ) = Indiv;
	}

	max_log		= Z( 1, Exch );
	avg_log		= Z( 1, Exch );
	exc_log		= Z( 1, Exch );

	return	{ Group, Group_Eleat, Indiv, max_log, avg_log, exc_log };
}

/*------------------------------------------------------*/
/*		生物集団発生ルーチン			*/
/*							*/
/*	設定された 生物数・染色体数・染色体長さ		*/
/*	に基づいて、ランダムに遺伝子を決定し		*/
/*	生物集団を形成する。				*/
/*							*/
/*------------------------------------------------------*/

Func List make_Indiv( Group, Num_Ind )
	List	Group;
	Integer	Num_Ind;
{
	List	Indiv;
	Matrix	Gene;
	Integer	G_Len;
	Integer	G_Num;
	Integer	i, j, k;

	Indiv = Group( 1, List );
	Gene = Indiv( 1, Matrix );
	G_Len = Cols( Gene );
	G_Num = Rows( Gene );

	for( i = 1 ; i <= Num_Ind ; i++ ){
		for( j = 1 ; j <= G_Num ; j++ ){
			for( k = 1 ; k <= G_Len ; k++ ){
				Gene( j, k ) = Integer( rand() );
			}
		}
		Indiv( 1 ) = Gene;
		Group( i ) = Indiv;
	}
	return	Group;
}

/*------------------------------------------------------*/
/*              突然変異ルーチン                        */
/*                                                      */
/*      設定された 突然変異率に基づいて 各生物の        */
/*      遺伝子を反転する。                              */
/*------------------------------------------------------*/
        
Func List mutate( Group, Num_Ind, R_Mut )
	List	Group;
	Integer	Num_Ind;
	Real	R_Mut;
{
	List	Indiv;
	Matrix	Gene;
	Integer	G_Len;
	Integer	G_Num;

        Integer i, j, k;  
        Real    rate; 
        Real    stack;

	Indiv = Group( 1, List );
	Gene  = Indiv( 1, Matrix );
	G_Len = Cols( Gene );
	G_Num = Rows( Gene );

        for( i=1 ; i <= Num_Ind ; i++ ){
                Indiv = Group( i, List );
                Gene  = Indiv( 1, Matrix );
                for( j=1 ; j<=G_Num ; j++ ){
                        for( k=1 ; k<=G_Len ; k++ ){
                                rate = rand();
                                if( rate < R_Mut ){
                                        stack = abs( Gene( j, k ) - 1.0 );
                                        Gene( j, k ) = stack;
                                }
                        }
                }
                Indiv( 1 ) = Gene;
                Group( i ) = Indiv;
        }
	return Group;
}

/*------------------------------------------------------*/
/*              遺伝子交差ルーチン                      */
/*                                                      */
/*      設定された 遺伝子交差率に基づいて ２つの        */
/*      生物の 遺伝子をランダムな位置で交換する。       */
/*------------------------------------------------------*/
                        
Func List crossover( Group, Num_Ind, R_Cros )
	List	Group;
	Integer	Num_Ind;
	Real	R_Cros;
{
	List	Indiv;
	Matrix	Gene;
	Integer	G_Len;
	Integer	G_Num;

        Integer i, j, k, l;
        Real    rate;
        Integer cut_pos;
        Matrix  stack;
        Matrix  gene1, gene2;
        List    indiv1, indiv2;

	Indiv = Group( 1, List );
	Gene  = Indiv( 1, Matrix );
	G_Len = Cols( Gene );
	G_Num = Rows( Gene );

        for( i=1 ; i<=(Num_Ind-1) ; i++ ){
                for( j=(i+1) ; j<Num_Ind ; j++ ){

                        indiv1 = Group( i, List );
                        gene1  = indiv1( 1, Matrix );
                         
                        indiv2 = Group( j, List );
                        gene2  = indiv2( 1, Matrix );
        
                        rate = rand();
                        if( rate < R_Cros ){
                                cut_pos = Integer( rand()*( G_Len -1 ) + 1 );
        
                                stack = gene1( :, cut_pos:G_Len );
                                gene1( :, cut_pos:G_Len ) = gene2( :, cut_pos:G_Len );
                                gene2( :, cut_pos:G_Len ) = stack;
                        }
        
                        indiv1( 1 ) = gene1;
                        Group( i ) = indiv1;
         
                        indiv2( 1 ) = gene2;
                        Group( j ) = indiv1;
                }
        }
	return Group;
}

/*------------------------------------------------------*/
/*              並べ変えルーチン                        */
/*                                                      */
/*      エリートと生物集団を併せたものを適応度に        */
/*      基づいて並べ変え、適応度の高いものをエリート    */
/*      に戻す。                                        */
/*      また、最後の１行を加えると生物集団にエリート    */
/*      を加え、収束を速めることもできる                */
/*------------------------------------------------------*/
        
Func List sorting( Group, Group_Eleat, Num_Ind, Num_Eleat )
	List	Group;
	List	Group_Eleat;
	Integer	Num_Ind;
	Integer	Num_Eleat;
{
        Integer i,j;
        List    indiv1, indiv2;
        Real    fit1, fit2;
        List    All_Group;
        
        All_Group = makelist( Num_Ind + Num_Eleat );
        
        for( i = 1 ; i <= Num_Eleat ; i++ ){
                indiv1 = Group_Eleat( i, List );
                All_Group( i ) = indiv1;
        }
                        
        for( i = 1 ; i <= Num_Ind ; i++ ){
                indiv1 = Group( i, List );
                All_Group( Num_Eleat + i ) = indiv1;
        }

        for( i = 1 ; i <= ( Num_Ind + Num_Eleat - 1 ) ; i++ ){
                for( j = i ; j <= ( Num_Ind + Num_Eleat ) ; j++ ){

                        indiv1 = All_Group( i, List );
                        fit1 = indiv1( 3, Real );

                        indiv2 = All_Group( j, List );
                        fit2 = indiv2( 3, Real );
                        
                        if( fit1 < fit2 ){
                                All_Group( i ) = indiv2;
                                All_Group( j ) = indiv1;
                        }  
                }
        }
        Group_Eleat( 1 : Num_Eleat ) = All_Group( 1 : Num_Eleat );
        
	return	{ Group, Group_Eleat };
}
                
/*----------------------------------------------*/
/*              データログルーチン              */
/*                                              */
/*      現在の世代における適応度の最大値・      */
/*      平均値を記録するルーチン                */
/*----------------------------------------------*/  
         
Func Array log_max( now, Group_Eleat, max_log )
        Integer         now;
	List		Group_Eleat;
	Array		max_log;
{
	List		indiv;
        Integer         i;
        Real            max_now;

        indiv = Group_Eleat( 1, List );
        max_now = indiv( 3, Real );
                         
        max_log( 1, now ) = max_now;

	return	max_log;
}

Func Array log_avg( now, Group, Num_Ind, avg_log )
        Integer         now;
	List		Group;
	Integer		Num_Ind;
	Array		avg_log;
{
	List		indiv;
        Integer         i;
	Real		Fit;
        Real            avg_now;

        avg_now = 0.0;

        for( i=1 ; i<= Num_Ind ; i++ ){
                indiv = Group( 1, List );
                avg_now = avg_now + indiv( 3, Real );
        }
        avg_now = avg_now / ( Num_Ind * 1.0 );

        avg_log( 1, now ) = avg_now;

	return	avg_log;
}
                
Func Array log_exc( now, exc_log )
	Integer	now;
	Array	exc_log;
{
	exc_log( 1, now ) = now;
	return	exc_log;
}

/*----------------------------------------------*/
/*              表現型計算ルーチン              */
/*                                              */
/*      各個体の遺伝子から表現型の計算          */
/*      を行うルーチン                          */
/*----------------------------------------------*/
 
Func List calc_phenotype( Group, Num_Ind )
	List	Group;
	Integer	Num_Ind;
{
	List	Indiv;
	Matrix	Gene;
	Matrix	Pheno;
	Integer	G_Len;
	Integer	G_Num;         

        Integer i, j, k;
        Matrix  gene_vec;
        Real    phenotype;


	Indiv = Group( 1, List );
	Gene  = Indiv( 1, Matrix );
	G_Len = Cols( Gene );
	G_Num = Rows( Gene );

        for( i=1 ; i<=Num_Ind ; i++ ){

                Indiv   = Group( i, List );
                Gene    = Indiv( 1, Matrix );
                Pheno   = Indiv( 2, Matrix );

                for( j=1 ; j<=G_Num ; j++ ){
                        gene_vec = Gene( j, : );
                        phenotype = 0.0;
        
                        for( k=1 ; k<=G_Len ; k++ ){
                                phenotype = phenotype + 2.0^( - k ) * gene_vec( 1, k );
                        } 
                        Pheno( 1, j ) = phenotype;
                }
                Indiv( 2 ) = Pheno;
                Group( i ) = Indiv;
        }
	return	Group;
}
         
/*----------------------------------------------*/
/*              適応度計算ルーチン              */
/*                                              */
/*      各個体の表現型から適応度の計算          */
/*      を行うルーチン。                        */
/*      Real    fitnessfunction() が必要        */
/*----------------------------------------------*/
        
Func List calc_fitness( Group, Num_Ind )
	List	Group;
	Integer	Num_Ind;
{  
	List	Indiv;
	Matrix	Pheno;
	Real	Fit;        

        Integer i,j;
        Real    fitnessfunction();      /*      適応度関数	*/
        
        for( i=1 ; i<=Num_Ind ; i++ ){
                Indiv = Group( i, List );
                Pheno   = Indiv( 2, Matrix );
                Fit   = fitnessfunction( Pheno );
                Indiv( 3 ) = Fit;
                Group( i ) = Indiv;
        }
	return	Group;
}
                        
/*----------------------------------------------*/
/*              淘汰・繁殖ルーチン              */
/*                                              */
/*      各個体の適応度に応じて確率的に          */
/*      次世代の個体を選択するルーチン          */
/*----------------------------------------------*/

Func List reproduction( Group, Num_Ind )
	List	Group;
	Integer	Num_Ind;
{
	List	Indiv;
        List    New_Group;
        Matrix  P;
        Integer i, j;
        Real    fitness;
        Real    sum_fitness;
        Real    rate;
        Integer check; 
        
	List	make_Indiv();

        P = Z( Num_Ind, 1 );
        sum_fitness = 0.0;
        
        New_Group = makelist( Num_Ind ); 
                
        for( i = 1 ; i <= Num_Ind ; i++ ){   
                Indiv = Group( i, List );
		fitness = Indiv( 3, Real );
		if( fitness < 0.0 ){
			error( "\nFitness is less than 0.\nCan't reproduct New generation !!\n" );
			exit(-1);
		}else{
	                sum_fitness = sum_fitness + fitness;
		}
        }
        
	if( sum_fitness <= 0.0 ){
		warning( "\nSum of fitness = 0.0, Create Individuals again \n");
		Group = make_Indiv( Group, Num_Ind );
	}else{
	        for( i = 1 ; i <= Num_Ind ; i++ ){
	                Indiv = Group( i, List );
	                fitness = Indiv( 3, Real );
	                P( i, 1 ) = fitness / sum_fitness;
	        }
	        for( i = 1 ; i <= Num_Ind ; i++ ){
	                rate = rand();
	                check = 0;
	                j = 1;
	                while( check == 0 ){
	                        if( rate <= P( j, 1 ) ){
	                                New_Group( i ) = Group( j, List );
	                                check = 1;
	                        }else{
	                                rate = rate - P( j, 1 );
	                        }
	                        j++;
	                }
	        }
	}
                 
	return	New_Group;
}