gringimadre

module F2 where
import Data.List

dnaLetter = "ACGT"
data SeqTypes = Dna | Protein deriving (Show, Eq)
data MolSeq = MolSeq String String SeqTypes deriving (Show)

string2seq :: String -> String -> MolSeq
string2seq namn sekvens = MolSeq namn sekvens (string2seqHelper namn sekvens)
string2seqHelper :: String -> String -> SeqTypes
string2seqHelper namn sekvens
   | (null sekvens) = Dna
   | (head sekvens) `elem` dnaLetter = string2seqHelper namn (tail sekvens)
   | otherwise = Protein

seqName :: MolSeq -> String
seqName (MolSeq namn sekvens typ) = namn

seqSequence :: MolSeq -> String
seqSequence (MolSeq namn sekvens typ) = sekvens

seqType :: MolSeq -> SeqTypes
seqType (MolSeq namn sekvens typ) = typ

seqLength :: MolSeq -> Double
seqLength (MolSeq namn sekvens typ) = sum [ 1 | _ <- sekvens] :: Double

getAlpha :: MolSeq -> MolSeq -> Double
getAlpha (MolSeq namn1 sekvens1 typ1) (MolSeq namn2 sekvens2 typ2) = sum [ if x == y then 0 else 1 | (x,y) <- (zip sekvens1 sekvens2)] :: Double
checkValidity :: MolSeq -> MolSeq -> Bool
checkValidity (MolSeq namn1 sekvens1 typ1) (MolSeq namn2 sekvens2 typ2) = if typ1 == typ2 then True else False
eDistanceDNA :: MolSeq -> MolSeq -> Double
eDistanceDNA (MolSeq namn1 sekvens1 typ1) (MolSeq namn2 sekvens2 typ2)
  | alpha > 0.74 = 3.3
  | otherwise = (-3/4)*log(1-(4*alpha/3))
  where alpha = gA/(seqLength (MolSeq namn1 sekvens1 typ1))
        gA = getAlpha (MolSeq namn1 sekvens1 typ1) (MolSeq namn2 sekvens2 typ2)
eDistanceProtein :: MolSeq -> MolSeq -> Double
eDistanceProtein (MolSeq namn1 sekvens1 typ1) (MolSeq namn2 sekvens2 typ2)
  | alpha > 0.94 = 3.7
  | otherwise = (-19/20)*log(1-(20*alpha/19))
  where alpha = gA/(seqLength (MolSeq namn1 sekvens1 typ1))
        gA = getAlpha (MolSeq namn1 sekvens1 typ1) (MolSeq namn2 sekvens2 typ2)
seqDistance :: MolSeq -> MolSeq -> Double
seqDistance molseq1 molseq2 =
  if (checkValidity molseq1 molseq2)
  then (if (seqType molseq1) == Dna then eDistanceDNA molseq1 molseq2 else eDistanceProtein molseq1 molseq2)
  else error "FU BRO"

-------------------UPPGIFT2-----------------------------------
-------------------UPPGIFT2-----------------------------------
-------------------UPPGIFT2-----------------------------------

data Profile = Profile [[Double]] SeqTypes Int String deriving (Show)

molseqs2profile :: String -> [MolSeq] -> Profile
molseqs2profile s mArr = Profile (getM mArr) (getSeqUpp3 mArr) (length mArr) s
-- getM xs = [ getMHelper x | x <- xs ]
-- getMHelper [] = []
-- getMHelper (x:xs) = snd x : getMHelper xs
-- getSeqUpp2 profMatrix = string2seqHelper "notImportnat" [ fst x | x <- (concat profMatrix)]
-- getkukiLength molseqArr = (length (head (makeProfileMatrix molseqArr)))
-- getM :: [MolSeq] -> [[Int]]
getM molseqArr = [ getMHelper (fromIntegral (length molseqArr)) x | x <- xs ]
--getM molseqArr = [ getMHelper x | x <- xs ]
  where xs = (makeProfileMatrix molseqArr)
-- getMHelper :: Int -> [[(Char, Int)]] -> [[Int]]
getMHelper mydiv [] = []
getMHelper mydiv (x:xs) = (fromIntegral (snd x)) / (fromIntegral mydiv) : getMHelper mydiv xs
getSeqUpp3 molseqArr = string2seqHelper "notImportnat" [ fst x | x <- (concat profMatrix)]
  where profMatrix = makeProfileMatrix molseqArr


nucleotides = "ACGT"
aminoacids = sort "ARNDCEQGHILKMFPSTWYVX"
makeProfileMatrix :: [MolSeq] -> [[(Char, Int)]]
makeProfileMatrix [] = error "Empty sequence list"
makeProfileMatrix sl = res
  where
    t = seqType (head sl)
    defaults =
      if (t == Dna) then
        zip nucleotides (replicate (length nucleotides) 0) -- Rad (i)
      else
        zip aminoacids (replicate (length aminoacids) 0)   -- Rad (ii)
    strs = map seqSequence sl                              -- Rad (iii)
    tmp1 = map (map (\x -> ((head x), (length x))) . group . sort)
               (transpose strs)                            -- Rad (iv)
    equalFst a b = (fst a) == (fst b)
    res = map sort (map (\l -> unionBy equalFst l defaults) tmp1)