/* 

  Markov process in Picat.

  Cf 
  - ppl_hiring_young_or_old.pi
  - my Gamble model gamble_markov_process.rkt

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import ppl_distributions, ppl_utils.
import util.
% import ordset.

main => go.


/*

Transitions = [[0.1,0.6,0.3],[0.333333,0.333333,0.333333],[0.2,0.2,0.6]]
Matrix
[0.1,0.6,0.3]
[0.333333,0.333333,0.333333]
[0.2,0.2,0.6]

[len1 = 3,len2 = 3]
init = [1,0,0]
Random:
[3,2,3,3,3,1,2,3,2,2,2,3,2,2,3,3,3,3,1,3,1,2,1,3,1,2,1,1,3,1]

PDF:
1 = 0.222223
2 = 0.333334
3 = 0.444443

CDF:
1 = 0.222223
2 = 0.555557
3 = 1.0

Quantile:
0.000001 = 1
0.00001 = 1
0.001 = 1
0.01 = 1
0.025 = 1
0.05 = 1
0.1 = 1
0.25 = 2
0.5 = 2
0.75 = 3
0.84 = 3
0.9 = 3
0.975 = 3
0.99 = 3
0.999 = 3
0.99999 = 3
0.999999 = 3

Markov chains:
start = 1 = [1,3,3,2,2,2,2,3,2,3,3,3,3,3,2,2,3,3,1,3,3,1,2,2,3,2,1,1,2,1]
start = 2 = [2,3,3,3,2,2,2,3,3,3,3,1,2,3,3,3,2,3,3,2,3,1,2,3,3,2,1,3,2,3]
start = 3 = [3,3,3,3,3,3,2,3,3,3,3,2,3,3,1,3,3,3,3,2,3,2,2,1,2,1,2,2,1,3]

Stationary:
stationary1 = [0.222222,0.333333,0.444444]
1: 0.22222222222227
2: 0.33333333333349
3: 0.44444444444424

stationary2 = [0.222222,0.333333,0.444444]
1: 0.22222222222227
2: 0.33333333333349
3: 0.44444444444424


*/
go ?=>
  Transitions = [[1/10,  6/10,  3/10],
                 [1/3 ,  1/3 ,  1/3 ],
                 [2/10, 2/10,  6/10]],

  % Transitions := [[1,0,0],
  %                 [0,1,0],
  %                 [0,0,1]],

  % Transitions := [[1/3,1/3,1/3],
  %                 [  0,1/2,1/2],
  %                 [  0,  0,   1]],

  % Transitions := [[0,1,0],
  %                 [0,0,1],
  %                 [1,0,0]],

  % Transitions := random_transition_matrix(30,1),
  
  println(transitions=Transitions),
  show_matrix(Transitions),
  println([len1=Transitions.len,len2=Transitions[1].len]),


  Len = Transitions.len,
  Init = one_hot(Len,1),
  println(init=Init),
  println("Random:"),
  println(discrete_markov_process_dist_n(Transitions,Init,10,30)),
  nl,
  println("PDF:"),
  foreach(I in 1..Len)
    println(I=discrete_markov_process_dist_pdf(Transitions,Init,10,I))
  end,
  nl,
  println("CDF:"),
  foreach(I in 1..Len)
    println(I=discrete_markov_process_dist_cdf(Transitions,Init,10,I))
  end,
  nl,  
  println("Quantile:"),
  foreach(Q in quantile_qs())
    println(Q=discrete_markov_process_dist_quantile(Transitions,Init,10,Q))
  end,
  nl,
  println("Markov chains:"),
  foreach(I in 1..Len)
    println(start=I=markov_chain(Transitions,one_hot(Len,I),30))
  end,
  nl,
  println("Stationary:"),
  Stationary1 = stationary_dist(Transitions),
  println(stationary1=Stationary1),
  foreach({I,S} in zip(1..Len,Stationary1))
    printf("%w: %.14f\n",I,S)
  end,
  nl,
  Stationary2 = stationary_dist1(Transitions,1.0e-16),
  println(stationary2=Stationary2),  
  foreach({I,S} in zip(1..Len,Stationary1))
    printf("%w: %.14f\n",I,S)
  end,
  nl,
  nl.