feat(product/MarkovChains): steady state + ergodicity + reducibility

.vscode/settings.json

@@ -6,6 +6,7 @@
     "boundingbox",
     "dismissable",
     "gsap",
+    "lusolve",
     "Pluggables",
     "preorder",
     "Spacebar",

client/src/graphs/plugins/characteristics/tarjans.ts

@@ -1,21 +1,15 @@
-// TypeScript program to find strongly connected
-// components in a given directed graph using Tarjan's algorithm (single DFS)
-
-// This class represents a directed graph using adjacency list representation
 class TarjanGraph {
   private V: number; // Number of vertices
   private adj: number[][]; // Adjacency list
   private Time: number; // Timer to keep track of discovery time
   private SCCs: number[][] = []; // To store SCCs
 
-  // Constructor
   constructor(v: number) {
     this.V = v;
     this.adj = new Array(v).fill(0).map(() => []);
     this.Time = 0;
   }
 
-  // Function to add an edge into the graph
   addEdge(v: number, w: number): void {
     this.adj[v].push(w);
   }

client/src/products/binary-trees/ui/useNodeColor.ts → client/src/graphs/themes/helpers/useNodeColor.ts

@@ -1,18 +1,17 @@
 import { useTheme } from '@graph/themes/useTheme';
 import type { GNode, Graph } from '@graph/types';
+
 import type { Color } from '@utils/colors';
 
 import type { MaybeRef } from 'vue';
 
 type ColorMap = Map<GNode['id'], Color>;
-type ColorGetter = (nodeId: GNode['id']) => Color;
-
-const DEFAULT_USETHEME_ID = 'node-colorer';
+type ColorGetter = (nodeId: GNode['id']) => Color | undefined;
 
 export const useNodeColor = (
   graph: Graph,
   mapOrGetter: MaybeRef<ColorMap> | ColorGetter,
-  themeId = DEFAULT_USETHEME_ID,
+  themeId: string,
 ) => {
   const get = (nodeId: GNode['id']) => {
     if (typeof mapOrGetter === 'function') return mapOrGetter(nodeId);

client/src/products/binary-trees/ui/useNodeLabel.ts → client/src/graphs/themes/helpers/useNodeLabel.ts

@@ -1,6 +1,5 @@
 import { useTheme } from '@graph/themes/useTheme';
 import type { GNode, Graph } from '@graph/types';
-
 import type { MaybeRef } from 'vue';
 
 /**
@@ -9,14 +8,12 @@ import type { MaybeRef } from 'vue';
 type LabelLike = string | number | boolean;
 
 type LabelMap = Map<GNode['id'], LabelLike>;
-type LabelGetter = (nodeId: GNode['id']) => LabelLike;
-
-const DEFAULT_USETHEME_ID = 'node-labeller';
+type LabelGetter = (nodeId: GNode['id']) => LabelLike | undefined;
 
 export const useNodeLabel = (
   graph: Graph,
   mapOrGetter: MaybeRef<LabelMap> | LabelGetter,
-  themeId = DEFAULT_USETHEME_ID,
+  themeId: string,
 ) => {
   const get = (nodeId: GNode['id']) => {
     if (typeof mapOrGetter === 'function') return mapOrGetter(nodeId);

client/src/products/binary-trees/ui/useBalanceFactorLabels.ts

@@ -3,8 +3,8 @@ import colors from '@utils/colors';
 import type { Color } from '@utils/colors';
 
 import type { TreeControls } from '../useTree';
-import { useNodeColor } from './useNodeColor';
-import { useNodeLabel } from './useNodeLabel';
+import { useNodeColor } from '../../../graphs/themes/helpers/useNodeColor';
+import { useNodeLabel } from '../../../graphs/themes/helpers/useNodeLabel';
 
 export const useBalanceFactorLabels = (graph: Graph, tree: TreeControls) => {
   const { nodeIdToBalanceFactor: nodeToBf } = tree;
@@ -21,8 +21,8 @@ export const useBalanceFactorLabels = (graph: Graph, tree: TreeControls) => {
     return MAP_COLOR[nodeToBf.value.get(nodeId) ?? 0] ?? UNBALANCED_COLOR;
   };
 
-  const { label, unlabel } = useNodeLabel(graph, nodeToBf);
-  const { color, uncolor } = useNodeColor(graph, colorGetter);
+  const { label, unlabel } = useNodeLabel(graph, nodeToBf, 'balance-factor-text');
+  const { color, uncolor } = useNodeColor(graph, colorGetter, 'balance-factor-color');
 
   const activate = () => {
     label();

client/src/products/binary-trees/ui/useHeightLabels.ts

@@ -3,8 +3,8 @@ import colors from '@utils/colors';
 
 import type { TreeControls } from '../useTree';
 import { numberToColor } from './numberToColor';
-import { useNodeColor } from './useNodeColor';
-import { useNodeLabel } from './useNodeLabel';
+import { useNodeColor } from '../../../graphs/themes/helpers/useNodeColor';
+import { useNodeLabel } from '../../../graphs/themes/helpers/useNodeLabel';
 
 export const useHeightLabels = (graph: Graph, tree: TreeControls) => {
   const { nodeIdToHeight } = tree;
@@ -17,8 +17,8 @@ export const useHeightLabels = (graph: Graph, tree: TreeControls) => {
   const colorGetter = (nodeId: GNode['id']) =>
     mapColor(nodeIdToHeight.value.get(nodeId) ?? 0);
 
-  const { label, unlabel } = useNodeLabel(graph, nodeIdToHeight);
-  const { color, uncolor } = useNodeColor(graph, colorGetter);
+  const { label, unlabel } = useNodeLabel(graph, nodeIdToHeight, 'height-text');
+  const { color, uncolor } = useNodeColor(graph, colorGetter, 'height-color');
 
   const activate = () => {
     label();

client/src/products/markov-chains/markov/definitions.ts

@@ -1,15 +1,14 @@
-/**
- * definitions for markov chain terms
- */
 export default {
   valid:
-    'A markov chain is valid if all states have an outgoing probability of 1.',
+    'For every state, the sum of the outgoing transition probabilities is 1.',
   periodic:
-    'A markov chain is said to be periodic if the greatest common divisor of the lengths of all possible cycles is greater than one.',
+    'A Markov chain state is periodic if the greatest common divisor of the lengths of all possible return cycles to that state is greater than 1. A chain is periodic if all its states are periodic.',
   absorbing:
-    'an absorbing markov chain is a chain with at least one absorbing state, which is a state that once entered cannot be left.',
+    'Contains at least one absorbing state, which is a state that, once entered, cannot be left (its self-transition probability is 1).',
   communicatingClasses:
-    'A communicating class is a subset of states in a Markov chain such that any state in the subset can be reached from any other state in the subset.',
+    'A subset of states such that each state in the subset can be reached from every other state in the subset.',
   steadyState:
-    'A steady state is a state in which the system is in equilibrium and the properties of the system do not change over time.',
+    'The probability distribution over states that does not change over time when the chain evolves.',
+  irreducible: 'All states belong to a single communicating class',
+  ergotic: 'Irreducible and Aperiodic'
 } as const;

client/src/products/markov-chains/markov/useMarkovChain.ts

@@ -1,10 +1,10 @@
 import type { Graph } from '@graph/types';
-import { reduceSet } from '@utils/sets';
+import { mergeSetArrayIntoSet } from '@utils/sets';
 
 import { computed } from 'vue';
 
 import { useMarkovClasses } from './useMarkovClasses';
-import { useMarkovNodeWeights } from './useMarkovNodeWeights';
+import { useNodeIdToOutboundWeight } from './useMarkovNodeWeights';
 import { useMarkovPeriodicity } from './useMarkovPeriodicity';
 import { useMarkovSteadyState } from './useMarkovSteadyState';
 
@@ -19,30 +19,55 @@ export const useMarkovChain = (graph: Graph) => {
     nodeIdToTransientClassIndex,
   } = useMarkovClasses(graph);
 
-  const recurrentStates = computed(() => reduceSet(recurrentClasses.value));
-  const transientStates = computed(() => reduceSet(transientClasses.value));
+  const recurrentStates = computed(() => mergeSetArrayIntoSet(recurrentClasses.value));
+  const transientStates = computed(() => mergeSetArrayIntoSet(transientClasses.value));
 
   const { isPeriodic, recurrentClassPeriods } = useMarkovPeriodicity(
     graph,
     recurrentClasses,
   );
 
-  // TODO check with a pro to see if this is correct.
-  const isAbsorbing = computed(() => {
-    if (recurrentClassPeriods.value.length === 0) return false;
-    return recurrentClasses.value.every((recurrentClass) => {
-      return recurrentClass.size === 1;
-    });
+  const absorbingStates = computed(() => {
+    return mergeSetArrayIntoSet(recurrentClasses.value.filter((rc) => rc.size === 1))
+  })
+
+  const isChainAbsorbing = computed(() => {
+    return absorbingStates.value.size > 0
   });
 
   const communicatingClasses = computed(() => {
     return graph.characteristics.stronglyConnectedComponents.value;
   });
 
-  const { nodeIdToOutgoingWeight, illegalNodeIds } =
-    useMarkovNodeWeights(graph);
+  const isIrreducible = computed(() => {
+    return recurrentClasses.value.length === 1 && transientClasses.value.length === 0
+  })
+
+  const isErgodic = computed(() => {
+    return !isPeriodic.value && isIrreducible.value
+  })
 
-  const steadyState = useMarkovSteadyState();
+  const nodeIdToOutgoingWeight = useNodeIdToOutboundWeight(graph);
+
+  const invalidStates = computed(() => {
+    const invalidStatesArr = graph.nodes.value
+      .map((node) => node.id)
+      .filter((nodeId) => {
+        const outgoingWeight = nodeIdToOutgoingWeight.value.get(nodeId)!
+        return outgoingWeight.valueOf() !== 1
+      })
+    return new Set(invalidStatesArr)
+  });
+
+  const isChainValid = computed(() => invalidStates.value.size === 0)
+
+  const steadyState = useMarkovSteadyState(graph);
+  const uniqueSteadyState = computed(() => {
+    if (!isChainValid.value) return { type: 'error-invalid' } as const;
+    if (recurrentClasses.value.length > 1) return { type: 'error-not-unique' } as const;
+    if (isPeriodic.value) return { type: 'error-no-convergence' } as const
+    return { type: 'success', data: steadyState.value } as const;
+  })
 
   return {
     communicatingClasses,
@@ -57,13 +82,17 @@ export const useMarkovChain = (graph: Graph) => {
     nodeIdToTransientClassIndex,
 
     isPeriodic,
-    isAbsorbing,
+    isChainAbsorbing,
+    absorbingStates,
+    isIrreducible,
 
-    // TODO implement a correct version of steady state
-    steadyState,
+    uniqueSteadyState,
 
     nodeIdToOutgoingWeight,
-    illegalNodeIds,
+
+    invalidStates,
+    isChainValid,
+    isErgodic,
   };
 };
 

client/src/products/markov-chains/markov/useMarkovClasses.ts

@@ -1,7 +1,6 @@
 import type { GNode, Graph } from '@graph/types';
 
 import { computed } from 'vue';
-import type { Ref } from 'vue';
 
 import type { ComponentAdjacencyMap } from './useComponentAdjacencyMap';
 

client/src/products/markov-chains/markov/useMarkovNodeWeights.ts

@@ -9,13 +9,12 @@ import { computed } from 'vue';
 export type NodeIdToOutgoingWeight = Map<GNode['id'], Fraction>;
 
 /**
- * reactive outgoing weights of nodes including computing illegal nodes/states
- * for a markov chain
+ * maps node ids to the sum of their outgoing edge weights
  */
-export const useMarkovNodeWeights = (graph: Graph) => {
+export const useNodeIdToOutboundWeight = (graph: Graph) => {
   const { getOutboundEdges, getEdgeWeightFraction } = graph.helpers;
 
-  const nodeIdToOutgoingWeight = computed(() => {
+  return computed(() => {
     return graph.nodes.value.reduce<NodeIdToOutgoingWeight>((acc, node) => {
       const outgoingEdges = getOutboundEdges(node.id);
       const weights = outgoingEdges.map((edge) =>
@@ -28,27 +27,4 @@ export const useMarkovNodeWeights = (graph: Graph) => {
       return acc;
     }, new Map());
   });
-
-  const illegalNodeIds = computed(() => {
-    return new Set(
-      graph.nodes.value
-        .filter(
-          (node) => nodeIdToOutgoingWeight.value.get(node.id)?.valueOf() !== 1,
-        )
-        .map((node) => node.id),
-    );
-  });
-
-  return {
-    /**
-     * maps node ids to the sum of their outgoing edge weights
-     */
-    nodeIdToOutgoingWeight,
-    /**
-     * set of node ids with outgoing edge weights not equal to 1 (within tolerance)
-     */
-    illegalNodeIds,
-  };
 };
-
-export type MarkovNodeWeights = ReturnType<typeof useMarkovNodeWeights>;

client/src/products/markov-chains/markov/useMarkovPeriodicity.ts

@@ -101,7 +101,6 @@ export const useMarkovPeriodicity = (
   const { adjacencyList } = graph.adjacencyList;
 
   const recurrentClassPeriods = computed(() => {
-    // console.log(getPeriod(adjacencyList.value, recurrentClasses.value[0]));
     const res = recurrentClasses.value.map((recurrentClass) =>
       getPeriod(adjacencyList.value, recurrentClass),
     );

client/src/products/markov-chains/markov/useMarkovSteadyState.ts

@@ -1,11 +1,75 @@
 import { computed } from 'vue';
+import type { Graph } from '@graph/types';
+import type { TransitionMatrix } from '@graph/useTransitionMatrix';
+import { divide, equal, fraction, Fraction, identity, Matrix, multiply, subtract, transpose } from 'mathjs';
 
-/**
- * TODO implement!
- * reactive unique steady state of a markov chain
- */
-export const useMarkovSteadyState = () => {
-  return computed(() => undefined);
-};
+// TODO add tests!!! rref function is untrusted AI output :(
+const rref = (matrix: TransitionMatrix<Fraction>) => {
+  const A = matrix.map(row => row.map(cell => fraction(cell))); // deep copy
+  let lead = 0;
+  const rowCount = A.length;
+  const colCount = A[0].length;
+
+  for (let r = 0; r < rowCount; r++) {
+    if (lead >= colCount) return A;
+    let i = r;
+    while (equal(A[i][lead], 0)) {
+      i++;
+      if (i === rowCount) {
+        i = r;
+        lead++;
+        if (lead === colCount) return A;
+      }
+    }
+
+    // Swap rows
+    [A[i], A[r]] = [A[r], A[i]];
+
+    // Normalize row
+    const val = A[r][lead];
+    for (let j = 0; j < colCount; j++) {
+      A[r][j] = divide(A[r][j], val) as Fraction;
+    }
+
+    // Eliminate other rows
+    for (let i2 = 0; i2 < rowCount; i2++) {
+      if (i2 !== r) {
+        const val2 = A[i2][lead];
+        for (let j = 0; j < colCount; j++) {
+          A[i2][j] = subtract(
+            A[i2][j],
+            multiply(val2, A[r][j]) as Fraction,
+          );
+        }
+      }
+    }
+    lead++;
+  }
+  return A;
+}
+
+const getSteadyStateVector = (transitionMatrix: TransitionMatrix<Fraction>) => {
+  if (transitionMatrix.length === 0) return []
+
+  const identityMatrix = identity(transitionMatrix.length) as Matrix<Fraction>
+  const transposedTransMatrix = transpose(transitionMatrix)
+
+  const subtractedMatrix = subtract(transposedTransMatrix, identityMatrix).valueOf() as TransitionMatrix<Fraction>
+  subtractedMatrix[subtractedMatrix.length - 1] = Array(subtractedMatrix.length).fill(fraction(1));
+
+  const b = Array(subtractedMatrix.length - 1)
+    .fill(fraction(0))
+    .concat([fraction(1)]);
+
+  const augmented: TransitionMatrix<Fraction> = subtractedMatrix.map((row, i) => [...row, b[i]]);
+
+  return rref(augmented).map((row) => row.at(-1)!)
+}
+
+export const useMarkovSteadyState = (graph: Graph) => computed(() => {
+  const matrix = graph.transitionMatrix.fracTransitionMatrix.value;
+  if (matrix.length === 0) return []
+  return getSteadyStateVector(graph.transitionMatrix.fracTransitionMatrix.value);
+});
 
 export type MarkovSteadyState = ReturnType<typeof useMarkovSteadyState>;