diff --git a/BioInf/RNAwolf/Bulge.hs b/BioInf/RNAwolf/Bulge.hs
new file mode 100644
--- /dev/null
+++ b/BioInf/RNAwolf/Bulge.hs
@@ -0,0 +1,143 @@
+{-# LANGUAGE RecordWildCards #-}
+
+module BioInf.RNAwolf.Bulge where
+
+import qualified Data.Vector.Unboxed as VU
+
+import Biobase.Primary
+import Biobase.Secondary
+import Data.PrimitiveArray
+import Data.PrimitiveArray.Ix
+
+import BioInf.Params
+import BioInf.RNAwolf.Types
+
+
+
+-- * Outer part of a bulge
+
+-- | The outer closing pair of a bulge loop (one unpaired region).
+
+fBulgeOuter :: BaseF (NBulgLoop -> ExtFeatures (VU.Vector (PairIdx,Double)))
+fBulgeOuter Params{..} inp (NBulgLoop nBulgLoop) i j ct eI eJ
+  | i<0 || j>n = error $ "fBulgeOuter: " ++ show (i,j)
+  | otherwise = VU.singleton s
+  where
+    s =
+      ( (i,j)
+      , nBulgLoop ! (i,j)
+      + bulgeClose ! ((nI,nJ),(ct,eI,eJ))
+      )
+    nI = inp VU.! i
+    nJ = inp VU.! j
+    n = VU.length inp -1
+{-# INLINE fBulgeOuter #-}
+
+-- | Outer part of a normal bulge
+
+btBulgeOuter
+  :: Params
+  -> Primary
+  -> EStem
+  -> NBulgLoop
+  -> NBT
+  -> ExtBT
+btBulgeOuter ps inp (EStem eStem) nBulgLoop btBULoop i j ct eI eJ d =
+  [ (ij:x,z)
+  | i>=0,i<j,j<=n
+  , ((k,l),enext) <- VU.toList $ fBulgeOuter ps inp nBulgLoop i j ct eI eJ
+  , let d' = newD d ehere enext
+--  , trc' ("btESbulge",ij,d') $ testD d'
+  , testD d'
+  , (x,z) <- btBULoop k l d'
+  ] where
+      ij = ((i,j),(ct,eI,eJ))
+      ehere = eStem!ij
+      n = VU.length inp -1
+
+
+
+-- * Loop part of a bulge
+
+-- | The loop-part of bulges. Increases speed by 2x
+
+fBulgeLoop :: BaseF (NBulg -> Features (VU.Vector (PairIdx,Double)))
+fBulgeLoop Params{..} inp (NBulg nBulg) i j
+  | j-i<=6 = VU.empty
+  | otherwise = VU.map f kls
+  where
+    f (k,l) =
+      ( (k,l)
+      , nBulg ! (k,l)
+      + bulgeLength ! (max (k-i-1) (j-l-1))
+      )
+    kls = ks VU.++ ls
+    ks = VU.fromList [ (k,l)
+                     | k <- takeWhile (\k -> k-i-1<=maxLength) [i+2 .. j-4], let l = j-1
+                     ]
+    ls = VU.fromList [ (k,l)
+                     | let k = i+1, l <- takeWhile (\l -> j-l-1<=maxLength) [j-2,j-3 .. i+4]
+                     ]
+{-# INLINE fBulgeLoop #-}
+
+-- | Backtrack the bulge loop part.
+
+btBulgeLoop
+  :: Params
+  -> Primary
+  -> NBulgLoop
+  -> NBulg
+  -> NBT
+  -> NBT
+btBulgeLoop ps inp (NBulgLoop nBulgLoop) nBulg btBU i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ((k,l),enext) <- VU.toList $ fBulgeLoop ps inp nBulg i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btBU k l d'
+  , testD z
+  ] where
+      ehere = nBulgLoop!(i,j)
+      n = VU.length inp -1
+
+
+-- * Inner part of a bulge
+
+-- | Inner part of a bulge to speed up bulge calculations
+
+fBulgeInner :: BaseF (EStem -> Features (VU.Vector (ExtPairIdx,Double)))
+fBulgeInner Params{..} inp (EStem eStem) i j
+  | j-i<2 = VU.empty
+  | otherwise = VU.map f kls
+  where
+    f ijExt@((i,j),(ctIJ,eI,eJ)) =
+      ( ijExt
+      , eStem ! ijExt
+      + bulgeClose ! ((nJ,nI),(ctIJ,eJ,eI))
+      ) where nI = inp VU.! i
+              nJ = inp VU.! j
+    kls = VU.fromList [ ((i,j),(ctIJ,eI,eJ))
+                      | eI<-wsh, eJ<-wsh, ctIJ<-citr
+                      ]
+{-# INLINE fBulgeInner #-}
+
+-- | Backtrack the inner part of a bulge
+
+btBulgeInner
+  :: Params
+  -> Primary
+  -> NBulg
+  -> EStem
+  -> ExtBT
+  -> NBT
+btBulgeInner ps inp (NBulg nBulg) eStem btES i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ((_,(eI,eJ,ct)),enext) <- VU.toList $ fBulgeInner ps inp eStem i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btES i j eI eJ ct d'
+  ] where
+      n = VU.length inp -1
+      ehere = nBulg!(i,j)
diff --git a/BioInf/RNAwolf/Extern.hs b/BioInf/RNAwolf/Extern.hs
new file mode 100644
--- /dev/null
+++ b/BioInf/RNAwolf/Extern.hs
@@ -0,0 +1,172 @@
+{-# LANGUAGE RecordWildCards #-}
+
+-- | External loops are complete substructures, of which zero to many sit on
+-- the chain of nucleotides.
+
+module BioInf.RNAwolf.Extern where
+
+import qualified Data.Vector.Unboxed as VU
+
+import Biobase.Primary
+import Biobase.Secondary
+import Data.PrimitiveArray
+import Data.PrimitiveArray.Ix
+
+import BioInf.Params
+import BioInf.RNAwolf.Types
+
+
+
+
+-- * Unpaired left nucleotide
+
+-- | An external loop with an unpaired nucleotide to the left
+
+fLeftUnpaired :: BaseF (NExtn -> Features (VU.Vector (PairIdx,Double)))
+fLeftUnpaired Params{..} inp (NExtn nExtn) i j
+  | i<0 || j>n || i>=j = error $ "Extern.fLeftUnpaired: " ++ show (i,j)
+  | otherwise = VU.singleton s
+  where
+    s =
+      ( (i+1,j)
+      , nExtn ! (i+1,j)
+      )
+    n = VU.length inp -1
+{-# INLINE fLeftUnpaired #-}
+
+-- | Backtracking a structure with an unpaired nucleotide to the left.
+
+btLeftUnpaired
+  :: Params
+  -> Primary
+  -> NExtn
+  -> NBT
+  -> NBT
+btLeftUnpaired ps inp nExtn btE i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ((k,l),enext) <- VU.toList $ fLeftUnpaired ps inp nExtn i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btE k l d'
+  , testD z
+  ] where
+      ehere = unNExtn nExtn !(i,j)
+      n = VU.length inp -1
+
+
+
+-- * exactly one stem with indices (i,k), i<k<=j
+
+-- | Energy for exactly one stem at (i,k)
+
+fStem :: BaseF (NStem -> Features (VU.Vector (PairIdx,Double)))
+fStem Params{..} inp (NStem nStem) i j
+  | i<0 || j>n || i>=j = error $ "Extern.fStem: " ++ show (i,j)
+  | otherwise = VU.map f $ VU.enumFromN (i+1) (j-i)
+  where
+    f k = ( (i,k)
+          , nStem !(i,k)
+          )
+    n = VU.length inp -1
+{-# INLINE fStem #-}
+
+-- | Backtrack one stem with right index k.
+
+btStem
+  :: Params
+  -> Primary
+  -> NExtn
+  -> NStem
+  -> NBT
+  -> NBT
+btStem ps inp nExtn nStem btNStem i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ((_,k),enext) <- VU.toList $ fStem ps inp nStem i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btNStem i k d'
+  , testD z
+  ] where
+      ehere = unNExtn nExtn !(i,j)
+      n = VU.length inp -1
+
+
+
+-- * The neutral element for energy. This is an unpaired stretch between (i,j)
+
+-- | This one is important as otherwise, some stretches of nucleotides will
+-- always have to be paired. (Obviously, I forgot to add this one for a
+-- time...)
+
+fOne :: BaseF (Features (VU.Vector (PairIdx,Double)))
+fOne Params{..} inp i j
+  | i<0 || j>n || i>j = error $ "Extern.fOne: " ++ show (i,j)
+  | otherwise = VU.singleton s
+  where
+    s = ( (i,j)
+        , 0
+        )
+    n = VU.length inp -1
+{-# INLINE fOne #-}
+
+btOne
+  :: Params
+  -> Primary
+  -> NExtn
+  -> NBT
+btOne ps inp nExtn i j d =
+  [ (x,z)
+  | i>=0,i<=j,j<=n
+  , (_,enext) <- VU.toList $ fOne ps inp i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , let x = []
+  , let z = d'
+  , testD z
+  ] where
+      ehere = unNExtn nExtn !(i,j)
+      n = VU.length inp -1
+
+
+
+-- * Two or more stems in the external structure
+
+-- | External structures with more than one stem have a NStem on the left and
+-- an external NExtn structure on the right.
+
+fStems :: BaseF (NStem -> NExtn -> Features (VU.Vector (Int,Double)))
+fStems Params{..} inp (NStem nStem) (NExtn nExtn) i j
+  | i<0 || j>n || i>j = error $ "Extern.fStems: " ++ show (i,j)
+  | otherwise = VU.map f $ VU.enumFromN (i+1) (j-i-1)
+  where
+    f k = ( k
+          , nStem !(i,k) + nExtn !(k+1,j)
+          )
+    n = VU.length inp -1
+{-# INLINE fStems #-}
+
+-- | Backtracking of an external structure with more than one stem
+
+btStems
+  :: Params
+  -> Primary
+  -> NStem
+  -> NExtn
+  -> NBT
+  -> NBT
+  -> NBT
+btStems ps inp nStem nExtn btNS btE i j d =
+  [ (x++y,z)
+  | i>=0,i<j,j<=n
+  , (k,enext) <- VU.toList $ fStems ps inp nStem nExtn i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z') <- btNS i k d'
+  , testD z'
+  , (y,z) <- btE (k+1) j z'
+  , testD z
+  ] where
+      ehere = unNExtn nExtn !(i,j)
+      n = VU.length inp -1
diff --git a/BioInf/RNAwolf/Hairpin.hs b/BioInf/RNAwolf/Hairpin.hs
new file mode 100644
--- /dev/null
+++ b/BioInf/RNAwolf/Hairpin.hs
@@ -0,0 +1,59 @@
+{-# LANGUAGE RecordWildCards #-}
+
+module BioInf.RNAwolf.Hairpin where
+
+import qualified Data.Vector.Unboxed as VU
+
+import Biobase.Primary
+import Biobase.Secondary
+import Data.PrimitiveArray
+import Data.PrimitiveArray.Ix
+
+import BioInf.Params
+import BioInf.RNAwolf.Types
+
+-- | A hairpin is a number of 0 or more unpaired nucleotides, enclosed by the
+-- nucleotides (i,j) which pair.
+--
+-- TODO should we allow loops with more than 30 unpaired nucleotides?
+--
+-- TODO should we allow hairpins with no unpaired nucleotides in the pin? They
+-- do occur, but only under special circumstances which we should model
+-- differently...
+--
+-- TODO re-allow IMI
+
+fHairpin :: [Int] -> BaseF (ExtFeatures (VU.Vector (ExtPairIdx,Double)))
+fHairpin imi Params{..} inp i j ct eI eJ
+--  | checkIMI imi    = VU.singleton (k,interMolInit)
+  | j-i<3           = VU.empty
+  | j-i-1>maxLength = VU.empty
+  | otherwise       = VU.singleton (k,v)
+  where
+    k = ((i,j),(ct,eI,eJ))
+    v = 0
+      + hairpinLength ! (j-i-1)
+      + hairpinClose  ! (((inp VU.! i,inp VU.! j),(ct,eI,eJ)),inp VU.! (i+1), inp VU.! (j-1))
+    checkIMI [] = False
+    checkIMI (x:xs) = i<x && j>x || checkIMI xs
+--      + pairDistance  ! (j-i-1)
+{-# INLINE fHairpin #-}
+
+-- | Backtracking hairpins.
+
+btHairpin
+  :: Params
+  -> Primary
+  -> EStem
+  -> ExtBT
+btHairpin ps inp (EStem eStem) i j ct eI eJ d =
+  [ ([ij],d')
+  | i>=0,j-i>=3,j<=n
+  , (_,enext) <- VU.toList $ fHairpin imi ps inp i j ct eI eJ
+  , let d' = newD d ehere enext
+  , testD d'
+  ] where
+      ij = ((i,j),(ct,eI,eJ))
+      ehere = eStem!ij
+      n = VU.length inp -1
+      imi = map fst . filter ((==nIMI).snd) $ zip [0..] (VU.toList inp)
diff --git a/BioInf/RNAwolf/Interior.hs b/BioInf/RNAwolf/Interior.hs
new file mode 100644
--- /dev/null
+++ b/BioInf/RNAwolf/Interior.hs
@@ -0,0 +1,196 @@
+
+{-# LANGUAGE RecordWildCards #-}
+
+module BioInf.RNAwolf.Interior where
+
+import qualified Data.Vector.Unboxed as VU
+
+import Biobase.Primary
+import Biobase.Secondary
+import Data.PrimitiveArray
+import Data.PrimitiveArray.Ix
+
+import BioInf.Params
+import BioInf.RNAwolf.Types
+
+
+
+-- * Outer part
+
+-- | The outer part of an interior loop. Given a certain basepair type, add the
+-- cost from the unpaired part.
+
+fInteriorOuter :: BaseF (NInteLoop -> ExtFeatures (VU.Vector (PairIdx,Double)))
+fInteriorOuter Params{..} inp (NInteLoop nInteLoop) i j ct eI eJ
+  | j-i<4     = VU.empty
+  | otherwise = VU.map f $ VU.singleton (i,j)
+  where
+    f (k,l) = ( (i,j)
+              , nInteLoop ! (i,j)
+              + ijSc
+--              + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0
+              )
+    ijSc = interiorClose ! (((nI,nJ),(ct,eI,eJ)),nIp1,nJm1)
+    nI   = inp VU.! i
+    nJ   = inp VU.! j
+    nIp1 = inp VU.! (i+1)
+    nJm1 = inp VU.! (j-1)
+{-# INLINE fInteriorOuter #-}
+
+-- |
+
+btInteriorOuter
+  :: Params
+  -> Primary
+  -> EStem
+  -> NInteLoop
+  -> NBT -- recursive backtracking function for loops
+  -> ExtBT
+btInteriorOuter ps inp (EStem eStem) nInteLoop btILoop i j ct eI eJ d =
+  [ (ij:x,z) -- interior loop
+  | i>=0,i<j,j<=n
+  , ((k,l),enext) <- VU.toList $ fInteriorOuter ps inp nInteLoop i j ct eI eJ
+  , i<k && l<j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btILoop k l d'
+  , testD z
+  ] where
+      ij = ((i,j),(ct,eI,eJ))
+      ehere = eStem!ij
+      n = VU.length inp -1
+
+
+
+-- * Loop part
+
+-- | Performs the interior loop calculations between (i,j) "outer" and (k,l)
+-- "inner" part. The score based on the unpaired nucleotides is independent of
+-- both, the outer and the inner basepair type.
+--
+-- NOTE / TODO -- fusion enabled for this function (due to it taking 50% of the
+-- time), full fusion is still dependent on other factors and needs to be
+-- checked (in particular, we still have allocation events)
+
+fInteriorLoop :: BaseF (NInte -> Features (VU.Vector (PairIdx,Double)))
+fInteriorLoop Params{..} inp (NInte nInte) i j
+  | j-i<4 = VU.empty
+  | otherwise = VU.map f kls
+  where
+    f (k,l) = ( (k,l)
+              , nInte ! (k,l)
+              + interiorLength ! (lenI+lenJ)
+              + interiorAsym ! (abs $ lenI - lenJ)
+              ) where lenI = k-i-1; lenJ = j-l-1
+    kls  = VU.map (\(dI,dJ) -> (i+dI,j-dJ)) $ fInteriorKLs i j
+{-# INLINE fInteriorLoop #-}
+
+-- | Backtrack the unpaired loop region
+
+btInteriorLoop
+  :: Params
+  -> Primary
+  -> NInteLoop
+  -> NInte
+  -> NBT
+  -> NBT
+btInteriorLoop ps inp (NInteLoop nInteLoop) nInte btIL i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ((k,l),enext) <- VU.toList $ fInteriorLoop ps inp nInte i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btIL k l d'
+  , testD z
+  ] where
+      ehere = nInteLoop!(i,j)
+      n = VU.length inp -1
+
+
+
+-- * Inner part
+
+-- | This opens up an interior loop. For each index (i,j) we minimize over all
+-- possible basepair types.
+
+fInteriorInner :: BaseF (EStem -> Features (VU.Vector (ExtPairIdx,Double)))
+fInteriorInner Params{..} inp (EStem eStem) i j
+  | j-i<2 = VU.empty
+  | i==0 || j+1==VU.length inp = VU.empty
+  | otherwise = VU.map f kls
+  where
+    f ijExt@((i,j),(ctIJ,eI,eJ)) =
+      ( ijExt
+      , eStem ! ijExt
+      + interiorClose ! (((nJ,nI),(ctIJ,eJ,eI)),nJp1,nIm1)
+      ) where nI = inp VU.! i
+              nJ = inp VU.! j
+              nIm1 = inp VU.! (i-1)
+              nJp1 = inp VU.! (j+1)
+    kls = VU.fromList [ ((i,j),(ctIJ,eI,eJ))
+                      | eI<-wsh, eJ<-wsh, ctIJ<-citr
+                      ]
+{-# INLINE fInteriorInner #-}
+
+-- | Backtrack from an NInte result to the corresponding EStem parts
+
+btInteriorInner
+  :: Params
+  -> Primary
+  -> NInte
+  -> EStem
+  -> ExtBT
+  -> NBT
+btInteriorInner ps inp (NInte nInte) eStem btES i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ((_,(eI,eJ,ct)),enext) <- VU.toList $ fInteriorInner ps inp eStem i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btES i j eI eJ ct d'
+  , testD z
+  ] where
+      n = VU.length inp -1
+      ehere = nInte!(i,j)
+
+
+
+-- * Helper functions
+
+-- | Since backtracking interior loops is mostly selfcontained, we encapsulate
+-- the above three functions -- which we can't do easily with the forward
+-- calculations as they actually have to save on runtime.
+
+{-
+btInteriorComplete
+  :: Params
+  -> Primary
+  -> EStem
+  -> NInteLoop
+  -> NInte
+  -> ExtBT
+  -> ExtBT
+btInteriorComplete ps pri eStem nInteLoop nInte btExtStem i j ct eI eJ d =
+  btInteriorOuter ps pri eStem btiloop i j ct eI eJ d
+  where btiloop i j d = btInteriorLoop  ps pri nInteLoop nInte btinner i j d
+        btinner = btInteriorInner ps pri nInte     eStem btExtStem
+-}
+
+-- | Given the outer indices (i,j), produces delta_i and delta_j so that
+-- i+delta_i and j-delta_j are the inner indices. 'fInteriorKLs' should fuse
+-- and should make sure that l-k>=4 is always true (maxd). Furthermore the
+-- maximal unpaired length of both sides combined is determined by 'maxLength'.
+--
+-- TODO better name than 'maxLength'
+
+fInteriorKLs :: Int -> Int -> VU.Vector (Int,Int)
+fInteriorKLs i j = didjs where
+  didjs = VU.unfoldr mkDs (4,2)
+  mkDs (d,s)
+    | d>maxd = Nothing
+    | s>=d-2 = Just ((d-s,s),(d+1,2))
+    | otherwise = Just ((d-s,s),(d,s+1))
+  {-# INLINE mkDs #-}
+  maxd = min maxLength (j-i-4)
+{-# INLINE fInteriorKLs #-}
+
diff --git a/BioInf/RNAwolf/Multibranched.hs b/BioInf/RNAwolf/Multibranched.hs
new file mode 100644
--- /dev/null
+++ b/BioInf/RNAwolf/Multibranched.hs
@@ -0,0 +1,350 @@
+{-# LANGUAGE RecordWildCards #-}
+
+-- | Functions for handling non-triplet multibranched loops.
+--
+-- TODO We can do the loop-splitting thing again to speed up multibranched
+-- closing by x10.
+
+module BioInf.RNAwolf.Multibranched where
+
+import qualified Data.Vector.Unboxed as VU
+
+import Biobase.Primary
+import Biobase.Secondary
+import Data.PrimitiveArray
+import Data.PrimitiveArray.Ix
+
+import BioInf.Params
+import BioInf.RNAwolf.Types
+
+
+
+-- * An unpaired nucleotide to the right of an NMbr structure.
+
+-- | Energy for having the rightmost nucleotide (at j) unpaired in NMBr.
+
+fUnpairedRight :: BaseF (NMbr -> Features (VU.Vector (PairIdx,Double)))
+fUnpairedRight Params{..} inp (NMbr nMbr) i j
+  | i<0 || j>n = error $ "Multibranched.fUnpairedRight: " ++ show (i,j)
+  | i==j       = VU.empty
+  | otherwise  = VU.singleton s
+  where
+    s = ( (i,j-1)
+        , nMbr !(i,j-1)
+        )
+    n = VU.length inp -1
+
+-- | Backtrack in NMbr if the nucleotide at j is unpaired.
+
+btUnpairedRight
+  :: Params
+  -> Primary
+  -> NMbr
+  -> NBT
+  -> NBT
+btUnpairedRight ps inp nMbr btM i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ((_,k),enext) <- VU.toList $ fUnpairedRight ps inp nMbr i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btM i k d'
+  , testD z
+  ] where
+      ehere = unNMbr nMbr !(i,j)
+      n = VU.length inp -1
+
+
+-- * An unpaired nucleotide to the right of an NMbr1 structure.
+
+-- | Energy for having the rightmost nucleotide (at j) unpaired in NMBr1.
+
+fUnpairedRight1 :: BaseF (NMbr1 -> Features (VU.Vector (PairIdx,Double)))
+fUnpairedRight1 Params{..} inp (NMbr1 nMbr1) i j
+  | i<0 || j>n = error $ "Multibranched.fUnpairedRight: " ++ show (i,j)
+  | i==j       = VU.empty
+  | otherwise  = VU.singleton s
+  where
+    s = ( (i,j-1)
+        , nMbr1 !(i,j-1)
+        )
+    n = VU.length inp -1
+
+-- | Backtrack NMbr1 if the nucleotide at j is unpaired.
+
+btUnpairedRight1
+  :: Params
+  -> Primary
+  -> NMbr1
+  -> NBT
+  -> NBT
+btUnpairedRight1 ps inp nMbr1 btM1 i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ((_,k),enext) <- VU.toList $ fUnpairedRight1 ps inp nMbr1 i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btM1 i k d'
+  , testD z
+  ] where
+      ehere = unNMbr1 nMbr1 !(i,j)
+      n = VU.length inp -1
+
+
+
+-- * EStem to Helix start
+
+-- | A multibranched helix (except the closing one). (i,j) are closed by a
+-- basepair. Backtracking into the EStem reveals the type of pairing.
+
+fMlHelix :: BaseF (EStem -> Features (VU.Vector (ExtPairIdx,Double)))
+fMlHelix Params{..} inp (EStem eStem) i j
+  | i==0 || j==VU.length inp -1 = VU.empty -- TODO not required ?!
+  | otherwise = VU.map f exts
+  where
+    f ext@(ct,eI,eJ) =
+      ( ijExt
+      , eStem ! ijExt
+      + mbClose ! (((nJ,nI),(ct,eJ,eI)),inp VU.! (j+1) ,inp VU.! (i-1))
+      + multiHelix
+      ) where
+          ijExt = ((i,j),(ct,eI,eJ))
+          nI = inp VU.! i
+          nJ = inp VU.! j
+    exts = VU.fromList [ (ct,eI,eJ)
+                       | j-i>=2, i>0, j+1<VU.length inp
+                       , eI<-wsh, eJ<-wsh, ct<-citr
+                       ]
+{-# INLINE fMlHelix #-}
+
+-- | Backtracks from (i,j) in NMult into the extended-pairing EStem.
+
+btMlHelix
+  :: Params
+  -> Primary
+  -> NMult
+  -> EStem
+  -> ExtBT
+  -> NBT
+btMlHelix ps inp (NMult nMult) eStem btES i j d =
+  [ (x,z)
+  | i>0,i<j,j<n -- correct boundaries since we access elements at (i-1) and (j+1)
+  , ((_,(eI,eJ,ct)),enext) <- VU.toList $ fMlHelix ps inp eStem i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btES i j eI eJ ct d'
+  , testD z
+  ] where
+      ehere = nMult !(i,j)
+      n = VU.length inp -1
+
+
+
+-- * Closes a multibranched loop and produces an extended structure at
+-- ((i,j),ext)
+
+-- | Closes a multibranch loop.
+--
+-- TODO make completely triplet compliant
+
+fMlClose :: BaseF (NMultLoop -> ExtFeatures (VU.Vector (PairIdx,Double)))
+fMlClose Params{..} inp (NMultLoop nMultLoop) i j ct eI eJ
+  | i>=j      = VU.empty
+  | otherwise = VU.singleton s
+  where
+    s = ( (i,j)
+        , nMultLoop !(i,j)
+        + mlc
+        )
+    mlc = 0
+        + multiBranched
+        + multiHelix
+        + mbClose ! ( ((inp VU.! i, inp VU.! j),(ct,eI,eJ))
+                    , inp VU.! (i+1)
+                    , inp VU.! (j-1)
+                    )
+  --      + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0
+{-# INLINE fMlClose #-}
+
+-- | Backtrack from and extended annotation (ij,ext) into the helper table
+-- NMultLoop.
+
+btMlClose
+  :: Params
+  -> Primary
+  -> EStem
+  -> NMultLoop
+  -> NBT
+  -> ExtBT
+btMlClose ps inp (EStem eStem) nMultLoop btMultLoop i j ct eI eJ d =
+  [ (ij:x,z)
+  | i>=0,i<j,j<=n
+  , (_,enext) <- VU.toList $ fMlClose ps inp nMultLoop i j ct eI eJ
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btMultLoop i j d'
+  , testD z
+  ] where
+      ij = ((i,j),(ct,eI,eJ))
+      ehere = eStem!ij
+      n = VU.length inp -1
+
+
+
+-- * Multibranched loop helper table
+
+-- | Multibranched loop helper function that combines "at least one stem" with
+-- "exactly one stem" but does not add the closing energy from (i,j).
+
+fMlLoop :: BaseF (NMbr -> NMbr1 -> Features (VU.Vector (Int,Double)))
+fMlLoop Params{..} inp (NMbr nMbr) (NMbr1 nMbr1) i j = VU.map f ks where
+  f k = ( k
+        , nMbr!(i+1,k)
+        + nMbr1!(k+1,j-1)
+        )
+  ks  = VU.enumFromN (i+2) (j-i-3) -- == [i+2 .. j-2]
+{-# INLINE fMlLoop #-}
+
+-- | Backtracking the multibranched loop.
+
+btMlLoop
+  :: Params
+  -> Primary
+  -> NMultLoop
+  -> NMbr
+  -> NMbr1
+  -> NBT
+  -> NBT
+  -> NBT
+btMlLoop ps inp (NMultLoop nMultLoop) nMbr nMbr1 btM btM1 i j d =
+  [ (x++y,z)
+  | i>=0,i<j,j<=n
+  , (k,enext) <- VU.toList $ fMlLoop ps inp nMbr nMbr1 i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , i+1<k
+  , (x,z') <- btM (i+1) k d'
+  , testD z'
+  , k+1<j-1
+  , (y,z)  <- btM1 (k+1) (j-1) z'
+  , testD z
+  ] where
+      ehere = nMultLoop !(i,j)
+      n = VU.length inp -1
+
+
+
+-- * Backtracking of a multibranched stem with unpaired nucleotides to the
+-- left.
+
+-- | Backtrack a single stem in NMbr, where the stem has zero or more unpaired
+-- nucleotides to the left.
+
+fMlStem :: BaseF (NMult -> Features (VU.Vector (Int,Double)))
+fMlStem Params{..} inp (NMult nMult) i j = VU.map f ks where
+  f k = ( k
+        , nMult!(k,j)
+        )
+  ks  = VU.enumFromN i (j-i-1) -- == [i..j-2]
+{-# INLINE fMlStem #-}
+
+-- | Backtrack by trying to find a multilooped helix.
+
+btMlStem
+  :: Params
+  -> Primary
+  -> NMbr
+  -> NMult
+  -> NBT
+  -> NBT
+btMlStem ps inp (NMbr nMbr) nMult btMH i j d =
+  [ (x,z) -- stem at (k,j)
+  | i>=0,i<j,j<=n
+  , (k,enext) <- VU.toList $ fMlStem ps inp nMult i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btMH k j d'
+  , testD z
+  ] where
+      ehere = nMbr!(i,j)
+      n = VU.length inp -1
+
+
+
+-- * Backtracking of at least two stems by finding one or more stems to the
+-- left and exactly one stem to the right.
+
+-- | Add a stem to a multibranch table containing already at least one stem.
+
+fMlStems :: BaseF (NMbr -> NMult -> Features (VU.Vector (Int,Double)))
+fMlStems Params{..} inp (NMbr nMbr) (NMult nMult) i j = VU.map f ks where
+  f k = ( k
+        , nMbr!(i,k)
+        + nMult!(k+1,j)
+        )
+  ks  = VU.enumFromN (i+2) (j-i-4) -- == [i+2..j-3]
+{-# INLINE fMlStems #-}
+
+-- | Backtrack by finding the splitting index between an NMbr composite
+-- structure and a single multibranched stem NMult (which can contain unpaired
+-- nucleotides to the left).
+
+btMlStems
+  :: Params
+  -> Primary
+  -> NMbr
+  -> NMult
+  -> NBT
+  -> NBT
+  -> NBT
+btMlStems ps inp nMbr nMult btM btMH i j d =
+  [ (x++y,z) -- nMbr ++ nStem
+  | i>=0,i<j,j<=n
+  , (k,enext) <- VU.toList $ fMlStems ps inp nMbr nMult i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z') <- btM i k d'
+  , testD z'
+  , (y,z) <- btMH (k+1) j z'
+  , testD z
+  ] where
+      ehere = unNMbr nMbr !(i,j)
+      n = VU.length inp -1
+
+
+
+-- * Backtrack a single stem in NMbr1. This stem is closed at (i,j).
+
+-- | Add a single stem to a multibranch table containing zero stems already.
+--
+-- TODO this would be equal to mlHelix, unify!
+
+fMl1Stem :: BaseF (NMult -> Features (VU.Vector ((Int,Int),Double)))
+fMl1Stem Params{..} inp (NMult nMult) i j = VU.singleton s where
+  s = ( (i,j)
+      , nMult!(i,j)
+      )
+{-# INLINE fMl1Stem #-}
+
+-- | Backtrack a single stem closed at (i,j) for NMbr1. Takes the route through
+-- NMult which solves for the exact pairtype.
+
+btMl1Stem
+  :: Params
+  -> Primary
+  -> NMbr1
+  -> NMult
+  -> NBT
+  -> NBT
+btMl1Stem ps inp (NMbr1 nMbr1) nMult btMH i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ((k,l),enext) <- VU.toList $ fMl1Stem ps inp nMult i j
+  , i==k && j==l
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btMH i j d'
+  , testD z
+  ] where
+      ehere = nMbr1!(i,j)
+      n = VU.length inp -1
diff --git a/BioInf/RNAwolf/Stem.hs b/BioInf/RNAwolf/Stem.hs
new file mode 100644
--- /dev/null
+++ b/BioInf/RNAwolf/Stem.hs
@@ -0,0 +1,101 @@
+
+{-# LANGUAGE RecordWildCards #-}
+
+module BioInf.RNAwolf.Stem where
+
+import qualified Data.Vector.Unboxed as VU
+
+import Biobase.Primary
+import Biobase.Secondary
+import Data.PrimitiveArray
+import Data.PrimitiveArray.Ix
+
+import BioInf.Params
+import BioInf.RNAwolf.Types
+
+
+
+-- * Collapses an extended stem into a normal stem, or allows going back from a
+-- normal stem to an extended stem.
+
+-- | A normal stem is created by taking the minimum over all possible basepairs
+-- of the extended stem.
+
+fNstem :: BaseF (EStem -> Features (VU.Vector (ExtPairIdx,Double)))
+fNstem Params{..} inp (EStem eStem) i j
+  | i<0 || j>n = error $ "Stem.fNstem: " ++ show (i,j)
+  | j-i<=2      = VU.empty
+  | otherwise   = VU.map f $ VU.fromList [ (ct,eI,eJ) | ct<-citr, eI<-wsh, eJ<-wsh ]
+  where
+    f ext = ( idx
+            , eStem !idx
+            ) where idx = ((i,j),ext)
+    n = VU.length inp -1
+{-# INLINE fNstem #-}
+
+-- | Backtrack from a normal stem back into the extended stem.
+
+btNstem
+  :: Params
+  -> Primary
+  -> NStem
+  -> EStem
+  -> ExtBT
+  -> NBT
+btNstem ps inp nStem eStem btES i j d =
+  [ (x,z)
+  | i>=0,i<j,j<=n
+  , ct <- citr, eI <- wsh, eJ <- wsh
+  , ((_,(ct,eI,eJ)),enext) <- VU.toList $ fNstem ps inp eStem i j
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btES i j ct eI eJ d'
+  , testD z
+  ] where
+     ehere = unNStem nStem !(i,j)
+     n = VU.length inp -1
+
+
+
+-- * Stacking with extended basepair information
+
+-- | A stem is extended by another pair. The score contribution is dependent on
+-- the previous pair. Note that for score lookup purposes, the inner pair is
+-- switched.
+
+fStem :: BaseF (EStem -> ExtFeatures (VU.Vector (ExtPairIdx,Double)))
+fStem Params{..} inp (EStem eStem) i j ct eI eJ
+  | j-i<3     = VU.empty
+  | otherwise = VU.map f $ eKLs
+  where
+    f (ctKL,eK,eL) = ( ((k,l),(ctKL,eK,eL))
+                     , eStem ! ((k,l),(ctKL,eK,eL))
+                     + stem ! (ijExt,((inp VU.! l,inp VU.! k),(ctKL,eL,eK)))
+--                     + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0
+                     )
+    ijExt = ((inp VU.! i, inp VU.! j),(ct,eI,eJ))
+    k = i+1
+    l = j-1
+    eKLs = VU.fromList [ (ctKL,eK,eL) | eK<-wsh, eL<-wsh, ctKL<-citr ]
+    {-# NOINLINE eKLs #-} -- speed-up for small input sizes
+{-# INLINE fStem #-}
+
+-- | Stem backtracking.
+
+btStem
+  :: Params
+  -> Primary
+  -> EStem
+  -> ExtBT -- recursive backtracking function
+  -> ExtBT
+btStem ps inp eStem btES i j ct eI eJ d =
+  [ (ij:x,z)
+  | i>=0,i<j,j<=n
+  , (((k,l),(ctKL,eK,eL)),enext) <- VU.toList $ fStem ps inp eStem i j ct eI eJ
+  , let d' = newD d ehere enext
+  , testD d'
+  , (x,z) <- btES k l ctKL eK eL d'
+  ] where
+      ij = ((i,j),(ct,eI,eJ))
+      ehere = unEStem eStem !ij
+      n = VU.length inp -1
diff --git a/BioInf/RNAwolf/TripletBulge.hs b/BioInf/RNAwolf/TripletBulge.hs
new file mode 100644
--- /dev/null
+++ b/BioInf/RNAwolf/TripletBulge.hs
@@ -0,0 +1,17 @@
+
+{-# LANGUAGE RecordWildCards #-}
+
+module BioInf.RNAwolf.TripletBulge where
+
+import qualified Data.Vector.Unboxed as VU
+
+import Biobase.Primary
+import Biobase.Secondary
+import Data.PrimitiveArray
+import Data.PrimitiveArray.Ix
+
+import BioInf.Params
+import BioInf.RNAwolf.Types
+
+
+
diff --git a/BioInf/RNAwolf/TripletStem.hs b/BioInf/RNAwolf/TripletStem.hs
new file mode 100644
--- /dev/null
+++ b/BioInf/RNAwolf/TripletStem.hs
@@ -0,0 +1,64 @@
+
+{-# LANGUAGE RecordWildCards #-}
+
+module BioInf.RNAwolf.TripletStem where
+
+import qualified Data.Vector.Unboxed as VU
+
+import Biobase.Primary
+import Biobase.Secondary
+import Data.PrimitiveArray
+import Data.PrimitiveArray.Ix
+
+import BioInf.Params
+import BioInf.RNAwolf.Types
+
+
+{-
+      let vTStem    = 999999 -- minimumVU $ fTripletStem ps inp eStem      i j ct eI eJ
+-}
+
+{-
+-- | Triplet stems have a shared nucleotide.
+--   _C
+-- A-_B where A is paired with both B and C.
+--
+-- If A,B is paired with cWW, then A,C can only use the S or H edge.
+
+fTripletStem :: BaseF (ExtTable -> ExtFeatures (VU.Vector (ExtPairIdx,Double)))
+fTripletStem Params{..} inp eStem i j ct eI eJ
+  | j-i<2     = VU.empty
+  | otherwise = {- VU.map fI iShared VU.++ -} VU.map fJ jShared
+  where
+    fI klExt@((k,l),(ctKL,eK,eL))
+      = ( klExt
+        , eStem ! klExt
+--        + stemTriplet ! ( ((nI,nL),(ctKL,eK,eL)), ((nI,nJ),(ct,eI,eJ)) )
+--        + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0
+        ) where nL = inp VU.! l
+    fJ klExt@((k,l),(ctKL,eK,eL))
+      = ( klExt
+        , eStem ! klExt
+        + stemTriplet ! ( ((nI,nJ),((ct,eI,eJ))), ((nJ,nK),(ctKL,eJ,eK)) ) --l==j
+--        + if j-i-1<=maxDistance then pairDistance ! (j-i-1) else 0
+        ) where nK = inp VU.! k
+    nI = inp VU.! i
+    nJ = inp VU.! j
+    iShared = VU.fromList [ ((k,l),(du,eK,eL))
+                          | let k=i
+                          , let l=j-1
+                          , eK<-[wc,sugar,hoogsteen]
+                          , eK/=eI
+                          , eL<-[wc,sugar,hoogsteen]
+                          , du<-[cis,trans]
+                          ]
+    jShared = VU.fromList [ ((k,l),(du,eK,eL))
+                          | let k=i+1
+                          , let l=j
+                          , eK<-[wc,sugar,hoogsteen]
+                          , eL<-[wc,sugar,hoogsteen]
+                          , du<-[cis,trans]
+                          , eL/=eJ
+                          ]
+{- INLINE fTripletStem #-}
+-}
diff --git a/BioInf/RNAwolf/Types.hs b/BioInf/RNAwolf/Types.hs
new file mode 100644
--- /dev/null
+++ b/BioInf/RNAwolf/Types.hs
@@ -0,0 +1,65 @@
+
+module BioInf.RNAwolf.Types where
+
+import Biobase.Primary
+import Biobase.Secondary
+import Data.PrimitiveArray
+
+import BioInf.Params
+
+
+
+newD d here next = d - (next - here)
+testD d = d>=0
+
+-- | Should really go into BiobaseXNA
+
+wsh = [wc,sugar,hoogsteen]
+
+-- | Should really go into BiobaseXNA
+
+citr = [cis,trans]
+
+
+
+
+type ExtBT = Int -> Int -> CTisomerism -> Edge -> Edge -> Double -> BTAnswer
+type NBT = Int -> Int -> Double -> BTAnswer
+
+type BTAnswer = [([ExtPairIdx],Double)]
+
+type Table    = PrimArray PairIdx    Double
+type ExtTable = PrimArray ExtPairIdx Double
+
+type BaseF a = Params -> Primary -> a
+type ExtFeatures a = Int -> Int -> CTisomerism -> Edge -> Edge -> a
+type Features a = Int -> Int -> a
+
+type Tables = ( EStem
+              , NStem
+              , NInte
+              , NInteLoop
+              , NBulg
+              , NBulgLoop
+              , NMult
+              , NMbr
+              , NMbr1
+              , NMultLoop
+              , NExtn
+              )
+
+-- ** Newtype wrappers for all tables.
+--
+-- NOTE Don't ever not newtype-wrap or you will hurt your brain.
+
+newtype EStem = EStem {unEStem :: ExtTable}
+newtype NStem = NStem {unNStem :: Table}
+newtype NInte = NInte {unNInte :: Table}
+newtype NInteLoop = NInteLoop {unNInteLoop :: Table}
+newtype NMult = NMult {unNMult :: Table}
+newtype NBulg = NBulg {unNBulg :: Table}
+newtype NBulgLoop = NBulgLoop {unBulgLoop :: Table}
+newtype NMbr  = NMbr  {unNMbr :: Table}
+newtype NMbr1 = NMbr1 {unNMbr1 :: Table}
+newtype NExtn = NExtn {unNExtn :: Table}
+newtype NMultLoop = NMultLoop {unMultLoop :: Table}
diff --git a/RNAwolf.cabal b/RNAwolf.cabal
--- a/RNAwolf.cabal
+++ b/RNAwolf.cabal
@@ -1,5 +1,5 @@
 name:           RNAwolf
-version:        0.3.0.0
+version:        0.3.0.1
 author:         Christian Hoener zu Siederdissen, Stephan H Bernhart, Peter F Stadler, Ivo L Hofacker
 copyright:      Christian Hoener zu Siederdissen, 2010-2011
 homepage:       http://www.tbi.univie.ac.at/software/rnawolf/
@@ -73,6 +73,16 @@
     BioInf.Params.Import
     BioInf.PassiveAggressive
     BioInf.RNAwolf
+    BioInf.RNAwolf.Bulge
+    BioInf.RNAwolf.Extern
+    BioInf.RNAwolf.Hairpin
+    BioInf.RNAwolf.Interior
+    BioInf.RNAwolf.Multibranched
+    BioInf.RNAwolf.Stem
+    BioInf.RNAwolf.TripletBulge
+    BioInf.RNAwolf.TripletStem
+    BioInf.RNAwolf.Types
+
   ghc-options:
     -O2
   if flag(llvm)