#!/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.5 2007/02/07 15:23:10 rmd Exp $
# $Log: Exptools.pm,v $
# Revision 1.5  2007/02/07  15:23:10  rmd
# Changed time format in add_notes so that it is more resistant to different
# locales.
#
# Revision 1.4  2006/10/02  15:29:58  rmd
# Changes to better handle the case of a read being in the trace server but not
# with the trace type and seq type requested.  In this case we will now the next
# method in search_order, if available, and if that doesn't work readExpFile
# will return 1.
#
# Revision 1.3  2006/08/29  09:08:01  jkb
# Major overhaul to support Internal Trace Server
#
# Revision 1.1  2006/03/13 12:29:40  jkb
# Imported my local hack of Rob's Exptools.pm package. Ultimately this
# needs merging back with badger-bin, but it's currently not backwards
# compatible.
#
# The changes involve white-space handling. All entries held in memory
# no longer have newline on the end of the records. Multi-line entries
# no longer store the leading additional white space either as this gets
# added back when the multi-line record is written.
#
# Also it now passes perl -w in "use strict" mode.
#
# 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;
use POSIX qw(strftime);

# TraceServer support only valid with Perl 5.8.0 onwards
if ($] >= 5.008) {
    require TraceServer;
}

# -----------------------------------------------------------------------------
# Main interface
sub new {
    my ($type, $opts) = @_;

    my $self = {expOrder     => [],
		expNentries  => {},
		expEntry     => {},
		iterators    => {},
		search_order => ["ts", "repos"],
	        sections     => [],
	        warn         => 1,
	        seq_type     => undef,
	        trace_type   => undef,
	        ts           => undef};

    foreach (keys %$opts) {
	if (exists($self->{$_})) {
	    # Quick and hacky config support for now.
	    if (ref($opts->{$_}) eq "ARRAY") {
		$self->{$_} = [@{$opts->{$_}}]
	    } else {
		$self->{$_} = $opts->{$_};
	    }
	}
    }

    return bless $self, $type;
}

#
# $sections is an array reference; if present it contains a list of which
# segments of an experiment file to populate, otherwise the entire thing
# will be loaded. Valid segment strings are grouped by experiment file line
# identifiers as follows:
#
# Section name     Exp line IDs
# -------------------------------------------------
# sequence	   SQ
# confidence	   AV
# date		   DT
# trace		   MT, LN, LT
# chemistry	   PR, DR, CH
# source	   LG, SV, SF, CN, CV, CF, TN, ST
# clips		   SL, SR, CL, CR, CS
# notes		   NT
# tags		   TG
# status 	   PS
#
sub readExpFile {
    my ($self, $filename, $options) = @_;

    confess "Wrong type" unless ref $self;

    # Populate the options array with defaults from $self too.
    $options = {} unless defined($options);
    foreach (("search_order",
	      "warn",
	      "seq_type",
	      "trace_type",
	      "sections",
	      "ts")) {
	if (!exists($options->{$_}) && exists($self->{$_})) {
	    $options->{$_} = $self->{$_};
	}
    }
    $self->{options} = $options;

    # Coerce the sections array into a hash for quick lookup.
    if (defined($options->{sections})) {
	foreach (@{$options->{sections}}) {
	    $self->{sec}{$_} = 1;
	}
    }

    # Loop through the search_order trying to read the experiment file.
    foreach my $order (@{$options->{search_order}}) {
	if ($order eq "ts") {
	    if ($] >= 5.008 && $self->read_in_trace_server($options,
							   $filename)) {
		my $r = $self->readExpFileFromTS($options->{ts}, $filename);
		if (0 == $r) { return 0; }
	    }
	} elsif ($order eq "repos") {
	    if (open(EXP, "<$filename")) {
		my $r = $self->readExpFileFromRepos(\*EXP);
		close(EXP);
		if (0 == $r) { return 0; }
	    }
	}
    }

    warn "readExpFile: cannot get $filename\n" if $options->{warn};
    return 1;
}

sub writeExpFile {
    my ($self, $file) = @_;
    confess "Wrong type" unless ref $self;

    if ($self->{from_ts}) {
	print "Unimplemented yet\n";
    } else {
	open(EXP, ">$file") || die "Could not write $file";
	$self->writeExpFileToRepos(\*EXP);
	close(EXP);
    }
}

# -----------------------------------------------------------------------------
# Code to read/write from the repository, ie real Experiment files.
sub readExpFileFromRepos {
    my ($self, $file) = @_;
    $self->{expOrder} = [];
    $self->{expNentries} = {};
    $self->{expEntry} = {};
    $self->{from_ts} = 0;
    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"
	    if $self->{options}{warn};
	return 1;
    }

#    foreach my $id (keys %{$self->{expEntry}}) {
#	foreach (@{$self->{expEntry}->{$id}}) {
#	    print "$id='$_'\n";
#	}
#    }
#
#    $self->writeExpFile(\*STDOUT);

    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 = defined($self->{expNentries}->{$id})
		? $self->{expNentries}->{$id} : 0;
	    $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);
	    chomp($value);
	    my $n = $self->{expNentries}->{$id};
	    $n = defined($self->{expNentries}->{$id})
		? $self->{expNentries}->{$id} : 0;
	    if ($id eq $multiline && $value =~ s/^     //) {
		$self->{expEntry}->{$id}->[$n-1] .= "\n$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" || $id eq "NT") {
		    $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 writeExpFileToRepos {
    my ($self, $file) = @_;
    confess "Wrong type" unless ref $self;
    my $data = $self->encodeExpFile();
    print $file $data;
}


# -----------------------------------------------------------------------------
# Code to read/write from the internal trace server. This fakes up the
# Experiment file structures. We cache a single global TraceServer connection
# for the duration of this module's lifetime.

my $global_ts;
sub read_in_trace_server {
    my ($self, $options, $filename) = @_;

    if (!defined($options->{ts})) {
	if (!defined($global_ts)) {
	    eval {$global_ts = TraceServer->new()};
	}
	$options->{ts} = $global_ts;
	return 0 if ($@);
    }

    $filename =~ s:.*/::;
    return $options->{ts}->read_exists($filename);
}

sub readExpFileFromTS {
    my ($self, $ts, $filename) = @_;
    $self->{expOrder} = [];
    $self->{expNentries} = {};
    $self->{expEntry} = {};
    $self->{from_ts} = 1;

    $filename =~ s:.*/::;
    return $self->extract_sequence($ts, $filename);
}

# Extract a single sequence
sub extract_sequence {
    my ($exp, $ts, $name) = @_;

    # Fetch the read, sequence, confidence and types.
    my ($read, $index) = $ts->get_read_by_name($name);
    my ($trace_type, $seq_type) =
	(defined($exp->{options}{trace_type}) &&
	 defined($exp->{options}{seq_type}))
	 ? ($exp->{options}{trace_type}, $exp->{options}{seq_type})
	 : $read->get_seq_type_by_index();

    my ($seq, $conf);
    if (!exists($exp->{sec}) || $exp->{sec}{sequence}
	                     || $exp->{sec}{clips}
	                     || $exp->{sec}{notes}
	                     || $exp->{sec}{tags}
	                     || $exp->{sec}{status}) {
	$seq = $read->get_sequence_by_type($trace_type, $seq_type);
	unless ($seq) {
	    return 1;
	}
    }
    if (!exists($exp->{sec}) || $exp->{sec}{confidence}) {
	$conf = $read->get_confidence_by_type($trace_type, $seq_type);
    }

    # Cache read object for any subsequent manipulation
    $exp->{read} = $read;

    # Create a new Experiment file and copy basic read info
    add_exp_details($exp, $read);
    if (!exists($exp->{sec}) || $exp->{sec}{source}) {
	add_clone_info($exp, $read);
    }
    add_seq_conf($exp, $read, $seq, $conf);
    if (!exists($exp->{sec}) || $exp->{sec}{clips}) {
	add_clips($exp, $seq);
    }
    if (!exists($exp->{sec}) || $exp->{sec}{notes}) {
	add_notes($exp, $read);
    }
    if (!exists($exp->{sec}) || $exp->{sec}{tags}) {
	add_tags($exp, $seq);
    }
    if (!exists($exp->{sec}) || $exp->{sec}{status}) {
	add_status($exp, $read);
    }

    # $exp->addEntry("OP", `whoami`);
    
    # No basecaller info yet?
    # $exp->addEntry("BC", "phred");

    return 0;
}

# Adds EXP field info from a Read.
# Basic stuff like ID, EN, MN, DT, MT, LN, LT, PR, DR, CH
sub add_exp_details {
    my ($exp, $read) = @_;

    # Identities
    $exp->addEntry("ID", $read->get_name());
    $exp->addEntry("EN", $read->get_name());
    $exp->addEntry("MN", $read->get_name());
    if (!exists($exp->{sec}) || $exp->{sec}{date}) {
	my $date = $read->get_run_datetime();
	if (defined($date)) {
	    my $dtstr = strftime("%Y-%m-%d", localtime($date));
	    $dtstr =~ s/-0/-/g;
	    $exp->addEntry("DT", $dtstr);
	}
    }
    if (!exists($exp->{sec}) || $exp->{sec}{trace}) {
	$exp->addEntry("MT", "SCF");
	$exp->addEntry("LN", $read->get_name() . "SCF");
	$exp->addEntry("LT", "SCF");
    }

    # Primer type/direction/chemistry
    if (!exists($exp->{sec}) || $exp->{sec}{chemistry}) {
	# Chemistry code
	my $ch = $read->get_chemistry();
	if (defined($ch)) {
	    if ($ch =~ /^TRM/) {
		$exp->addEntry("CH", "1");
	    } elsif ($ch =~ /^PRI/) {
		$exp->addEntry("CH", "0");
	    }
	}

	# Primer + direction
	my $pr = 1;
	$pr += 2 if ($read->get_universal_primer() == 0);
	if ($read->get_direction() &&
	    $read->get_direction() == TraceServer::TSR_REVERSE()) {
	    $pr++;
	    $exp->addEntry("PR", $pr);
	    $exp->addEntry("DR", "-");
	} else {
	    $exp->addEntry("PR", $pr);
	    $exp->addEntry("DR", "+");
	}
    }
}

# Adds clone information to an EXP file.
sub add_clone_info {
    my ($exp, $read) = @_;

    my ($lg, $tn);
    my $src = $read->get_dnasource();

    if (defined($src)) {
	if ($src->get_type() eq "clone") {
	    my $size   = $src->get_library_size();
	    my $stddev = $src->get_library_stddev();
	    my $vec    = $src->get_library_vector();
	    $tn = $src->get_name();
	    $tn =~ s/.*:://;
	    if ($size) {
		$exp->addEntry("SI", ($size-$stddev)."..".($size+$stddev));
	    }
	    $lg = $src->get_library_name();
	    $lg =~ s/.*:://;
	    $exp->addEntry("LG", $lg);
	    $exp->addEntry("SV", $vec);
	    $exp->addEntry("SF", "BADGER/lib/seqs/staden/".$vec.".seq");
	} else {
	    $tn = $src->get_name();
	}

	# top-most dnasource of type clone is our "Clone" line.
	my $clone = undef;
	for ($src = $src->get_parent(); $src; $src = $src->get_parent()) {
	    if ($src->get_type() eq "clone") {
		$clone = $src;

		# Also set ligation if it was unset before due to the immediate
		# source being PCR (for example).
		if (!defined($lg)) {
		    $lg = $src->get_library_name();
		    $lg =~ s/.*:://;
		    $exp->addEntry("LG", $lg);
		}
	    }
	}

	if (defined($clone) && $clone->get_type() eq "clone") {
	    my $cn = $clone->get_name();
	    $cn =~ s/.*CLONE:://;
	    $exp->addEntry("CN", $cn);
	    my $vec    = $clone->get_library_vector();
	    if (defined($vec)) {
		$exp->addEntry("CV", $vec);
		$exp->addEntry("CF", "BADGER/lib/seqs/staden/".$vec.".seq");
	    }
	}

	$exp->addEntry("TN", $tn);
	$exp->addEntry("ST", "2");
    }

    # Vector primer/cut sites are currently missing.
    # 
    # $exp->addEntry("SP", "-46");
    # $exp->addEntry("SC", "271");
}

# Adds sequence and confidence info
sub add_seq_conf {
    my ($exp, $read, $seq, $conf) = @_;

    # Sequence
    if ((!exists($exp->{sec}) || $exp->{sec}{sequence}) && defined($seq)) {
	my $dna = $seq->get_dna();
	$dna =~ s/N/-/g;
	$exp->addEntry("SQ", $dna);
    }

    # Confidence
    if ((!exists($exp->{sec}) || $exp->{sec}{confidence}) && defined($conf)) {
	my $pcon = $conf->get_phred();
	my $len = length($pcon);
	my $j;
	my $str = "";
	for ($j = 0; $j < $len; $j += 20) {
	    my @tmp = unpack("C*", substr($pcon, $j, 20));
	    $str .= "@tmp\n";
	}
	$exp->addEntry("AV", "$str");
    }
}

# Adds clipping info; SL/SR, CL/CR/CS, QL/QL
sub add_clips {
    my ($exp, $seq) = @_;

    return unless (defined($seq));

    my $nclips = $seq->get_num_clips();
    my $len = length($seq->get_dna());
    for (my $i = 0; $i < $nclips; $i++) {
	my $clip = $seq->get_clip($i);
	if ($clip->get_type() eq "CVEC") {
	    if ($clip->get_end() == 0) {
		$exp->addEntry("CS", "1.." . ($clip->get_start()-1));
	    } elsif ($clip->get_start() == 0) {
		$exp->addEntry("CS", ($clip->get_end()+1) . "..$len");
	    } else {
		warn "Unable to compute CS line";
	    }
	} else {
	    my $type = "Q";
	    $type = "S" if ($clip->get_type() eq "SVEC");
	    if ($clip->get_start() != 0) {
		$exp->addEntry($type . "L", $clip->get_start()-1);
	    }
	    if ($clip->get_end() != 0) {
		$exp->addEntry($type . "R", $clip->get_end()+1);
	    }
	}
    }
}

# Adds notes (NT).
sub add_notes {
    my ($exp, $read) = @_;

    my $nnotes = $read->get_num_notes();
    for (my $i = 0; $i < $nnotes; $i++) {
	my $note = $read->get_note($i);
	my $time = $note->get_time();
	my $tstr = strftime("%a %b %e %H:%M:%S %Y %Z",localtime($time));
	my $data = $note->get_type() . " "
	    . "ctime=$tstr ($time)\n"
	    . "mtime=$tstr ($time)\n"
	    . "from=reading " . $read->get_name() . "\n"
	    . "comment=" . $note->get_text();
	$exp->addEntry("NT", $data);
    }
}

# Adds tags (TG).
sub add_tags {
    my ($exp, $seq) = @_;

    return unless (defined($seq));

    my $ntags = $seq->get_num_tags();
    for (my $i = 0; $i < $ntags; $i++) {
	my $tag = $seq->get_tag($i);
	my $data = $tag->get_type() . " "
	    . $tag->get_direction() . " "
	    . $tag->get_start() . ".." . $tag->get_end();
	my $text = $tag->get_text();
	if ($text ne "") {
	    $data .= ($text =~ /^\n/s) ? $text : " $text";
	}
	$exp->addEntry("TG", $data);
    }
}

# Add processing status (PS).
sub add_status {
    my ($exp, $read) = @_;
    my $PS = $read->get_status_description();
    if (defined($PS)) {
	if ($PS eq "PASS") {
	    undef($PS);
	} else {
	    $PS =~ s/^FAIL //;
	}
    }
    $exp->addEntry("PS", $PS) if defined($PS);
}

# -----------------------------------------------------------------------------
# Experiment file object manipulation methods
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;
    my $n = $self->{expNentries}->{$id};
    $n = 0 unless defined($n);
    $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;
    my $n = $self->{expNentries}->{$id};
    if ($n) {
	$self->{expEntry}->{$id}->[$n-1] = $value;
    } else {
	$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->{expNentries}->{$id};
    @{$self->{expOrder}} = grep { $_ ne $id; } @{$self->{expOrder}};
}

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;
}


# Serialises the Exp object into an experiment file as a single string.
sub encodeExpFile {
    my ($self) = @_;
    my $id = $self->getEntry('ID');

    my $buffer = "";

    $buffer .= "ID   $id\n" if $id;

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

    return $buffer;
}

1;


__END__

=head1 NAME

Exptools - Perl module for reading Experiment files.

=head1 SYNOPSIS

 use Exptools;
 $exp = Exptools->new();
 $exp = Exptools->new($option_hash);

 $err = $exp->readExpFile($filename);
 $err = $exp->readExpFile($filename, $option_hash);
 $exp->writeExpFile($filename);

 $item = $exp->getFirstEntry($line_type);
 $item = $exp->getNextEntry($line_type);
 $item = $exp->getEntry($line_type);
 $exp->addEntry($line_type, $value);
 $exp->replaceEntry($line_type, $value);
 $exp->killEntry($line_type);

 ($left, $right) = $exp->determineActiveRegion();
 $text = $exp->encodeExpFile();

=head1 DESCRIPTION

The Exptools module provides a unified interface to accessing
Experiment File contents from either the local repository disk or the
internal trace server.

When creating a new Exptools instance an optional hash reference may
be passed in containing parameters to control the behaviour of the
module. If the hash reference is passed over not all elements need to
be defined. The defaults are equivalent to:

 my $exp = Exptools->new({search_order => ["ts", "repos"],
                          warn         => 1,
                          trace_type   => undef,
			  seq_type     => undef,
			  sections     => []});

The B<search_order> is an array reference containing one or more
strings matching "ts" for Internal Trace Server and "repos" for the
local Repository disks. It's purpose is to act as a resolver
indicating where to search for data and in what order.

B<warn> controls whether to emit warnings or act in silence.

B<trace_type> and B<seq_type> only have relevance when accessing the
Internal Trace Server. They control which trace and sequence versions
are fetched out of the server. (On disk we only have one experiment
file so such control has traditionally been achieved through using
multiple directories with the same sequence names in each.)

B<sections> is an array reference containing a symbolic list of
sections to extract when accessing the Internal Trace Server. It is
ignored for direct repository file reading. The purpose is to speed up
database accessing as if the intention of the caller is to do little
more than convert an experiment file into a fasta file (for example)
then the database queries to populate ligation, clone name, date,
etc. are all wasted time. Valid symbols in sections and their
associated experiment file line types are:

=over

=item sequence

SQ

=item confidence

AV

=item date

DT

=item trace

MT, LN, LT

=item chemistry

PR, DR, CH

=item source

LG, SV, SF, CN, CV, CF, TN, ST

=item clips

SL, SR, CL, CR, CS

=item notes

NT

=item tags

TG

=item status

PS

=back

To read an experiment file use the B<readExpFile> function with the
filename as the first argument. As an optional second argument the
options hash reference may be passed in again. This will act to modify
the default options for this readExpFile call only.

For example:

  my $exp = Exptools->new({search_order => ["ts", "repos"]});
  $exp->readExpFile($filename, {search_order => ["repos"],
			        sections     => ["sequence"]});

=head1 TODO

Add documentation for the other functions.

Add support (or document lack of) for writing back to the Trace Server.

=cut
