###################################
# A wrap to run the founder       #
# simulation for a number of      #
# trials and calculate the        #
# average output per generation   #
# (effectively sampling as        #
# many genes as trials            #
# Robin Allaby UoW 2018           #
###################################

$N = $ARGV[0]; #population size
$Nb = $ARGV[1]; # founders
$a = $ARGV[2]; # no. alleles
$f = $ARGV[3]; #no. founder events
$m = $ARGV[4]; #model of gene frequency generation
$trials = $ARGV[5];
$ms = $ARGV[6]; # mating strategy % inbreeding
$deltams = $ARGV[7]; # rate of change of mating strategy
$mu = $ARGV[8];
$fast = $ARGV[9];
$starthet = $ARGV[10];
$startms = $ARGV[11];





$x = 0;
while ($x < $trials) {
$trial_number = $x+1;
print "trial $trial_number\n";
@argv = ();
@argv = ("perl founderv4.pl $N $Nb $a $f $m $ms $deltams $mu $fast $starthet $startms");
system(@argv) == 0 or die "Cannot run founder.pl\n";
open (IN, "founder.out") or die "Cannot open output file\n";
@frequencies = (); @frequencies = <IN>;
push (@allfreqs, [@frequencies]);
close IN;
$x++;
}

# now we have a matrix of all simulations
print "calculating means...\n";
$x = 0; $y = 0;

@avfreq = ();
@sd = ();

while ($y < $f) {
@collect = ();
	$x = 0;
	while ($x < $trials) {
		push (@collect, $allfreqs[$x][$y]); $x++;
	}
	$collect = join ("\t", @collect);
	$mean = find_mean ($collect);
	$sd = find_standard_deviation ($collect);
	push (@avfreq, $mean); push (@sd, $sd);
	$y++;
}

open (OUT, ">mean_founder.out");
$k = 0;
foreach $mean (@avfreq) {
	print OUT "$mean\t$sd[$k]\n"; $k++;
}
close OUT;


#####################
# SUBROUTINES       #
#####################

sub find_mean {
                my $list = $_[0]; my @list = ();
                my $x = 0; my $l = (); my $t = (); my $mean = ();
                @list = split (/\t/, $list);
                $l = scalar @list;
                if ($l == 0) { return $l;}
                foreach $t (@list) { $x = $x+$t;
		#print TEST "t: $t\n";
		}
		if ($x == 0) {return $x;}
        $mean = $x/$l;
        return $mean;

}

sub find_variance {
        my $list = $_[0]; my @list = (); my $xsquares = ();
        my $mean = (); my $x = (); my $xsquare = (); my @xsquares = ();
        my $meanofsquares = (); my $var = ();

        $mean = find_mean ($list);
        @list = split (/\t/, $list);
        foreach $x (@list) {
                $xsquare = $x*$x;
                push (@xsquares, $xsquare)
        }
        $xsquares = join ("\t", @xsquares);
        $meanofsquares = find_mean ($xsquares);
        $var = $meanofsquares - ($mean*$mean);
        return $var;
}

sub find_standard_deviation {
                my $list = $_[0]; my $var = (); my $sd = ();

                $var = find_variance ($list);
                # strange exceptions can occur which return a negative variance
                # probably due to perl's handling of floating point values
                # very close to zero (I guess), so here is a dodge
                if ($var < 0) { return 0;}
                $sd = sqrt ($var);
        return $sd;
}