/* 

  Growth in yeast culture in Picat.

  From A First Course in Mathematical Modeling, 4th edition,
  page 10ff

  Δp(n) = p(n+1) - p(n) = k * (capacity - p(n)) * p(n)
  p(n+1) = p(n) + k * (capacity - p(n))*p(n) 
  capacity ~ 665
  k        ~ 0.00082


  Cf my Gamble model gamble_growth_in_yeast_culture.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.

/*

  * With an observe_abc/3 scale of 1/100 (100 accepted samples)
    4.1s

    Num accepted samples: 99  Total samples: 88812  (0.001%)
    Num accepted samples: 100  Total samples: 89070  (0.001%)

    var : k
    mean = 0.00333136

    var : capacity
    mean = 1788.96

    var : v
    mean = 4.99377


    Summary:
    [orig = 9.6,post = 11.7261,diff = 2.12611]
    [orig = 18.3,post = 20.6325,diff = 2.33249]
    [orig = 29.0,post = 31.2551,diff = 2.25512]
    [orig = 47.2,post = 48.6684,diff = 1.46841]
    [orig = 71.1,post = 72.2751,diff = 1.1751]
    [orig = 119.1,post = 120.667,diff = 1.56747]
    [orig = 174.6,post = 175.359,diff = 0.759106]
    [orig = 257.3,post = 257.598,diff = 0.298275]
    [orig = 350.7,post = 351.549,diff = 0.848643]
    [orig = 441.0,post = 441.72,diff = 0.719673]
    [orig = 513.3,post = 513.363,diff = 0.0630627]
    [orig = 559.7,post = 559.815,diff = 0.11471]
    [orig = 594.8,post = 594.277,diff = -0.523112]
    [orig = 629.4,post = 629.512,diff = 0.112186]
    [orig = 640.8,post = 640.825,diff = 0.0251489]
    [orig = 651.1,post = 651.201,diff = 0.101409]
    [orig = 655.9,post = 655.402,diff = -0.498149]
    [orig = 659.6,post = 659.277,diff = -0.323228]
    sum_diffs = 12.6224


  * With an observe_abc/3 scale of 1/400 (100 accepted samples)
    28.9s

    Num accepted samples: 98  Total samples: 637321  (0.000%)
    Num accepted samples: 99  Total samples: 643644  (0.000%)
    Num accepted samples: 100  Total samples: 654404  (0.000%)

    var : k
    mean = 0.00181904

    var : capacity
    mean = 1929.51

    var : v
    mean = 4.98605


    Summary:
    [orig = 9.6,post = 10.9584,diff = 1.35838]
    [orig = 18.3,post = 20.3022,diff = 2.00222]
    [orig = 29.0,post = 30.555,diff = 1.55501]
    [orig = 47.2,post = 48.8618,diff = 1.66182]
    [orig = 71.1,post = 71.7277,diff = 0.627687]
    [orig = 119.1,post = 120.563,diff = 1.46312]
    [orig = 174.6,post = 175.899,diff = 1.29933]
    [orig = 257.3,post = 258.383,diff = 1.08329]
    [orig = 350.7,post = 350.85,diff = 0.150041]
    [orig = 441.0,post = 441.205,diff = 0.205324]
    [orig = 513.3,post = 513.867,diff = 0.567262]
    [orig = 559.7,post = 559.243,diff = -0.457132]
    [orig = 594.8,post = 594.479,diff = -0.321162]
    [orig = 629.4,post = 628.747,diff = -0.653388]
    [orig = 640.8,post = 640.367,diff = -0.432632]
    [orig = 651.1,post = 650.889,diff = -0.211083]
    [orig = 655.9,post = 655.393,diff = -0.507225]
    [orig = 659.6,post = 659.799,diff = 0.199223]
    sum_diffs = 9.59008

  * With an observe_abc scale of 1/400 and 1000 accepted samples
    3min37s

    Num accepted samples: 998  Total samples: 6884491  (0.000%)
    Num accepted samples: 999  Total samples: 6886165  (0.000%)
    Num accepted samples: 1000  Total samples: 6894531  (0.000%)

    var : k
    mean = 0.00170913

    var : capacity
    mean = 2097.79

    var : v
    mean = 4.78281


    Post summary:
    [orig = 9.6,post = 9.7751,diff = 0.175096]
    [orig = 18.3,post = 18.475,diff = 0.174983]
    [orig = 29.0,post = 29.2003,diff = 0.200321]
    [orig = 47.2,post = 47.1524,diff = -0.0476171]
    [orig = 71.1,post = 71.4761,diff = 0.376107]
    [orig = 119.1,post = 119.272,diff = 0.17151]
    [orig = 174.6,post = 174.774,diff = 0.174184]
    [orig = 257.3,post = 257.226,diff = -0.0738421]
    [orig = 350.7,post = 350.67,diff = -0.0296914]
    [orig = 441.0,post = 441.235,diff = 0.234696]
    [orig = 513.3,post = 513.406,diff = 0.105655]
    [orig = 559.7,post = 559.52,diff = -0.179692]
    [orig = 594.8,post = 594.805,diff = 0.0052886]
    [orig = 629.4,post = 629.22,diff = -0.180126]
    [orig = 640.8,post = 640.981,diff = 0.180759]
    [orig = 651.1,post = 651.233,diff = 0.133457]
    [orig = 655.9,post = 655.662,diff = -0.237734]
    [orig = 659.6,post = 659.493,diff = -0.107414]
    sum_diffs = 1.07594    

*/
go ?=>
  % In hours
  Yeast = [9.6,18.3,29.0,47.2,71.1,119.1,174.6,257.3,350.7,441.0,
           513.3,559.7,594.8,629.4,640.8,651.1,655.9,659.6,661.8],
  println([len=Yeast.len,min=Yeast.min,mean=Yeast.mean,max=Yeast.max,stdev=Yeast.stdev]),
  reset_store,
  run_model(10_000,$model(Yeast),[% show_probs_trunc,
                                  mean,
                                  % show_percentiles,show_histogram,
                                  % show_hpd_intervals,hpd_intervals=[0.84],
                                  min_accepted_samples=1000,show_accepted_samples=true,
                                  presentation=["k","capacity","v"]
                          ]),
  nl,
  % Adjust the size for comparison with the generated posterior in the model.
  post_summary(Yeast[1..Yeast.len-1]), % Use with add_post/2
  nl.
go => true.

model(Yeast) =>

  % YeastMean = Yeast.mean,
  N = Yeast.len,

  K = uniform(0,1),
  Capacity = uniform(0,10000),
  V = normal_dist(5,1),
  observe(V > 0), % truncated

  Y = [normal_dist(Yeast[T-1] + K * (Capacity - Yeast[T-1]),V)  : T in 2..N],

  % We have too reduce the allowed error quite much
  % so get a neat result. Also, reducing the number of minimum accepted samples
  % speeds it up.
  ErrorScale = 1/400, % 1/100, %
  observe_abc(Yeast[1..N-1],Y,ErrorScale),
  
  % observe_abc_adj(Yeast[1..N-1],Y,1), % Very experimental

  
  if observed_ok then
    add("k",K),
    add("capacity",Capacity),
    add("v",V),
    add_post(Yeast[1..N-1],Y) % This is used with the post_summary/1 above
  end.