/* 

  Markov chains in Picat.

  From Hamdy Taha "Operations Research" (8th edition), page 649ff.
  Fertilizer example.


  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 =>
  Cost =  [100.0, 125.0, 160.0],
  Mat = [[0.3, 0.6,  0.1],
         [0.1, 0.6,  0.3],
         [0.05, 0.4, 0.55]],

  % the transition matrix page 650
  % Mat = [[0.35, 0.6, 0.05],
  %        [0.3 , 0.6, 0.1 ],
  %        [0.25, 0.4, 0.35]],


  N = Cost.length,

  %% Not used since it involves nonlinear *
  % MeanFirstReturnTime = new_list(N),
  % MeanFirstReturnTime :: 0.0..1.0,

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

  % TotCost #= sum([Cost[I]*P[I] : I in 1..N]),
  TotCost #= sum([T : I in 1..N, T #= Cost[I]*P[I]]),

  steady_state_prob(Mat, P),
  % get_mean_first_return_time(P, MeanFirstReturnTime),

  % solve(P ++ MeanFirstReturnTime),
  solve($[min(TotCost)], P),

  println(p=P),
  println(totCost=TotCost),
  % println(meanFirstReturnTime=MeanFirstReturnTime),

  nl.


%
% Calculates the steady state probablity of a transition matrix m
%
steady_state_prob(M, Prob) =>
  Len = M.length,
  foreach(I in 1..Len)
     Prob[I] #= sum([Prob[J]* M[J,I] : J in 1..Len])
  end,
  sum(Prob) #= 1.0.


%
% calculate the mean first return time from a steady state probability array
%
% Gives error: nonlinear constraint,*
get_mean_first_return_time(Prob, Mfrt) =>
   foreach(I in 1..Prob.length)
     Mfrt[I] * Prob[I] #= 1.0
   end.