/* 

  Changepoint detection in Picat.

  From https://nishanthu.github.io/articles/ChangePointAnalysis.html

  Note: This is not easy for Picat PPL. We have to set a quite high value
        for the acceptable interval in
        S = [ cond(abs(Y[I]-Ys[I]) <= 3.0,1,0) : I in 1..N].sum,

  This is a port of my Gamble model gamble_changepoint.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.

main => go.

/*
  The change point is position 13:

  For the original (R generated data) with NumRuns = 100_000

  ys_len = 33
  var : mu1
  Probabilities (truncated):
  2.479500198048375: 0.0096153846153846
  2.447513537649881: 0.0096153846153846
  2.315999894736822: 0.0096153846153846
  2.190430345848288: 0.0096153846153846
  .........
  0.007489296974852: 0.0096153846153846
  0.004954568974014: 0.0096153846153846
  0.003207456341353: 0.0096153846153846
  0.001634187067241: 0.0096153846153846
  mean = 0.881959

  var : mu2
  Probabilities (truncated):
  2.594887529184876: 0.0096153846153846
  2.59114247014478: 0.0096153846153846
  2.590477661419693: 0.0096153846153846
  2.589387711195376: 0.0096153846153846
  .........
  1.950473232318729: 0.0096153846153846
  1.900622601710653: 0.0096153846153846
  1.837519713039057: 0.0096153846153846
  1.669302996371576: 0.0096153846153846
  mean = 2.36126

  var : cp
  Probabilities (truncated):
  14: 0.2115384615384615
  13: 0.1442307692307692
  12: 0.1442307692307692
  10: 0.0961538461538462
  .........
  18: 0.0096153846153846
  17: 0.0096153846153846
  8: 0.0096153846153846
  2: 0.0096153846153846
  mean = 12.5



  * Using min_accepted_samples=1000

  ...
  Num accepted samples: 997  Total samples: 1014197  (0.001%)
  Num accepted samples: 998  Total samples: 1014423  (0.001%)
  Num accepted samples: 999  Total samples: 1015934  (0.001%)
  Num accepted samples: 1000  Total samples: 1018015  (0.001%)
  var : mu1
  Probabilities (truncated):
  2.548775555655669: 0.0010000000000000
  2.543106970891154: 0.0010000000000000
  2.524649793793332: 0.0010000000000000
  2.524366092828493: 0.0010000000000000
  .........
  0.00555788186716: 0.0010000000000000
  0.005489988911137: 0.0010000000000000
  0.003027786641816: 0.0010000000000000
  0.001398106584301: 0.0010000000000000
  mean = 0.944525

  var : mu2
  Probabilities (truncated):
  2.596772147334582: 0.0010000000000000
  2.595995324845883: 0.0010000000000000
  2.595954777479061: 0.0010000000000000
  2.595918912468468: 0.0010000000000000
  .........
  1.632666607496248: 0.0010000000000000
  1.600443098718068: 0.0010000000000000
  1.265094447118136: 0.0010000000000000
  0.991761564904515: 0.0010000000000000
  mean = 2.3638

  var : cp
  Probabilities (truncated):
  14: 0.1870000000000000
  13: 0.1650000000000000
  12: 0.1090000000000000
  11: 0.0910000000000000
  .........
  28: 0.0010000000000000
  26: 0.0010000000000000
  25: 0.0010000000000000
  22: 0.0010000000000000
  mean = 12.316

*/
go ?=>
  % This was generated by R:
  %  > c(rnorm(13, mean=0, sd=0.5), rnorm(20, mean=4, sd=0.5))
  % I.e. the changepoint is at 13
  Ys = [0.90685440,0.38219559,0.13508408,-0.96054590,0.78553871,0.06791336,-0.09750806,
        -0.85431057,0.33478539,0.76388689,0.52911522,-0.41226421,1.53302740,
        5.01118838,4.75164473,3.79149212,3.97544567,4.33985561,4.52511631,4.21723600,
        4.22733898,3.80925885,3.65504091,3.25267865,3.85339052,3.40136529,4.21462117,
        3.94033196,4.00900359,4.44148942,4.23449516,4.52032000,4.41335270],
  % Picat> Ys=normal_dist_n(0,0.5,13) ++ normal_dist_n(4,0.5,20)
  % Ys = [-0.56171842476207,-0.426284918711003,0.635156161799009,-0.159162666097351,-0.361197417932874,0.291740884621344,-0.776407593769485,0.19604689500553,0.147032998014698,0.241381612561416,-0.15536817530333,-0.182741429736374,-0.20098528889876,3.301351952648498,3.49744705747228,4.137318723033274,3.733036385398599,3.813803399756961,3.209041041127601,4.374991813129283,4.373639193787166,3.561698242471097,4.267851694025974,3.466886751431267,3.332424815035318,3.921145949217203,3.644655004618383,3.705261941607204,3.377084712028713,3.768849554800923,4.589789657718582,3.923162060177829,3.335307515593086],

  println(ys_len=Ys.len),
  reset_store(),
  run_model(100_000,$model(Ys),[show_probs_trunc,mean,
                                show_hpd_intervals,hpd_intervals=[0.94],
                                min_accepted_samples=1000,show_accepted_samples=true
                                ]),
  show_store_lengths(), nl,
  nl.
go => true.
  

model(Ys) =>
  N = Ys.len,
  Mu = Ys.mean,

  Mu1 = uniform(0,Mu),
  Mu2 = uniform(0,Mu),

  CP = random_integer1(N),

  Mu0 = [ cond(I < CP,Mu1,Mu2) : I in 1..N],

  Y = [ normal_dist(Mu0[I],1) : I in 1..N],
  
  % Observe the valyes in Y
  % We have to increase the accepted limit quite much
  % to make it work.
  S = [ cond(abs(Y[I]-Ys[I]) <= 3.0,1,0) : I in 1..N].sum,
  
  observe(S == N),
  if observed_ok then
    add("mu1",Mu1),
    add("mu2",Mu2),
    add("cp",CP)
  end.