#!/usr/statlocal/bin/perl -wT

use CGI;
use CGI::Carp 'fatalsToBrowser';

use strict;
use warnings;
use POSIX;

use constant TOOBIG => 1.0e6;

sub htmldie(@)
{
    my $query = shift;
    print "<p><h1>CGI Error</h1><p>" . join(' ', @_ ) . "<p>\n";
    exit;
}

sub print_head($;$)
{
    my $query = shift;
    my $title = (@_ > 0) ? shift : "";

    $query->print( $query->header, "\n" );
    $query->print( $query->start_html( $title ) );
    $query->print( $query->h1( $title ), "\n" );
}

sub rnd_pc($)
{
    my $x = shift;

    return( sprintf( "%5.2f%%", 100*$x) );
}

sub print_form($)
{
    my $query = shift;

    $query->print( $query->p( <<END_OF_INTRO ) );
Welcome to the Low-Tech Hardy-Weinberg Simulator. Here you can
explore the models we discussed in class without shuffling any cards.

END_OF_INTRO

    $query->print( $query->startform, "\n" );

    $query->print( $query->p( "Set your initial population for each genotype. (Total will be made even. No more than a million, please.)\n", $query->br,
                   "AA ", $query->textfield( -name=>'nAA', -default=>"100",
                             -size=>12, -maxlength=>128 ), "\n", $query->br,
                   "Aa&nbsp;", $query->textfield( -name=>'nAa', -default=>"200",
                             -size=>12, -maxlength=>128 ), "\n", $query->br,
                   "aa&nbsp;&nbsp;", $query->textfield( -name=>'naa', -default=>"100",
                             -size=>12, -maxlength=>128 ), "\n", $query->br ) );

    $query->print( $query->p( "How many generations to simulate? &nbsp;&nbsp; ",
                   $query->textfield( -name=>'m', -default=>"10",
                             -size=>12, -maxlength=>128 ), "\n", $query->br ) );

    $query->print( $query->p( "Survival probabilities for offspring of each genotype.\n", $query->br,
                   "AA ", $query->textfield( -name=>'sAA', -default=>"1",
                             -size=>12, -maxlength=>128 ), "\n", $query->br,
                   "Aa&nbsp;", $query->textfield( -name=>'sAa', -default=>"1",
                             -size=>12, -maxlength=>128 ), "\n", $query->br,
                   "aa&nbsp;&nbsp;", $query->textfield( -name=>'saa', -default=>"1",
                             -size=>12, -maxlength=>128 ), "\n", $query->br ) );

    $query->print( $query->p( "Emigration probabilities for offspring of each genotype.\n", $query->br,
                   "AA ", $query->textfield( -name=>'eAA', -default=>"0",
                             -size=>12, -maxlength=>128 ), "\n", $query->br,
                   "Aa&nbsp;", $query->textfield( -name=>'eAa', -default=>"0",
                             -size=>12, -maxlength=>128 ), "\n", $query->br,
                   "aa&nbsp;&nbsp;", $query->textfield( -name=>'eaa', -default=>"0",
                             -size=>12, -maxlength=>128 ), "\n", $query->br ) );

    $query->print( $query->p( "Average number of immigrants per generation, by genotype.\n", $query->br,
                   "AA ", $query->textfield( -name=>'iAA', -default=>"0",
                             -size=>12, -maxlength=>128 ), "\n", $query->br,
                   "Aa&nbsp;", $query->textfield( -name=>'iAa', -default=>"0",
                             -size=>12, -maxlength=>128 ), "\n", $query->br,
                   "aa&nbsp;&nbsp;", $query->textfield( -name=>'iaa', -default=>"0",
                             -size=>12, -maxlength=>128 ), "\n", $query->br ) );

    $query->print( $query->p( "Mutation probabilities for each allele.\n", $query->br,
                   "A ", $query->textfield( -name=>'mA', -default=>"0",
                             -size=>12, -maxlength=>128 ), "\n", $query->br,
                   "a&nbsp;", $query->textfield( -name=>'ma', -default=>"0",
                             -size=>12, -maxlength=>128 ), "\n", $query->br ) );

    $query->print( $query->hidden( 'mode', 'yes' ), "\n" );

    $query->print( $query->submit( 'Let the mating begin' ), $query->defaults( 'Reset Fields' ), $query->p, "\n" );
 
    $query->print( $query->endform, "\n", $query->hr, $query->p, "\n" );
}

sub mate($$)
{
    my $male = shift;    # both in number of a alleles
    my $female = shift;

    $male /= 2;
    $female /= 2;

    my ($x, $y) = (rand, rand);

    return( 0 + (($x > $male) ? 0 : 1 ) + (($y > $female) ? 0 : 1 ) );
}

sub rpois($)
{
    my $lambda = shift;

    my $z = exp(-$lambda);
    my $F = $z;

    my $u = rand;

    if( $u <= $F )
    {
        return( 0 );
    }
    else
    {
        my $k = 0;
        
        while( $u > $F )
        {
            ++$k;
            $z *= $lambda/$k;
            $F += $z;
        }

        return( $k );
    }
}

sub one_generation($$$$$)
{
    my $n = shift;  # number
    my $s = shift;  # survival
    my $e = shift;  # emigration
    my $i = shift;  # immigration
    my $m = shift;  # mutation

    my @N = @{$n};
    my @Nprime = (0,0,0);
    my $total = $N[0] + $N[1] + $N[2]; # Assume even because it's set that way.

    #my $query = shift;
    #$query->print( sprintf("%3d %3d %3d %4d", @N, $total ), $query->br );
    #$query->print( "s: ", join( ", ", @{$s} ), $query->br );
    #$query->print( "e: ", join( ", ", @{$e} ), $query->br );
    #$query->print( "i: ", join( ", ", @{$i} ), $query->br );
    #$query->print( "m: ", join( ", ", @{$m} ), $query->br );

    foreach (1..2)
    {
        while( $total > 0 )
        {
            my @pair = (0,0);

            for( my $j = 0; $j < 2; ++$j )
            {
                my $X = 1 + POSIX::floor( rand($total) );

                if( $X <= $N[0] )
                {
                    $pair[$j] = 0;
                    --$N[0];
                }
                elsif( $X <= $N[0] + $N[1] )
                {
                    $pair[$j] = 1;
                    --$N[1];
                }
                elsif( $X <= $N[0] + $N[1] + $N[2] )
                {
                    $pair[$j] = 2;
                    --$N[2];
                }

                --$total;
            }

            my $offspring = mate($pair[0], $pair[1]);

            if( (rand() <= $s->[$offspring]) &&  (rand() > $e->[$offspring]) )
            {
                if( $m->[0] > 0 || $m->[1] > 0 ) # mutation
                {
                    if( $offspring == 0 )
                    {
                        $offspring = ((rand() <= $m->[0]) ? 1 : 0) + ((rand() <= $m->[0]) ? 1 : 0);
                    }
                    elsif( $offspring == 1 )
                    {
                        $offspring = ((rand() <= $m->[0]) ? 1 : 0) + ((rand() <= $m->[1]) ? 0 : 1);
                    }
                    elsif( $offspring == 2 )
                    {
                        $offspring = ((rand() <= $m->[1]) ? 0 : 1) + ((rand() <= $m->[1]) ? 0 : 1);
                    }
                }
                
                $Nprime[$offspring] += 1;
            }
        }

        # mate again to keep the population the same
        @N = @{$n};
        $total = $N[0] + $N[1] + $N[2];
    }

    foreach my $im (0..2)
    {
        if( $i->[$im] > 0 )
        {
            $Nprime[$im] += rpois( $i->[$im] );
        }
    }

    return(@Nprime);
}

sub detaint_prop(@)
{
    my $query = shift;

    my @res = ();

    foreach (@_)
    {
        my $foo = $query->param($_);
           $foo =~ m/\A\s*(0|1|0?\.\d+)\s*\Z/;

        push @res, $1;
    }

    return( @res );
}

sub detaint_int(@)
{
    my $query = shift;

    my @res = ();

    foreach (@_)
    {
        my $foo = $query->param($_);
           $foo =~ m/(\d+)/;

        push @res, $1;
    }

    return( @res );
}

sub table_vals($)
{
    my $n = shift;
    my $N = $n->[0] + $n->[1] + $n->[2];
    $N = -1 if( $N == 0 );
    my @p = map {rnd_pc($_/$N)} @{$n};
    my @t;


    if( $N > 0 )
    {
        @t = ($n->[0], $n->[1], $n->[2], $N, @p, rnd_pc((2*$n->[0] + $n->[1])/(2*$N)), rnd_pc((2*$n->[2] + $n->[1])/(2*$N)) );
    }
    else
    {
        @t = (@$n, $N, undef, undef, undef, undef, undef );
    }

    return( @t );
}

sub print_results($)
{
    my $query = shift;

    $query->print( $query->p("And the results are in!\n"), $query->br );

    # untaint variables

    my @n = detaint_int( $query, 'nAA', 'nAa', 'naa' );
    my $m = (detaint_int( $query, 'm' ))[0];
    my @s = detaint_prop( $query, 'sAA', 'sAa', 'saa' );
    my @e = detaint_prop( $query, 'eAA', 'eAa', 'eaa' );
    my @i = detaint_int( $query, 'iAA', 'iAa', 'iaa' );
    my @m = detaint_prop( $query, 'mA', 'ma' );

    if( $m > 250 )
    {
        htmldie( "No more than 250 generations allowed, sorry." );
    }

    if( $n[0] + $n[1] + $n[2] > 1e6 )
    {
        htmldie( "No more than a total population of one million, sorry." );
    }

    my @Results = ($query->Tr( {-align=>'center'}, $query->td( [0, table_vals([@n])] ) ),);

    # Run sim

    for( my $j = 0; $j < $m; ++$j )
    {
        @n = one_generation( [@n], [@s], [@e], [@i], [@m] );

        push @Results, $query->Tr( {-align=>'center'}, $query->td( [$j+1,table_vals([@n])] ) );
    }
        

    # Display results

    $query->print( $query->table( {-cellpadding=>2},
                                 $query->Tr( {-align=>"CENTER"}, $query->th('Generation'), $query->th( {-align=>"center", -colspan=>4}, '&nbsp;&nbsp;&nbsp;&nbsp;Genotype Counts&nbsp;&nbsp;&nbsp;&nbsp;' ),
                                             $query->th( {-align=>"center", -colspan=>3}, '&nbsp;&nbsp;&nbsp;&nbsp;Genotype Frequencies&nbsp;&nbsp;' ),
                                             $query->th( {-align=>"center", -colspan=>2}, '&nbsp;&nbsp;&nbsp;&nbsp;Gene (Allele) Frequencies&nbsp;&nbsp;&nbsp;&nbsp;' ) ), 
                                 $query->Tr( {-align=>"CENTER"}, $query->th( ['', 'AA', 'Aa', 'aa', 'Total', 'AA', 'Aa', 'aa', 'A', 'a' ] ) ),
                                 @Results ), "\n" );

    $query->delete( 'mode' );
    $query->print( $query->br, $query->hr, $query->a( {href=>$query->self_url}, "New Simulation" ), "\n" );

}

sub print_foot($)
{
    my $query = shift;

    print $query->hr,"\n";
    print $query->a({href=>"http://www.stat.cmu.edu/~genovese/"},"Home"), "\n";
    print $query->end_html,"\n";
}


my $q = new CGI;
my $title = q/Hardy-Weinberg Simulator/;

print_head( $q, $title );

if( not defined( $q->param( 'mode' ) ) )
{
    print_form( $q );
}
else
{
    print_results( $q );
}

print_foot($q);



# Local Variables:
# mode: perl
# End: