#!/usr/local/bin/perl
# -*-Perl-*- 
#######################################################################
# This software has been created by Genome Research Limited (GRL).    # 
# GRL hereby grants permission to use, copy, modify and distribute    # 
# this software and its documentation for non-commercial purposes     # 
# without fee at the user's own risk on the basis set out below.      #
# GRL neither undertakes nor accepts any duty whether contractual or  # 
# otherwise in connection with the software, its use or the use of    # 
# any derivative, and makes no representations or warranties, express #
# or implied, concerning the software, its suitability, fitness for   #
# a particular purpose or non-infringement.                           #
# In no event shall the authors of the software or GRL be responsible # 
# or liable for any loss or damage whatsoever arising in any way      # 
# directly or indirectly out of the use of this software or its       # 
# derivatives, even if advised of the possibility of such damage.     #
# Our software can be freely distributed under the conditions set out # 
# above, and must contain this copyright notice.                      #
#######################################################################
# 
# $Id: Exptools.pm,v 1.2 1997/08/19 16:07:03 rmd Exp $
# $Log: Exptools.pm,v $
# Revision 1.2  1997/08/19  16:07:03  rmd
# Fixed iteration error
#
# Revision 1.1  1997/06/09  15:44:48  rmd
# Initial revision
#
#

package Exptools;

use Carp;
use strict;
use Uncovered;

sub new {
    my $type = shift;
    my $self = {expOrder => [], expNentries => {}, expEntry => {}, iterators => {}};
    return bless $self, $type;
}

sub readExpFile {
    my ($self, $file, $dontWarn) = @_;
    confess "Wrong type" unless ref $self;
    $self->{expOrder} = [];
    $self->{expNentries} = {};
    $self->{expEntry} = {};
    my $id;
    my $firstLine = <$file>;
    if (($id) = $firstLine =~ /^ID\s+(\S+)/) {
	# experiment file
	$self->addEntry('ID', $id);
	$self->_readExpFormat($file);
    } elsif (/^;/) {
	my ($lt, $ln) = /;\s*\d+\s*\d+\s*\d+(....)(\S+)\s*$/;
	$self->addEntry('LT', $lt);
	$self->addEntry('LN', $ln);
	$self->_readStadenFormat($file);
    } elsif (/^>(\S)/) {
	$self->addEntry('ID', $1);
	$self->_readFastaFormat($file);
    } elsif ($self->_readPlainSeq($file, $_)) {
	warn "readExpFile: unknown format - cannot read $file\n" unless $dontWarn;
	return 1;
    }
    return 0;
}

sub _readExpFormat {
    my ($self, $file) = @_;
    confess "Wrong type" unless ref $self;
    my $multiline = "";
    while (<$file>) {
	my ($id) = /^(..)/;
	$id = $multiline if ($id =~ /^\s+$/ && $multiline ne "");
	if ($id eq "SQ") {
	    # Special treatment for sequences
	    push(@{$self->{expOrder}}, $id);
	    my $n = $self->{expNentries}->{$id};
	    $self->{expEntry}->{$id} = [] unless $self->{expEntry}->{$id};
	    $self->{expEntry}->{$id}->[$n] = "";
	    $multiline = "";
	    while (<$file>) {
		last if /^\/\//;
		$_ =~ s/\s//g; # Kill white space
####################$_ =~ s/\-/N/g; # Replace '-' by 'N' 
		$self->{expEntry}->{$id}->[$n] .= $_;
	    }
	    $self->{expNentries}->{$id}++;
	} else {
	    my $value = substr($_, 5);
	    my $n = $self->{expNentries}->{$id};
	    if ($id eq $multiline && $value =~ /^     /) {
		$self->{expEntry}->{$id}->[$n-1] .= $value;
	    } else {
		push(@{$self->{expOrder}}, $id); #Keep a track of line order
		$multiline = "";
		$self->{expEntry}->{$id}->[$n] = $value;
		$self->{expNentries}->{$id}++;
		if ($id eq "TG" || $id eq "TC" || $id eq "ON" || $id eq "AV") { $multiline = $id; }
	    }
	}
    }
}

sub _readStadenFormat {
    my ($self, $file) = @_;
    confess "Wrong type" unless ref $self;
    my ($left, $right, $seq) = ("", "", "");
    while (<$file>) {
	if (/^;</) {
	    $left .= substr($_,2);
	} elsif (/^;>/) {
	    $right .= substr($_,2);
	} elsif (/^;;/) {
	    #;;ALUS      1    146      0.873
	    my ($type,$from,$len,$comment) = /^;;(....)\s+(\d+)\s+(\d+)\s*(.*)$/;
	    my $to = $from + $len - 1;
	    
	    # FIXME YUK YIKES BUG HACK
	    # BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG
	    # In staden format, tag positions are relative to the first used base (less one)
	    # and must not fall outside the reading.
	    # Currently this is not supported
	    # BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG
	    
	    if (defined($comment)) {
		$self->addEntry("TG","$type + $from..$to\n     $comment");
	    } else {
		$self->addEntry("TG","$type + $from..$to");
	    }
	} elsif (/^;/) {
	    # Comment
	    $self->addEntry("CC",substr($_,1));
	} else {
	    $seq .= $_;
	}
    }
    # remove white space
    $left =~ s/\s//g;  
    my $leftlen = length($left);
    $right =~ s/\s//g; 
    my $rightlen = length($right);
    $seq =~ s/\s//g;   
    my $seqlen = length($seq);
    #Now dump out sequence...
    $seq = $left . $seq . $right;
    $self->addEntry("QL",$leftlen);
    $self->addEntry("QR",$leftlen+$seqlen+1);
    $self->addEntry("SQ",$seq);    
}

sub _readFastaFormat {
    my ($self, $file) = @_;
    confess "Wrong type" unless ref $self;
    my $seq = "";
    while(<$file>) {
	chomp();
	s/\s*//;
	$seq .= $_;
    }
    $self->addEntry('SQ', $seq) if $seq;
}

sub _readPlainSeq {
    my ($self, $file, $firstLine) = @_;
    confess "Wrong type" unless ref $self;
    return 1 unless $firstLine =~ /^[\s\*ACGT\-acgtNn]+$/;
    $firstLine =~ s/\s*//g;
    my $seq = $firstLine;
    while (<$file>) {
	return 1 unless /^[\s\*ACGT\-acgtNn]+$/;
	s/\s*//g;
	$seq .= $_;
    }
    $self->addEntry('SQ', $seq);
}

sub getFirstEntry {
    my ($self, $id) = @_;
    confess "Wrong type" unless ref $self;
    unless ($self->{expNentries}->{$id}) { return undef;}
    $self->{iterators}->{$id} = 1;
    return $self->{expEntry}->{$id}->[0];
}

sub getNextEntry {
    my ($self, $id) = @_;
    confess "Wrong type" unless ref $self;
    my $i = $self->{iterators}->{$id};
    unless ($i) { return $self->getFirstEntry($id);}
    my $n = $self->{expNentries}->{$id};
    unless ($n && $i < $n) { return undef;}
    return $self->{expEntry}->{$id}->[$self->{iterators}->{$id}++];
}

sub getEntry {
    my ($self, $id) = @_;
    confess "Wrong type" unless ref $self;
    my $n = $self->{expNentries}->{$id};
    return undef unless $n;
    return $self->{expEntry}->{$id}->[$n - 1];
}

sub addEntry {
    my ($self, $id, $value) = @_;
    confess "Wrong type" unless ref $self;
    $value .= "\n" unless $value =~ /\n$/;
    my $n = $self->{expNentries}->{$id};
    $self->{expEntry}->{$id} = [] unless $self->{expEntry}->{$id};
    $self->{expEntry}->{$id}->[$n] = $value;
    $self->{expNentries}->{$id}++;
    push(@{$self->{expOrder}}, $id);
}

sub replaceEntry {
    my ($self, $id, $value) = @_;
    confess "Wrong type" unless ref $self;
    $self->killEntry($id);
    $self->addEntry($id, $value);
}

sub killEntry {
    my ($self, $id) = @_;
    confess "Wrong type" unless ref $self;
    $self->{expEntry}->{$id} = [] if $self->{expEntry}->{$id};
    $self->{expNentries}->{$id} = 0 if $self->{expNentry}->{$id};
    $self->{expOrder} = grep { $_ ne $id; } @{$self->{expOrder}};
}

sub writeExpFile {
    my ($self, $file) = @_;
    confess "Wrong type" unless ref $self;
    my $id = $self->getEntry('ID');
    print $file "ID   $id" if $id;
    my %count;
    foreach $id (@{$self->{expOrder}}) {
	if ($id eq "ID") {
	    # skip
	} elsif ($id eq 'SQ') {
	    print $file "$id\n";
	    my $seq = $self->{expEntry}->{$id}->[$count{$id}];
	    $seq =~ s/\s//g;
	    my $len = length($seq);
	    my $j;
	    for ($j = 0; $j < $len;) {
		print $file "    ";
		my $k;
		for ($k = 0; $k < 6 && $j < $len; $k++, $j += 10) {
		    print $file " ",substr($seq, $j, 10);
		}
		print $file "\n";
	    }
	    print $file "//\n";
	} elsif ($id eq "AV") {
	    my $av_count = 1;
	    foreach (split(/\n/, $self->{expEntry}->{$id}->[$count{$id}])) {
		if ($av_count) {
		    print $file "$id   $_\n";
		    $av_count = 0;
                } else {
		    print $file "$id        $_\n";
                }	
	    }
	} else {
	    foreach ( split(/\n/,$self->{expEntry}->{$id}->[$count{$id}])) {
		print $file "$id   $_\n";
	    }
	}
	$count{$id}++;
    }
}

sub determineActiveRegion {
    my ($self) = @_;
    confess "Wrong type" unless ref $self;
    my $SQ = $self->getEntry('SQ');
    warn "determineActiveRegion - No Sequence\n" unless $SQ;
    my $SQlen = ($SQ)?length($SQ):4096;
    my $SL = $self->getEntry('SL');
    my $SR = $self->getEntry('SR');
    my $QL = $self->getEntry('QL');
    my $QR = $self->getEntry('QR');
    $SL = 0 unless $SL;
    $SR = $SQlen + 1 unless $SR;
    $QL = 0 unless $QL;
    $QR = $SQlen + 1 unless $QR;

    # Find where the good bits are
    my $l = ($SL > $QL) ? $SL : $QL;
    my $r = ($SR < $QR) ? $SR : $QR;
    
    my $CS = $self->getEntry('CS');
    if ($CS) {
	my ($CSfrom, $CSto) = $CS =~ /^(\d+)..(\d+)$/;
        #/*
        # * Six posibilities:
        # *           l                r
        # * 1  -----  >                <
        # * 2  -------|------->        <
        # * 3  -------|<---------------|-------
        # * 4         >    <--------   |
        # * 5         >    <-----------|-------
        # * 6         >                <  -----
        # *
        # * key: [-] CS [>] new left c.o [<] new right c.o [|] orig c.o
        # *
        # * Coded for understanding not efficiency!!!
        # */
        if ($CSfrom <= $l && $CSto <= $l) { #/*1*/
            $l = $l; #/* null statement */
        } elsif ($CSfrom <= $l+1 && $CSto < $r) { #/*2*/
            $l = $CSto;
        } elsif ($CSfrom <= $l+1 && $CSto >= $r) { #/*3*/
            $r = $l+1;
        } elsif ($CSfrom < $r && $CSto < $r) { #/*4*/
            $r = $CSfrom;
        } elsif ($CSfrom < $r && $CSto >= $r) { #/*5*/
            $r = $CSfrom;
        } elsif ($CSfrom >= $r && $CSto >= $r) { #/*6*/
            $l = $l; #/* null statement */
	}
    }
    if ($l>$r) { $l = $r-1; } #/* just a quick consistency check */
    return ($l+1, $r-1);
}

############################################################
# ($left, $right) = $expFile->myDetermineActiveRegion($mode)
#
# This is a rewrite of determineActiveRegion() using &cover()
#
# Mode:
#     -c  - mask cloning vector (CS)
#     -C  - mask everything but the cloning vector
#     -ql - mask poor quality regions (1..QL)
#     -qr - mask poor quality regions (QR..END)
#     -sl - mask sequencing vector (1..SL)
#     -sr - mask sequencing vector (SR..END)
#     -xl - mask using quality revert Squirrels final clip points
#     -xr - mask using quality revert Squirrels final clip points
#
############################################################

sub myDetermineActiveRegion {
    my ($self, $switches) = @_;
    confess "Wrong type" unless ref $self;
    my $SQ = $self->getEntry('SQ');
    warn "determineActiveRegion - No sequence\n" unless $SQ;
    my $SQlen = ($SQ) ? length($SQ) : 4096;
    my $cover = new Uncover(1, $SQlen);
    my $SL = ($switches =~ /-sl/) ? $self->getEntry('SL') : undef;
    $cover->cover(0, $SL) if defined($SL);
    my $SR = ($switches =~ /-sr/) ? $self->getEntry('SR') : undef;
    $cover->cover($SR, $SQlen + 1) if defined($SR);
    my $QL = ($switches =~ /-ql/) ? $self->getEntry('QL') : undef;
    $cover->cover(0, $QL) if defined($QL);
    my $QR = ($switches =~ /-qr/) ? $self->getEntry('QR') : undef;
    $cover->cover($QR, $SQlen + 1) if defined($QR);
    my $CS = ($switches =~ /(?:-c|-C)/) ? $self->getEntry('CS') : undef;
    if (defined($CS)) {
	my ($CSfrom,$CSto) = $CS =~ /^(\d+)..(\d+)$/;
	if ($switches =~ /-c/) {
	    $cover->cover($CSfrom,$CSto);
	} elsif ($switches =~ /-C/) {
	  $cover->cover(0,$CSfrom-1);
	  $cover->cover($CSto+1,$SQlen+1);  
	}
    }
    # Precheck
    # We should ignore overall clips
    if ($switches =~ /-x/) {
	my $overall = 0;
	foreach (@{$self->{expEntry}->{CC}}) {
	    if (/Set final clip points for gap/) {
		$overall = 1;
		last;
	    }
	}
	my $exQL = ($switches =~ /-xl/) ? $self->{expEntry}->{QL}->[$self->{expNentries}->{QL} - 1 - $overall] : undef;
	$cover->cover(0, $exQL) if defined($exQL);
	my $exQR = ($switches =~ /-xr/) ? $self->{expEntry}->{QR}->[$self->{expNentries}->{QR} - 1 - $overall] : undef;
	$cover->cover($exQR, $SQlen + 1) if defined($exQR);
    }
    return $cover;
}

1;
