# perl fasta.pl 25000000

# The Computer Language Benchmarks game
# http://benchmarksgame.alioth.debian.org/
#
# contributed by Barry Walsh
# port of fasta.rb #6
#
# MCE version by Mario Roy

use strict;
use warnings;
use feature 'say';

use MCE;
use MCE::Candy;
use MCE::Shared;

use constant IM => 139968;
use constant IA => 3877;
use constant IC => 29573;

my $LAST = MCE::Shared->scalar(42);

my $alu =
    'GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG' .
    'GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA' .
    'CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT' .
    'ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA' .
    'GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG' .
    'AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC' .
    'AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA';

my $iub = [
    [ 'a', 0.27 ], [ 'c', 0.12 ], [ 'g', 0.12 ],
    [ 't', 0.27 ], [ 'B', 0.02 ], [ 'D', 0.02 ],
    [ 'H', 0.02 ], [ 'K', 0.02 ], [ 'M', 0.02 ],
    [ 'N', 0.02 ], [ 'R', 0.02 ], [ 'S', 0.02 ],
    [ 'V', 0.02 ], [ 'W', 0.02 ], [ 'Y', 0.02 ]
];

my $homosapiens = [
    [ 'a', 0.3029549426680 ],
    [ 'c', 0.1979883004921 ],
    [ 'g', 0.1975473066391 ],
    [ 't', 0.3015094502008 ]
];

sub make_repeat_fasta {
    my ($src, $n) = @_;
    my $width = qr/(.{1,60})/;
    my $l     = length $src;
    my $s     = $src x (($n / $l) + 1);
    substr($s, $n, $l) = '';

    while ($s =~ m/$width/g) { say $1 }
}

sub make_random_fasta {
    my ($table, $n) = @_;
    my $rand   = undef;
    my $width  = 60;
    my $prob   = 0.0;
    my $output = '';
    my ($c1, $c2, $last);

    $_->[1] = ($prob += $_->[1]) for @$table;

    $c1  = '$rand = ($last = ($last * IA + IC) % IM) / IM;';
    $c1 .= "\$output .= '$_->[0]', next if $_->[1] > \$rand;\n" for @$table;

    my $code1 = q{
        my ($mce, $seq, $chunk_id) = @_;

        # process code-snippet orderly between workers
        MCE->relay_recv;

        my $last = $LAST->get;
        my $temp = $last;

        # pre-compute $LAST for the next worker
        for (1 .. ($seq->[1] - $seq->[0] + 1) * $width) {
            $temp = ($temp * IA + IC) % IM;
        }

        $LAST->set($temp);
        MCE->relay;

        # process code-snippet in parallel
        for ($seq->[0] .. $seq->[1]) {
            for (1..$width) { !C! }
            $output .= "\n";
        }

        # gather output orderly
        MCE->gather($chunk_id, $output);

        $output = '';
    };

    $code1 =~ s/!C!/$c1/g;

    MCE->new(
        bounds_only => 1,
        chunk_size  => 2000,
        init_relay  => 0,
        max_workers => 4, ## MCE::Util->get_ncpu || 4,
        sequence    => [ 1, ($n / $width) ],
        gather      => MCE::Candy::out_iter_fh(\*STDOUT),
        user_func   => sub { eval $code1; },
        use_threads => 0,
    )->run;

    $last = $LAST->get;

    $c2  = '$rand = ($last = ($last * IA + IC) % IM) / IM;';
    $c2 .= "print('$_->[0]'), next if $_->[1] > \$rand;\n" for @$table;

    my $code2 = q{
        if ($n % $width != 0) {
            for (1 .. $n % $width) { !C! }
            print "\n";
        }
    };

    $code2 =~ s/!C!/$c2/g;
    eval $code2;

    $LAST->set($last);
}

my $n = $ARGV[0] || 27;

say ">ONE Homo sapiens alu";
make_repeat_fasta($alu, $n*2);

say ">TWO IUB ambiguity codes";
make_random_fasta($iub, $n*3);

say ">THREE Homo sapiens frequency";
make_random_fasta($homosapiens, $n*5);