#! /usr/bin/perl -w

# Author : HqD
# Version : 0.0.1

use strict;
#use Switch;
use POSIX qw(ceil floor);

&alphabetializing();

sub alphabetializing {
	my $flag = 0;
	my $in_file = ""; my $out_file = "";
	my $letter = 5;

	foreach my $arg (@ARGV){
		if ($arg =~ /^-[i]/) {$flag = 1; next;}
		elsif ($arg =~ /^-[k]/) {$flag = 3; next;}
		
		if ($flag==1) {$in_file = $arg;next;}
		elsif ($flag==3) {$letter = convert($arg);next;}
	}
	
	print "Running multiple marks in one file ...\n";
	multiple_mark($in_file, $letter);
}

sub multiple_mark {
	my ($in_file, $letter) = @_;
	open(IN, $in_file ) or die "$!\n";
	my $N_regs = 0; my $N_marks = 0; my @mark_names = undef;
	my @coors = undef; my @signals = undef;
	
	my @str = <IN>; close(IN);
	$N_regs = scalar(@str)-1;
	
	my $ind_reg = 0;
	open(IN, $in_file ) or die "$!\n";
	while (<IN>) {
		chomp;
		my @line=split /\s+/,$_;
		if ($N_marks == 0) {
			$N_marks = @line; $N_marks = $N_marks-1;
			for (my $i=1; $i<=$N_marks; $i++) {$mark_names[$i-1] = $line[$i];}
			print "@mark_names\n"; 

			for (my $i=0; $i<$N_marks; $i++) {
				my @temp = undef;
				for (my $j=0; $j<$N_regs; $j++) {$temp[$j] = 0;}
				$signals[$i] = [ @temp ];
			}
			print $N_regs, " regions & ", $N_marks, " histone marks\n";
		}
		else {
				$coors[$ind_reg] = $line[0];
				for (my $i=1; $i<=$N_marks; $i++) {
					if ($line[$i] ne "NA") { $signals[$i-1][$ind_reg] = sprintf("%.15f", $line[$i]);}
					else {$signals[$i-1][$ind_reg] = 1000;}
				}
				$ind_reg++;
		}
	}
	close(IN);
	
	for (my $mark=0; $mark<$N_marks; $mark++) {
		my @index = undef; my @region = undef; my $N_temp = 0;
		for (my $i=0; $i<$N_regs; $i++) {
			if ($signals[$mark][$i]<1000) {
				$index[$N_temp] = $i; $region[$N_temp] = $signals[$mark][$i];
				$N_temp++;
			}
		}
		quicksort( \@region, \@index);
		my $out_file = $mark_names[$mark].".out";
		open(OUT, ">$out_file") or die "$out_file: $!\n";
		my $i=0; my $value = -1000; my $N_region_new = 0;
		#print "@region\n";
		#print "@index\n"; exit;
		while($i<$N_temp) {
			if ($region[$i] != $value) {	# for the special case of DNAmeth where a lot of regions have the same signal
				$value = $region[$i];
				$N_region_new++;
			} 
			$i++;
		}
		print OUT $N_region_new, "\t", $letter, "\n"; 
		$i=0; $value = -1;
		while($i<$N_temp) {
			if ($region[$i] != $value) {
				print OUT $region[$i], "\n";
				$value = $region[$i];
			} 
			$i++;
		}
		close(OUT);	
		system("./epi_letter -f ".$out_file." -o ".substr($out_file,0,index($out_file,".")).".cutoff");
		system("rm ".$out_file);
	}
}


sub convert {
  my $strNumber = $_[0];
  $strNumber = sprintf("%04d", $strNumber);
  my $number = $strNumber -1 + 1;
  return $number;
}

sub partition {
    my ( $array, $first, $last, $index_arr ) = @_;
	
    my $i = $first;
    my $j = $last - 1;
    my $pivot = $array->[ $last ];
	
SCAN: {
	do {
		# $first <= $i <= $j <= $last - 1
		# Point 1.
		
		# Move $i as far as possible.
		while ( $array->[ $i ] < $pivot ) {  
			$i++;
			last SCAN if $j < $i;
		}
		
		# Move $j as far as possible.
		while ( $array->[ $j ] > $pivot ) {
			$j--;
			last SCAN if $j < $i;
		}
		
		# $i and $j did not cross over, so swap a low and a high value.
		@$array[ $j, $i ] = @$array[ $i, $j ];
		@$index_arr [ $j, $i] = @$index_arr [ $i, $j ];
	} while ( --$j >= ++$i );
}
    # $first - 1 <= $j < $i <= $last
    # Point 2.
	
	# Swap the pivot with the first larger
	# element (if there is one).
    if ( $i < $last ) {
        @$array[ $last, $i ] = @$array[ $i, $last ];
		@$index_arr[ $last, $i ] = @$index_arr[ $i, $last ];
        ++$i;
    }
	
    # Point 3.
	
    return ( $i, $j );   # The new bounds exclude the middle.
}

sub quicksort_recurse {
    my ( $array, $first, $last, $index_arr ) = @_;
	
    if ( $last > $first ) {
        my ( $first_of_last, $last_of_first ) =
		partition( $array, $first, $last, $index_arr );
		
        local $^W = 0; # Silence deep recursion warning.
        quicksort_recurse($array, $first,         $last_of_first, $index_arr);
        quicksort_recurse($array, $first_of_last, $last, $index_arr);
    }
}

sub quicksort {
	# The recursive version is bad with BIG lists
	# because the function call stack gets REALLY deep.
	
	quicksort_recurse($_[ 0 ], 0, $#{ $_[ 0 ]} , $_[ 1 ]);
}