Algorithmic Information Theory

1
#include <iostream>
2
#include <string>
3
#include <unordered_map>
4
#include <cmath>
5
#include <iomanip>
6
7
// Estimate compressed size (in bits) using a simple educational LZ78 scheme.
8
// This is not an optimized compressor; it is self-contained for teaching.
9
// Encoding model: emit pairs (index_of_phrase, next_char). Index uses
10
// ceil(log2(dict_size + 1)) bits to reference any existing dictionary entry,
11
// and next_char uses 8 bits. We also emit a final pair if a phrase remains.
12
13
static size_t lz78CompressedBits(const std::string &s) {
14
    std::unordered_map<std::string, int> dict; // phrase -> index
15
    dict.reserve(s.size() * 2);
16
    dict.max_load_factor(0.7f);
17
18
    // Index 0 is the empty phrase ""
19
    dict[""] = 0;
20
    std::string w; // current phrase
21
    size_t bits = 0;
22
23
    auto indexBitWidth = [&](int dictSize) -> int {
24
        // Need to address indices in [0..dictSize-1]; use at least 1 bit
25
        if (dictSize <= 1) return 1; // only empty phrase exists
26
        return (int)std::ceil(std::log2((double)dictSize));
27
    };
28
29
    for (char c : s) {
30
        std::string wc = w;
31
        wc.push_back(c);
32
        if (dict.find(wc) != dict.end()) {
33
            // Extend current phrase
34
            w.swap(wc); // w = wc
35
        } else {
36
            // Emit (index(w), c)
37
            int idx = dict[w];
38
            int idxBits = indexBitWidth((int)dict.size());
39
            bits += (size_t)idxBits + 8; // index + next_char
40
            // Add new phrase wc to dictionary with next index
41
            dict[wc] = (int)dict.size();
42
            // Reset current phrase
43
            w.clear();
44
        }
45
    }
46
47
    // If there's a remaining phrase w, emit (index(w), special_end_char)
48
    if (!w.empty()) {
49
        int idxBits = indexBitWidth((int)dict.size());
50
        bits += (size_t)idxBits + 8; // pay a char-sized marker cost
51
    }
52
53
    return bits;
54
}
55
56
static double ncd(const std::string &x, const std::string &y) {
57
    size_t Cx = lz78CompressedBits(x);
58
    size_t Cy = lz78CompressedBits(y);
59
    std::string xy = x + y;
60
    size_t Cxy = lz78CompressedBits(xy);
61
62
    size_t maxC = std::max(Cx, Cy);
63
    size_t minC = std::min(Cx, Cy);
64
    if (maxC == 0) return 0.0; // both empty strings
65
    return (double)(Cxy - minC) / (double)maxC;
66
}
67
68
int main() {
69
    std::string a = std::string(10, 'A') + std::string(10, 'B'); // AAAAAAAAAABBBBBBBBBB
70
    std::string b;
71
    for (int i = 0; i < 20; ++i) b += (i % 2 == 0 ? 'A' : 'B'); // ABABAB... (20)
72
    std::string c = "QWERTYUIOPASDFGHJKLZ"; // looks irregular
73
74
    size_t Ca = lz78CompressedBits(a);
75
    size_t Cb = lz78CompressedBits(b);
76
    size_t Cc = lz78CompressedBits(c);
77
78
    std::cout << "Estimated LZ78 bits (upper bounds on K up to constants):\n";
79
    std::cout << "C(a) = " << Ca << " bits\n";
80
    std::cout << "C(b) = " << Cb << " bits\n";
81
    std::cout << "C(c) = " << Cc << " bits\n\n";
82
83
    std::cout << std::fixed << std::setprecision(4);
84
    std::cout << "NCD(a,b) = " << ncd(a, b) << "\n";
85
    std::cout << "NCD(a,c) = " << ncd(a, c) << "\n";
86
    std::cout << "NCD(b,c) = " << ncd(b, c) << "\n";
87
88
    return 0;
89
}
90

1
#include <iostream>
2
#include <string>
3
#include <cmath>
4
#include <random>
5
#include <iomanip>
6
7
// Compute prequential KT code length (in bits) for a binary string under i.i.d. Bernoulli.
8
// Probability of next 1 at time t is (n1 + 1/2) / (t + 1); similarly for 0.
9
static double ktBernoulliCodeLenBits(const std::string &s) {
10
    int n0 = 0, n1 = 0; // counts so far
11
    double L = 0.0;
12
    for (size_t t = 0; t < s.size(); ++t) {
13
        int bit = (s[t] == '1') ? 1 : 0;
14
        double p = (bit == 1) ? (n1 + 0.5) / (t + 1.0) : (n0 + 0.5) / (t + 1.0);
15
        L += -std::log2(p);
16
        if (bit) ++n1; else ++n0;
17
    }
18
    return L;
19
}
20
21
// Compute prequential KT code length (in bits) for a binary string under a 1st-order Markov model.
22
// We maintain two KT predictors, one for context 0 and one for context 1.
23
static double ktMarkov1CodeLenBits(const std::string &s) {
24
    if (s.empty()) return 0.0;
25
    // Cost of first symbol with KT Bernoulli from scratch: p = 1/2 => 1 bit
26
    double L = 1.0;
27
    int c[2][2] = {{0,0},{0,0}}; // c[prev][next]
28
    int tot[2] = {0,0};          // totals per prev
29
30
    for (size_t i = 1; i < s.size(); ++i) {
31
        int prev = (s[i-1] == '1') ? 1 : 0;
32
        int bit  = (s[i]   == '1') ? 1 : 0;
33
        double p = (c[prev][bit] + 0.5) / (tot[prev] + 1.0);
34
        L += -std::log2(p);
35
        c[prev][bit]++; tot[prev]++;
36
    }
37
    return L;
38
}
39
40
struct MDLResult {
41
    std::string chosen;
42
    double L_bernoulli;
43
    double L_markov1;
44
};
45
46
static MDLResult mdlSelect(const std::string &s) {
47
    const int k_bern = 1; // number of parameters (asymptotically)
48
    const int k_mark = 2; // p(1|0), p(1|1)
49
    int n = (int)s.size();
50
    double L_D_bern = ktBernoulliCodeLenBits(s);
51
    double L_D_mark = ktMarkov1CodeLenBits(s);
52
    // Two-part MDL: L(M) ~ 0.5 * k * log2(n) + const; add 1 bit to identify the model
53
    double L_M_bern = 1.0 + 0.5 * k_bern * (n > 0 ? std::log2((double)n) : 0.0);
54
    double L_M_mark = 1.0 + 0.5 * k_mark * (n > 0 ? std::log2((double)n) : 0.0);
55
56
    double L_bern = L_M_bern + L_D_bern;
57
    double L_mark = L_M_mark + L_D_mark;
58
59
    MDLResult r;
60
    r.L_bernoulli = L_bern;
61
    r.L_markov1 = L_mark;
62
    r.chosen = (L_bern <= L_mark) ? "Bernoulli(iid)" : "Markov(1st-order)";
63
    return r;
64
}
65
66
int main() {
67
    // Example 1: Alternating pattern should favor Markov(1)
68
    std::string alt;
69
    for (int i = 0; i < 200; ++i) alt.push_back((i % 2) ? '1' : '0');
70
    // Example 2: Random-looking pattern should favor Bernoulli
71
    std::mt19937 rng(12345);
72
    std::bernoulli_distribution coin(0.5);
73
    std::string rnd;
74
    for (int i = 0; i < 200; ++i) rnd.push_back(coin(rng) ? '1' : '0');
75
76
    auto R1 = mdlSelect(alt);
77
    auto R2 = mdlSelect(rnd);
78
79
    std::cout << std::fixed << std::setprecision(3);
80
    std::cout << "ALT chosen: " << R1.chosen
81
              << ", L_bern=" << R1.L_bernoulli
82
              << ", L_mark=" << R1.L_markov1 << "\n";
83
    std::cout << "RND chosen: " << R2.chosen
84
              << ", L_bern=" << R2.L_bernoulli
85
              << ", L_mark=" << R2.L_markov1 << "\n";
86
87
    return 0;
88
}
89