// M.Sc. model
// converted into C++

// mscm0001.cpp
// variation of mscmodel that outputs a range of altruisms for given
// lambda, alpha and beta.

#include <iostream.h>
#include <math.h>
#include <fstream.h>

void initialise_sub_populations();
void increment_reality_function();
void input_parameters();
void loop();
void calculate_benefits();
void find_total_benefit();
void calculate_new_populations();
void find_relative_populations();
void display_populations();
void display_to_file();
void metaloop();


long double benefit[10], alpha, beta;
long double lambda, reality = 0.5;
long double number[10], sensitivity[10];
long double relative_number[10], total_number;
long double total_benefit, altruism;

int i, c, trials;

ofstream output_file("mscout.txt");

void main()
{
input_parameters();

output_file << "Lambda = " << lambda;
output_file << "   Alpha = " << alpha;
output_file << "	 Beta = " << beta;
output_file << "Trials: " << trials << endl;

metaloop();
output_file.close();
}

void metaloop()
{
for (int m=0; m<21; m++)
{
	altruism = (double)m / 20.;
	initialise_sub_populations();
	loop();
	display_to_file();
   }
}


void display_to_file()
{

//output_file << "Lambda = " << lambda;
//output_file << "   Alpha = " << alpha;
//output_file << "	 Beta = " << beta;
output_file << "   Altruism = " << altruism;

//output_file << "Trials: " << trials;
output_file << "  Total population: " << total_number << endl;

//for (i=0; i<10; i++)
//{
//	output_file << "Sensitivity: " << sensitivity[i];
//	output_file << "  Number: " << number[i];
//	output_file << "  Relative: " << relative_number[i] << endl;
//
//	}

}

void input_parameters()
{
cout << "Enter value of lambda (0 - 4): ";
cin >> lambda;

cout << "Enter number of trials: ";
cin >> trials;

//cout << "Enter degree of altruism (1 = altruistic, 0 = selfish): ";
//cin >> altruism;

cout << "Enter selection parameters: " << endl;
cout << "alpha: ";
cin >> alpha;
cout << " beta: ";
cin >> beta;
}


void initialise_sub_populations()
{
for (i=0; i<10; i++)
	{
	number[i]=1.;
	sensitivity[i]=0.1* (long double)i + 0.05;
	cout << "Number " << number[i] << " Sensitivity: " << sensitivity[i] << endl;
	};
}


void loop()
{
for (c=0; c<trials; c++)
{
	calculate_benefits();
	find_total_benefit();
	calculate_new_populations();
	find_relative_populations();
   cout << c << endl;
	//display_populations();
	increment_reality_function();
	}
}


void increment_reality_function()
{
	reality = lambda * reality * (1. - reality);
}


void calculate_benefits()
{
for (i=0; i<10; i++)
{
	long double square = reality - sensitivity[i];
   square = -1. * beta * square * square;
	benefit[i]=alpha * expl(square);
	}
}

void find_total_benefit()
{
total_benefit = 0.;
total_number = 0.;

for (i=0; i<10; i++)
{
	total_benefit += benefit[i]*number[i];
	total_number += number[i];
	}

total_benefit = altruism * total_benefit/total_number;
}


void calculate_new_populations()
{
total_number = 0.;

for (i=0; i<10; i++)
{
	number[i] = number[i] * (total_benefit + (benefit[i] * (1. - altruism)));
   total_number += number[i];
	}
}


void display_populations()
{
cout << "Trial: " << c;
cout << "  'Reality' function: " << reality;
cout << "  Total population: " << total_number << endl;

for (i=0; i<10; i++)
{
	cout << "Sensitivity: " << sensitivity[i];
	cout << "  Number: " << number[i];
	cout << "  Relative: " << relative_number[i] << endl;
	}

cout << endl;
}


void find_relative_populations()
{
for (i=0; i<10; i++)
	relative_number[i] = number[i]/total_number;
}