{-# LANGUAGE InstanceSigs  #-}
{-# LANGUAGE TupleSections #-}
module Bio.Chain.Alignment.Algorithms where

import           Control.Lens             (Index, IxValue)
import qualified Data.Array.Unboxed       as A (bounds, range)
import           Data.List                (maximumBy)

import           Bio.Chain                hiding ((!))
import           Bio.Chain.Alignment.Type
import           Control.Monad            (forM_)
import           Data.Ord                 (comparing)

import           Control.Monad.ST         (ST)
import           Data.Array.Base          (readArray, writeArray)
import           Data.Array.ST            (MArray (..), STUArray, newArray,
                                           runSTUArray)
import           Data.Array.Unboxed       (Ix (..), UArray, (!))


-- Alignnment methods

newtype EditDistance e1 e2       = EditDistance        (e1 -> e2 -> Bool)
data GlobalAlignment a e1 e2     = GlobalAlignment     (Scoring e1 e2) a
data LocalAlignment a e1 e2      = LocalAlignment      (Scoring e1 e2) a
data SemiglobalAlignment a e1 e2 = SemiglobalAlignment (Scoring e1 e2) a

-- Common functions

-- | Lift simple substitution function to a ChainLike collection
--
{-# SPECIALISE substitute :: (Char -> Char -> Int) -> Chain Int Char -> Chain Int Char -> Int -> Int -> Int #-}
{-# INLINE substitute #-}
substitute :: (Alignable m, Alignable m') => (IxValue m -> IxValue m' -> Int) -> m -> m' -> Index m -> Index m' -> Int
substitute :: (IxValue m -> IxValue m' -> Int)
-> m -> m' -> Index m -> Index m' -> Int
substitute IxValue m -> IxValue m' -> Int
f m
s m'
t Index m
i Index m'
j = IxValue m -> IxValue m' -> Int
f (m
s m -> Index m -> IxValue m
forall m. ChainLike m => m -> Index m -> IxValue m
`unsafeRead` (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i)) (m'
t m' -> Index m' -> IxValue m'
forall m. ChainLike m => m -> Index m -> IxValue m
`unsafeRead` (Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j))

-- | Simple substitution function for edit distance
--
substituteED :: EditDistance e1 e2 -> (e1 -> e2 -> Int)
substituteED :: EditDistance e1 e2 -> e1 -> e2 -> Int
substituteED (EditDistance e1 -> e2 -> Bool
genericEq) e1
x e2
y = if e1
x e1 -> e2 -> Bool
`genericEq` e2
y then Int
1 else Int
0

-- | Default traceback stop condition.
--
{-# SPECIALISE defStop :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> Int -> Int -> Bool #-}
{-# INLINE defStop #-}
defStop :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool
defStop :: Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool
defStop Matrix m m'
_ m
s m'
t Index m
i Index m'
j = let (Index m
lowerS, Index m
_) = m -> (Index m, Index m)
forall m. ChainLike m => m -> (Index m, Index m)
bounds m
s
                        (Index m'
lowerT, Index m'
_) = m' -> (Index m', Index m')
forall m. ChainLike m => m -> (Index m, Index m)
bounds m'
t
                    in  Index m
i Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS Bool -> Bool -> Bool
&& Index m'
j Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT

-- | Traceback stop condition for the local alignment.
--
{-# SPECIALISE localStop :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> Int -> Int -> Bool #-}
{-# INLINE localStop #-}
localStop :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool
localStop :: Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool
localStop Matrix m m'
m' m
s m'
t Index m
i Index m'
j = let (Index m
lowerS, Index m
_) = m -> (Index m, Index m)
forall m. ChainLike m => m -> (Index m, Index m)
bounds m
s
                           (Index m'
lowerT, Index m'
_) = m' -> (Index m', Index m')
forall m. ChainLike m => m -> (Index m, Index m)
bounds m'
t
                       in  Index m
i Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS Bool -> Bool -> Bool
|| Index m'
j Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT Bool -> Bool -> Bool
|| Matrix m m'
m' Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Match) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0

{-# INLINE move #-}
move
  :: (Alignable m, Alignable m', IsGap g)
  => g
  -> Move m m'
move :: g -> Move m m'
move g
g
  | g -> Bool
forall a. IsGap a => a -> Bool
isAffine g
g = g -> Move m m'
forall m m' g.
(Alignable m, Alignable m', IsGap g) =>
g -> Move m m'
moveAffine g
g
  | Bool
otherwise = g -> Move m m'
forall m m' g.
(Alignable m, Alignable m', IsGap g) =>
g -> Move m m'
moveSimple g
g

{-# INLINE moveSimple #-}
moveSimple
  :: (Alignable m, Alignable m', IsGap g)
  => g
  -> Move m m'
moveSimple :: g -> Move m m'
moveSimple g
g Matrix m m'
m m
s m'
t Index m
i Index m'
j EditOp
_
  | Index m
i Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS = (EditOp
Match, Index m
lowerS, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, Index m' -> Operation (Index m) (Index m')
forall i j. j -> Operation i j
INSERT (Index m' -> Operation (Index m) (Index m'))
-> Index m' -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j)
  | Index m'
j Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT = (EditOp
Match, Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m'
lowerT, Index m -> Operation (Index m) (Index m')
forall i j. i -> Operation i j
DELETE (Index m -> Operation (Index m) (Index m'))
-> Index m -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i)
  | Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m'
j, EditOp
Match) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Match) = (EditOp
Match, Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m'
j, Index m -> Operation (Index m) (Index m')
forall i j. i -> Operation i j
DELETE (Index m -> Operation (Index m) (Index m'))
-> Index m -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i)
  | Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, EditOp
Match) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Match) = (EditOp
Match, Index m
i, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, Index m' -> Operation (Index m) (Index m')
forall i j. j -> Operation i j
INSERT (Index m' -> Operation (Index m) (Index m'))
-> Index m' -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j)
  | Bool
otherwise = (EditOp
Match, Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, Index m -> Index m' -> Operation (Index m) (Index m')
forall i j. i -> j -> Operation i j
MATCH (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i) (Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j))
    where
      (Index m'
lowerT, Index m'
_) = m' -> (Index m', Index m')
forall m. ChainLike m => m -> (Index m, Index m)
bounds m'
t
      (Index m
lowerS, Index m
_) = m -> (Index m, Index m)
forall m. ChainLike m => m -> (Index m, Index m)
bounds m
s

{-# INLINE moveAffine #-}
-- | Move function for affine alignment traceback. Implements a "Manhattan grid".
--
-- See here: <http://www.csbio.unc.edu/mcmillan/Comp555S16/Lecture14.html>
-- or file @doc/Affine_Alignment.pdf@ in the repository.
moveAffine
  :: (Alignable m, Alignable m', IsGap g)
  => g
  -> Move m m'
moveAffine :: g -> Move m m'
moveAffine g
g Matrix m m'
m m
s m'
t Index m
i Index m'
j EditOp
prevOp
  | Index m
i Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS = (EditOp
Insert, Index m
lowerS, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, Index m' -> Operation (Index m) (Index m')
forall i j. j -> Operation i j
INSERT (Index m' -> Operation (Index m) (Index m'))
-> Index m' -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j)
  | Index m'
j Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT = (EditOp
Delete, Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m'
lowerT, Index m -> Operation (Index m) (Index m')
forall i j. i -> Operation i j
DELETE (Index m -> Operation (Index m) (Index m'))
-> Index m -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i)
  | Bool
otherwise =
    case EditOp
prevOp of
      EditOp
Delete
        | Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Delete) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m'
j, EditOp
Delete) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
deleteCostExtend g
g -> (EditOp
Delete, Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m'
j, Index m -> Operation (Index m) (Index m')
forall i j. i -> Operation i j
DELETE (Index m -> Operation (Index m) (Index m'))
-> Index m -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i)
        | Bool
otherwise -> (EditOp
Match, Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m'
j, Index m -> Operation (Index m) (Index m')
forall i j. i -> Operation i j
DELETE (Index m -> Operation (Index m) (Index m'))
-> Index m -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i)
      EditOp
Insert
        | Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Insert) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, EditOp
Insert) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
insertCostExtend g
g -> (EditOp
Insert, Index m
i, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, Index m' -> Operation (Index m) (Index m')
forall i j. j -> Operation i j
INSERT (Index m' -> Operation (Index m) (Index m'))
-> Index m' -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j)
        | Bool
otherwise -> (EditOp
Match, Index m
i, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, Index m' -> Operation (Index m) (Index m')
forall i j. j -> Operation i j
INSERT (Index m' -> Operation (Index m) (Index m'))
-> Index m' -> Operation (Index m) (Index m')
forall a b. (a -> b) -> a -> b
$ Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j)
      EditOp
Match
        | Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Match) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Delete) -> g -> Move m m'
forall m m' g.
(Alignable m, Alignable m', IsGap g) =>
g -> Move m m'
move g
g Matrix m m'
m m
s m'
t Index m
i Index m'
j EditOp
Delete
        | Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Match) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Insert) -> g -> Move m m'
forall m m' g.
(Alignable m, Alignable m', IsGap g) =>
g -> Move m m'
move g
g Matrix m m'
m m
s m'
t Index m
i Index m'
j EditOp
Insert
        | Bool
otherwise -> (EditOp
Match, Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, Index m -> Index m' -> Operation (Index m) (Index m')
forall i j. i -> j -> Operation i j
MATCH (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i) (Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j))
    where
      (Index m'
lowerT, Index m'
_) = m' -> (Index m', Index m')
forall m. ChainLike m => m -> (Index m, Index m)
bounds m'
t
      (Index m
lowerS, Index m
_) = m -> (Index m, Index m)
forall m. ChainLike m => m -> (Index m, Index m)
bounds m
s

-- | Default condition of moving diagonally in traceback.
--
{-# SPECIALISE defDiag :: (Char -> Char -> Int) -> Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> Int -> Int -> Bool #-}
{-# INLINE defDiag #-}
defDiag :: (Alignable m, Alignable m') => (IxValue m -> IxValue m' -> Int) -> Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool
defDiag :: (IxValue m -> IxValue m' -> Int)
-> Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool
defDiag IxValue m -> IxValue m' -> Int
sub' Matrix m m'
m m
s m'
t Index m
i Index m'
j = let sub :: Index m -> Index m' -> Int
sub = (IxValue m -> IxValue m' -> Int)
-> m -> m' -> Index m -> Index m' -> Int
forall m m'.
(Alignable m, Alignable m') =>
(IxValue m -> IxValue m' -> Int)
-> m -> m' -> Index m -> Index m' -> Int
substitute IxValue m -> IxValue m' -> Int
sub' m
s m'
t
                         in  Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
i, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
j, EditOp
Match) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Index m -> Index m' -> Int
sub Index m
i Index m'
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! (Index m
i, Index m'
j, EditOp
Match)

-- | Default start condition for traceback.
--
{-# SPECIALISE defStart :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> (Int, Int) #-}
{-# INLINE defStart #-}
defStart :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> (Index m, Index m')
defStart :: Matrix m m' -> m -> m' -> (Index m, Index m')
defStart Matrix m m'
m m
_ m'
_ = let ((Index m
_, Index m'
_, EditOp
_), (Index m
upperS, Index m'
upperT, EditOp
_)) = Matrix m m'
-> ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> (i, i)
A.bounds Matrix m m'
m in (Index m
upperS, Index m'
upperT)

-- | Default start condition for traceback in local alignment.
--
{-# SPECIALISE localStart :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> (Int, Int) #-}
{-# INLINE localStart #-}
localStart :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> (Index m, Index m')
localStart :: Matrix m m' -> m -> m' -> (Index m, Index m')
localStart Matrix m m'
m m
_ m'
_ = let ((Index m
lowerS, Index m'
lowerT, EditOp
_), (Index m
upperS, Index m'
upperT, EditOp
_)) = Matrix m m'
-> ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> (i, i)
A.bounds Matrix m m'
m
                       range' :: [(Index m, Index m', EditOp)]
range' = ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
-> [(Index m, Index m', EditOp)]
forall a. Ix a => (a, a) -> [a]
A.range ((Index m
lowerS, Index m'
lowerT, EditOp
Match), (Index m
upperS, Index m'
upperT, EditOp
Match))
                   in  (\(Index m
a, Index m'
b, EditOp
_) -> (Index m
a, Index m'
b)) ((Index m, Index m', EditOp) -> (Index m, Index m'))
-> (Index m, Index m', EditOp) -> (Index m, Index m')
forall a b. (a -> b) -> a -> b
$ ((Index m, Index m', EditOp)
 -> (Index m, Index m', EditOp) -> Ordering)
-> [(Index m, Index m', EditOp)] -> (Index m, Index m', EditOp)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (((Index m, Index m', EditOp) -> Int)
-> (Index m, Index m', EditOp)
-> (Index m, Index m', EditOp)
-> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
!)) [(Index m, Index m', EditOp)]
range'

-- | Default start condition for traceback in semiglobal alignment.
--
{-# SPECIALISE semiStart :: Matrix (Chain Int Char) (Chain Int Char) -> Chain Int Char -> Chain Int Char -> (Int, Int) #-}
{-# INLINE semiStart #-}
semiStart :: (Alignable m, Alignable m') => Matrix m m' -> m -> m' -> (Index m, Index m')
semiStart :: Matrix m m' -> m -> m' -> (Index m, Index m')
semiStart Matrix m m'
m m
_ m'
_ = let ((Index m
lowerS, Index m'
lowerT, EditOp
_), (Index m
upperS, Index m'
upperT, EditOp
_)) = Matrix m m'
-> ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> (i, i)
A.bounds Matrix m m'
m
                      lastCol :: [(Index m, Index m', EditOp)]
lastCol = (, Index m'
upperT, EditOp
Match) (Index m -> (Index m, Index m', EditOp))
-> [Index m] -> [(Index m, Index m', EditOp)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Index m
lowerS .. Index m
upperS]
                      lastRow :: [(Index m, Index m', EditOp)]
lastRow = (Index m
upperS, , EditOp
Match) (Index m' -> (Index m, Index m', EditOp))
-> [Index m'] -> [(Index m, Index m', EditOp)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Index m'
lowerT .. Index m'
upperT]
                  in  (\(Index m
a, Index m'
b, EditOp
_) -> (Index m
a, Index m'
b)) ((Index m, Index m', EditOp) -> (Index m, Index m'))
-> (Index m, Index m', EditOp) -> (Index m, Index m')
forall a b. (a -> b) -> a -> b
$ ((Index m, Index m', EditOp)
 -> (Index m, Index m', EditOp) -> Ordering)
-> [(Index m, Index m', EditOp)] -> (Index m, Index m', EditOp)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (((Index m, Index m', EditOp) -> Int)
-> (Index m, Index m', EditOp)
-> (Index m, Index m', EditOp)
-> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Matrix m m'
m Matrix m m' -> (Index m, Index m', EditOp) -> Int
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
!)) ([(Index m, Index m', EditOp)] -> (Index m, Index m', EditOp))
-> [(Index m, Index m', EditOp)] -> (Index m, Index m', EditOp)
forall a b. (a -> b) -> a -> b
$ [(Index m, Index m', EditOp)]
lastCol [(Index m, Index m', EditOp)]
-> [(Index m, Index m', EditOp)] -> [(Index m, Index m', EditOp)]
forall a. [a] -> [a] -> [a]
++ [(Index m, Index m', EditOp)]
lastRow

-- Alignment algorithm instances

-------------------
  --
  --                        About affine gaps
  --
  -- There are three matrices used in all the algorithms below:
  -- 1) One stores the resulting scores for each prefix pair;
  -- 2) One stores insertion costs in the first sequence for each prefix pair;
  -- 3) One stores insertion costs in the second sequence for each prefix pair.
  --
  -- Matrices 2 and 3 are used in affine penalty calculation:
  -- Let `gapOpen` and `gapExtend` be affine gap penalties (for opening and extending a gap correspondingly).
  -- M2[i, j] = gapOpen   if neither of sequences has insertion or deletion at prefix (i, j);
  -- M2[i, j] = gapExtend if sequence1 has insertion or, equivalently, sequence2 has deletion at prefix (i, j);
  --
  -- gap penalty in the first sequence for the prefix (i, j) is gapOpen + M2[i, j]
  -- So, if there are no gaps in the sequence before, the penalty will be `gapOpen`.
  -- Otherwise it will be `gapExtend`.
  --
  -- So, M2 holds insertion costs for the first sequence, and M3 holds insertion costs for the second sequence.
  --
  -- The resulting score is the same as in plain gap penalty:
  -- the biggest one between substitution, insertion and deletion scores.
  --
-------------------

instance IsGap g => SequenceAlignment (GlobalAlignment g) where

    -- Conditions of traceback are described below
    --
    {-# INLINE cond #-}
    cond :: GlobalAlignment g (IxValue m) (IxValue m') -> Conditions m m'
cond (GlobalAlignment Scoring (IxValue m) (IxValue m')
_ g
gap) = Stop m m' -> Move m m' -> Conditions m m'
forall m m'. Stop m m' -> Move m m' -> Conditions m m'
Conditions Stop m m'
forall m m'.
(Alignable m, Alignable m') =>
Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool
defStop (g -> Move m m'
forall m m' g.
(Alignable m, Alignable m', IsGap g) =>
g -> Move m m'
move g
gap)

    -- Start from bottom right corner
    --
    {-# INLINE traceStart #-}
    traceStart :: GlobalAlignment g (IxValue m) (IxValue m')
-> Matrix m m' -> m -> m' -> (Index m, Index m')
traceStart = (Matrix m m' -> m -> m' -> (Index m, Index m'))
-> GlobalAlignment g (IxValue m) (IxValue m')
-> Matrix m m'
-> m
-> m'
-> (Index m, Index m')
forall a b. a -> b -> a
const Matrix m m' -> m -> m' -> (Index m, Index m')
forall m m'.
(Alignable m, Alignable m') =>
Matrix m m' -> m -> m' -> (Index m, Index m')
defStart

    scoreMatrix :: forall m m' . (Alignable m, Alignable m')
                => GlobalAlignment g (IxValue m) (IxValue m')
                -> m
                -> m'
                -> Matrix m m'
    scoreMatrix :: GlobalAlignment g (IxValue m) (IxValue m')
-> m -> m' -> Matrix m m'
scoreMatrix (GlobalAlignment Scoring (IxValue m) (IxValue m')
subC g
g) m
s m'
t | g -> Bool
forall a. IsGap a => a -> Bool
isAffine g
g = Matrix m m'
uMatrixAffine
                                             | Bool
otherwise  = Matrix m m'
uMatrixSimple
      where
        uMatrixSimple :: UArray (Index m, Index m', EditOp) Int
        uMatrixSimple :: Matrix m m'
uMatrixSimple = (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall i e. (forall s. ST s (STUArray s i e)) -> UArray i e
runSTUArray ((forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
 -> Matrix m m')
-> (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall a b. (a -> b) -> a -> b
$ do
          STUArray s (Index m, Index m', EditOp) Int
matrix <- ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
-> Int -> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Index m
lowerS, Index m'
lowerT, EditOp
Match), (Index m
nilS, Index m'
nilT, EditOp
Match)) Int
0 :: ST s (STUArray s (Index m, Index m', EditOp) Int)

          [Index m] -> (Index m -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m
lowerS .. Index m
nilS] ((Index m -> ST s ()) -> ST s ())
-> (Index m -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m
ixS ->
            [Index m'] -> (Index m' -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m'
lowerT .. Index m'
nilT] ((Index m' -> ST s ()) -> ST s ())
-> (Index m' -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m'
ixT ->

              -- Next cell = max (d_i-1,j + gap, d_i,j-1 + gap, d_i-1,j-1 + s(i,j))
              --
              if | Index m
ixS Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS -> STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ (g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Index m', Index m') -> Index m' -> Int
forall a. Ix a => (a, a) -> a -> Int
index (Index m'
lowerT, Index m'
nilT) Index m'
ixT
                 | Index m'
ixT Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT -> STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ (g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Index m, Index m) -> Index m -> Int
forall a. Ix a => (a, a) -> a -> Int
index (Index m
lowerS, Index m
nilS) Index m
ixS
                 | Bool
otherwise -> do
                   Int
predDiag <- STUArray s (Index m, Index m', EditOp) Int
matrix STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
ixS, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
ixT, EditOp
Match)
                   Int
predS    <- STUArray s (Index m, Index m', EditOp) Int
matrix STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
ixS,      Index m'
ixT, EditOp
Match)
                   Int
predT    <- STUArray s (Index m, Index m', EditOp) Int
matrix STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (     Index m
ixS, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
ixT, EditOp
Match)
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [ Int
predDiag Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Index m -> Index m' -> Int
sub Index m
ixS Index m'
ixT
                                                                 , Int
predS Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g
                                                                 , Int
predT Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g
                                                                 ]
          STUArray s (Index m, Index m', EditOp) Int
-> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (f :: * -> *) a. Applicative f => a -> f a
pure STUArray s (Index m, Index m', EditOp) Int
matrix

        uMatrixAffine :: UArray (Index m, Index m', EditOp) Int
        uMatrixAffine :: Matrix m m'
uMatrixAffine = (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall i e. (forall s. ST s (STUArray s i e)) -> UArray i e
runSTUArray ((forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
 -> Matrix m m')
-> (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall a b. (a -> b) -> a -> b
$ do
          STUArray s (Index m, Index m', EditOp) Int
matrix <- ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
-> Int -> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Index m
lowerS, Index m'
lowerT, EditOp
Insert), (Index m
nilS, Index m'
nilT, EditOp
Match)) Int
0 :: ST s (STUArray s (Index m, Index m', EditOp) Int)
          [Index m] -> (Index m -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m
lowerS .. Index m
nilS] ((Index m -> ST s ()) -> ST s ())
-> (Index m -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m
ixS ->
            [Index m'] -> (Index m' -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m'
lowerT .. Index m'
nilT] ((Index m' -> ST s ()) -> ST s ())
-> (Index m' -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m'
ixT ->

              -- Next cell = max (d_i-1,j + gap, d_i,j-1 + gap, d_i-1,j-1 + s(i,j))
              -- Matrices with gap costs are also filled as follows:
              -- gepMatrix[i, j] <- gapExtend if one of strings has gap at this position else gapOpen
              --
              if | Index m
ixS Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS Bool -> Bool -> Bool
&& Index m'
ixT Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT -> do
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) Int
0
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Insert) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Delete) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g
                 | Index m
ixS Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS -> do
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (g -> Int
forall a. IsGap a => a -> Int
insertCostExtend g
g) Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int -> Int
forall a. Enum a => a -> a
pred ((Index m', Index m') -> Index m' -> Int
forall a. Ix a => (a, a) -> a -> Int
index (Index m'
lowerT, Index m'
nilT) Index m'
ixT)
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Insert) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
insertCostExtend g
g
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Delete) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g
                 | Index m'
ixT Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT -> do
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (g -> Int
forall a. IsGap a => a -> Int
deleteCostExtend g
g) Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int -> Int
forall a. Enum a => a -> a
pred ((Index m, Index m) -> Index m -> Int
forall a. Ix a => (a, a) -> a -> Int
index (Index m
lowerS, Index m
nilS) Index m
ixS)
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Delete) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
deleteCostExtend g
g
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Insert) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g
                 | Bool
otherwise -> g
-> STUArray s (Index m, Index m', EditOp) Int
-> Bool
-> (Index m -> Index m' -> Int)
-> Index m
-> Index m'
-> ST s ()
forall g ix ix' s.
(IsGap g, Ix ix, Ix ix', Enum ix, Enum ix') =>
g
-> STUArray s (ix, ix', EditOp) Int
-> Bool
-> (ix -> ix' -> Int)
-> ix
-> ix'
-> ST s ()
fillInnerMatrix g
g STUArray s (Index m, Index m', EditOp) Int
matrix Bool
False Index m -> Index m' -> Int
sub Index m
ixS Index m'
ixT
          STUArray s (Index m, Index m', EditOp) Int
-> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (f :: * -> *) a. Applicative f => a -> f a
pure STUArray s (Index m, Index m', EditOp) Int
matrix

        (Index m
lowerS, Index m
upperS) = m -> (Index m, Index m)
forall m. ChainLike m => m -> (Index m, Index m)
bounds m
s
        (Index m'
lowerT, Index m'
upperT) = m' -> (Index m', Index m')
forall m. ChainLike m => m -> (Index m, Index m)
bounds m'
t
        nilS :: Index m
nilS = Index m -> Index m
forall a. Enum a => a -> a
succ Index m
upperS
        nilT :: Index m'
nilT = Index m' -> Index m'
forall a. Enum a => a -> a
succ Index m'
upperT

        sub :: Index m -> Index m' -> Int
        sub :: Index m -> Index m' -> Int
sub = Scoring (IxValue m) (IxValue m')
-> m -> m' -> Index m -> Index m' -> Int
forall m m'.
(Alignable m, Alignable m') =>
(IxValue m -> IxValue m' -> Int)
-> m -> m' -> Index m -> Index m' -> Int
substitute Scoring (IxValue m) (IxValue m')
subC m
s m'
t

instance IsGap g => SequenceAlignment (LocalAlignment g) where

    -- Conditions of traceback are described below
    --
    {-# INLINE cond #-}
    cond :: LocalAlignment g (IxValue m) (IxValue m') -> Conditions m m'
cond (LocalAlignment Scoring (IxValue m) (IxValue m')
_ g
gap) = Stop m m' -> Move m m' -> Conditions m m'
forall m m'. Stop m m' -> Move m m' -> Conditions m m'
Conditions Stop m m'
forall m m'.
(Alignable m, Alignable m') =>
Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool
localStop (g -> Move m m'
forall m m' g.
(Alignable m, Alignable m', IsGap g) =>
g -> Move m m'
move g
gap)

    -- Start from bottom right corner
    --
    {-# INLINE traceStart #-}
    traceStart :: LocalAlignment g (IxValue m) (IxValue m')
-> Matrix m m' -> m -> m' -> (Index m, Index m')
traceStart = (Matrix m m' -> m -> m' -> (Index m, Index m'))
-> LocalAlignment g (IxValue m) (IxValue m')
-> Matrix m m'
-> m
-> m'
-> (Index m, Index m')
forall a b. a -> b -> a
const Matrix m m' -> m -> m' -> (Index m, Index m')
forall m m'.
(Alignable m, Alignable m') =>
Matrix m m' -> m -> m' -> (Index m, Index m')
localStart

    scoreMatrix :: forall m m' . (Alignable m, Alignable m')
                => LocalAlignment g (IxValue m) (IxValue m')
                -> m
                -> m'
                -> Matrix m m'
    scoreMatrix :: LocalAlignment g (IxValue m) (IxValue m') -> m -> m' -> Matrix m m'
scoreMatrix (LocalAlignment Scoring (IxValue m) (IxValue m')
subC g
g) m
s m'
t | g -> Bool
forall a. IsGap a => a -> Bool
isAffine g
g = Matrix m m'
uMatrixAffine
                                            | Bool
otherwise  = Matrix m m'
uMatrixSimple
      where
        uMatrixSimple :: UArray (Index m, Index m', EditOp) Int
        uMatrixSimple :: Matrix m m'
uMatrixSimple = (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall i e. (forall s. ST s (STUArray s i e)) -> UArray i e
runSTUArray ((forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
 -> Matrix m m')
-> (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall a b. (a -> b) -> a -> b
$ do
          STUArray s (Index m, Index m', EditOp) Int
matrix <- ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
-> Int -> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Index m
lowerS, Index m'
lowerT, EditOp
Match), (Index m
nilS, Index m'
nilT, EditOp
Match)) Int
0  :: ST s (STUArray s (Index m, Index m', EditOp) Int)
          [Index m] -> (Index m -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m
lowerS .. Index m
nilS] ((Index m -> ST s ()) -> ST s ())
-> (Index m -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m
ixS ->
            [Index m'] -> (Index m' -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m'
lowerT .. Index m'
nilT] ((Index m' -> ST s ()) -> ST s ())
-> (Index m' -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m'
ixT ->

              -- Next cell = max (d_i-1,j + gap, d_i,j-1 + gap, d_i-1,j-1 + s(i,j))
              --
              if | Index m
ixS Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS -> STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) Int
0
                 | Index m'
ixT Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT -> STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) Int
0
                 | Bool
otherwise -> do
                   Int
predDiag <- STUArray s (Index m, Index m', EditOp) Int
matrix STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
ixS, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
ixT, EditOp
Match)
                   Int
predS    <- STUArray s (Index m, Index m', EditOp) Int
matrix STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
ixS,      Index m'
ixT, EditOp
Match)
                   Int
predT    <- STUArray s (Index m, Index m', EditOp) Int
matrix STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (     Index m
ixS, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
ixT, EditOp
Match)
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [ Int
predDiag Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Index m -> Index m' -> Int
sub Index m
ixS Index m'
ixT
                                                                 , Int
predS Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g
                                                                 , Int
predT Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g
                                                                 , Int
0
                                                                 ]
          STUArray s (Index m, Index m', EditOp) Int
-> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (f :: * -> *) a. Applicative f => a -> f a
pure STUArray s (Index m, Index m', EditOp) Int
matrix

        uMatrixAffine :: UArray (Index m, Index m', EditOp) Int
        uMatrixAffine :: Matrix m m'
uMatrixAffine = (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall i e. (forall s. ST s (STUArray s i e)) -> UArray i e
runSTUArray ((forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
 -> Matrix m m')
-> (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall a b. (a -> b) -> a -> b
$ do
          STUArray s (Index m, Index m', EditOp) Int
matrix <- ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
-> Int -> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Index m
lowerS, Index m'
lowerT, EditOp
Insert), (Index m
nilS, Index m'
nilT, EditOp
Match)) Int
0 :: ST s (STUArray s (Index m, Index m', EditOp) Int)
          [Index m] -> (Index m -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m
lowerS .. Index m
nilS] ((Index m -> ST s ()) -> ST s ())
-> (Index m -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m
ixS ->
            [Index m'] -> (Index m' -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m'
lowerT .. Index m'
nilT] ((Index m' -> ST s ()) -> ST s ())
-> (Index m' -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m'
ixT ->

              -- Next cell = max (d_i-1,j + gap, d_i,j-1 + gap, d_i-1,j-1 + s(i,j))
              -- Matrices with gap costs are also filled as follows:
              -- gepMatrix[i, j] <- gapExtend if one of strings has gap at this position else gapOpen
              --
              if | Index m
ixS Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS Bool -> Bool -> Bool
|| Index m'
ixT Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT -> do
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match)  Int
0
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Insert) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Delete) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g
                 | Bool
otherwise -> g
-> STUArray s (Index m, Index m', EditOp) Int
-> Bool
-> (Index m -> Index m' -> Int)
-> Index m
-> Index m'
-> ST s ()
forall g ix ix' s.
(IsGap g, Ix ix, Ix ix', Enum ix, Enum ix') =>
g
-> STUArray s (ix, ix', EditOp) Int
-> Bool
-> (ix -> ix' -> Int)
-> ix
-> ix'
-> ST s ()
fillInnerMatrix g
g STUArray s (Index m, Index m', EditOp) Int
matrix Bool
True Index m -> Index m' -> Int
sub Index m
ixS Index m'
ixT
          STUArray s (Index m, Index m', EditOp) Int
-> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (f :: * -> *) a. Applicative f => a -> f a
pure STUArray s (Index m, Index m', EditOp) Int
matrix

        (Index m
lowerS, Index m
upperS) = m -> (Index m, Index m)
forall m. ChainLike m => m -> (Index m, Index m)
bounds m
s
        (Index m'
lowerT, Index m'
upperT) = m' -> (Index m', Index m')
forall m. ChainLike m => m -> (Index m, Index m)
bounds m'
t
        nilS :: Index m
nilS = Index m -> Index m
forall a. Enum a => a -> a
succ Index m
upperS
        nilT :: Index m'
nilT = Index m' -> Index m'
forall a. Enum a => a -> a
succ Index m'
upperT

        sub :: Index m -> Index m' -> Int
        sub :: Index m -> Index m' -> Int
sub = Scoring (IxValue m) (IxValue m')
-> m -> m' -> Index m -> Index m' -> Int
forall m m'.
(Alignable m, Alignable m') =>
(IxValue m -> IxValue m' -> Int)
-> m -> m' -> Index m -> Index m' -> Int
substitute Scoring (IxValue m) (IxValue m')
subC m
s m'
t

instance IsGap g => SequenceAlignment (SemiglobalAlignment g) where

    -- The alignment is semiglobal, so we have to perform some additional operations
    --
    {-# INLINE semi #-}
    semi :: SemiglobalAlignment g e1 e2 -> Bool
semi = Bool -> SemiglobalAlignment g e1 e2 -> Bool
forall a b. a -> b -> a
const Bool
True

    -- Conditions of traceback are described below
    --
    {-# INLINE cond #-}
    cond :: SemiglobalAlignment g (IxValue m) (IxValue m') -> Conditions m m'
cond (SemiglobalAlignment Scoring (IxValue m) (IxValue m')
_ g
gap) = Stop m m' -> Move m m' -> Conditions m m'
forall m m'. Stop m m' -> Move m m' -> Conditions m m'
Conditions Stop m m'
forall m m'.
(Alignable m, Alignable m') =>
Matrix m m' -> m -> m' -> Index m -> Index m' -> Bool
defStop (g -> Move m m'
forall m m' g.
(Alignable m, Alignable m', IsGap g) =>
g -> Move m m'
move g
gap)

    -- Start from bottom right corner
    --
    {-# INLINE traceStart #-}
    traceStart :: SemiglobalAlignment g (IxValue m) (IxValue m')
-> Matrix m m' -> m -> m' -> (Index m, Index m')
traceStart = (Matrix m m' -> m -> m' -> (Index m, Index m'))
-> SemiglobalAlignment g (IxValue m) (IxValue m')
-> Matrix m m'
-> m
-> m'
-> (Index m, Index m')
forall a b. a -> b -> a
const Matrix m m' -> m -> m' -> (Index m, Index m')
forall m m'.
(Alignable m, Alignable m') =>
Matrix m m' -> m -> m' -> (Index m, Index m')
semiStart

    scoreMatrix :: forall m m' . (Alignable m, Alignable m')
                => SemiglobalAlignment g (IxValue m) (IxValue m')
                -> m
                -> m'
                -> Matrix m m'
    scoreMatrix :: SemiglobalAlignment g (IxValue m) (IxValue m')
-> m -> m' -> Matrix m m'
scoreMatrix (SemiglobalAlignment Scoring (IxValue m) (IxValue m')
subC g
g) m
s m'
t | g -> Bool
forall a. IsGap a => a -> Bool
isAffine g
g = Matrix m m'
uMatrixAffine
                                                 | Bool
otherwise  = Matrix m m'
uMatrixSimple
      where
        uMatrixSimple :: UArray (Index m, Index m', EditOp) Int
        uMatrixSimple :: Matrix m m'
uMatrixSimple = (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall i e. (forall s. ST s (STUArray s i e)) -> UArray i e
runSTUArray ((forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
 -> Matrix m m')
-> (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall a b. (a -> b) -> a -> b
$ do
          STUArray s (Index m, Index m', EditOp) Int
matrix <- ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
-> Int -> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Index m
lowerS, Index m'
lowerT, EditOp
Match), (Index m
nilS, Index m'
nilT, EditOp
Match)) Int
0  :: ST s (STUArray s (Index m, Index m', EditOp) Int)
          [Index m] -> (Index m -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m
lowerS .. Index m
nilS] ((Index m -> ST s ()) -> ST s ())
-> (Index m -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m
ixS ->
            [Index m'] -> (Index m' -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m'
lowerT .. Index m'
nilT] ((Index m' -> ST s ()) -> ST s ())
-> (Index m' -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m'
ixT ->

              -- Next cell = max (d_i-1,j + gap, d_i,j-1 + gap, d_i-1,j-1 + s(i,j))
              --
              if | Index m
ixS Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS -> STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) Int
0
                 | Index m'
ixT Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT -> STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) Int
0
                 | Bool
otherwise -> do
                   Int
predDiag <- STUArray s (Index m, Index m', EditOp) Int
matrix STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
ixS, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
ixT, EditOp
Match)
                   Int
predS    <- STUArray s (Index m, Index m', EditOp) Int
matrix STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (Index m -> Index m
forall a. Enum a => a -> a
pred Index m
ixS,      Index m'
ixT, EditOp
Match)
                   Int
predT    <- STUArray s (Index m, Index m', EditOp) Int
matrix STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (     Index m
ixS, Index m' -> Index m'
forall a. Enum a => a -> a
pred Index m'
ixT, EditOp
Match)
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [ Int
predDiag Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Index m -> Index m' -> Int
sub Index m
ixS Index m'
ixT
                                                                 , Int
predS Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g
                                                                 , Int
predT Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g
                                                                 ]
          STUArray s (Index m, Index m', EditOp) Int
-> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (f :: * -> *) a. Applicative f => a -> f a
pure STUArray s (Index m, Index m', EditOp) Int
matrix

        uMatrixAffine :: UArray (Index m, Index m', EditOp) Int
        uMatrixAffine :: Matrix m m'
uMatrixAffine = (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall i e. (forall s. ST s (STUArray s i e)) -> UArray i e
runSTUArray ((forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
 -> Matrix m m')
-> (forall s. ST s (STUArray s (Index m, Index m', EditOp) Int))
-> Matrix m m'
forall a b. (a -> b) -> a -> b
$ do
          STUArray s (Index m, Index m', EditOp) Int
matrix <- ((Index m, Index m', EditOp), (Index m, Index m', EditOp))
-> Int -> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> e -> m (a i e)
newArray ((Index m
lowerS, Index m'
lowerT, EditOp
Insert), (Index m
nilS, Index m'
nilT, EditOp
Match)) Int
0 :: ST s (STUArray s (Index m, Index m', EditOp) Int)
          [Index m] -> (Index m -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m
lowerS .. Index m
nilS] ((Index m -> ST s ()) -> ST s ())
-> (Index m -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m
ixS ->
            [Index m'] -> (Index m' -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Index m'
lowerT .. Index m'
nilT] ((Index m' -> ST s ()) -> ST s ())
-> (Index m' -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Index m'
ixT ->
              if | Index m
ixS Index m -> Index m -> Bool
forall a. Eq a => a -> a -> Bool
== Index m
lowerS Bool -> Bool -> Bool
|| Index m'
ixT Index m' -> Index m' -> Bool
forall a. Eq a => a -> a -> Bool
== Index m'
lowerT -> do
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Match)  Int
0
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Insert) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g
                   STUArray s (Index m, Index m', EditOp) Int
-> (Index m, Index m', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (Index m, Index m', EditOp) Int
matrix (Index m
ixS, Index m'
ixT, EditOp
Delete) (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g
                 | Bool
otherwise -> g
-> STUArray s (Index m, Index m', EditOp) Int
-> Bool
-> (Index m -> Index m' -> Int)
-> Index m
-> Index m'
-> ST s ()
forall g ix ix' s.
(IsGap g, Ix ix, Ix ix', Enum ix, Enum ix') =>
g
-> STUArray s (ix, ix', EditOp) Int
-> Bool
-> (ix -> ix' -> Int)
-> ix
-> ix'
-> ST s ()
fillInnerMatrix g
g STUArray s (Index m, Index m', EditOp) Int
matrix Bool
False Index m -> Index m' -> Int
sub Index m
ixS Index m'
ixT
          STUArray s (Index m, Index m', EditOp) Int
-> ST s (STUArray s (Index m, Index m', EditOp) Int)
forall (f :: * -> *) a. Applicative f => a -> f a
pure STUArray s (Index m, Index m', EditOp) Int
matrix

        (Index m
lowerS, Index m
upperS) = m -> (Index m, Index m)
forall m. ChainLike m => m -> (Index m, Index m)
bounds m
s
        (Index m'
lowerT, Index m'
upperT) = m' -> (Index m', Index m')
forall m. ChainLike m => m -> (Index m, Index m)
bounds m'
t
        nilS :: Index m
nilS = Index m -> Index m
forall a. Enum a => a -> a
succ Index m
upperS
        nilT :: Index m'
nilT = Index m' -> Index m'
forall a. Enum a => a -> a
succ Index m'
upperT

        sub :: Index m -> Index m' -> Int
        sub :: Index m -> Index m' -> Int
sub = Scoring (IxValue m) (IxValue m')
-> m -> m' -> Index m -> Index m' -> Int
forall m m'.
(Alignable m, Alignable m') =>
(IxValue m -> IxValue m' -> Int)
-> m -> m' -> Index m -> Index m' -> Int
substitute Scoring (IxValue m) (IxValue m')
subC m
s m'
t

{-# INLINE fillInnerMatrix #-}
fillInnerMatrix
  :: (IsGap g, Ix ix, Ix ix', Enum ix, Enum ix')
  => g
  -> STUArray s (ix, ix', EditOp) Int
  -> Bool               -- ^ Is this local alignment?
  -> (ix -> ix' -> Int) -- ^ Substitution function
  -> ix
  -> ix'
  -> ST s ()
fillInnerMatrix :: g
-> STUArray s (ix, ix', EditOp) Int
-> Bool
-> (ix -> ix' -> Int)
-> ix
-> ix'
-> ST s ()
fillInnerMatrix g
g STUArray s (ix, ix', EditOp) Int
matrix Bool
isLocal ix -> ix' -> Int
sub ix
ixS ix'
ixT = do
    Int
predDiag <- STUArray s (ix, ix', EditOp) Int
matrix STUArray s (ix, ix', EditOp) Int -> (ix, ix', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (ix -> ix
forall a. Enum a => a -> a
pred ix
ixS, ix' -> ix'
forall a. Enum a => a -> a
pred ix'
ixT, EditOp
Match)
    Int
predS    <- STUArray s (ix, ix', EditOp) Int
matrix STUArray s (ix, ix', EditOp) Int -> (ix, ix', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (ix -> ix
forall a. Enum a => a -> a
pred ix
ixS,      ix'
ixT, EditOp
Match)
    Int
predT    <- STUArray s (ix, ix', EditOp) Int
matrix STUArray s (ix, ix', EditOp) Int -> (ix, ix', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (     ix
ixS, ix' -> ix'
forall a. Enum a => a -> a
pred ix'
ixT, EditOp
Match)

    Int
delCost  <- STUArray s (ix, ix', EditOp) Int
matrix STUArray s (ix, ix', EditOp) Int -> (ix, ix', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (ix -> ix
forall a. Enum a => a -> a
pred ix
ixS,      ix'
ixT, EditOp
Delete)
    Int
insCost  <- STUArray s (ix, ix', EditOp) Int
matrix STUArray s (ix, ix', EditOp) Int -> (ix, ix', EditOp) -> ST s Int
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
`readArray` (     ix
ixS, ix' -> ix'
forall a. Enum a => a -> a
pred ix'
ixT, EditOp
Insert)

    let
      updInsCost :: Int
updInsCost = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Int
insCost Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
insertCostExtend g
g) (Int
predT Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
insertCostOpen g
g)
      updDelCost :: Int
updDelCost = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Int
delCost Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
deleteCostExtend g
g) (Int
predS Int -> Int -> Int
forall a. Num a => a -> a -> a
+ g -> Int
forall a. IsGap a => a -> Int
deleteCostOpen g
g)
      repCost :: Int
repCost = Int
predDiag Int -> Int -> Int
forall a. Num a => a -> a -> a
+ ix -> ix' -> Int
sub ix
ixS ix'
ixT
      maxCost :: Int
maxCost = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
repCost (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
updInsCost Int
updDelCost

      finalCost :: Int
finalCost = if Bool
isLocal then Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 Int
maxCost else Int
maxCost

    STUArray s (ix, ix', EditOp) Int
-> (ix, ix', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (ix, ix', EditOp) Int
matrix (ix
ixS, ix'
ixT, EditOp
Delete) Int
updDelCost
    STUArray s (ix, ix', EditOp) Int
-> (ix, ix', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (ix, ix', EditOp) Int
matrix (ix
ixS, ix'
ixT, EditOp
Insert) Int
updInsCost
    STUArray s (ix, ix', EditOp) Int
-> (ix, ix', EditOp) -> Int -> ST s ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray STUArray s (ix, ix', EditOp) Int
matrix (ix
ixS, ix'
ixT, EditOp
Match) Int
finalCost