//my stupid expectation formula requires iterating over all sequences in {0,...,n-1}^k.
//the following two functions let us do that by iterating instead over {0,...,n^k -1},
//which is significantly easier.


//this takes a number a between 0 and base^power - 1 (inclusive) and converts it into
//the corresponding natural sequence embedded in {0,...,base}^power.
int[int] to_sequence(int a, int base, int power){
	int[int] seq;
	int leftovers = a;
	
	//doesn't matter for our purposes whether we preserve orientation, so we might not.
	for i from 1 to (power){
		seq[power-i] = leftovers/(base**(power - i));
		//debugging.
		//print((base**(power - i))+" goes into "+leftovers+" this many times: "+seq[power-i]);
		leftovers -= seq[power-i]* (base**(power - i));
	}
	return seq;
}

//turns those sequences into sets for cleaner use in foreach loops.
boolean[int] to_set(int[int] seq){
		boolean[int] set;
		foreach i in seq{
			set[seq[i]] = true;
		}
		return set;
}

//for any appropriate number X, "foreach i in X.to_sequence(base,power).to_set()" iterates
//i over precisely the unique digits found in the sequence corresponding to X.
//if you don't discard the multiplicity here, the iteration becomes quite a bit messier.
//Coding natural morphisms is tricky, huh?


//the formula includes the multiplicity function, so we write that here.
//this'll return the number of times a digit j appears in the sequence seq.
int mult(int[int] seq, int j){
	int count;
	foreach i in seq{
		if(seq[i] == j){
			count++;
		}
	}
	return count;
}

//phi is the underlying probability distribution. I guess we'll normalize it in this function
//for ease of use (so any proportional distribution could be input).
//r will be the rejection rate array with positions corresponding to entries in phi.
//k will be queue length, since this can't be deduced from r and phi.
float expectation(float[int] phi, float[int] r, int k){
	//deduce size of outcome set from input rather than requiring that input.
	int n = count(phi);
	//normalize phi.
	float phisum = 0.0;
	foreach i in phi{
		phisum += phi[i];
	}
	foreach i in phi{
		phi[i] /= phisum;
	}
	//so now phi is a probability distribution unless some sick fart put negative/complex values in.
	
	//so now we have to do this. https://i.imgur.com/DsYTyrd.png
	//the top and bottom are nearly identical, so we'll be computing these simultaneously
	//while leaving out that factor of #a from the summands on the bottom.
	float top_outer_sum = 0.0;
	float bottom_outer_sum = 0.0;
	float inner_sum;
	float indexed_product;
	int[int] index_seq;
	boolean[int] index_set;
	for index from 0 to ((n**k) -1){
		clear(index_seq);
		clear(index_set);
		inner_sum = 0.0;
		indexed_product = 1.0;
		index_seq = to_sequence(index,n,k);
		index_set = to_set(index_seq);
		
		//make the inner sum.
		foreach j in index_set{
			inner_sum += phi[j]*r[j];
		}
		//do the 1- thing.
		inner_sum = 1 - inner_sum;
		
		//now make the indexed product.
		
		foreach j in index_set{
			indexed_product *= ( ( phi[j]*(1-r[j]) ) ** mult(index_seq,j) ) / ( 1 - r[j] );
		}
		
		//now we increment the bottom by those two things multiplied together.
		bottom_outer_sum += inner_sum * indexed_product ;
		
		//and increment the top by the same thing but with a factor of #a.
		top_outer_sum += inner_sum * indexed_product * count(index_set);
		
	}
	
	//smush the top and bottom together and return it.
	return (top_outer_sum / bottom_outer_sum );
}

/*

//check that this works by looking at dice.
float[int] fee = {1,1,1,1,1,1};
float[int] arrr = {0,0,0,0,0,0};
float[int] herman_cain = {0.999,0.999,0.999,0.999,0.999,0.999};
print(expectation(fee , arrr ,1));
//the above should print the average number of unique die roll outcomes occurring
//in the last 1 die rolls. If it's not 1.0, we're very sad.
//dear lord it returns exactly 1.0.

debugging stuff
print("let\'s write "+122+" in base 5 with 4 digits");
int[int] boop = 122.to_sequence(5,4);
foreach i in boop{
print(i+": "+boop[i]);
}

print("let\'s write "+69+" in base 3 with 11 digits");
boop = 69.to_sequence(3,11);
foreach i in boop{
print(i+": "+boop[i]);
}
*/


//let's make data for some very pretty graphs.
float[float] data;
float[int] fee = {1,1,1,1,1};
float[int] arrr = {0,0,0,0,0};

//graphing how much rejection rate matters to steady state frequency.
for i from 1 to 8{
	arrr[4] = 0.0;
	for j from 0 to 999{
		arrr[4] += 0.001;
		data[arrr[4]] = expectation(fee,arrr,i)/(i.to_float());
	}
	//when i is 1, this should generate a text file with a column of 1000 1.0's.
	//only letting script do that for debugging.
	map_to_file(data,"plot of normalized expected variety in uniform phi with rejection on one outcome for n=5 k="+i+".txt");
	clear(data);
	print("done with making graph number "+i);
}