/* 

  Student interviews in Picat.

  From SPPL model student-interviews.pynb

  Here is the output from the SPPL model of gpa_prior and gpa_posterior for gpa_student(1):

        gpa                   gpa_prior (cumsum < gpa)  gpa_posterior (cumsum < gpa)
        -----------------------------------------------------------------------------
        0.0                   0.0                    0.0
        0.21052631578947367   2.9321229365490116e-08 3.938526167438616e-08
        0.42105263157894735   3.399923850680484e-06  4.566892092513705e-06
        0.631578947368421     5.235987400911553e-05  7.03315441989632e-05
        0.8421052631578947    0.0003515966799506029  0.0004722764885157954
        1.0526315789473684    0.0014935106089171875  0.0020061336928425886
        1.263157894736842     0.004734542792041406   0.0063595971522462345
        1.4736842105263157    0.012227795632612991   0.01642478644698973
        1.6842105263157894    0.027098290512422067   0.03639933542542442
        1.894736842105263     0.05325350529777253    0.07153189980839067
        2.1052631578947367    0.09485332168120815    0.12741017262720986
        2.3157894736842106    0.15542018995700252    0.20876562761532186
        2.526315789473684     0.23665218641802727    0.31787917797401755
        2.7368421052631575    0.3371086961787463     0.45281574132555596
        2.9473684210526314    0.45107070187375614    0.6058933411522036
        3.1578947368421053    0.5680350963535515     0.7630038505974432
        3.3684210526315788    0.6734850651989842     0.9046477961757059
        3.578947368421052     0.7517864026032522     0.983997710206192
        3.789473684210526     0.7922926319251767     0.9970552758034202
        4.0                   0.8                    0.9995398186972839

        
  The mean value of the prior is about: 3.050505050505051
  The mean value of the posterior is about: 2.8080808080808084

  (I am not sure how to get an exact value of the mean in SPPL. 
   For continous distributions, SPPL requires that one work with 
   intervals.)


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

main => go.

/*
  With the observations:
    * the num recruiters was > 30
    * student 1 got 1 offer

  var : num recruiters
  Probabilities:
  31: 0.3146067415730337
  32: 0.2471910112359550
  33: 0.1573033707865168
  34: 0.1011235955056180
  35: 0.0674157303370786
  36: 0.0561797752808989
  37: 0.0337078651685393
  38: 0.0224719101123595
  mean = 32.7753

  var : student 1 gpa
  Probabilities (truncated):
  4.0: 0.1910112359550562
  3.948711928746043: 0.0112359550561798
  3.927603922891366: 0.0112359550561798
  3.828394200772268: 0.0112359550561798
  .........
  2.025023943263929: 0.0112359550561798
  1.920772354951433: 0.0112359550561798
  1.80016802081324: 0.0112359550561798
  1.625747928164307: 0.0112359550561798
  mean = 3.24138

  var : student 1 interviews
  Probabilities (truncated):
  12: 0.1348314606741573
  9: 0.1348314606741573
  14: 0.1235955056179775
  11: 0.1123595505617977
  .........
  18: 0.0224719101123595
  8: 0.0224719101123595
  6: 0.0224719101123595
  5: 0.0224719101123595
  mean = 11.9551

  var : student 1 offers
  Probabilities:
  1: 1.0000000000000000
  mean = 1.0

  var : student 1 perfect
  Probabilities:
  false: 0.8089887640449438
  true: 0.1910112359550562
  mean = [false = 0.808989,true = 0.191011]

  var : student 2 gpa
  Probabilities (truncated):
  4.0: 0.1460674157303371
  3.745046320921631: 0.0112359550561798
  3.629652581677798: 0.0112359550561798
  3.6221684361338: 0.0112359550561798
  .........
  1.840351681255545: 0.0112359550561798
  1.637474342039162: 0.0112359550561798
  1.469359448945245: 0.0112359550561798
  1.461037291355876: 0.0112359550561798
  mean = 2.96963

  var : student 2 interviews
  Probabilities (truncated):
  17: 0.1685393258426966
  16: 0.1685393258426966
  18: 0.1460674157303371
  13: 0.1011235955056180
  .........
  22: 0.0112359550561798
  11: 0.0112359550561798
  10: 0.0112359550561798
  9: 0.0112359550561798
  mean = 16.2921

  var : student 2 offers
  Probabilities (truncated):
  8: 0.1910112359550562
  7: 0.1797752808988764
  6: 0.1797752808988764
  5: 0.1460674157303371
  .........
  3: 0.0337078651685393
  2: 0.0224719101123595
  13: 0.0112359550561798
  12: 0.0112359550561798
  mean = 6.85393

  var : student 2 perfect
  Probabilities:
  false: 0.8539325842696629
  true: 0.1460674157303371
  mean = [false = 0.853933,true = 0.146067]

  var : student gpa prior
  Probabilities (truncated):
  4: 0.1910112359550562
  3.794524864943851: 0.0112359550561798
  3.699832828136565: 0.0112359550561798
  3.643047347302603: 0.0112359550561798
  .........
  1.572949686087224: 0.0112359550561798
  1.559160713111676: 0.0112359550561798
  1.479778615463971: 0.0112359550561798
  1.132941225002102: 0.0112359550561798
  mean = 3.01237

  var : student perfect prior
  Probabilities:
  false: 0.8089887640449438
  true: 0.1910112359550562
  mean = [false = 0.808989,true = 0.191011]

*/
go ?=>
  reset_store,
  run_model(100_000,$model,[show_probs_trunc,mean]),
  nl,
  show_store_lengths,
  % fail,
  nl.
go => true.
  
% Note: The SPPL model has beta(a=7,b=3,scale=4),
% (it's not directly supported in Gamble or in WebPPL)
beta_scale_dist(A,B,Scale) = beta_dist(A,B) * Scale.
model() =>
  NumStudents = 2,

  % Did student get a perfect score?
  StudentPerfect = flip_n(0.2,NumStudents),

  % GPA
  StudentGPA = [condt(StudentPerfect[S],4.0,beta_scale_dist(7,3,4)) : S in 1..NumStudents],

  % Number of recruiters, between 10 and 40
  NumRecruiters = poisson_dist(25),
  observe(NumRecruiters >= 10, NumRecruiters <= 40),

  % Number of interviews
  StudentInterviews = [cases([[StudentGPA[S] == 4, binomial_dist(NumRecruiters,0.9)],
                              [(StudentGPA[S] > 3.5, StudentGPA[S] < 4),
                                                  binomial_dist(NumRecruiters,0.4)],
                              [true, binomial_dist(NumRecruiters,0.5)]
                             ])
                       : S in 1..NumStudents],

  % How many offers did a student get?
  StudentOffers = [ cond(StudentInterviews[S] == 0, 0, binomial_dist(StudentInterviews[S],0.4)) :
                    S in 1..NumStudents],

  % Observations

  observe(NumRecruiters > 30),
  
  % Student 1
  observe(StudentOffers[1] == 1, NumRecruiters > 30),

  % Student 2
  % observe(StudentOffers[2] == 13),

  % Priors
  StudentPerfectPrior = flip(0.2),
  StudentGPAPrior     = cond(not StudentPerfectPrior, 4*beta_dist(7,3),4),

  if observed_ok then
    add("num recruiters",NumRecruiters),
    add("student perfect prior",StudentPerfectPrior),
    add("student gpa prior",StudentGPAPrior),
    
    add("student 1 perfect",StudentPerfect[1]),
    add("student 1 gpa",StudentGPA[1]),
    add("student 1 interviews",StudentInterviews[1]),
    add("student 1 offers",StudentOffers[1]),    
    
    add("student 2 perfect",StudentPerfect[2]),
    add("student 2 gpa",StudentGPA[2]),
    add("student 2 interviews",StudentInterviews[2]),
    add("student 2 offers",StudentOffers[2]),        
  end.