#!/bin/perl
#
# Simulate AE research
#
print STDERR "Days until normal airframe availability [1000] ";
( $days = <> ) or die "Error reading airframe delay";
chomp ($days);
if ( $days == 0 ) {
	$days = 1000;
};

print STDERR "Factory maximum size [30] ";
( $maxsize = <> ) or die "Error reading factory max size";
chomp ( $maxsize );
if ( $maxsize == 0 ) {
	$maxsize = 30;
};

print STDERR "How much already repaired [0] ";
( $cursize = <> ) or die "Error reading current factory size";
chomp ( $cursize );
if ( $cursize == 0 ) {
	$cursize = 0;
};

print STDERR "Number of facilities like this one [1] ";
( $facilities = <> ) or die "Error reading number of facilities";
chomp ( $facilities );
if ( $facilities == 0 ) {
	$facilities = 1;
};

# The 2d @sizedist array is indexed first by day number and then by size
# $sizedist[i] gives the distribution of factory sizes on the i'th day (0 = start of first day, $days = end of last day )
# $sizedist[i,j] gives the probability that the factory is size j on the i'th day.
#
# Initialize factory size at start of the first day...
for ( $size=0; $size<=$maxsize; $size++ ) {
	if ( $size == $cursize ) {
		$sizedist[0][$size] = 1;
	} else {
		$sizedist[0][$size] = 0;
	};
};

# Initialize the cumulative research array...
$research[0] = 0;

# Walk distribution forward each day of the period... day 0 is the first day, day $days-1 is the last day.
print STDERR "Number of days = ", $days, "\n";

# Save the originally specified number of days to arrival
$nominal_days = $days;

for ( $day=0; $day<$nominal_days; $day++ ) {
	# $p is the probability of a factory size increase on this day
	if ( $day >= $days ) {
		# We're past modified expected due date.  Assume automatic increase by +1 size
		$p = 1;
	} else {		
		$p = $maxsize / ( $days - $day );
		if ( $p > 1 ) {
			$p = 1;
		};
	};
	# generate size distribution for next day;

	# Size zero factory that does not grow
	$sizedist[$day+1][0] = $sizedist[$day][0] * ( 1-$p );

	# Size n factory that either grew from n-1 or didn't grow from n.
	for ( $size=1; $size<=$maxsize; $size++ ) {
		$sizedist[$day+1][$size] = $sizedist[$day][$size-1] * $p +
					   $sizedist[$day][$size] * (1-$p);
	};

	# If the factory was already max size, it can't grow further.
	$sizedist[$day+1][$maxsize] = $sizedist[$day+1][$maxsize] + $sizedist[$day][$maxsize] * $p;

	# The code gets a little funny here due to an observation by kmitahj
	# http://www.matrixgames.com/forums/tm.asp?m=3022096&mpage=2&key=&#3095091
	#
	# If the suite of factories has generated enough research to advance
	# availability by one or more months then that means that repairs to
	# the not-yet-fully-repaired factories will be accelerated because of
	# the updated availability date.
	#
	# This effect is random and is not modelled exactly.  However, by
	# adding up the expected value of the research generated by all the
	# _other_ factories devoted to this airframe, we get a plausible
	# estimate of the modified due date that affects the repair rate of
	# _this_ factory.
	#
	# Note that when modelling a single factory, this tweak will have
	# no effect whatsoever.
	#
	# I'm not going to swear by the correcness of the (n-1)/n factor on
	# this correction, but it feels like it should unbias the resulting estimator.
	# It's pretty much down in the noise anyway.

	# Update cumulative expected value of research to date

	# First, figure the number of research points per factory per turn
	# complete repairs.
	$p = ($maxsize+1)/30;
	if ( $p > 1 ) {
		$p = 1;
	};

	$today_research = $p * $sizedist[$day][$maxsize] * $facilities;
	$research[$day+1] = $research[$day] + $today_research;
	$other_factory_research = $other_factory_research + $today_research * ($facilities-1) / $facilities;

	if ( $other_factory_research >= 100 ) {
		# Availability date for purposes of factory growth moves up by 30 days
		$days = $days - 30;
		$other_factory_research = $other_factory_research - 100;
	};
		
};


#Debugging -- print the detailed distribution of factory sizes on a day by day basis
#for ( $day=0; $day<=$nominal_days; $day++ ) {
#	print "Size distro for day ", $day, " is:";
#	for ( $size=0; $size<=$maxsize; $size++ ) {
#		print " ", $sizedist[$day][$size];
#	};
#	print "\n";
#};

#Probability of full repairs
$days = $nominal_days;
for ( $day=1; $day<=$days; $day++ ) {
	# Print at 5 percentage break points, but fudge the 100% level down to the 99% threshold
	$percentage = int($sizedist[$day][$maxsize]*100);
	if ( $percentage >= 99 ) {
		$percentage = 100;
	};
	$yesterday = int($sizedist[$day-1][$maxsize]*100);
	if ( $yesterday >= 99 ) {
		$yesterday = 100;
	};
	if ( int($percentage/5) != int($yesterday/5) ) {
		print "Day ", $day, " probability of full repairs = ", $percentage, "%\n";
	};
};

#Production-ready table -- !!do not trust this part, not tested!!
$months_advanced = 0;
for ( $day=0; $day<$days; $day++ ) {
	if ( $day/30 == int($day/30) && $research[$day] >= 1 ) {
		print "Cumulative expected research by month ", int($day/30), " is: ", int($research[$day]*100)/100, "\n";
	};
	if ( $research[$day] - $months_advanced*100 >= 100 ) {
		$days = $days - 30;
		print "Expected advance by one month at day number ", $day, "\n";
		$months_advanced++;
	};
};
print "Expected modified availability at day ", $day, "\n";
print "This represents an acceleration by ", $months_advanced, " months or ", $nominal_days - $day, " days\n";