<commit_before>{-# LANGUAGE OverloadedStrings #-}

import Test.HUnit
import Data.Binary
import qualified Data.Text as T

import TestFileUtils
import CladeModel
import PepModel
import Alignment

-- Numeric testing of PepModels:

small_prob = 0.0001 :: Double -- need to be specific, else logBase complains
scale_factor = 1000 :: Double

test_aln1 = [
    AlnRow "my-OTU" "VIKD-" 1,
    AlnRow "my-OTU" "VARE-" 1,
    AlnRow "my-OTU" "VARNX" 1,
    AlnRow "my-OTU" "VIK--" 1,
    AlnRow "my-OTU" "VIHQY" 1
    ]

aln1Mod = alnToPepModel small_prob scale_factor "my-OTU" test_aln1 

-- These tests are calques from TestNucModel (this is also why the test numbers
-- are discontinuous), but while for nucleotides it is practical to test all
-- five possibilities (ACTG-), for the 20 amino acids we won't go all that way.
-- We test the residues that are found at a given position, plus perhaps one (or
-- the gap) that is not.
--
-- There are 5 out of 5 'V's in the first column, hence 5/5 ... 1
test_1 = "A@1(aln1)" ~: round (scale_factor * (logBase 10 (5/5)))
    ~?= pepScoreOf aln1Mod 'V' 1
-- There are 0 out of 5 'C's in the first column: use small_prob
test_2 = "C@1(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= pepScoreOf aln1Mod 'C' 1
-- 'D' is for 'dash'
test_5 = "D@1(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= (pepScoreOf aln1Mod '-' 1)

test_6  = "A@2(aln1)" ~: (round (scale_factor * (logBase 10 (2/5))))
    ~?= (pepScoreOf aln1Mod 'A' 2)
test_8  = "G@2(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= (pepScoreOf aln1Mod 'G' 2)
test_9  = "T@2(aln1)" ~: (round (scale_factor * (logBase 10 (3/5))))
    ~?= (pepScoreOf aln1Mod 'I' 2)

test_11 = "A@3(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'H' 3)
test_12 = "C@3(aln1)" ~: (round (scale_factor * (logBase 10 (2/5))))
    ~?= (pepScoreOf aln1Mod 'R' 3)
test_13 = "G@3(aln1)" ~: (round (scale_factor * (logBase 10 (2/5))))
    ~?= (pepScoreOf aln1Mod 'K' 3)
test_14 = "T@3(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= (pepScoreOf aln1Mod 'T' 3)

test_16 = "A@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'D' 4)
test_17 = "C@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'E' 4)
test_18 = "G@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'N' 4)
test_19 = "T@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'Q' 4)
test_20 = "D@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod '-' 4)

test_21 = "A@5(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= (pepScoreOf aln1Mod 'A' 5)
-- we _always_ divide by the number of sequences in the aln, even when non-{ATGC-} residues occur (like N)
test_24 = "T@5(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'Y' 5)
test_25 = "D@5(aln1)" ~: (round (scale_factor * (logBase 10 (3/5))))
    ~?= (pepScoreOf aln1Mod '-' 5)

-- Test Binary functions, i.e. storage as binary file to disk, and reading from
-- it, then check every single value in the model, just as above. These cases
-- are all under the same test, because I don't want to write and read the file
-- every time.

test_26 = TestCase
    (do
        removeIfExists "aln1Mod.bmod"
        encodeFile "aln1Mod.bmod" aln1Mod
        aln2Mod <- decodeFile "aln1Mod.bmod"
        assertEqual "store-read" aln1Mod aln2Mod
        round (scale_factor * (logBase 10 (5/5)))
            @=? pepScoreOf aln1Mod 'V' 1
        -- There are 0 out of 5 'C's in the first column: use small_prob
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod 'C' 1
        -- 'D' is for 'dash'
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod '-' 1
        
        round (scale_factor * (logBase 10 (2/5)))
            @=? pepScoreOf aln1Mod 'A' 2
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod 'G' 2
        round (scale_factor * (logBase 10 (3/5)))
            @=? pepScoreOf aln1Mod 'I' 2
        
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'H' 3
        round (scale_factor * (logBase 10 (2/5)))
            @=? pepScoreOf aln1Mod 'R' 3
        round (scale_factor * (logBase 10 (2/5)))
            @=? pepScoreOf aln1Mod 'K' 3
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod 'T' 3
        
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'D' 4
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'E' 4
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'N' 4
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'Q' 4
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod '-' 4
        
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod 'A' 5
        -- we _always_ divide by the number of sequences in the aln, even when non-{ATGC-} residues occur (like N)
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'Y' 5
        round (scale_factor * (logBase 10 (3/5)))
            @=? pepScoreOf aln1Mod '-' 5
        

        )

-- Computes the expected score from the list of positional frequencies. Short
-- name to save some space.
--
es ps = sum $ map (round . (scale_factor *) . (logBase 10)) ps

test_27 = "AAAAA" ~: (pepScoreSeq aln1Mod "AAAAA") ~?= (es [5/5, 2/5, 1/5, 1/5, small_prob])
test_28 = "-----" ~: (pepScoreSeq aln1Mod "-----") ~?= (es [small_prob, small_prob, small_prob, 1/5, 3/5])
test_29 = "ATCG-" ~: (pepScoreSeq aln1Mod "ATCG-") ~?= (es [5/5, 3/5, 2/5, 1/5, 3/5])
test_30 = "ATCGN" ~: (pepScoreSeq aln1Mod "ATCGN") ~?= (es [5/5, 3/5, 2/5, 1/5, small_prob])
                             
-- Test model length

test_31 = "modLength" ~: (pepModLength aln1Mod) ~?= 5

-- Tests for empty alignments
--
aln2Mod = alnToPepModel small_prob scale_factor "anonymous" []<commit_msg>added numerical tests of PepModel<commit_after>{-# LANGUAGE OverloadedStrings #-}

import Test.HUnit
import Data.Binary
import qualified Data.Text as T

import TestFileUtils
import CladeModel
import PepModel
import Alignment

-- Numeric testing of PepModels:

small_prob = 0.0001 :: Double -- need to be specific, else logBase complains
scale_factor = 1000 :: Double

test_aln1 = [
    AlnRow "my-OTU" "VIKD-" 1,
    AlnRow "my-OTU" "VARE-" 1,
    AlnRow "my-OTU" "VARNX" 1,
    AlnRow "my-OTU" "VIK--" 1,
    AlnRow "my-OTU" "VIHQY" 1
    ]

aln1Mod = alnToPepModel small_prob scale_factor "my-OTU" test_aln1 

-- These tests are calques from TestNucModel (this is also why the test numbers
-- are discontinuous), but while for nucleotides it is practical to test all
-- five possibilities (ACTG-), for the 20 amino acids we won't go all that way.
-- We test the residues that are found at a given position, plus perhaps one (or
-- the gap) that is not.
--
-- There are 5 out of 5 'V's in the first column, hence 5/5 ... 1
test_1 = "A@1(aln1)" ~: round (scale_factor * (logBase 10 (5/5)))
    ~?= pepScoreOf aln1Mod 'V' 1
-- There are 0 out of 5 'C's in the first column: use small_prob
test_2 = "C@1(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= pepScoreOf aln1Mod 'C' 1
-- 'D' is for 'dash'
test_5 = "D@1(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= (pepScoreOf aln1Mod '-' 1)

test_6  = "A@2(aln1)" ~: (round (scale_factor * (logBase 10 (2/5))))
    ~?= (pepScoreOf aln1Mod 'A' 2)
test_8  = "G@2(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= (pepScoreOf aln1Mod 'G' 2)
test_9  = "T@2(aln1)" ~: (round (scale_factor * (logBase 10 (3/5))))
    ~?= (pepScoreOf aln1Mod 'I' 2)

test_11 = "A@3(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'H' 3)
test_12 = "C@3(aln1)" ~: (round (scale_factor * (logBase 10 (2/5))))
    ~?= (pepScoreOf aln1Mod 'R' 3)
test_13 = "G@3(aln1)" ~: (round (scale_factor * (logBase 10 (2/5))))
    ~?= (pepScoreOf aln1Mod 'K' 3)
test_14 = "T@3(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= (pepScoreOf aln1Mod 'T' 3)

test_16 = "A@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'D' 4)
test_17 = "C@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'E' 4)
test_18 = "G@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'N' 4)
test_19 = "T@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'Q' 4)
test_20 = "D@4(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod '-' 4)

test_21 = "A@5(aln1)" ~: (round (scale_factor * (logBase 10 small_prob)))
    ~?= (pepScoreOf aln1Mod 'A' 5)
-- we _always_ divide by the number of sequences in the aln, even when non-{ATGC-} residues occur (like N)
test_24 = "T@5(aln1)" ~: (round (scale_factor * (logBase 10 (1/5))))
    ~?= (pepScoreOf aln1Mod 'Y' 5)
test_25 = "D@5(aln1)" ~: (round (scale_factor * (logBase 10 (3/5))))
    ~?= (pepScoreOf aln1Mod '-' 5)

-- Test Binary functions, i.e. storage as binary file to disk, and reading from
-- it, then check every single value in the model, just as above. These cases
-- are all under the same test, because I don't want to write and read the file
-- every time.

test_26 = TestCase
    (do
        removeIfExists "aln1Mod.bmod"
        encodeFile "aln1Mod.bmod" aln1Mod
        aln2Mod <- decodeFile "aln1Mod.bmod"
        assertEqual "store-read" aln1Mod aln2Mod
        round (scale_factor * (logBase 10 (5/5)))
            @=? pepScoreOf aln1Mod 'V' 1
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod 'C' 1
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod '-' 1
        
        round (scale_factor * (logBase 10 (2/5)))
            @=? pepScoreOf aln1Mod 'A' 2
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod 'G' 2
        round (scale_factor * (logBase 10 (3/5)))
            @=? pepScoreOf aln1Mod 'I' 2
        
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'H' 3
        round (scale_factor * (logBase 10 (2/5)))
            @=? pepScoreOf aln1Mod 'R' 3
        round (scale_factor * (logBase 10 (2/5)))
            @=? pepScoreOf aln1Mod 'K' 3
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod 'T' 3
        
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'D' 4
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'E' 4
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'N' 4
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'Q' 4
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod '-' 4
        
        round (scale_factor * (logBase 10 small_prob))
            @=? pepScoreOf aln1Mod 'A' 5
        round (scale_factor * (logBase 10 (1/5)))
            @=? pepScoreOf aln1Mod 'Y' 5
        round (scale_factor * (logBase 10 (3/5)))
            @=? pepScoreOf aln1Mod '-' 5
        )

-- Computes the expected score from the list of positional frequencies. Short
-- name to save some space.
--
es ps = sum $ map (round . (scale_factor *) . (logBase 10)) ps

test_27 = "AAAAA" ~: (pepScoreSeq aln1Mod "AAAAA") ~?= (es [small_prob, 2/5, small_prob, small_prob, small_prob])
test_28 = "-----" ~: (pepScoreSeq aln1Mod "-----") ~?= (es [small_prob, small_prob, small_prob, 1/5, 3/5])
test_29 = "ATCG-" ~: (pepScoreSeq aln1Mod "ATCG-") ~?= (es [5/5, 3/5, 2/5, 1/5, 3/5])
test_30 = "ATCGN" ~: (pepScoreSeq aln1Mod "ATCGN") ~?= (es [5/5, 3/5, 2/5, 1/5, small_prob])
                             
-- Test model length

test_31 = "modLength" ~: (pepModLength aln1Mod) ~?= 5

-- Tests for empty alignments
--
aln2Mod = alnToPepModel small_prob scale_factor "anonymous" []
{'fraction_non_alphanumeric': 0.1568451896467524, 'fraction_numerical': 0.055998697704704545, 'mean_word_length': 2.8134698944754812, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2931460', 'n_tokens_mistral': 5221, 'n_tokens_neox': 4837, 'n_words': 1379}