#!/usr/bin/perl 

use strict;
use warnings;


my %tag_counts;
my %sort;
my %end2tag;
my %final_count;
my @final_tags = ();
my @tags = ();
my $input = "";
my $output = "";
my $routput = "";
my $total_count = 0; 
my $coesive = 0;
my $total_di = 0;
my @all_tags = ();
my $dup = 0;

my $version = 0.52;
my $vname = "counts2ditags.pl";

if (!$ARGV[0] || $ARGV[0] eq "-h") {
	print "\n\t-c <integer> \t:: join tags assuming cohesive ends of # bases [default=$coesive]\n",
	"\t-n <integer> \t:: create # ditags using reposition on the selection process [default=$total_di]\n",
	"\t-d <real> \t:: mean of duplicated ditags artificially created [default=$dup]\n",
	"\t-i <FILE> \t:: input file name[default=STDIN]\n",
	"\t-o <FILE> \t:: output file name[default=STDOUT]\n",
	"\t-r <FILE> \t:: report file (countaining the new count) output name [default=STDERR]\n",
	"\t-h      \t:: display help information\n",
	"\t-v      \t:: display current version\n",
	"\n\t$vname Version $version\n\n";
	exit;
}

if ($ARGV[0] eq "-v") {
	print ("\n\t$vname Version $version\n\n");
	exit;
}



#
#first read command line to collect arguments
#
while (my $arg=shift(@ARGV)) {
	#
	#output file
	#
	if ($arg eq "-o") {
		if ($arg = shift(@ARGV)) {
			$output = $arg;
		}
		else {
			print STDERR "ERROR: OUTPUT FILE NOT SPECIFIED!\n";
			exit;
		}
	}

	#
	#input file
	#
	elsif ($arg eq "-i") {
		if ($arg = shift(@ARGV)) {
			$input = $arg;
		}
		else {
			print STDERR "ERROR: INPUT FILE NOT SPECIFIED!\n";
			exit;
		}
	}
	#
	#report output file
	#
	elsif ($arg eq "-r") {
		if ($arg = shift(@ARGV)) {
			$routput = $arg;
		}
		else {
			print STDERR "ERROR: REPORT OUTPUT FILE NOT SPECIFIED!\n";
			exit;
		}
	}
	#
	#sets coesive end generation
	#
	elsif ($arg eq "-c") {
		if ($arg = shift(@ARGV)) {
			$coesive = $arg;
			if ($coesive < 0) {
				print STDERR "ERROR: INVALID COHESIVE END SIZE\n";
				exit;
			}
		}
	}
	#
	#set the number of ditags to create
	#
	elsif ($arg eq "-n") {
		if ($arg = shift(@ARGV)) {
			$total_di = $arg;
			if ($total_di <= 0) {
				print STDERR "ERROR: INVALID NUMBER OF DITAGS TO CREATE\n";
				exit;
			}
		}
		else {
			print STDERR "ERROR: INVALID COHESIVE END SIZE\n";
			exit;
		}
	}
	#
	#set the number of duplicated ditags created
	#
	elsif ($arg eq "-d") {
		if ($arg = shift(@ARGV)) {
			$dup = $arg;
			if ($dup <= 0) {
				print STDERR "ERROR: INVALID NUMBER OF DITAGS TO DUPLICATE\n";
				exit;
			}
		}
		else {
			$dup = 0;
		}
	}
}


if ($input) {
	open INPUT, "<$input" or die "COULD NOT OPEN $input : $!\n";
	$input = "<INPUT>";
}
else {
	$input = "<>";
}

if ($output) {
	open OUTPUT, ">$output" or die "COULD NOT OPEN $output : $!\n";
}

if ($routput) {
	open ROUTPUT, ">$routput" or die "COULD NOT OPEN $routput : $!\n";
}

#
#get the tag counts and tags
#

my $total = 0;

while (my $entry = eval ($input)){
	chomp ($entry); 
	my $tag;
	my $count;
	if ($entry =~ /=/) {
		($tag, $count) = split(/=/, $entry);
	}
	else {
		($tag, $count) = split(/ /, $entry);
	}
	$tag = &clean ($tag);
	$count = &clean ($count);
	if (!$tag || !$count) {
		next;
	}
	$tag_counts{$tag} = $count;
	$total += $count;
}

close INPUT;

foreach my $tag (keys (%tag_counts)) {
	while ($tag_counts{$tag}) {
		$tag_counts{$tag}--;
		push(@all_tags,rand()."=>".$tag);
	}
}

#
#sort the ditags by the given random number
#

my @sorted = sort (@all_tags);

foreach my $entry (@sorted) {
	my ($junk, $tag) = split ("=>", $entry);
	push (@final_tags, $tag);
}


my $togo = $total_di;
my $tolerance = 100;
if ($coesive) {
	while (@final_tags) {
		my $tag = &getTag();
		my $last2 = substr ($tag, -$coesive);
		my $last2reverse = &reverseComplement ($last2);
		if (my $tag2 = shift (@{$end2tag{$last2reverse}})){
			my $ditag = $tag2 . substr (&reverseComplement ($tag), $coesive);
			my $n = &poisson ($dup);
			for (my $i = 0; $i < $n + 1; $i++) {
				&printout ($ditag, "\n");
			}
			if (defined $final_count {$tag}) {
				$final_count {$tag}++;
			}
			else {
				$final_count {$tag} = 1;
			}
			if (defined $final_count {$tag2}) {
				$final_count {$tag2}++;
			}
			else {
				$final_count {$tag2} = 1;
			}
			if ($total_di) {
				$togo--;
				if ($togo <= 0) {
					@final_tags = ();
				}
				
			}
			$tolerance = 100;

		}
		else {
			push (@{$end2tag{$last2}},$tag);
		}
		if ($total_di && $tolerance-- == 0) {
			@final_tags = ();
			&printerr ("Couldn't find any more tags to pair up\n");
		}
	}

	while (my ($tag, $count) = each(%final_count)) {
		&printerr ("$tag=$count\n");
	}
}

else {
	while (@final_tags) {
		my $first = &getTag;
		my $second = &getTag;

		if ($second){
			my $ditag = $first . &reverseComplement ($second);
			my $n = &poisson ($dup);
			for (my $i = 0; $i < $n + 1; $i++) {
				&printout ($ditag, "\n");
			}
			if (defined $final_count {$first}) {
				$final_count {$first}++;
			}
			else {
				$final_count {$first} = 1;
			}
			if (defined $final_count {$second}) {
				$final_count {$second}++;
			}
			else {
				$final_count {$second} = 1;
			}

		}
		if ($total_di) {
			$togo--;
			if ($togo <= 0) {
				@final_tags = ();
			}
		}

	}
	if ($total_di) {
		while (my ($tag, $count) = each(%final_count)) {
			&printerr ("$tag=$count\n");
		}
	}
}

#
#return a  tag from the bin, randomly or not, as setted
#
sub getTag {
	if ($total_di) {
		my $tag = int (rand ($total));
		return $final_tags[$tag];
	}
	else {
		return shift (@final_tags);
	}
}

#
#return the reverse complement of the DNA sequence
#
sub reverseComplement {
	my $string = shift();
	my $base;
	my @split;
	@split = split (//, $string);
	$string = "";
	while (@split) {
		$base = pop (@split);
		if ($base eq "a") {
			$string .= "t";
		}
		elsif ($base eq "t") {
			$string .= "a";
		}
		elsif ($base eq "c") {
			$string .= "g";
		}
		else {
			$string .= "c";
		}
	}

	return $string;
}


sub poisson {
    my $lambda = shift;
    my $k = 0;
    my $p = 1;
    my $e = 2.71828183;
    my $L = $e**(-$lambda);

    while ($p >= $L) {
        $k++;
        $p = $p * rand();
    }

    return $k - 1;
}

sub printout {
	if ($output) {
		while (my $in = shift()) {
			print OUTPUT $in;
		}
	}
	else {
		while (my $in = shift()) {
			print $in;
		}
	}
}

sub printerr {
	if ($routput) {
		while (my $in = shift()) {
			print ROUTPUT $in;
		}
	}
	else {
		while (my $in = shift()) {
			print STDERR $in;
		}
	}
}

sub clean {
	my $string = shift;

	if ($string) {
		$string =~ s/^\s+//;
		$string =~ s/\s+$//;
	}

	return $string;
}

