// stats.ccc  version to support queue2.cc, from GPQUICK-2.1 chrome.cc
// W.Langdon cs.ucl.ac.uk $Revision: 1.6 $

//Modifications (reverse order):
// WBL 13 Oct 1994  Display new meaning of adf1

// WBL 13 Sep 1994  Add displaying adf1

// WBL 13 Sep 1994  Make compilable by solaris gcc

// WBL 25 Aug 1994  Remove BestFitness check

// WBL 22 Aug 1994  Add queue2.h
//                  Take hits value from nfitness, rather than Chrome
//                  Add display of mem_used.

// WBL 15 Aug 1994  Add display of hits and total_hits

// WBL	6 Jul 1994  Fix bug whereby may not set last few values of nfit

// WBL  1 Jul 1994  New file based upon chrome.cxx Revision: 1.9
//                  Add pTraceDynamic.
//		    Display primative frequency, covarience and correllation

// WBL 28 Jun 1994  Make nfit take average of cube rather than linear average
//                  Tidy displays in DisplayDStats

// WBL 27 Jun 1994  Correct display of gencount in DisplayStats
//                  Add class stats and make DisplayStats use it.
//                  Complete impementation of DisplayDStats

// WBL 24 Jun 1994  Direct crossover to a particular tree
//                  Protect sqrt from rounding errors

// WBL 14 Jun 1994  Add support for pPopSeed, pTestSeed and PTrace

// WBL 14 Jun 1994  Add ChromeParams::Save() and ChromeParams::Load()
//                  based upon gpcplus3.0 cpv.cc.
//                  Add Pop::DisplayDStats()
//		    Add Pop::Write()

// WBL 16 May 1994  Add Pop::DisplayStats(). Fix MULTREE CrossOver bug

// WBL 12 May 1994  fix % by zero error in rnd() on SETNODE if varnum==0,
//		    Probably only arises if you have the wrong compiler switch

// WBL  5 May 1994  Make pMaxAge = 0 disable age limit,
//		    Add support for P-M random numbers

// GPQUICK
// C++ prefix stack implementation
// Divides GP system into classes:

#include <math.h>	//for DisplayStats
#include <assert.h>
#include "pch.h"
#include "chrome.h"
#include "queue2.h"
#include "suit.h"

class stats {						//WBL
	int cases = 0;
	float min = FLT_MAX;
	float max = -FLT_MAX;
	float sum = 0.0;
	float sqr = 0.0;

float	Mean() { return sum / cases; };
float	Variance();
public:
float GetMean() {return sum / cases;};
void 	include ( float in );
int	Cases() { return cases; };
void	display ( char* name, ostream & fout,
		  stats*,     stats*         );
};//end class stats

float stats::Variance()					//WBL
{ 	float v = (sqr - cases*Mean()*Mean() )/cases;
	if ( v >= 0 ) return v; else return 0.0; //protect against rounding
};

void stats::include( float in )  			//WBL
{	cases++;
	sum += in;
    	sqr += in * in;
	if ( max < in ) max = in;
	if ( min > in ) min = in;
}//end stats::include

#ifdef INTERFACE
void stats::display (char * objectname, ostream & fout = cout,stats* own_stats = NULL, stats* co_stats = NULL)
{
if (cases > 0)
//cout << "got here";
	{
	float variance = Variance();
	float SD = sqrt(variance);
	char varc[100];
	char SDc[100];
	sprintf(varc,"%f",variance);
	sprintf(SDc,"%f",SD);
	SUIT_setText(SUIT_getChild(SUIT_name(objectname),0),LABEL,varc);
	SUIT_setText(SUIT_getChild(SUIT_name(objectname),1),LABEL,SDc);
	}
}
#else

void stats::display ( char* name, ostream & fout = cout,
                    stats* own_stats = NULL, stats* co_stats = NULL )
                                                      //WBL
{
if ( cases > 0 )
    {	float mean     = Mean();
	float variance = Variance();
	fout << "Mean " << name << "= " << mean;
	fout << "\tvar = " << variance;
	fout << "\tSD = "<<sqrt(variance);
	if ( own_stats != NULL && co_stats != NULL)
	     { 	fout << "\tcovar ";
//
//covariance = 	sum (x(i)*y(i))
//              ---------------  -  mean x * mean y
//               num cases
//
//In our case sum (x(i)*y(i)) = sum because we count y times. 
//The mean x is number of x (ie cases) divided by the popsize
//
		float covar;

		covar = (float)sum - (float) cases * co_stats->Mean();
		covar = covar/(float) own_stats->Cases();
		fout << covar;

// Correlleation co-efficient =             covariance
//				------------------------------------
//				sqrt( own variance * others variance
//
// Own variance includes cases where chrome contains none of the
// opcode, nb including those not in cases
//
// own variance = sqr / popsize - ( sum / popsize )**2
//
		{
		float d = sqrt( own_stats->Variance() * co_stats->Variance() );
		if ( d > 0 )			//protect rounding error
		      { fout << "\tcorel ";
			fout << covar / d;
		      }
		} 
	     }
	fout << "\tmin = " << min;
	fout << "\tmax = " << max <<endl;
    }//endif
}//end stats::display

#endif
void Pop::DisplayStats(ostream & fout) //WBL
{
	time_t now;
//	ChromeParams* params;
//	Problem* problem;
	time (&now);
#ifndef INTERFACE

	fout << ctime(&now);
	fout << "popsize =\t"<< popsize  <<endl; 
	fout << "individuals =\t" <<gencount+1  <<endl;
	fout << "BestFitness =\t" <<*BestFitness  <<endl;
	fout << "BestMember =\t" <<BestMember  <<endl;
	fout << "initeval =\t" <<initeval  <<endl;
	fout << "nexteval =\t" <<nexteval  <<endl;
#endif
	stats fitness;
	stats length;
	stats hits [empty+1];
	stats total_hits;
	stats memory;
	stats adf1;

	for (UINT i=0;i<popsize;i++)
	      { fitness.include (pop[i]->nfitness->fvalue);
		length.include (pop[i]->tree_starts[NUM_TREES-1]
                          + pop[i]->SubLen(pop[i]->tree_starts[NUM_TREES-1]));
		int total = 0;
                                                             //force right type
		ptrmyfitnessvalue temp = ptrmyfitnessvalue (pop[i]->nfitness);
		for (int j=0; j<=empty; j++)
			{ total += temp->hits[j];
			  hits[j].include (temp->hits[j]);
		        }
		total_hits.include (total);
		memory.include (temp->mem_used);
		adf1.include(temp->adf1);
	      }
#ifdef INTERFACE
	hits[0].display("treestack1",fout);
	hits[1].display("treestack2",fout);
	hits[2].display("treestack3",fout);
	hits[3].display("treestack4",fout);
	hits[4].display("treestack5",fout);
	hits[5].display("treestack6",fout);
#else                                                      
	fitness.display( "fitness", fout );
	length.display( "length ", fout );

	for (int j=0; j<=empty; j++) 
		{ problem->WriteTreeName(j, fout);

		  fout << "\t";
		  hits[j].display("", fout);
	        }
#endif
#ifdef INTERFACE
	SUIT_setText(SUIT_name("tl1"),LABEL,problem->GetTreeName(0));
	SUIT_setText(SUIT_name("tl2"),LABEL,problem->GetTreeName(1));
	SUIT_setText(SUIT_name("tl3"),LABEL,problem->GetTreeName(2));
	SUIT_setText(SUIT_name("tl4"),LABEL,problem->GetTreeName(3));
	SUIT_setText(SUIT_name("tl5"),LABEL,problem->GetTreeName(4));
	SUIT_setText(SUIT_name("tl6"),LABEL,problem->GetTreeName(5));
#endif

#ifndef INTERFACE
	fout << "total\t";
	total_hits.display( "", fout );
	fout << "memory\t";
	memory.display( "", fout );
	fout << "adf1ver\t";
	adf1.display( "", fout );
#endif

}//end Pop::DisplayStats

#ifdef INTERFACE
double Pop::UpdateInterface(float max, double increment,double x)
{
	float best_fitness;
	float mean_fitness;
	stats fitness;
	for (UINT i =0;i<popsize;i++)
	{
		
		fitness.include (pop[i]->nfitness->fvalue);
	}
	mean_fitness = fitness.GetMean();
	best_fitness = *BestFitness;
	double normalized_mean = mean_fitness/max;
	double normalized_best = best_fitness/max;
	SUIT_viewport graph1;
        SUIT_viewport graph2; /* The active Viewport */
        rectangle       recvView1;
        rectangle       recvView2;
        GP_pushGraphicsState();

                        graph1 = SUIT_getViewport(SUIT_name("Best Fitness"), VIEWPORT);
                        recvView1.bottom_left.x = graph1.bottom_left.x + 25;
                        recvView1.bottom_left.y = graph1.bottom_left.y + 25;
                        recvView1.top_right.x   = graph1.top_right.x  - 20;
                        recvView1.top_right.y   = graph1.top_right.y - 20;
//                        GP_setViewport(recvView1);

                        graph2= SUIT_getViewport(SUIT_name("Average Fitness"), VIEWPORT);
                        recvView2.bottom_left.x = graph2.bottom_left.x + 25 ;
                        recvView2.bottom_left.y = graph2.bottom_left.y + 25;
                        recvView2.top_right.x   = graph2.top_right.x - 20 ;
                        recvView2.top_right.y   = graph2.top_right.y - 20;
//                        GP_setViewport(recvView);
			GP_setViewport(recvView1);
			GP_setColor(GP_defColor("black",TRUE));
			GP_lineCoord(x,previous_mean,(x+increment),normalized_mean);
			previous_mean = normalized_mean;
			GP_setViewport(recvView2);
			GP_setColor(GP_defColor("black",TRUE));
			GP_lineCoord(x,previous_best,(x+increment),normalized_best);
			previous_best = normalized_best;
			
	GP_popGraphicsState();
	x = x + increment;
//	cout << x << "Thst was x" << endl;
	return x;
}
#endif
		
void Pop::DisplayDStats(ostream & fout ) //WBL
//
// Convert chromosome nfitness into a rank, and then into an averag
// rank and then into an approximate chance of selection. Cf roulette
// wheel selection.
//    Gather and display loads of statistics by tree/opcode
//
{
Chrome* rank [popsize];
float nrank [popsize];   //normalised rank
float nfit [popsize];    //approx equivelent roulette wheel
				   
for (UINT i = 0; i < popsize; i++ ) rank [i] = pop[i];
UINT j;
for ( i = 0; i < popsize; i++ )
      { BOOL swapflag = FALSE;
	for (j = 1; j < popsize; j++ )
	      { if ( rank[j-1]->nfitness->IsBetter(rank[j]->nfitness) )
		      { Chrome* r  = rank [j-1];
			rank [j-1] = rank [j];
			rank [j]   = r;
			swapflag   = TRUE;
                      }
	      }//end loop j
	if (swapflag == FALSE) break;
      }//end loop i

//Nolonger true WBL 25/8/94
//assert( BestFitness->fvalue == rank[popsize-1]->nfitness->fvalue ); //check

assert (params->params[pTournSize] == 4);
UINT k = 0;
float P = popsize;
float sum_normalised_fitness = 1.0/(P*P*P); //large population approximation
					    //((rank/p)**3)
UINT l;
for ( j = 1; j < popsize; j++ )
      { float J = j;
	float K = k;
	if ( rank[j]->nfitness->IsBetter(rank[k]->nfitness))
	      { 
		nrank [ k ] = (J*(J-1)-(K-1)*K)/(2*(J-K));
		nfit  [ k ] = sum_normalised_fitness/(J-K);
		for ( l = k+1; l<j; l++ )
		      { nrank[ l ] = nrank [ k ];
			nfit [ l ] = nfit  [ k ];
		      }
		k = j;
		sum_normalised_fitness = 0.0;
	      }//end j is not same rank as k
	sum_normalised_fitness +=((J + 1.0)/P) * ((J + 1.0)/P) * ((J + 1.0)/P);
      }

for ( l = k; l<popsize; l++ )
	      { float J = popsize; //NB include j
		float K = k;
		nrank [ l ] = (J*(J-1)-(K-1)*K)/(2*(J-K));
		nfit  [ l ] = sum_normalised_fitness/(J-K);
	      }


stats length [NUM_TREES];
struct x
{ stats frequency; stats length; stats children;
  stats fitness;   stats nrank;  stats nfit;     }; 
x summary;
x counts [NUM_TREES][problem->funccount];

for ( i = 0; i < popsize; i++ )
{
//summary.frequency JUNK
//summary.length    JUNK
summary.children.	include( rank[i]->num_children );
summary.fitness.	include( rank[i]->nfitness->fvalue );
summary.nrank.		include( nrank [i] );
summary.nfit.		include( nfit [i] );

int len;
for ( int t = 0; t<NUM_TREES; t++)
      { len = rank[i]->SubLen(rank[i]->tree_starts[t]);
	length[t].  include ( len );

	int f [problem->funccount];
	memset(f,0,sizeof(f));

	{for ( int n = rank[i]->tree_starts[t];
              n < (rank[i]->tree_starts[t] + len);
	      n++ 				)
	{
		f [rank[i]->expr[n].op]++;

counts[t][rank[i]->expr[n].op].length.  include ( len );
counts[t][rank[i]->expr[n].op].children.include ( rank[i]->num_children );
counts[t][rank[i]->expr[n].op].fitness. include ( rank[i]->nfitness->fvalue );
counts[t][rank[i]->expr[n].op].nrank.   include ( nrank[i] );
counts[t][rank[i]->expr[n].op].nfit.    include ( nfit[i] );

	}};//end for each op in this tree

	{for ( int op = 0; op < problem->funccount; op++ )
	      { counts[t][op].frequency.include ( f[op] ); }}

      };//end for each tree
};//end for whole population

summary.children.display ("tree child"  , fout );
summary.fitness. display ("tree fitness", fout );
summary.nrank.   display ("tree nrank"  , fout );
summary.nfit.    display ("tree nfit"   , fout );

for ( int t = 0; t<NUM_TREES; t++)
{
fout << endl;
problem->WriteTreeName(t, fout);
fout << "\t";
length[t].display ("tree size", fout );
for ( int op = 0; op < problem->funccount; op++)
{
	problem->WriteTreeName(t, fout);
	fout << "\t" << problem->funclist[op]->getname() << "\t";
	fout << "occurs " << counts[t][op].length.Cases() << " times\t";
	counts[t][op].frequency.  display ("frequency ", fout );


	if ( counts[t][op].length.Cases() >0 )
	{
	problem->WriteTreeName(t, fout);
	fout << "\t" << problem->funclist[op]->getname() << "\t";

	counts[t][op].length.  display ("length ", fout, 
		&counts[t][op].frequency, &length[t] );
	
	problem->WriteTreeName(t, fout);
	fout << "\t" << problem->funclist[op]->getname() << "\t";

	counts[t][op].children.display ("childs ", fout,  
		&counts[t][op].frequency, &summary.children );

	problem->WriteTreeName(t, fout);
	fout << "\t" << problem->funclist[op]->getname() << "\t";

	counts[t][op].fitness. display ("fitness", fout,  
		&counts[t][op].frequency, &summary.fitness );

	problem->WriteTreeName(t, fout);
	fout << "\t" << problem->funclist[op]->getname() << "\t";

	counts[t][op].nrank.   display ("nrank  ", fout,  
		&counts[t][op].frequency, &summary.nrank );

	problem->WriteTreeName(t, fout);
	fout << "\t" << problem->funclist[op]->getname() << "\t";

	counts[t][op].nfit.    display ("nfit   ", fout,
		&counts[t][op].frequency, &summary.nfit );
	}//end if op code is used in tree
};//end for each op code
};//end for each tree

}//end Pop::DisplayDStats