# package main; # still; repackage ?

use strict;

BEGIN{ 
  # warn "# FIXME: correct this perl for ref/scaffold with multi-ref input ";
}

use constant SAMEBASE1 => 10; # for _sameloc, slop allowed in loca == locb

our ($VERSION)= "1.1";

our %EQlegend;
our ($debug);

# dang our, do we get main values here?
our $BINSIZE ;
our $GENEBINSIZE ; # use only where also checking gene ID
our $NEARDIST; #?

# sort gene exon set by ref>start>stop
sub _exlocsort {
  return $a->[2] cmp $b->[2] or $a->[3] <=> $b->[3] or $a->[4] <=> $b->[4];
}


sub printgenes_gff {
  my($genes, $geneloc, @lastid)=@_;
  my %skipid= map { $_,1 } @lastid; # may be partial
  
  foreach my $geneid (sort keys %$genes) {
    next if($skipid{$geneid});
    my($gref,$gb,$ge,$gstrand, $geq, $galt)= @{ $geneloc->{$geneid} };
    
    my $matchID= "td_$geneid"; # FIXME: option
    (my $geqtype= $EQlegend{$geq}) =~ s/\W/_/g;
    my $gat= "ID=$matchID"; 
    if($galt) { 
      $galt=~s/[\+\*\@]//g; # really should compress this to single hashed id and match w/ others
      my @galt= split ",",$galt;
      foreach (@galt) { 
        $_ = split_exonid($_,1);  # drop .exnum:location
        # s/\:[\d\-]+$//; s/\.[\d]+$//;  # drop .exnum:location
        }
      my %galt= map{$_,1} @galt; $galt= join".",sort keys %galt; # put in std order
      $gat.= ";altid=$galt"; 
      }
    print join("\t",$gref,"tandy.".$geqtype,"match", $gb, $ge,".",$gstrand,".",$gat),"\n";

    my @exons= @{ $genes->{$geneid} };
    foreach my $ex (@exons) {
      my($exnum, $xid, $ref, $b,$e, $p, $eq, $xalt)= @$ex;
      (my $xeqtype= $EQlegend{$eq}) =~ s/\W/_/g;
      my $xat= "Parent=$matchID;xid=$xid;ix=$exnum"; ##;eq=$eq
      if($galt) { #? or xalt 
        $xat.= ";altpar=$galt";  
        }
      print join("\t",$gref,"tandy.".$xeqtype,"match_part",$b,$e,$p,$gstrand,".",$xat),"\n";
      }
      
    ## clear hash
    delete $genes->{$geneid}; ## unless( grep {$_ eq $geneid} @lastid );
    delete $geneloc->{$geneid}; # this too
  }
    
}

sub printgenes2_gff {
  my($genes, $geneloc, $ibin)=@_;
  # my %skipid= map { $_,1 } @lastid; # may be partial
  
  ## for this algo, need to condense here or before the same-location gene/exon items == altids of same match
  ## cant use bin as it may cover 2+ gene locs
  
  return unless(defined $ibin);
  foreach my $geneid (sort keys %$genes) {
    # next if($skipid{$geneid});
    next unless exists $geneloc->{$geneid}{$ibin};
    
    my($gref,$gb,$ge,$gstrand, $geq, $galt)= @{ $geneloc->{$geneid}{$ibin} };
    
    my $matchID= "td_$geneid"; # FIXME: option
    (my $geqtype= $EQlegend{$geq}) =~ s/\W/_/g;
    my $gat= "ID=$matchID;tdclass=$geqtype"; 
    if($galt) { 
      $galt=~s/[\+\*\@]//g; # really should compress this to single hashed id and match w/ others
      my @galt=split ",",$galt;
      foreach (@galt) { 
        $_ = split_exonid($_,1);  # drop .exnum:location
        # s/\:[\d\-]+$//; s/\.[\d]+$//;  # drop .exnum:location
        }
      my %galt= map{$_,1}@galt; $galt= join".",sort keys %galt; # put in std order
      $gat.= ";altid=$galt"; 
      }
    print join("\t",$gref,"tandy.".$geqtype,"match", $gb, $ge,".",$gstrand,".",$gat),"\n";

    my($l_ref, $l_b, $l_e)=(0)x3;
    my @exons= @{ $genes->{$geneid}{$ibin} };
    foreach my $ex (sort _exlocsort @exons) {
      my($exnum, $xid, $ref, $b,$e, $p, $eq, $xalt)= @$ex;
      
      next if($ref eq $l_ref && $b eq $l_b && $e eq $l_e);
      ($l_ref, $l_b, $l_e)= ($ref,$b,$e);
      ## check for dups here; sort by location
      
      (my $xeqtype= $EQlegend{$eq}) =~ s/\W/_/g;
      my $xat= "Parent=$matchID;tdclass=$xeqtype;xid=$xid;ix=$exnum"; ##;eq=$eq
      if($galt) { #? or xalt 
        $xat.= ";altpar=$galt";  
        }
      print join("\t",$gref,"tandy.".$xeqtype,"match_part",$b,$e,$p,$gstrand,".",$xat),"\n";
      }
      
    ## clear hash
    delete $genes->{$geneid}{$ibin}; ## unless( grep {$_ eq $geneid} @lastid );
    delete $geneloc->{$geneid}{$ibin}; # this too
  }
    
}



=item tandem algo2

    -- not quite good enough; 
    foreach my $ib (@bins)  {
        store $genes{$geneid}{$ib} , $geneloc{$geneid}{$ib}
        }
    
    my %bins= map{ $_,1 } @bins; ## need to track back before current @bins
    foreach my $ib ( sort{$a<=>$b} keys %bin_save ) {
      next if( exists $bins{$ib} || $ib > $b - $GENEBINSIZE );
      printgenes2_gff( \%genes, \%geneloc, $ib); 
      delete $bin_save{$ib};
      }
    map { $bin_save{$_}++; }  @bins;

   ## finally:
  foreach my $ib ( sort{$a<=>$b} keys %bin_save )  {
    printgenes2_gff( \%genes, \%geneloc, $ib);
    }

=cut

=item tandem algo3

  1. find interesting ++ marked exon/gene
  
  2. collect all exons of ++marked in region i - 100 .. i + 200 or such
      save as geneid->@exons; this should include exons of tandem dupls.
      
  3. add in same/almost same exon locations for more geneids to add common id
      redo step 2. to collect new geneid exons in region
      
  4. segregate exons to gene models in tandem region (how??)
      - using exon.nums from ids, break at change in exon.num sequence?  1,2,3,4,5 | 2,3,4,5
      - using max intron size?  max/ave of step b/n base exons?
  
  5. mark as done all gene/exons in region and move on (step 1)


=cut


sub tandem_gff {
  my($exon_table, $genomefa, $queryfa, $location, $idmap)=@_;

  tandy_gffhead($exon_table, $genomefa, $queryfa, $location, $idmap);

  my (%genes, %geneloc, %bin_save, $atid, $lastid, $lastid2, $lastid3);
  my $nx= scalar(@$exon_table);
  foreach my $ix (0..$nx-1) {
  
    my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt)= @{ $exon_table->[$ix] };
    
    next unless($eq == -2 || $eq == 0  || $eq == 4); 

    my @bins= (int($tb/$GENEBINSIZE) .. int($te/$GENEBINSIZE));
    my %bins= map{ $_,1 } @bins; ## need to track back before current @bins
    foreach my $ib (@bins) {
    
    ## FIXME: now can have multiple +++IDs in alt
    while($alt =~ m/(\+{2,})([^,+]+)/g) 
    # if($alt =~ m/(\+{2,})([^,+]+)/) 
    
    {
      my ($pluss, $xid)= ($1, $2);

      ## xid may have location: ID=Dappu1_FM5_96539.2:169647-169896 ; drop it
      my($geneid,$exnum)= ($xid,1);
      if($xid =~ m/^(.+)\.(\d+)/) { $geneid=$1; $exnum=$2; }
      
      my $commonid= $idmap->{$geneid}{common} || $geneid; #??
      
      unless(exists $genes{$commonid}{$ib}) { $genes{$commonid}{$ib}=[]; }
      push( @{ $genes{$commonid}{$ib}}, [$exnum, $xid, $ref,$tb,$te, $p, $eq, $alt]);
        ## FIXME: ^ with common ids, can have many exons at same location ; test 


      my $sameloc= 0;
      if(exists $genes{$commonid}{$ib}) { 
        my $lgene= ${ $genes{$commonid}{$ib}} [-1];
        $sameloc= _sameloc( $ref, $tb, $te, 0,  @{$lgene}[2,3,4], 0);
      } else {
        $genes{$commonid}{$ib}=[]; 
      }
      push( @{ $genes{$commonid}{$ib}}, [$exnum, $xid, $ref,$tb,$te, $p, $eq, $alt])
        unless($sameloc); 

      unless(exists $geneloc{$commonid}{$ib}) { $geneloc{$commonid}{$ib}=[$ref,$tb,$te,$strand,$eq, $alt]; }
      else {
        $geneloc{$commonid}{$ib}[1]= $tb if( $tb<$geneloc{$commonid}{$ib}[1]);
        $geneloc{$commonid}{$ib}[2]= $te if( $te>$geneloc{$commonid}{$ib}[2]);
        }

      # backtrack in exontab for other exons: $ix-1,-2,-3...
      if($pluss eq '++') {
        my $need=1;
        for (my $jx= $ix-1; $jx>$ix-30; $jx--) {
          my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @{$exon_table->[$jx]};
          if ( ($jeq<=0 || $jeq>=4)  
            && ($jqid =~ /$geneid\.(\d+)/ || $jalt =~ /$geneid\.(\d+)/ )
            # || ($jqid =~ /^$commonid\.(\d+)/ || $jalt =~ /^$commonid\.(\d+)/ ) #?? FIXME
             ){
            my($jexnum, $jxid);
            $jexnum= $1; $jxid= $geneid.".".$jexnum;
            push( @{ $genes{$commonid}{$ib}}, [$jexnum,$jxid,$jref,$jb,$je, $jp, $jeq, $jalt]);
            $geneloc{$commonid}{$ib}[1]= $jb if( $jb<$geneloc{$commonid}{$ib}[1]);
            $geneloc{$commonid}{$ib}[2]= $je if( $je>$geneloc{$commonid}{$ib}[2]);
            --$need; last if($need<1);
            }
          }
        }
 
      #? change this to test region, location-bin ?? input table is sorted by location
#       $lastid3= $lastid2; $lastid2= $lastid; $lastid= $geneid;
#       
#       printgenes_gff(\%genes, \%geneloc, $lastid, $lastid2, $lastid3) 
#         if($lastid3 && $lastid ne $lastid2 && $lastid2 ne $lastid3);
      ## think this is bad; storing only by geneid in hash causes mutli-loc matches
      #   to be joined in one geneloc .. bad idea; delete geneloc{geneid} cure?
      
      }
    }

    ## this is after 1st $ib (@bins) loop; check if we can flush prior locations
    my $minb= int( ($tb - $GENEBINSIZE)/$GENEBINSIZE);
    foreach my $ib ( sort{$a<=>$b} keys %bin_save ) {
      next if( exists $bins{$ib} || $ib > $minb );
      printgenes2_gff( \%genes, \%geneloc, $ib); 
      delete $bin_save{$ib};
      }
    map { $bin_save{$_}++; }  @bins;

  }
  
  ## finally flush all remaining genes:
#  printgenes_gff(\%genes, \%geneloc);
  foreach my $ib ( sort{$a<=>$b} keys %bin_save ) {
    printgenes2_gff( \%genes, \%geneloc, $ib);
    }


}


1;