# ***** BEGIN LICENSE BLOCK *****
# Version: MPL 2.0
#
# The contents of this file are subject to the Mozilla Public License Version
# 2.0 (the "License"); you may not use this file except in compliance with
# the License. You may obtain a copy of the License at
# http://www.mozilla.org/MPL/2.0/
#
# Software distributed under the License is distributed on an "AS IS" basis,
# WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
# for the specific language governing rights and limitations under the
# License.
#
# The Original Code is OMA standalone.
#
# The Initial Developer of the Original Code is CBRG Research Group; 
# ETH Zurich; Switzerland.
# Portions created by the Initial Developer are Copyright (C) 2005-2013
# the Initial Developer. All Rights Reserved.
#
# Contributor(s):
#   Christophe Dessimoz <cdessimoz@inf.ethz.ch>
#   Adrian Altenhoff <adrian.altenhoff@inf.ethz.ch>
#   Stefan Zoller <stefan.zoller@inf.ethz.ch>
#   Adrian Schneider <adrian.schneider@inf.ethz.ch>
#   Alexander Roth <alexander.roth@inf.ethz.ch>
#   Gaston Gonnet <gonnet@inf.ethz.ch>
#
# ***** END LICENSE BLOCK *****
IdentifyPosition := proc(fAugustus:string,fContigLengths:string)
    # Read and store contig lengths into a table
    cLengths := table();
    OpenReading(fContigLengths);
    t1 := ReadRawLine();
    while t1 <> 'EOF' do
        t2 := SearchDelim('\t',t1);
        cLengths[t2[1]] :=parse(t2[2,1..-2]);
        t1 := ReadRawLine();
    od; 

    # Read and store gene positions into a table
    gRanges := table();
    OpenReading(fAugustus);
    t1 := ReadRawLine();
    while t1 <> 'EOF' do
        if t1[1] <> '#' then
            t2 := SearchDelim('\t',t1);
            if t2[3] ='gene' then 
                # positive strand
                if t2[7] = '+' then
                    gRanges[t2[9,1..-2]] := 
                        [t2[1],parse(t2[4])..parse(t2[5])];
                else
                    gRanges[t2[9,1..-2]] := 
                        [t2[1],parse(t2[5])..parse(t2[4])];
                fi;
            fi;
        fi;
        t1 := ReadRawLine();
    od; 
    
    return(gRanges,cLengths);
end:


DistBounds := proc(g,gRanges,cLengths)
    r := gRanges[g,2];
    l := cLengths[gRanges[g,1]];
    # positive strand
    if r[2] > r[1] then
        dbegin := r[1]-1;
        dend   := l-r[2];
    else
        dbegin := l-r[2];
        dend   := r[1]-1;
    fi;
    return([min(dbegin,dend),dbegin,dend]);
end: