use strict; use warnings;
use 5.010;

use List::Util qw(min max);

open my $ov_debug, '>', 'ov_debug.txt' or die $!;
open my $valid_fast, '<', 'all.Yc2h2.nrdb.fasta' or die $!;

sub overlap {
	return () if ( $_[0][0] > $_[1][-1]
				   or $_[0][-1] < $_[1][0]
				   or $_[1][0] > $_[0][-1]
				   or $_[1][-1] < $_[0][0] );
	my ($small_start, $small_end) = (max($_[0][0],$_[1][0]),min($_[0][-1],$_[1][-1]));
	my ($s1_len, $small_len, $s2_len) = ($_[0][-1]-$_[0][0], $small_end-$small_start, $_[1][-1]-$_[1][0]);
	my ($s1_lap, $s2_lap) = ($small_len/$s1_len, $small_len/$s2_len);
	my ($large_start, $large_end) = (min($_[0][0],$_[1][0]),max($_[0][-1],$_[1][-1]));
	say $ov_debug "---";
	say $ov_debug "s1@{$_[0]}\tlen$s1_len\tlap$s1_lap\t",
		map { $_ >= $_[0][0] && $_ <= $_[0][-1] ? '|' : ' ' } ($large_start .. $large_end);
	say $ov_debug "s2@{$_[1]}\tlen$s2_len\tlap$s2_lap\t",
		map { $_ >= $_[1][0] && $_ <= $_[1][-1] ? '|' : ' ' } ($large_start .. $large_end);
	say $ov_debug "$large_start - $small_start --> $small_len,v",(($s1_lap**2+$s2_lap**2)**(1/1.7))," <-- $small_end - $large_end";
	($large_start, (($s1_lap**2+$s2_lap**2)**(1/1.7)), $large_end)
}

sub in_cseq {
	for my $sl (@{$_[0]}) {
		for (@_[1..$#_]) {
			return 1 if $sl == $_;
		}
	}
	return;
}

sub insert_into_cseq {
	my $target;
	for my $cseq (@{$_[0]}) {
		if (in_cseq($cseq, @_[1..$#_])) {
			$target = $cseq;
			last;
		}
	}
	unless (defined $target) {
		push@{$_[0]},[@_[1..$#_]];
	}
	else {
		push @$target, grep { !in_cseq($target, $_) } @_[1..$#_];
	}
}

my ($mod, $mod_id, $id, $gene);
my (%gene_an);
my %valid_tbl;

say STDERR 'valid fasta...';
while (<$valid_fast>) {
	chomp;
	if (my ($res) = /^>(.*)/) {
		my @ids = split '[+]', $res;
		for (@ids) {
			$valid_tbl{$_}=undef;
		}
	}
}

use Data::Dumper;

#print Dumper \%valid_tbl;

say STDERR 'loading...';
while (<>) {
	chomp;
	if (my @res = /^Modell: (\d+)\.hmm_search (..) ID: (\d+) >> (.*?)\s+/) {
		($mod, $mod_id, $id, $gene) = @res;
	}
	elsif (my ($sid, $start, $end) = /^(\d+) S:(\d+) E:(\d+)$/) {
		push @{ $gene_an{$gene} }, [ $start, $mod_id, $id, $sid, $end ]
			if exists $valid_tbl{(sprintf '%.2s%04d%03d', $mod_id, $id, $sid)};
	}
	else {
		warn "unreco $ARGV ^$_\$\n";
	}
}

my $j=0; for $gene (sort keys %gene_an) {
	my $i=0; print STDERR $gene,'(',++$j,'/',scalar keys %gene_an,')', "~", ' 'x length("$i/@{[scalar@{$gene_an{$gene}}]}");
	my $can = [];

	for my $s1pos (@{ $gene_an{$gene} }) {
		print STDERR "\010 \010"x length("$i/@{[scalar@{$gene_an{$gene}}]}"),++$i,'/',scalar @{ $gene_an{$gene} };
		for my $s2pos (@{ $gene_an{$gene} }) {
			next if $s1pos == $s2pos;
			my ($start, $lap, $end) = overlap($s1pos, $s2pos);
			if (defined $lap && $lap >= 7/10) {
				say $ov_debug "MAP";
				insert_into_cseq($can, $s1pos, $s2pos);
				say $ov_debug "s1:@{$s1pos}, s2:@{$s2pos}, [$start, 'debug', $lap, $end]";
			}
			elsif (defined $lap) {
				say $ov_debug "NOM";
			}
		}
		insert_into_cseq($can, $s1pos);
	}
	print STDERR "\010 \010"x(length("$i/@{[scalar@{$gene_an{$gene}}]}")+length("$gene($j/@{[scalar keys %gene_an]})~"));

	say $gene;
	for my $cseq_exp (sort {$a->[0] <=> $b->[0] || $b->[-1] <=> $a->[-1] } map { [ (min map { $_->[0] } @$_), $_, (max map { $_->[-1] } @$_) ] } @{ $can }) {
		say "\t", $cseq_exp->[0], "\t", $cseq_exp->[-1];
		say "\t"x3, $_->[0], "\t", (join ".", @{$_}[1..$#$_-1]), "\t", $_->[-1], "\t",
			do { my $l=$_; map { $_ >= $l->[0] && $_ <= $l->[-1] ? '|' : ' ' } ($cseq_exp->[0] .. $cseq_exp->[-1]) }
				for @{$cseq_exp->[1]};
	}
	say '==========';
	print $ov_debug Dumper $can;
}


