Markov-Solver/main.cpp at main · secondplayerIs/Markov-Solver · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#include <iostream>
#include <vector>
#include <valarray>
using namespace std;
#define float_vec vector<float>


vector<float_vec > multiply(vector<float_vec > A, vector<float_vec > B, int N) {
    vector<float_vec > C(N, vector<float>(N, 0));
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
         for (int k = 0; k < N; ++k)
            C[i][j] += A[i][k] * B[k][j];
   return C;
}

vector<float_vec > power(vector<float_vec > M, int p, int n) {
    vector<float_vec > A(n, vector<float>(n, 0));
    for (int i = 0; i < n; ++i)
        A[i][i] = 1;
    while (p) {
      if (p % 2)
         A = multiply(A, M, n);
      M = multiply(M, M, n);
      p /= 2;
    }
    return A;
}

float get_probs(vector<float_vec > M, int States, int Final, int Source, int Time) {
    vector<float_vec > trans_mat = power(M, Time, States);
    return trans_mat[Final - 1][Source - 1];
}

void as_fraction(double number, int cycles = 10, double precision = 5e-4){
    int sign  = (number > 0) ? 1 : -1;
    number = number * sign;
    double new_number,whole_part;
    double decimal_part =  number - (int)number;
    int counter = 0;

    valarray<double> vec_1{double((int) number), 1}, vec_2{1,0}, temporary;

    while(decimal_part > precision & counter < cycles){
        new_number = 1 / decimal_part;
        whole_part = (int) new_number;

        temporary = vec_1;
        vec_1 = whole_part * vec_1 + vec_2;
        vec_2 = temporary;

        decimal_part = new_number - whole_part;
        counter += 1;
    }
    cout<<"p: "<< number <<"\tFraction: " << sign * vec_1[0]<<'/'<< vec_1[1]<<endl;
}

int main() {
    vector<float_vec > Trans{
      {0, 0.3333333, 0.33333333, 0.3333333},
      {0.3333333, 0, 0.3333333, 0.3333333},
      {0.3333333, 0.3333333, 0, 0.3333333},
      {0.3333333, 0.3333333, 0.3333333, 0}
   };

    int States = 4;
    int S = 1, F = 1, T = 4;
    float res = get_probs(Trans, States, F, S, T);
    cout << "Probability of reaching node " << F << " within time " << T << " after starting from node " << S << " is: " << res;
    cout << "\nAs a fraction, this is:\n";
    as_fraction(res);
    return 0;
}