    {-| Data type for a sequence location consiting of a contiguous range
    of positions on the sequence.

    Throughout, /sequence position/ refers to a 'Pos.Pos' which includes a
    strand.  An index into a sequence is referred to as an /offset/, and
    is generally of type 'Offset'.

     -}

    module Bio.Location.ContigLocation (
      -- * Sequence locations
      ContigLoc(..), fromStartEnd, fromPosLen

      -- * Locations and positions
      , bounds, startPos, endPos, posInto, posOutof, isWithin, overlaps

      -- * Extracting subsequences
      , seqData, seqDataPadded

      -- * Transforming locations
      , slide, extend

      -- * Displaying locations
      , display
                                       ) where

    import Prelude hiding (length)

    import Control.Monad.Error
    import qualified Data.ByteString.Lazy.Char8 as LBS

    import Bio.Sequence.SeqData

    import qualified Bio.Location.Position as Pos
    import Bio.Location.Strand

    -- | Contiguous sequence location defined by a span of sequence
    -- positions, lying on a specific strand of the sequence.
    data ContigLoc = ContigLoc { offset5 :: !Offset   -- ^ The offset of the 5\' end of the location, as a 0-based index
                               , length :: !Offset    -- ^ The length of the location
                               , strand :: !Strand    -- ^ The strand of the location
                               } deriving (Eq, Ord, Show)

    instance Stranded ContigLoc where
        revCompl (ContigLoc seq5 len str) = ContigLoc seq5 len $ revCompl str

    -- | Create a sequence location lying between 0-based starting and
    -- ending offsets.  When @start < end@, the location
    -- be on the forward strand, otherwise it will be on the
    -- reverse complement strand.

    fromStartEnd :: Offset -> Offset -> ContigLoc
    fromStartEnd start end
        | start < end = ContigLoc start (1 + end - start) Fwd
        | otherwise   = ContigLoc end   (1 + start - end) RevCompl

    -- | Create a sequence location from the sequence position of the
    -- start of the location and the length of the position.  The strand
    -- of the location, and the direction it extends from the starting
    -- position, are determined by the strand of the starting position.
    fromPosLen :: Pos.Pos -> Offset -> ContigLoc
    fromPosLen (Pos.Pos off5 Fwd)      len = ContigLoc off5               len Fwd
    fromPosLen (Pos.Pos off3 RevCompl) len = ContigLoc (off3 - (len - 1)) len RevCompl

    -- | The bounds of a sequence location.  This is a pair consisting of
    -- the lowest and highest sequence offsets covered by the region.  The
    -- bounds ignore the strand of the sequence location, and the first
    -- element of the pair will always be lower than the second.
    bounds :: ContigLoc -> (Offset, Offset)
    bounds (ContigLoc seq5 len _) = (seq5, seq5 + len - 1)

    -- | Sequence position of the start of the location.  This is the 5'
    -- end on the location strand, which will have a higher offset than
    -- 'endPos' if the location is on the 'RevCompl' strand.
    startPos :: ContigLoc -> Pos.Pos
    startPos (ContigLoc seq5 len str)
        = case str of
            Fwd      -> Pos.Pos seq5             str
            RevCompl -> Pos.Pos (seq5 + len - 1) str

    -- | Sequence position of the end of the location, as described in 'startPos'.
    endPos :: ContigLoc -> Pos.Pos
    endPos (ContigLoc seq5 len str)
        = case str of
            Fwd      -> Pos.Pos (seq5 + len - 1) str
            RevCompl -> Pos.Pos seq5             str

    -- | Returns a location resulting from sliding the original location
    -- along the sequence by a specified offset.  A positive offset will
    -- move the location away from the 5\' end of the forward stand of the
    -- sequence regardless of the strand of the location itself.  Thus,
    --
    -- > slide (revCompl cloc) off == revCompl (slide cloc off)
    slide :: Offset -> ContigLoc -> ContigLoc
    slide dpos (ContigLoc seq5 len str) = ContigLoc (seq5 + dpos) len str

    -- | Extract the nucleotide 'SeqData' for the sequence location.  If
    -- any part of the location lies outside the bounds of the sequence,
    -- an error results.
    seqData :: (Error e, MonadError e m) => SeqData -> ContigLoc -> m SeqData
    seqData sequ (ContigLoc seq5 len str)
        | seq5 < 0 = outOfBounds
        | otherwise = case LBS.take len $ LBS.drop seq5 sequ of
                        fwdseq | LBS.length fwdseq == len -> return $ stranded str fwdseq
                               | otherwise -> outOfBounds
        where outOfBounds = throwError $ strMsg $ "contig seq loc " ++ show (seq5, seq5 + len - 1) ++ " out of SeqData bounds"

    -- | As 'seqData', extract the nucleotide subsequence for the
    -- location.  Any positions in the location lying outside the bounds
    -- of the sequence are returned as @N@ rather than producing an error.
    seqDataPadded :: SeqData -> ContigLoc -> SeqData
    seqDataPadded sequ (ContigLoc seq5 len str) = stranded str fwdseq
        where fwdseq
                  | seq5 + len <= 0         = LBS.replicate len 'N'
                  | seq5 >= LBS.length sequ = LBS.replicate len 'N'
                  | seq5 < 0  = LBS.replicate (negate seq5) 'N' `LBS.append` takePadded (len + seq5) sequ
                  | otherwise = takePadded len $ LBS.drop seq5 sequ
              takePadded sublen subsequ
                  | sublen <= LBS.length subsequ = LBS.take sublen subsequ
                  | otherwise = subsequ `LBS.append` LBS.replicate (sublen - LBS.length subsequ) 'N'

    -- | Given a sequence position and a sequence location relative to the
    -- same sequence, compute a new position representing the original
    -- position relative to the subsequence defined by the location.  If
    -- the sequence position lies outside of the sequence location,
    -- @Nothing@ is returned; thus, the offset of the new position will
    -- always be in the range @[0, length cloc - 1]@.
    posInto :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos
    posInto (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)
        | pos < seq5 || pos >= seq5 + len = Nothing
        | otherwise = Just $ case cStrand of
                               Fwd      -> Pos.Pos (pos - seq5)              pStrand
                               RevCompl -> Pos.Pos (seq5 + len  - (pos + 1)) (revCompl pStrand)

    -- | Given a sequence location and a sequence position within that
    -- location, compute a new position representing the original position
    -- relative to the outer sequence.  If the sequence position lies
    -- outside the location, @Nothing@ is returned.
    --
    -- This function inverts 'posInto' when the sequence position lies
    -- within the position is actually within the location.
    posOutof :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos
    posOutof (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)
        | pos < 0 || pos >= len = Nothing
        | otherwise = Just $ case cStrand of
                               Fwd      -> Pos.Pos (pos + seq5)             pStrand
                               RevCompl -> Pos.Pos (seq5 + len - (pos + 1)) (revCompl pStrand)

    -- | Returns a sequence location produced by extending the original
    -- location on each end, based on a pair of (/5\' extension/, /3\'
    -- extension/).  The 5\' extension is applied to the 5\' end of the
    -- location on the location strand; if the location is on the
    -- 'RevCompl' strand, the 5\' end will have a higher offset than the
    -- 3\' end and this offset will increase by the amount of the 5\'
    -- extension.  Similarly, the 3\' extension is applied to the 3\' end
    -- of the location.
    extend :: (Offset, Offset) -> ContigLoc -> ContigLoc
    extend (ext5, ext3) (ContigLoc seq5 len str)
        = case str of
            Fwd -> ContigLoc (seq5 - ext5) (len + ext5 + ext3) str
            RevCompl -> ContigLoc (seq5 - ext3) (len + ext5 + ext3) str

    -- | Returns @True@ when  a sequence position lies within a sequence
    -- location on the same sequence, and occupies the same strand.
    isWithin :: Pos.Pos -> ContigLoc -> Bool
    isWithin (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)
        = (pos >= seq5) && (pos < seq5 + len) && (cStrand == pStrand)

    -- | Returns @True@ when two sequence locations overlap at any
    -- position.
    overlaps :: ContigLoc -> ContigLoc -> Bool
    overlaps contig1 contig2
        = case (bounds contig1, bounds contig2) of
            ((low1, high1),(low2, high2)) -> (strand contig1 == strand contig2)
                                             && (low1 <= high2) && (low2 <= high1)

    -- | Display a human-friendly, zero-based representation of a sequence location.
    display :: ContigLoc -> String
    display cloc = show (Pos.offset $ startPos cloc) ++ "to" ++ show (Pos.offset $ endPos cloc)