/* 

  Markov chains in Picat.

  From the (Swedish) book "Statistisk Dataanalys", page 299ff.
  
  Transition matrix of market shares for three products A, B, and C.
    
    From   To
           A    B     C
      A   0.7  0.1   0,2
      B   0.2  0.6   0.2
      C   0.4  0.1   0.5

  What is the stable state of the transitions?

  Also, this model solves the reversed problem: 
  given a stable state, generate a transition matrix.


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

*/


import mip.

main => go.

go =>
  transitions(1,Transitions),
  markov_chain(Transitions, X),
  println(X),
  nl.

go2 =>
  transitions(2,Transitions),
  markov_chain(Transitions, X),
  println(X),
  nl.


% reversed problem
go3 => 
  % This is the result from go/0
  X = [ 0.514285714285714,  0.2, 0.285714285714286],
  markov_chain_rev(X, Transitions),
  foreach(Row in Transitions) println(Row) end,
  nl.

% reversed problem 2
go4 =>
  % This is the result from go2/0
  X = [0.367182763904145,0.07548189012462,0.102060886308352,0.15336027885601,0.034415433900105,0.080532266138785,0.13700725112718,0.023360147677889,0.01100743321783,0.015591648745084],
  markov_chain_rev(X, Transitions),

  foreach(Row in Transitions) println(Row) end,
  nl.



%
% Given a transition matrix, yield stable state X
%
markov_chain(Transitions, X) => 

  N = Transitions.length,

  X = new_list(N),
  X :: 0.0..1.0,

  foreach(I in 1..N) 
     X[I] #= sum([Transitions[I,J]*X[J] : J in 1..N])
  end,
  
  sum(X) #= 1.0,

  solve(X).


transitions(1,Transitions) =>
 Transitions = 
 [[0.7,0.2,0.4],
  [0.1,0.6,0.1],
  [0.2,0.2,0.5]].

% Another problem from the book "Optimieringslara":
transitions(2,Transitions) =>
 Transitions = 
[[0.01782,0.54024,0.46916,0.85467,0.30838,0.62996,0.42815,0.22189,0.37009,0.74681],
 [0.04633,0.04427,0.14569,0.03971,0.50937,0.09334,0.01055,0.12776,0.36217,0.04528],
 [0.12289,0.22650,0.12335,0.02186,0.10853,0.14409,0.05414,0.01577,0.01709,0.03761],
 [0.28441,0.01805,0.09189,0.04428,0.00565,0.00538,0.16482,0.20185,0.15423,0.11391],
 [0.00980,0.07623,0.03248,0.00753,0.01028,0.00424,0.12144,0.11167,0.05904,0.00009],
 [0.15721,0.09004,0.07690,0.00448,0.03511,0.05270,0.00859,0.00057,0.00754,0.04810],
 [0.29881,0.00202,0.00490,0.00909,0.00343,0.00296,0.17159,0.05433,0.00924,0.00041],
 [0.02110,0.00004,0.02863,0.00255,0.01232,0.06338,0.02043,0.16868,0.00080,0.00129],
 [0.01843,0.00002,0.02181,0.00453,0.00290,0.00011,0.00819,0.00118,0.00312,0.00163],
 [0.02320,0.00259,0.00519,0.01130,0.00403,0.00384,0.01210,0.09630,0.01668,0.00487]].


%
% reversed problem: Given stable state X, yield a transition matrix.
%
markov_chain_rev(X, Transitions) =>

  N = length(X),

  Transitions = new_array(N,N),
  Transitions.vars() :: 0.0..1.0,

  % Transitions[1,1] #<= 0.99, % to make it more interesting...
  foreach(I in 1..N) 
     X[I] #= sum([Transitions[I,J]*X[J] : J in 1..N])
  end,
  
  sum(X) #= 1.0,


  solve(Transitions).