#!/usr/bin/env perl
use warnings;
use strict;
use POSIX;

my $f_chr_len = $ARGV[0]; # chromosome length file for whole genome (e.g., ./data/rg_human_hg19.len)
my $resolution = $ARGV[1]; # the resolution of interest: the bin size to evenly split a chromosome (e.g., 1000000)
my $dir_out = $ARGV[2]; # files in this folder with suffixs:"fends_GC", "indexes_mapped_uniquely", and "fends_mappability"

######################## get chr len
my %chr_lens;
open IN, "$f_chr_len";
while (my $line = <IN>) {
  chomp $line; $line =~ s/^\s+//; $line =~ s/\s+$//;
  my @items = split(/\s+/, $line);
  my $chr = $items[0];
  $chr =~ s/chr//;
  $chr_lens{$chr} = $items[1];
}
close IN;

for my $chr (keys %chr_lens) {

  my %fends_GC;
  my $f_GC = "$dir_out/$chr.fends_GC";
  if (! -e $f_GC) { next; }
  open IN, "$f_GC" or die $!;
  while (my $line = <IN>) {
    chomp $line; $line =~ s/^\s+//; $line =~ s/\s+$//;
    my @items = split(/\s+/, $line);
    my $key = $items[0] + 1;
    push(@{$fends_GC{$key}}, $items[1]);
    push(@{$fends_GC{$key}}, $items[2]);
  }
  close IN;
  my %fends_map;
  my $f_map = "$dir_out/$chr.fends_mappability";
  if (! -e $f_map) { next; }
  open IN, "$f_map" or die $!;
  while (my $line = <IN>) {
    chomp $line; $line =~ s/^\s+//; $line =~ s/\s+$//;
    my @items = split(/\s+/, $line);
    my $key = $items[0];
    push(@{$fends_map{$key}}, $items[1]);
    push(@{$fends_map{$key}}, $items[2]);
  }
  close IN;

  my $chr_len = $chr_lens{$chr};
  my $num_bins = ceil($chr_len / $resolution);
  # feature 1: density - number of cutting sites in each bin
  my (%bin_density, %bin_sum_GC, %bin_sum_map);
  for my $i (1..$num_bins) {
    $bin_density{$i} = 0;
  }
  for my $key (sort {$a<=>$b} keys %fends_GC) {
    my $tmp_bin = ceil($key / $resolution);
    $bin_density{$tmp_bin}++;
    my ($bin1, $bin2);
    $bin1 = ceil(($key - 200) / $resolution);
    $bin2 = ceil(($key - 1) / $resolution);
    if ($bin1 == $bin2) { push(@{$bin_sum_GC{$bin1}}, $fends_GC{$key}->[0]); }
    $bin1 = ceil(($key + 1) / $resolution);
    $bin2 = ceil(($key + 200) / $resolution);
    if ($bin1 == $bin2) { push(@{$bin_sum_GC{$bin1}}, $fends_GC{$key}->[1]); }
    $bin1 = ceil(($key - 500) / $resolution);
    $bin2 = ceil(($key - 1) / $resolution);
    if ($bin1 == $bin2) { push(@{$bin_sum_map{$bin1}}, $fends_map{$key}->[0]); }
    $bin1 = ceil(($key + 1) / $resolution);
    $bin2 = ceil(($key + 500) / $resolution);
    if ($bin1 == $bin2) { push(@{$bin_sum_map{$bin1}}, $fends_map{$key}->[1]); }
  }
  my (%bin_GC, %bin_map); 
  for my $key (sort {$a<=>$b} keys %bin_sum_GC) {
    my $n = @{$bin_sum_GC{$key}};
    my $sum = 0;
    for (my $i = 0; $i < $n; $i++) {
      $sum += $bin_sum_GC{$key}->[$i];
    }
    $bin_GC{$key} = sprintf "%0.3f", $sum / $n;
  }
  for my $key (sort {$a<=>$b} keys %bin_sum_map) {
    my $n = @{$bin_sum_map{$key}};
    my $sum = 0;
    for (my $i = 0; $i < $n; $i++) {
			if ($bin_sum_map{$key}->[$i] eq 'NA'){
      	$sum += 0;
			} else {
      	$sum += $bin_sum_map{$key}->[$i];
			}
    }
    $bin_map{$key} = sprintf "%0.3f", $sum / $n;
  }
  # output
  my $file_output = "$dir_out/$chr.bin_features_$resolution";
  open OUT, ">$file_output";
  for my $i (1..$num_bins) {
    my $bin_start = ($i - 1) * $resolution + 1;
    my $bin_end = $i * $resolution;
    if ($i == $num_bins) { $bin_end = $chr_len; }
    print OUT "$bin_start\t$bin_end\t$bin_density{$i}\t";
    if (exists $bin_GC{$i}) {
      print OUT "$bin_GC{$i}\t";
    } else {
      print OUT "0\t";
    }
    if (exists $bin_map{$i}) {
      print OUT "$bin_map{$i}\n";
    } else {
      print OUT "0\n";
    }
  }
  close OUT;
}

# End

