# 用吸收马克洛夫链解Doomsday Fuel问题

Riguz Lee 发布于 2022-06-22

# 问题描述

## Doomsday Fuel

Making fuel for the LAMBCHOP's reactor core is a tricky process because of the exotic matter involved. It starts as raw ore, then during processing, begins randomly changing between forms, eventually reaching a stable form. There may be multiple stable forms that a sample could ultimately reach, not all of which are useful as fuel.

Commander Lambda has tasked you to help the scientists increase fuel creation efficiency by predicting the end state of a given ore sample. You have carefully studied the different structures that the ore can take and which transitions it undergoes. It appears that, while random, the probability of each structure transforming is fixed. That is, each time the ore is in 1 state, it has the same probabilities of entering the next state (which might be the same state). You have recorded the observed transitions in a matrix. The others in the lab have hypothesized more exotic forms that the ore can become, but you haven't seen all of them.

Write a function solution(m) that takes an array of array of nonnegative ints representing how many times that state has gone to the next state and return an array of ints for each terminal state giving the exact probabilities of each terminal state, represented as the numerator for each state, then the denominator for all of them at the end and in simplest form. The matrix is at most 10 by 10. It is guaranteed that no matter which state the ore is in, there is a path from that state to a terminal state. That is, the processing will always eventually end in a stable state. The ore starts in state 0. The denominator will fit within a signed 32-bit integer during the calculation, as long as the fraction is simplified regularly.

For example, consider the matrix m:

[
[0,1,0,0,0,1],  # s0, the initial state, goes to s1 and s5 with equal probability
[4,0,0,3,2,0],  # s1 can become s0, s3, or s4, but with different probabilities
[0,0,0,0,0,0],  # s2 is terminal, and unreachable (never observed in practice)
[0,0,0,0,0,0],  # s3 is terminal
[0,0,0,0,0,0],  # s4 is terminal
[0,0,0,0,0,0],  # s5 is terminal
]

So, we can consider different paths to terminal states, such as:

• s0 -> s1 -> s3
• s0 -> s1 -> s0 -> s1 -> s0 -> s1 -> s4
• s0 -> s1 -> s0 -> s5

Tracing the probabilities of each, we find that

• s2 has probability 0
• s3 has probability 3/14
• s4 has probability 1/7
• s5 has probability 9/14

So, putting that together, and making a common denominator, gives an answer in the form of [s2.numerator, s3.numerator, s4.numerator, s5.numerator, denominator] which is [0, 3, 2, 9, 14].

## 问题分析

• 有一系列的状态，如s0, s1, ..., s5
• 一个状态以固定的概率迁移到另一个状态，采用矩阵的形式描述了不同状态间迁移的概率。注意矩阵里面给的是次数而不是概率，比如第一列数据 [0,1,0,0,0,1]，那么s0迁移到s1和s5的概率就是1/2, 1/2
• 终止状态不能迁移到其他状态

• s0 -> s1 -> s3，概率为$$\frac{1}{2} \times \frac{3}{9}=\frac{1}{6}$$
• s0 -> s1 -> s0 -> s1 -> s3, 概率为$$\frac{1}{2} \times \frac{4}{9} \times \frac{1}{2} \times \frac{3}{9}=\frac{2}{9} \times \frac{1}{6}$$
• s0 -> s1 -> s0 -> s1 -> s0 -> s1-> s3, 概率为$$\frac{1}{2} \times \frac{4}{9} \times \frac{1}{2} \times \frac{4}{9} \times \frac{1}{2} \times \frac{3}{9}=\left(\frac{2}{9} \right)^2 \times \frac{1}{6}$$
• ...

# 马尔科夫链

• 存在至少一个吸收态
• 从任意一个状态都可以经过有限步到达一个吸收态

• Q为txt的瞬时态转移概率
• R为txr瞬时态到吸收态的概率
• I为rxr单位矩阵

\begin{align} &N =(I_{t} - Q)^{-1} \\ &B = NR \end{align}

[0,1,0,0,0,1],  # s0, the initial state, goes to s1 and s5 with equal probability
[4,0,0,3,2,0],  # s1 can become s0, s3, or s4, but with different probabilities
[0,0,0,0,0,0],  # s2 is terminal, and unreachable (never observed in practice)
[0,0,0,0,0,0],  # s3 is terminal
[0,0,0,0,0,0],  # s4 is terminal
[0,0,0,0,0,0],  # s5 is terminal

s0, s1为瞬时态，s2, s3, s4, s5为吸收态。有： \begin{align} & Q=\begin{pmatrix} 0 & \frac{1}{2} \\ \frac{4}{9} & 0 \end{pmatrix} \\ & R=\begin{pmatrix} 0 & 0 & 0 & \frac{1}{2} \\ 0 & \frac{3}{9} & \frac{2}{9} & 0\end{pmatrix} \\ & I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \end{align}

Input:
Solution.solution({
{0, 2, 1, 0, 0},
{0, 0, 0, 3, 4},
{0, 0, 0, 0, 0},
{0, 0, 0, 0, 0},
{0, 0, 0, 0, 0}})
Output:
[7, 6, 8, 21]

\begin{align} & Q=\begin{pmatrix} 0 & \frac{2}{3} \\ 0 & 0 \end{pmatrix} \\ & R=\begin{pmatrix} \frac{1}{3} & 0 & 0 \\ 0 & \frac{3}{7} & \frac{4}{7}\end{pmatrix} \\ & I_{t} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\ & N = \begin{pmatrix} 1 & \frac{2}{3} \\ 0 & 1 \end{pmatrix} \\ & B = \begin{pmatrix} \frac{1}{3} & \frac{2}{7} & \frac{8}{21} \\ 0 & \frac{3}{7} & \frac{4}{7} \end{pmatrix} \end{align}

# 代码实现

## 分式计算

public class Fraction {
private final int numerator;
private final int denominator;

public Fraction(int numerator, int denominator) {
if (denominator == 0) {
throw new IllegalArgumentException("Denominator should not be zero");
}
this.numerator = numerator;
this.denominator = denominator;
}
}

### 最大公约数、最小公倍数计算

$gcd(a,b) = gcd(b, a\%b)$

public static int gcd(int m, int n) {
return m % n == 0 ? n : gcd(n, m % n);
}

public static int lcm(int m, int n) {
return m * n / gcd(m, n);
}

### 分式四则运算及约分

public Fraction reduce() {
int gcd = Math.gcd(numerator, denominator);
int n = numerator / gcd, d = denominator / gcd;
boolean swap = d < 0;
return new Fraction(swap ? -n : n, swap ? -d : d);
}

public Fraction negative() {
return new Fraction(-numerator, denominator);
}

public Fraction reciprocal() {
return new Fraction(denominator, numerator);
}

public Fraction convert(int common) {
int factor = common / denominator;
return new Fraction(numerator * factor, common);
}

int common = Math.lcm(denominator, another.denominator);
Fraction a = convert(common), b = another.convert(common);
return new Fraction(a.numerator + b.numerator,
common)
.reduce();
}

public Fraction subtract(Fraction another) {
}

public Fraction multiply(Fraction another) {
return new Fraction(numerator * another.numerator,
denominator * another.denominator)
.reduce();
}

public Fraction divide(Fraction another) {
return multiply(another.reciprocal());
}

## 矩阵运算

public class Matrix {
private final Fraction[][] matrix;
private final int rows;
private final int columns;

public Matrix(int rows, int columns) {
this.rows = rows;
this.columns = columns;
this.matrix = new Fraction[rows][columns];
}
}

### 矩阵减法

public Matrix subtract(Matrix another) {
if (this.rows != another.rows || this.columns != another.columns) {
throw new IllegalArgumentException("Unable to add two matrix with different size");
}
Fraction[][] result = new Fraction[rows][columns];
for (int i = 0; i < rows; i++) {
for (int j = 0; j < columns; j++) {
result[i][j] = matrix[i][j].subtract(another.matrix[i][j]);
}
}
return new Matrix(result, rows, columns);
}

### 矩阵乘法

public Matrix multiply(Matrix another) {
if (this.columns != another.rows) {
throw new IllegalArgumentException("Unable to add two matrix with different size");
}
Matrix r = new Matrix(rows, another.columns);
for (int i = 0; i < rows; i++) {
for (int j = 0; j < another.columns; j++) {
Fraction sum = new Fraction(0, 1);
for (int m = 0; m < columns; m++) {
}
r.matrix[i][j] = sum;
}
}
return r;
}

### 余子式及代数余子式

public Matrix complementMinor(int row, int column) {
Matrix m = new Matrix(rows - 1, columns - 1);
int rowOffset = 0;
for (int i = 0; i < rows; i++) {
if (i == row) {
continue;
}
int columnOffset = 0;
for (int j = 0; j < columns; j++) {
if (j == column) {
continue;
}
m.getMatrix()[rowOffset][columnOffset] = matrix[i][j];
columnOffset++;
}
rowOffset++;
}
return m;
}

### 行列式的值

public Fraction determinantValue() {
if (rows != columns || rows < 1) {
throw new IllegalArgumentException("Not supported");
}
if (rows == 1) {
return matrix[0][0].reduce();
} else if (rows == 2) {
return matrix[0][0].multiply(matrix[1][1]).subtract(matrix[0][1].multiply(matrix[1][0]));
} else {
Fraction sum = new Fraction(0, 1);
for (int i = 0; i < columns; i++) {
Fraction sign = Fraction.of(Math.pow(-1, i));
}
return sum;
}
}

### 矩阵的逆

public Matrix adjugateMatrix() {
Matrix m = new Matrix(rows, columns);
for (int i = 0; i < rows; i++) {
for (int j = 0; j < columns; j++) {
Fraction sign = Fraction.of(Math.pow(-1, (i + j)));
m.matrix[j][i] = sign.multiply(complementMinor(i, j).determinantValue());
}
}
return m;
}

public Matrix inverse() {
Fraction f = determinantValue();
Matrix m = new Matrix(rows, columns);
for (int i = 0; i < rows; i++) {
for (int j = 0; j < columns; j++) {
}
}
return m;
}

## 马尔科夫链计算

public static int[] solution(int[][] nums) {
int totalStates = nums.length, absorbingStates = 0;

Set<Integer> absorbingStateIds = new HashSet<>();
for (int i = 0; i < totalStates; i++) {
if (sum(nums[i]) == 0) {
}
}
absorbingStates = absorbingStateIds.size();
int[][] transformed = new int[totalStates][totalStates];

int [] indexMapping = new int[totalStates];
int offset = 0;
for (int i = 0; i < totalStates; i++) {
if(!absorbingStateIds.contains(i))
indexMapping[offset++] = i;
}
for(int id: absorbingStateIds){
indexMapping[offset++] = id;
}

for(int i = 0; i < totalStates; i++) {
for(int j = 0; j < totalStates; j++) {
transformed[i][j] = nums[indexMapping[i]][indexMapping[j]];
}
}
nums = transformed;
// ...
}

Matrix matrixQ = new Matrix(transientStates, transientStates),
matrixR = new Matrix(transientStates, absorbingStates),
matrixI = new Matrix(transientStates, transientStates);

for (int i = 0; i < transientStates; i++) {
int common = sum(nums[i]);
for (int j = 0; j < transientStates; j++) {
matrixQ.getMatrix()[i][j] = new Fraction(nums[i][j], common);
matrixI.getMatrix()[i][j] = i == j ? new Fraction(1, 1) :
new Fraction(0, 1);
}
}

for (int i = 0; i < transientStates; i++) {
int common = sum(nums[i]);
for (int j = 0; j < absorbingStates; j++) {
matrixR.getMatrix()[i][j] = new Fraction(nums[i][j + transientStates], common);
}
}
Matrix matrixN = matrixI.subtract(matrixQ).inverse();
Matrix matrixB = matrixN.multiply(matrixR);

### 结果约分

int[] result = new int[matrixB.getColumns() + 1];
int lcm = 0;
for (int i = 0; i < matrixB.getColumns(); i++) {
if (matrixB.getMatrix()[0][i].getNumerator() == 0) {
continue;
}
if (lcm == 0) {
lcm = matrixB.getMatrix()[0][i].getDenominator();
} else {
lcm = Math.lcm(lcm, matrixB.getMatrix()[0][i].getDenominator());
}
}
for (int i = 0; i < matrixB.getColumns(); i++) {
result[i] = matrixB.getMatrix()[0][i].getNumerator() * lcm /
matrixB.getMatrix()[0][i].getDenominator();
}
result[result.length - 1] = lcm;

### Test Case 3

foobar:~/doomsday-fuel solee.linux\$ verify Solution.java
Verifying solution...
Test 1 passed!
Test 2 passed!
Test 3 failed  [Hidden]
Test 4 passed! [Hidden]
Test 5 passed! [Hidden]
Test 6 passed! [Hidden]
Test 7 passed! [Hidden]
Test 8 passed! [Hidden]
Test 9 passed! [Hidden]
Test 10 passed! [Hidden]