#!/home/martink/bin/perl -d:DProf

=pod

=head1 NAME

bundlelinks - bundle neighbouring links together in an effort to reduce complexity of link structure

=head1 SYNOPSIS

  bundlelinks -links linkfile.txt -max_gap NUM
              {-min_bundle_membership NUM | -strict }
              {-min_bundle_size NUM}
              {-min_bundle_identity FRACTION}

  bundlelinks -links data/dog.vs.human.all.txt > data/dog.vs.human.bundles.txt

=head1 DESCRIPTION

The purpose of this script is to turn a large set of links into a
smaller one by merging neighbouring links.

=head1 OUTPUT

The script produces a new link file to STDOUT.

A tally is sent to STDERR that lists the following

  # total number of links read in
  num_links 39978
  # number of links in initial bundles (filtered for membership)
  num_links_in_init_bundles 39839
  # total number of initial bundles
  num_init_bundles 524
  # number of accepted bundles (filtered for both size and membership)
  num_passed_bundles 277 (52.86%)
  # number of links in accepted bundles
  num_links_in_passed_bundles 36732 (92.20%)

=head1 OPTIONS

=head2 -max_gap

Adjacent links are merged into bundles if their start/end coordinates
are sufficiently close. Given two links L1 and L2, they are merged
into a bundle if

  chr( start(L1) ) == chr( start(L2) )
  chr( end(L1) ) == chr( end(L2) )

  distance( start(L1), start(L2) ) <= MAX_GAP
  distance( end(L1), end(L2) ) <= MAX_GAP

If a link does not have any acceptable adjacent neighbouring links, it forms a single-link bundle.

=head2 -min_bundle_membership, -strict

These parameters filter bundles with few links. You can set the
minimum number of links required in a bundle for the bundle to be accepted (-min_bundle_membership).

The -strict option is equivalent to -min_bundle_membership 2.

=head2 -min_bundle_size

In addition to, or in place of, filtering bundles based on the number
of links they comprise, you can accept bundles based on the size of
the links in the bundle.

The minimum size parameter is applied independently to both ends of all links in a bundle.

=head2 -min_bundle_identity

This parameter filters bundles based on the bundle identity, which is defined as

 identity = size(merged links) / extent(merged links)

Both ends of the bundle are evaluated independently.

=head1 HISTORY

=over

=item * 28 Aug 2008

=item * 16 July 2008

Started and versioned.

=back 

=head1 BUGS

=head1 AUTHOR

Martin Krzywinski

=head1 CONTACT

  Martin Krzywinski
  Genome Sciences Centre
  Vancouver BC Canada
  www.bcgsc.ca
  martink@bcgsc.ca

=cut

################################################################
#
# Copyright 2002-2008 Martin Krzywinski
#
# This file is part of the Genome Sciences Centre Perl code base.
#
# This script is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This script is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this script; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#
################################################################

################################################################
#                           Martin Krzywinski (martink@bcgsc.ca)
#                                                           2008
################################################################

use strict;
use Config::General;
use Data::Dumper;
use File::Basename;
use FindBin;
use Getopt::Long;
use IO::File;
use Math::VecStat qw(sum min max average);
use Memoize;
use Devel::DProf;
use Pod::Usage;
use Set::IntSpan;
use Time::HiRes qw(gettimeofday tv_interval);
use lib "$FindBin::RealBin";
use lib "$FindBin::RealBin/../lib";
use lib "$FindBin::RealBin/lib";
use vars qw(%OPT %CONF);

################################################################
#
# *** YOUR MODULE IMPORTS HERE
#
################################################################

#memoize("get_num_links");

GetOptions(\%OPT,
	   "links=s",
	   "max_gap=i",
	   "max_gap_1=i",
	   "max_gap_2=i",
	   "min_bundle_membership=i",
	   "min_bundle_identity=f",
	   "min_bundle_size=i",
	   "strict",
	   "zreverse",
	   "zmult=f",
	   "zshift=f",
	   "configfile=s","help","man","debug+");

pod2usage() if $OPT{help};
pod2usage(-verbose=>2) if $OPT{man};
loadconfiguration($OPT{configfile});
populateconfiguration(); # copy command line options to config hash
validateconfiguration(); 
if($CONF{debug} > 1) {
  $Data::Dumper::Pad = "debug parameters";
  $Data::Dumper::Indent = 1;
  $Data::Dumper::Quotekeys = 0;
  $Data::Dumper::Terse = 1;
  print Dumper(\%CONF);
}

my ($links,$nlinks,$chrs) = parse_links($CONF{links});
printerr("num_links",$nlinks);

my $nlinks_in_bundles;
my $in_bundle;
my $bundles;
for my $c (@$chrs) {
  # ls - all links on this chromosome
  my @ls = sort { 
    ($a->{chr}[1] cmp $b->{chr}[1])
      ||
	($a->{set}[0]->min <=> $b->{set}[0]->min)
	  ||
	    ($a->{set}[1]->min <=> $b->{set}[1]->min) } @{$links->{$c}};
 LINKI:
  for my $i ( 0 .. @ls-1 ) {
    # first link
    my $li = $ls[$i];
    next if $in_bundle->{ $li->{id} }++;
    my $bundle = [ $li ];
  LINKJ:
    for my $j ( $i+1 .. @ls-1 ) {
      # next link
      my $lj = $ls[$j];
      # gap between starts (gap1) and ends (gap2) of the links li,lj
      my $gap1 = span_distance($lj->{set}[0],get_bundle_span($bundle,0));
      my $gap2 = span_distance($lj->{set}[1],get_bundle_span($bundle,1));
      $CONF{debug} && printdebug("linki",$li->{chr}[0],$li->{set}[0]->run_list,$li->{chr}[1],$li->{set}[1]->run_list);
      $CONF{debug} && printdebug("linkj",$lj->{chr}[0],$lj->{set}[0]->run_list,$lj->{chr}[1],$lj->{set}[1]->run_list);
      my $gap_pass_1 = 1;
      my $gap_pass_2 = 1;
      $gap_pass_1 = 0 if $CONF{max_gap_1} && $gap1 > $CONF{max_gap_1};
      $gap_pass_2 = 0 if $CONF{max_gap_2} && $gap2 > $CONF{max_gap_2};
      $gap_pass_1 = 0 if $CONF{max_gap}   && $gap1 > $CONF{max_gap};
      $gap_pass_2 = 0 if $CONF{max_gap}   && $gap2 > $CONF{max_gap};
      printdebug("gaps",$gap1,$gap2,"pass",$gap_pass_1,$gap_pass_2);
      if($gap_pass_1 && $gap_pass_2 && 
	 $li->{chr}[1] eq $lj->{chr}[1]) {
	# link li,lj are within bundle gap parameters - add lj to current bundle
	# if lj has not been previously bundled
	if(! $in_bundle->{ $lj->{id} }++) {
	  $CONF{debug} && printdebug("+bundle");
	  $CONF{debug} && printdebug(get_bundle_span($bundle,0)->run_list);
	  $CONF{debug} && printdebug(get_bundle_span($bundle,1)->run_list);
	  push @$bundle, $lj;
	}
      }
      elsif ( $gap_pass_1 ) {
	# 
	next LINKJ;
      }
      else {
	last LINKJ;
      }
    }
    if(@$bundle >= $CONF{min_bundle_membership}) {
      push @$bundles, $bundle;
      $nlinks_in_bundles += @$bundle;
    }
  }
}

printerr(sprintf("num_links_in_init_bundles %d",$nlinks_in_bundles));
printerr(sprintf("num_init_bundles %d",int(@$bundles)));

my $num_passed_bundles;
my $num_links_in_passed_bundles;
for my $bundle (@$bundles) {
  my $pos1 = Set::IntSpan->new();
  my $pos2 = Set::IntSpan->new();
  my $z;
  for my $lb (@$bundle) {
    map { $pos1->U($_) } $lb->{set}[0];
    map { $pos2->U($_) } $lb->{set}[1];
    $z += $lb->{set}[0]->cardinality;
    $z += $lb->{set}[1]->cardinality;
  }
  if($CONF{min_bundle_size}) {
    next if $pos1->cardinality < $CONF{min_bundle_size};
    next if $pos2->cardinality < $CONF{min_bundle_size};
  }
  if($CONF{min_bundle_identity}) {
    next if $pos1->cardinality / ($pos1->max-$pos1->min+1) < $CONF{min_bundle_identity};
    next if $pos2->cardinality / ($pos2->max-$pos2->min+1) < $CONF{min_bundle_identity};
  }
  $num_passed_bundles++;
  $num_links_in_passed_bundles += @$bundle;
  #printerr($pos1->cardinality,
  #	   $pos2->cardinality,
  #	   int(@$bundle));
  $z /= $CONF{max_gap} ? $CONF{max_gap} : average($CONF{max_gap_1},$CONF{max_gap_2});
  $z *= $CONF{zmult} if $CONF{zmult};
  $z += $CONF{zshift} if $CONF{zshift};
  $z *= -1 if $CONF{zreverse};
  printinfo($bundle->[0]{id}."bundle",
	    $bundle->[0]{chr}[0],
	    $pos1->min,
	    $pos1->max,
	    sprintf("z=%d,color=chr%s",$z,substr($bundle->[0]{chr}[1],2)),
	   );
  printinfo($bundle->[0]{id}."bundle",
	    $bundle->[0]{chr}[1],
	    $pos2->min,
	    $pos2->max,
	    sprintf("z=%d,color=chr%s",$z,substr($bundle->[0]{chr}[1],2)),
	   );
}
printerr(sprintf("num_passed_bundles %d (%.2f%%)",
		 $num_passed_bundles,
		 100*$num_passed_bundles/int(@$bundles)));
printerr(sprintf("num_links_in_passed_bundles %d (%.2f%%)",
		 $num_links_in_passed_bundles,
		 100*$num_links_in_passed_bundles/$nlinks_in_bundles));

sub get_bundle_span {
  my $bundle   = shift;
  my $spanside = shift;
  my $bundle_min;
  my $bundle_max;
  for my $link (@$bundle) {
    if(! defined $bundle_min || $link->{set}[$spanside]->min < $bundle_min) {
      $bundle_min = $link->{set}[$spanside]->min;
    }
    if(! defined $bundle_max || $link->{set}[$spanside]->max > $bundle_max) {
      $bundle_max = $link->{set}[$spanside]->max;
    }
  }
  return Set::IntSpan->new("$bundle_min-$bundle_max");
}

sub span_distance {
  my ($s1,$s2) = @_;
  if($s1->intersect($s2)->cardinality) {
    return 0;
  } else {
    if($s1->max < $s2->min) {
      return $s2->min - $s1->max + 1;
    } elsif ($s1->min > $s2->max) {
      return $s1->min - $s2->max + 1;
    } else {
      die "strange span arrangement";
    }
  }
}

sub parse_links {
  my $file = shift;
  open(F,$file);
  my $links;
  my $chash;
  my $nlinks;
  while(<F>) {
    chomp;
    my @tok1 = split;
    my $line2 = <F>;
    last unless $line2;
    chomp $line2;
    my @tok2 = split(" ",$line2);
    my ($c1,$c2) = $tok1[1],
    my $link = {
		#list=>[\@tok1,\@tok2],
		set=>[Set::IntSpan->new(sprintf("%d-%d",@tok1[2,3])),
		      Set::IntSpan->new(sprintf("%d-%d",@tok2[2,3]))],
		id=>$tok1[0],
#		size=>[ $tok1[3]-$tok1[2]+1,
#			$tok2[3]-$tok2[2]+1 ],
#		mid=>[ ($tok1[3]+$tok1[2])/2,
#		       ($tok2[3]+$tok2[2])/2 ],
		chr=>[$tok1[1],$tok2[1]]};
    $nlinks++;
    $chash->{ $tok1[1] }++;
    push @{$links->{$tok1[1]}}, $link;
  }
  return ($links,$nlinks,[sort keys %$chash]);
}

sub parse_anchors {
  my $file = shift;
  open(F,$file);
  my $anchors;
  while(<F>) {
    chomp;
    push @$anchors, $_;
  }
  return $anchors;
}

sub validateconfiguration {
  $CONF{min_bundle_membership} = 2 if ! $CONF{min_bundle_membership} && $CONF{strict};
  $CONF{min_bundle_membership} ||= 1;
}

################################################################
#
# *** DO NOT EDIT BELOW THIS LINE ***
#
################################################################

sub populateconfiguration {
  foreach my $key (keys %OPT) {
    $CONF{$key} = $OPT{$key};
  }

  # any configuration fields of the form __XXX__ are parsed and replaced with eval(XXX). The configuration
  # can therefore depend on itself.
  #
  # flag = 10
  # note = __2*$CONF{flag}__ # would become 2*10 = 20

  for my $key (keys %CONF) {
    my $value = $CONF{$key};
    while($value =~ /__([^_].+?)__/g) {
      my $source = "__" . $1 . "__";
      my $target = eval $1;
      $value =~ s/\Q$source\E/$target/g;
      #printinfo($source,$target,$value);
    }
    $CONF{$key} = $value;
  }

}

sub loadconfiguration {
  my $file = shift;
  my ($scriptname) = fileparse($0);
  if(-e $file && -r _) {
    # great the file exists
  } elsif (-e "/home/$ENV{LOGNAME}/.$scriptname.conf" && -r _) {
    $file = "/home/$ENV{LOGNAME}/.$scriptname.conf";
  } elsif (-e "$FindBin::RealBin/$scriptname.conf" && -r _) {
    $file = "$FindBin::RealBin/$scriptname.conf";
  } elsif (-e "$FindBin::RealBin/etc/$scriptname.conf" && -r _) {
    $file = "$FindBin::RealBin/etc/$scriptname.conf";
  } elsif (-e "$FindBin::RealBin/../etc/$scriptname.conf" && -r _) {
    $file = "$FindBin::RealBin/../etc/$scriptname.conf";
  } else {
    return undef;
  }
  $OPT{configfile} = $file;
  my $conf = new Config::General(-ConfigFile=>$file,
				 -AllowMultiOptions=>"yes",
				 -LowerCaseNames=>1,
				 -AutoTrue=>1);
  %CONF = $conf->getall;
}

sub printdebug {
  printinfo("debug",@_)  if $CONF{debug};
}

sub printinfo {
  printf("%s\n",join(" ",@_));
}

sub printerr {
  printf STDERR ("%s\n",join(" ",@_));
}