Markov experiment in Nim

markov.nim

import random, math, times

#randomize()

type
  TransitionMatrix* = object
    order: int
    tMatrix: seq[seq[float]]
    multiplier: seq[int]

  Chain* = object
    state: int
    lastState: seq[int]


func mkMultiplier*(probLen:int, order:int): seq[int] {.inline.}=
  var mp:seq[int]
  mp.add(probLen^(order - 1))
  for i in 1..<order:
    mp.add(mp[i - 1] div probLen)
  mp


func matrixOrder*(nRows: int, nCols: int): int {.inline.}= 
  int(pow(float(nRows), 1/nCols))


func initMatrix*(tmatrix: seq[seq[float]]): TransitionMatrix =
  var tm = TransitionMatrix()
  tm.tMatrix = tmatrix
  tm.order = matrixOrder(len(tmatrix), len(tmatrix[0]))
  tm.multiplier = mkMultiplier(len(tmatrix[0]), tm.order)
  tm


func initChain*(initial: seq[int]): Chain =
  var chain = Chain()
  chain.lastState = initial
  chain.state = initial[^1]
  chain


func rowIndex*(tm: TransitionMatrix, lastState: seq[int]): int {.inline.}= 
  var 
    rowIndex: int = 0 
  for i in 0..<len(lastState):
    rowIndex = rowIndex + (lastState[i] * tm.multiplier[i])  
  rowIndex


func setState*(chain: var Chain, tm: TransitionMatrix, state:int) {.inline.}=
  chain.state = state
  for i in 0..<len(chain.lastState):
    if i == len(chain.lastState)-1:
      chain.lastState[i] = state
    else:
      chain.lastState[i] = chain.lastState[i+1]


func nextState*(chain: var Chain, tm: TransitionMatrix, selector: float):int = 
  let 
    rowIndex = rowIndex(tm, chain.lastState)
    row = tm.tMatrix[rowIndex]
  var 
    probSum: float = 0.0
  for i in 0..<len(row):
    probSum += row[i]
    if probSum > selector:
      setState(chain, tm, i)
      break
  chain.state


#-----------------
let 
  state = @["A", "B", "C"]


echo "First Order"
let 
  pm1:seq[seq[float]] = @[
    @[0.4, 0.3, 0.3],
    @[0.6, 0.1, 0.3],
    @[0.1, 0.7, 0.2],
  ]
  m1 = initMatrix(pm1)
var 
  stream = initRand(7)
  c1 = initChain(@[0])
for i in 0..9:
  let rnd1 = stream.rand(1.0)
  echo state[nextState(c1, m1, rnd1)]

#--------------------
echo "\nSecond Order"
let 
  hist_2 = @[
    "AA", "AB", "AC", 
    "BA", "BB", "BC", 
    "CA", "CB", "CC"
  ]
  pm2 = @[
    @[0.4, 0.3, 0.3],
    @[0.6, 0.1, 0.3],
    @[0.1, 0.7, 0.2],
    
    @[0.4, 0.3, 0.3],
    @[0.6, 0.1, 0.3],
    @[0.1, 0.7, 0.2],
    
    @[0.4, 0.3, 0.3],
    @[0.6, 0.1, 0.3],
    @[0.1, 0.7, 0.2]
  ]
  m2 = initMatrix(pm2)
var
  nstream = initRand(7)
  c2 = initChain(@[0,0]) 
echo hist_2[0]
let time = cpuTime()
for i in 0..100000:
  let 
    rnd2 = nstream.rand(1.0)
    s = nextState(c2, m2, rnd2)
  #echo state[s]
  #echo hist_2[rowIndex(m2, c2.lastState)]
echo cpuTime() - time