1:- module(mcintyre,[
    2  mc_prob/2,
    3  mc_prob/3,
    4  mc_rejection_sample/5,
    5  mc_rejection_sample/4,
    6  mc_sample/3,
    7  mc_sample/4,
    8  mc_sample_arg/4,
    9  mc_sample_arg/5,
   10  mc_mh_sample/4,
   11  mc_mh_sample/5,
   12  mc_lw_sample/4,
   13  mc_gibbs_sample/5,
   14  mc_gibbs_sample/4,
   15  mc_gibbs_sample/3,
   16  mc_rejection_sample_arg/6,
   17  mc_rejection_sample_arg/5,
   18  mc_mh_sample_arg/5,
   19  mc_mh_sample_arg/6,
   20  mc_gibbs_sample_arg/4,
   21  mc_gibbs_sample_arg/5,
   22  mc_gibbs_sample_arg/6,
   23  mc_sample_arg_first/4,
   24  mc_sample_arg_first/5,
   25  mc_sample_arg_one/4,
   26  mc_sample_arg_one/5,
   27  mc_sample_arg_raw/4,
   28  mc_expectation/4,
   29  mc_mh_expectation/5,
   30  mc_mh_expectation/6,
   31  mc_gibbs_expectation/4,
   32  mc_gibbs_expectation/5,
   33  mc_gibbs_expectation/6,
   34  mc_rejection_expectation/5,
   35  set_mc/2,setting_mc/2,
   36  mc_load/1,mc_load_file/1,
   37  take_a_sample/5,
   38  sample_head/5,
   39  mc_lw_sample_arg/5,
   40  mc_lw_sample_arg_log/5,
   41  mc_lw_expectation/5,
   42  mc_particle_sample_arg/5,
   43  mc_particle_sample/4,
   44  mc_particle_expectation/5,
   45  gaussian/5,
   46  gaussian/4,
   47  add_prob/3,
   48  op(600,xfy,'::'),
   49  op(600,xfy,'~'),
   50  op(500,xfx,'~='),
   51  op(1200,xfy,':='),
   52  op(1150,fx,mcaction),
   53  ~= /2,
   54  swap/2,
   55  msw/2,
   56  set_sw/2
   57  ]).

mcintyre

This module performs reasoning over Logic Programs with Annotated Disjunctions and CP-Logic programs. It reads probabilistic program and computes the probability of queries using sampling.

See https://friguzzi.github.io/cplint/ for details.

Reexports cplint_util and clpr.

author
- Fabrizio Riguzzi
copyright
- Fabrizio Riguzzi */
license
- Artistic License 2.0 https://opensource.org/licenses/Artistic-2.0
   75:- reexport(library(cplint_util)).   76:- reexport(library(clpr)).   77
   78:-meta_predicate s(:,-).   79:-meta_predicate mc_prob(:,-).   80:-meta_predicate mc_prob(:,-,+).   81:-meta_predicate mc_sample(:,+,-,-,-).   82:-meta_predicate mc_rejection_sample(:,:,+,-,-,-).   83:-meta_predicate mc_rejection_sample(:,:,+,-,+).   84:-meta_predicate mc_rejection_sample(:,:,+,-).   85:-meta_predicate mc_mh_sample(:,:,+,-).   86:-meta_predicate mc_mh_sample(:,:,+,-,+).   87:-meta_predicate mc_mh_sample(:,:,+,+,+,-,-,-).   88:-meta_predicate mc_gibbs_sample(:,:,+,-,+).   89:-meta_predicate mc_gibbs_sample(:,+,-,+).   90:-meta_predicate mc_gibbs_sample(:,+,-).   91:-meta_predicate mc_sample(:,+,-).   92:-meta_predicate mc_sample(:,+,-,+).   93:-meta_predicate mc_sample_arg(:,+,?,-).   94:-meta_predicate mc_sample_arg(:,+,?,-,+).   95:-meta_predicate mc_rejection_sample_arg(:,:,+,?,-,+).   96:-meta_predicate mc_rejection_sample_arg(:,:,+,?,-).   97:-meta_predicate mc_mh_sample_arg0(:,:,+,+,+,?,-).   98:-meta_predicate mc_mh_sample_arg(:,:,+,?,-,+).   99:-meta_predicate mc_mh_sample_arg(:,:,+,?,-).  100:-meta_predicate mc_gibbs_sample_arg0(:,:,+,+,+,?,-).  101:-meta_predicate mc_gibbs_sample_arg0(:,+,+,+,?,-).  102:-meta_predicate mc_gibbs_sample_arg(:,:,+,?,-,+).  103:-meta_predicate mc_gibbs_sample_arg(:,+,?,-,+).  104:-meta_predicate mc_gibbs_sample_arg(:,+,?,-).  105:-meta_predicate mc_sample_arg_first(:,+,?,-).  106:-meta_predicate mc_sample_arg_first(:,+,?,-,+).  107:-meta_predicate mc_sample_arg_one(:,+,?,-,+).  108:-meta_predicate mc_sample_arg_one(:,+,?,-).  109:-meta_predicate mc_sample_arg_raw(:,+,?,-).  110:-meta_predicate mc_expectation(:,+,?,-).  111:-meta_predicate mc_mh_expectation(:,:,+,?,-).  112:-meta_predicate mc_mh_expectation(:,:,+,?,-,+).  113:-meta_predicate mc_gibbs_expectation(:,+,?,-).  114:-meta_predicate mc_gibbs_expectation(:,+,?,-,+).  115:-meta_predicate mc_gibbs_expectation(:,:,+,?,-,+).  116:-meta_predicate mc_rejection_expectation(:,:,+,?,-).  117:-meta_predicate montecarlo_cycle(-,-,:,-,-,-,-,-,-).  118:-meta_predicate montecarlo(-,-,-,:,-,-).  119:-meta_predicate initial_sample_cycle(:).  120:-meta_predicate gibbs_sample_cycle(:).  121:-meta_predicate initial_sample(:).  122:-meta_predicate lw_sample_bool(+,:,:,-).  123:-meta_predicate initial_sample_neg(:).  124
  125:-meta_predicate mc_lw_sample(:,:,+,-).  126:-meta_predicate mc_lw_sample_arg(:,:,+,?,-).  127:-meta_predicate mc_lw_sample_arg_log(:,:,+,?,-).  128:-meta_predicate mc_lw_expectation(:,:,+,?,-).  129:-meta_predicate mc_particle_expectation(:,:,+,?,-).  130:-meta_predicate mc_particle_sample_arg(:,:,+,?,-).  131:-meta_predicate mc_particle_sample(:,:,+,-).  132:-meta_predicate particle_sample_arg(:,:,+,?,-).  133:-meta_predicate particle_sample_first(+,+,:,:,?).  134:-meta_predicate particle_sample_arg_gl(:,:,+,+,+,-).  135:-meta_predicate particle_sample_first_gl(+,+,:,:,-,-).  136:-meta_predicate particle_sample(+,+,:,-).  137:-meta_predicate lw_sample_cycle(:).  138:-meta_predicate lw_sample_weight_cycle(:,-).  139:-meta_predicate ~=(:,-).  140:-meta_predicate msw(:,-).  141:-meta_predicate set_sw(:,+).  142
  143:-meta_predicate mh_sample_arg(+,+,+,:,:,?,+,-,+,-).  144:-meta_predicate mh_montecarlo(+,+,+,+,+,+,-,:,:,+,-).  145:-meta_predicate gibbs_montecarlo(+,+,+,:,-).  146:-meta_predicate gibbs_montecarlo(+,+,+,:,:,-).  147
  148:-meta_predicate mc_gibbs_sample(:,:,+,+,+,-,-,-).  149:-meta_predicate mc_gibbs_sample(:,+,+,+,-,-,-).  150:-meta_predicate gibbs_sample_arg(+,:,:,+,?,+,-).  151:-meta_predicate gibbs_sample_arg(+,:,+,?,+,-).  152
  153
  154:-meta_predicate rejection_montecarlo(+,+,+,:,:,-,-).  155:-meta_predicate set_mc(:,+).  156:-meta_predicate setting_mc(:,-).  157
  158:-use_module(library(lists)).  159:-use_module(library(apply)).  160:-use_module(library(assoc)).  161:-use_module(library(clpr)).  162:-use_module(library(clpfd)).  163:-use_module(library(matrix)).  164:-use_module(library(option)).  165:-use_module(library(predicate_options)).  166:-use_module(cplint_util).  167
  168:-predicate_options(mc_prob/3,3,[bar(-)]).  169:-predicate_options(mc_sample/4,4,[successes(-),failures(-),bar(-)]).  170:-predicate_options(mc_mh_sample/5,5,[successes(-),failures(-),mix(+),lag(+)]).  171:-predicate_options(mc_gibbs_sample/5,5,[successes(-),failures(-),mix(+),block(+)]).  172:-predicate_options(mc_gibbs_sample/4,4,[successes(-),failures(-),mix(+),block(+)]).  173:-predicate_options(mc_rejection_sample/5,5,[successes(-),failures(-)]).  174:-predicate_options(mc_sample_arg/5,5,[bar(-)]).  175:-predicate_options(mc_rejection_sample_arg/6,6,[bar(-)]).  176:-predicate_options(mc_mh_sample_arg/6,6,[mix(+),lag(+),bar(-)]).  177:-predicate_options(mc_gibbs_sample_arg/6,6,[mix(+),bar(-),block(+)]).  178:-predicate_options(mc_gibbs_sample_arg/5,5,[mix(+),bar(-),block(+)]).  179:-predicate_options(mc_sample_arg_first/5,5,[bar(-)]).  180:-predicate_options(mc_sample_arg_one/5,5,[bar(-)]).  181:-predicate_options(mc_mh_expectation/6,6,[mix(+),lag(+)]).  182:-predicate_options(mc_gibbs_expectation/6,6,[mix(+),block(+)]).  183:-predicate_options(mc_gibbs_expectation/5,5,[mix(+),block(+)]).  184:-predicate_options(histogram/3,3,[max(+),min(+),nbins(+)]).  185:-predicate_options(density/3,3,[max(+),min(+),nbins(+)]).  186:-predicate_options(density2d/3,3,[xmax(+),xmin(+),ymax(+),ymin(+),nbins(+)]).  187
  188:- style_check(-discontiguous).  189
  190
  191
  192:- thread_local v/3,rule_n/1,mc_input_mod/1,local_mc_setting/2.  193
  194/*:- multifile one/2,zero/2,and/4,or/4,bdd_not/3,init/3,init_bdd/2,init_test/1,
  195  end/1,end_bdd/1,end_test/0,ret_prob/3,em/9,randomize/1,
  196  get_var_n/5,add_var/5,equality/4.*/
  197%  remove/3.
  198
  199
  200/* k
  201 * -
  202 * This parameter shows the amount of items of the same type to consider at once.
  203 *
  204 * Default value:	500
  205 * Applies to:		bestfirst, bestk, montecarlo
  206 */
  207default_setting_mc(k, 500).
  208/* min_error
  209 * ---------
  210 * This parameter shows the threshold for the probability interval.
  211 *
  212 * Default value:	0.02
  213 * Applies to:		bestfirst, montecarlo
  214 */
  215default_setting_mc(min_error, 0.02).
  216
  217default_setting_mc(max_samples,5e4).
  218
  219
  220default_setting_mc(epsilon_parsing, 1e-5).
  221/* on, off */
  222
  223default_setting_mc(compiling,off).
  224
  225:-set_prolog_flag(unknown,warning).  226
  227default_setting_mc(depth_bound,false).  %if true, it limits the derivation of the example to the value of 'depth'
  228default_setting_mc(depth,2).
  229default_setting_mc(single_var,false). %false:1 variable for every grounding of a rule; true: 1 variable for rule (even if a rule has more groundings),simpler.
  230default_setting_mc(prism_memoization,false). %false: original prism semantics, true: semantics with memoization
 mc_load(++File:atom) is det
Loads File.lpad if it exists, otherwise loads File.cpl if it exists. /
  236mc_load(File):-
  237  must_be(atom,File),
  238  atomic_concat(File,'.lpad',FileLPAD),
  239  (exists_file(FileLPAD)->
  240    mc_load_file(FileLPAD)
  241  ;
  242    atomic_concat(File,'.cpl',FileCPL),
  243    (exists_file(FileCPL)->
  244      mc_load_file(FileCPL)
  245    )
  246  ).
 mc_load_file(++FileWithExtension:atom) is det
Loads FileWithExtension. /
  253mc_load_file(File):-
  254  must_be(atom,File),
  255  begin_lpad_pred,
  256  user:consult(File),
  257  end_lpad_pred.
 s(:Query:conjunction_of_literals, -Probability:float) is nondet
The predicate computes the probability of the ground query Query. If Query is not ground, it returns in backtracking all instantiations of Query together with their probabilities /
  266s(Mo:Goal,P):-
  267  Mo:local_mc_setting(min_error, MinError),
  268  Mo:local_mc_setting(k, K),
  269% Resetting the clocks...
  270% Performing resolution...
  271  copy_term(Goal,Goal1),
  272  numbervars(Goal1),
  273  save_samples(Mo,Goal1),
  274  montecarlo_cycle(0, 0, Mo:Goal, K, MinError, _Samples, _Lower, P, _Upper),
  275  !,
  276  erase_samples(Mo),
  277  restore_samples_delete_copy(Mo,Goal1).
  278
  279save_samples_tab(M,I,S):-
  280  M:sampled(R,Sub,V),
  281  assert(M:mem(I,S,R,Sub,V)),
  282  retract(M:sampled(R,Sub,V)),
  283  fail.
  284
  285save_samples_tab(M,I,S):-
  286  M:sampled_g(Sub,V),
  287  assert(M:mem(I,S,rw,Sub,V)),
  288  retract(M:sampled_g(Sub,V)),
  289  fail.
  290
  291save_samples_tab(M,I,S):-
  292  M:sampled_g(Sub),
  293  assert(M:mem(I,S,r,Sub,1)),
  294  retract(M:sampled_g(Sub)),
  295  fail.
  296
  297save_samples_tab(_M,_I,_S).
  298
  299save_samples(M,I,S):-
  300  M:sampled(R,Sub,V),
  301  assert(M:mem(I,S,R,Sub,V)),
  302  retract(M:sampled(R,Sub,V)),
  303  fail.
  304
  305save_samples(M,_I,_S):-
  306  retractall(M:sampled_g(_,_)),
  307  retractall(M:sampled_g(_)).
  308
  309save_samples(M,G):-
  310  forall(retract(M:sampled(R,Sub,V)),
  311         assertz(M:mem(G,R,Sub,V))).
  312
  313restore_samples(M,I,S):-
  314  M:mem(I,S,R,Sub,V),
  315  assert_samp(R,M,Sub,V),
  316  fail.
  317
  318restore_samples(_M,_I,_S).
  319
  320assert_samp(r,M,Sub,_V):-!,
  321  assertz(M:sampled_g(Sub)).
  322
  323assert_samp(rw,M,Sub,V):-!,
  324  assertz(M:sampled_g(Sub,V)).
  325
  326assert_samp(R,M,Sub,V):-
  327  assertz(M:sampled(R,Sub,V)).
  328
  329
  330restore_samples(M,G):-
  331  M:mem(G,R,Sub,V),
  332  assertz(M:sampled(R,Sub,V)),
  333  fail.
  334
  335restore_samples(_M,_G).
  336
  337
  338restore_samples_delete_copy(M,G):-
  339  retract(M:mem(G,R,Sub,V)),
  340  assertz(M:sampled(R,Sub,V)),
  341  fail.
  342
  343restore_samples_delete_copy(_M,_G).
  344
  345save_samples_copy(M,G):-
  346  M:sampled(R,Sub,V),
  347  assert(M:mem(G,R,Sub,V)),
  348  fail.
  349
  350save_samples_copy(_M,_G).
  351
  352delete_samples_copy(M,G):-
  353  retract(M:mem(G,_R,_Sub,_V)),
  354  fail.
  355
  356delete_samples_copy(_M,_G).
  357
  358count_samples(M,N):-
  359  findall(a,M:sampled(_Key,_Sub,_Val),L),
  360  length(L,N).
  361
  362
  363resample(_M,0):-!.
  364
  365resample(M,N):-
  366  findall(sampled(Key,Sub,Val),M:sampled(Key,Sub,Val),L),
  367  sample_one(L,S),
  368  retractall(M:S),
  369  N1 is N-1,
  370  resample(M,N1).
  371
  372
  373erase_samples(M):-
  374  retractall(M:sampled(_,_,_)),
  375  retractall(M:sampled_g(_,_)),
  376  retractall(M:sampled_g(_)).
  377
  378print_samples(M):-
  379  M:sampled(Key,Sub,Val),
  380  write(Key-(Sub,Val)),nl,
  381  fail.
  382
  383print_samples(_M):-
  384  write(end),nl.
  385
  386montecarlo_cycle(N0, S0, M:Goals, K, MinError, Samples, Lower, Prob, Upper):-!,
  387  montecarlo(K,N0, S0, M:Goals, N, S),
  388  P is S / N,
  389  D is N - S,
  390  Semi is 1.95996 * sqrt(P * (1 - P) / N),
  391  Int is 2 * Semi,
  392  M:local_mc_setting(max_samples,MaxSamples),
  393  /*   N * P > 5;   N * S / N > 5;   S > 5
  394  *   N (1 - P) > 5;   N (1 - S / N) > 5;   N (N - S) / N > 5;   N - S > 5
  395  */
  396  %format("Batch: samples ~d positive ~d interval ~f~n",[N,S,Int]),
  397% flush_output,
  398  (((S > 5, D > 5, (Int < MinError; Int =:= 0));
  399    ((Int < MinError; Int =:= 0),N>MaxSamples)) ->
  400    Samples is N,
  401    Lower is P - Semi,
  402    Prob is P,
  403    Upper is P + Semi %,
  404%   writeln(Semi)
  405   ;
  406     montecarlo_cycle(N, S, M:Goals, K, MinError, Samples, Lower, Prob, Upper)
  407   ).
  408
  409montecarlo(0,N,S , _Goals,N,S):-!.
  410
  411montecarlo(K1,Count, Success, M:Goals,N1,S1):-
  412  erase_samples(M),
  413  copy_term(Goals,Goals1),
  414  (M:Goals1->
  415    Valid=1
  416  ;
  417    Valid=0
  418  ),
  419  N is Count + 1,
  420  S is Success + Valid,
  421  %format("Sample ~d Valid ~d~n",[N,Valid]),
  422  %flush_output,
  423  K2 is K1-1,
  424  montecarlo(K2,N, S, M:Goals, N1,S1).
  425
  426
  427rejection_montecarlo(0,N,S , _Goals,_Ev,N,S):-!.
  428
  429rejection_montecarlo(K1,Count, Success, M:Goals,M:Ev,N1,S1):-
  430  erase_samples(M),
  431  copy_term(Ev,Ev1),
  432  (M:Ev1->
  433    copy_term(Goals,Goals1),
  434    (M:Goals1->
  435      Succ=1
  436    ;
  437      Succ=0
  438    ),
  439    N is Count + 1,
  440    S is Success + Succ,
  441  %format("Sample ~d Valid ~d~n",[N,Valid]),
  442  %flush_output,
  443    K2 is K1-1
  444  ;
  445    N = Count,
  446    S = Success,
  447    K2 = K1
  448  ),
  449  rejection_montecarlo(K2,N, S, M:Goals,M:Ev, N1,S1).
  450
  451mh_montecarlo(_L,K,_NC0,N,S,Succ0,Succ0, _Goals,_Ev,N,S):-
  452  K=<0,!.
  453
  454mh_montecarlo(L,K0,NC0,N0, S0,Succ0, SuccNew,M:Goal, M:Evidence, N, S):-
  455  resample(M,L),
  456  copy_term(Evidence,Ev1),
  457  (M:Ev1->
  458    copy_term(Goal,Goal1),
  459    (M:Goal1->
  460      Succ1=1
  461    ;
  462      Succ1=0
  463    ),
  464    count_samples(M,NC1),
  465    (accept(NC0,NC1)->
  466      Succ = Succ1,
  467      delete_samples_copy(M,Goal),
  468      save_samples_copy(M,Goal)
  469    ;
  470      Succ = Succ0,
  471      erase_samples(M),
  472      restore_samples(M,Goal)
  473    ),
  474    N1 is N0 + 1,
  475    S1 is S0 + Succ,
  476  %format("Sample ~d Valid ~d~n",[N,Valid]),
  477  %flush_output,
  478    K1 is K0-1
  479  ;
  480    NC1 = NC0,
  481    Succ = Succ0,
  482    N1 is N0 + 1,
  483    K1 is K0-1,
  484    S1 is S0 + Succ,
  485    erase_samples(M),
  486    restore_samples(M,Goal)
  487  ),
  488  mh_montecarlo(L,K1,NC1,N1, S1,Succ, SuccNew,M:Goal,M:Evidence, N,S).
  489
  490accept(NC1,NC2):-
  491  P is min(1,NC1/NC2),
  492  random(P0),
  493  P>P0.
 mc_prob(:Query:atom, -Probability:float, +Options:list) is det
The predicate computes the probability of the query Query If Query is not ground, it considers it as an existential query and returns the probability that there is a satisfying assignment to the query.

Options is a list of options, the following are recognised by mc_prob/3:

bar(-BarChart:dict)
BarChart is a dict for rendering with c3 as a bar chart with a bar for the probability of success and a bar for the probability of failure. /
  508mc_prob(M:Goal,P,Options):-
  509  must_be(nonvar,Goal),
  510  must_be(var,P),
  511  must_be(list,Options),
  512  s(M:Goal,P),
  513  option(bar(Chart),Options,no),
  514  (nonvar(Chart)->
  515    true
  516  ;
  517    bar(P,Chart)
  518  ).
 mc_prob(:Query:conjunction_of_literals, -Probability:float) is det
Equivalent to mc_prob/2 with an empty option list. /
  524mc_prob(M:Goal,P):-
  525  must_be(nonvar,Goal),
  526  must_be(var,P),
  527  mc_prob(M:Goal,P,[]).
 mc_sample(:Query:conjunction_of_literals, +Samples:int, -Probability:float, +Options:list) is det
The predicate samples Query a number of Samples times and returns the resulting Probability (Successes/Samples) If Query is not ground, it considers it as an existential query

Options is a list of options, the following are recognised by mc_sample/4:

successes(-Successes:int)
Number of successes
failures(-Failures:int)
Number of failures
bar(-BarChart:dict)
BarChart is a dict for rendering with c3 as a bar chart with a bar for the number of successes and a bar for the number of failures. /
  545mc_sample(M:Goal,S,P,Options):-
  546  must_be(nonvar,Goal),
  547  must_be(nonneg,S),
  548  must_be(var,P),
  549  must_be(list,Options),
  550  option(successes(T),Options,_T),
  551  option(failures(F),Options,_F),
  552  mc_sample(M:Goal,S,T,F,P),
  553  option(bar(Chart),Options,no),
  554  (nonvar(Chart)->
  555    true
  556  ;
  557    bar(T,F,Chart)
  558  ).
 mc_sample(:Query:conjunction_of_literals, +Samples:int, -Probability:float) is det
Equivalent to mc_sample/4 with an empty option list. /
  565mc_sample(M:Goal,S,P):-
  566  must_be(nonvar,Goal),
  567  must_be(nonneg,S),
  568  must_be(var,P),
  569  mc_sample(M:Goal,S,P,[]).
 mc_sample(:Query:conjunction_of_literals, +Samples:int, -Successes:int, -Failures:int, -Probability:float) is det
The predicate samples Query a number of Samples times and returns the number of Successes, of Failures and the Probability (Successes/Samples) If Query is not ground, it considers it as an existential query /
  579mc_sample(M:Goal,S,T,F,P):-
  580  must_be(nonvar,Goal),
  581  must_be(nonneg,S),
  582  must_be(var,T),
  583  must_be(var,F),
  584  must_be(var,P),
  585  copy_term(Goal,Goal1),
  586  numbervars(Goal1),
  587  save_samples(M,Goal1),
  588  montecarlo(S,0, 0, M:Goal, N, T),
  589  P is T / N,
  590  F is N - T,
  591  erase_samples(M),
  592  restore_samples_delete_copy(M,Goal1).
 mc_rejection_sample(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, -Probability:float, +Options:list) is det
The predicate samples Query a number of Samples times given that Evidence is true and returns the Probability of Query. It performs rejection sampling: if in a sample Evidence is false, the sample is discarded. If Query/Evidence are not ground, it considers them an existential queries.

Options is a list of options, the following are recognised by mc_rejection_sample/5:

successes(-Successes:int)
Number of succeses
failures(-Failures:int)
Number of failueres /
  610mc_rejection_sample(M:Goal,M:Evidence,S,P,Options):-
  611  must_be(nonvar,Goal),
  612  must_be(nonvar,Evidence),
  613  must_be(nonneg,S),
  614  must_be(var,P),
  615  must_be(list,Options),
  616  option(successes(T),Options,_T),
  617  option(failures(F),Options,_F),
  618  mc_rejection_sample(M:Goal,M:Evidence,S,T,F,P).
 mc_rejection_sample(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, -Probability:float) is det
Equivalent to mc_rejection_sample/5 with an empty option list. /
  625mc_rejection_sample(M:Goal,M:Evidence,S,P):-
  626  must_be(nonvar,Goal),
  627  must_be(nonvar,Evidence),
  628  must_be(nonneg,S),
  629  must_be(var,P),
  630  mc_rejection_sample(M:Goal,M:Evidence,S,P,[]).
 mc_rejection_sample(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, -Successes:int, -Failures:int, -Probability:float) is det
The predicate samples Query a number of Samples times given that Evidence is true and returns the number of Successes, of Failures and the Probability (Successes/Samples). It performs rejection sampling: if in a sample Evidence is false, the sample is discarded. If Query/Evidence are not ground, it considers them an existential queries. /
  642mc_rejection_sample(M:Goal,M:Evidence0,S,T,F,P):-
  643  must_be(nonvar,Goal),
  644  must_be(nonvar,Evidence0),
  645  must_be(nonneg,S),
  646  must_be(var,T),
  647  must_be(var,F),
  648  must_be(var,P),
  649  test_prism(M),
  650  deal_with_ev(Evidence0,M,Evidence,UpdatedClausesRefs,ClausesToReAdd),
  651  rejection_montecarlo(S,0, 0, M:Goal,M:Evidence, N, T),
  652  P is T / N,
  653  F is N - T,
  654  erase_samples(M),
  655  maplist(erase,UpdatedClausesRefs),
  656  maplist(M:assertz,ClausesToReAdd).
  657
  658deal_with_ev(Ev,M,EvGoal,UC,CA):-
  659  list2and(EvL,Ev),
  660  partition(ac,EvL,ActL,EvNoActL),
  661  deal_with_actions(ActL,M,UC,CA),
  662  list2and(EvNoActL,EvGoal).
  663
  664deal_with_actions(ActL,M,UC,CA):-
  665  empty_assoc(AP0),
  666  foldl(get_pred_const,ActL,AP0,AP),
  667  assoc_to_list(AP,LP),
  668  maplist(update_clauses(M),LP,UCL,CAL),
  669  partition(nac,ActL,_NActL,PActL),
  670  maplist(assert_actions(M),PActL,ActRefs),
  671  append([ActRefs|UCL],UC),
  672  append(CAL,CA).
  673
  674assert_actions(M,do(A),Ref):-
  675  M:assertz(A,Ref).
  676
  677update_clauses(M,P/0- _,[],LCA):-!,
  678  findall(Ref,M:clause(P,_B,Ref),UC),
  679  findall((P:-B),M:clause(P,B),LCA),
  680  maplist(erase,UC).
  681
  682update_clauses(M,P/A-Constants,UC,CA):-
  683  functor(G,P,A),
  684  G=..[_|Args],
  685  findall((G,B,Ref),M:clause(G,B,Ref),LC),
  686  maplist(get_const(Args),Constants,ConstraintsL),
  687  list2and(ConstraintsL,Constraints),
  688  maplist(add_cons(G,Constraints,M),LC,UC,CA).
  689
  690add_cons(G,C,M,(H,B,Ref),Ref1,(H:-B)):-
  691  copy_term((G,C),(G1,C1)),
  692  G1=H,
  693  erase(Ref),
  694  M:assertz((H:-(C1,B)),Ref1).
  695
  696
  697get_const(Args,Constants,Constraint):-
  698  maplist(constr,Args,Constants,ConstraintL),
  699  list2and(ConstraintL,Constraint).
  700
  701constr(V,C,dif(V,C)).
  702
  703get_pred_const(do(Do0),AP0,AP):-
  704  (Do0= (\+ Do)->
  705    true
  706  ;
  707    Do=Do0
  708  ),
  709  functor(Do,F,A),
  710  Do=..[_|Args],
  711  (get_assoc(F/A,AP0,V)->
  712    put_assoc(F/A,AP0,[Args|V],AP)
  713  ;
  714    put_assoc(F/A,AP0,[Args],AP)
  715  ).
  716
  717
  718ac(do(_)).
  719nac(do(\+ _)).
 mc_gibbs_sample(:Query:conjunction_of_literals, +Samples:int, -Probability:float, +Options:list) is det
The predicate samples Query a number of Mix+Samples (Mix is set with the options, default value 0) times. The first Mix (that is set with the options, default value 0) samples are discarded (mixing time). It performs Gibbs sampling: each sample is obtained from the previous one by resampling a variable given the values of the variables in its Markov blanket. If Query/Evidence are not ground, it considers them as existential queries.

Options is a list of options, the following are recognised by mc_gibbs_sample/4:

block +Block:int
Perform blocked Gibbs: Block variables are sampled together, default value 1
mix(+Mix:int)
The first Mix samples are discarded (mixing time), default value 0
successes(-Successes:int)
Number of succeses
failures(-Failures:int)
Number of failueres /
  741mc_gibbs_sample(M:Goal,S,P,Options):-
  742  must_be(nonvar,Goal),
  743  must_be(nonneg,S),
  744  must_be(var,P),
  745  must_be(list,Options),
  746  option(mix(Mix),Options,0),
  747  option(block(Block),Options,1),
  748  option(successes(T),Options,_T),
  749  option(failures(F),Options,_F),
  750  mc_gibbs_sample(M:Goal,S,Mix,Block,T,F,P).
 mc_gibbs_sample(:Query:conjunction_of_literals, +Samples:int, -Probability:float) is det
Equivalent to mc_gibbs_sample/4 with an empty option list. /
  757mc_gibbs_sample(M:Goal,S,P):-
  758  must_be(nonvar,Goal),
  759  must_be(nonneg,S),
  760  must_be(var,P),
  761  mc_gibbs_sample(M:Goal,S,P,[]).
  762
  763mc_gibbs_sample(M:Goal,S,Mix,Block,T,F,P):-
  764  copy_term(Goal,Goal1),
  765  (M:Goal1->
  766    Succ=1
  767  ;
  768    Succ=0
  769  ),
  770  (Mix=0->
  771    T1=Succ,
  772    S1 is S-1
  773  ;
  774    T1=0,
  775    S1=S,
  776    Mix1 is Mix-1,
  777    gibbs_montecarlo(Mix1,0,Block,M:Goal,_TMix)
  778  ),
  779  gibbs_montecarlo(S1,T1,Block,M:Goal,T),
  780  P is T / S,
  781  F is S - T,
  782  erase_samples(M).
  783
  784gibbs_montecarlo(K,T,_Block,_Goals,T):-
  785  K=<0,!.
  786
  787gibbs_montecarlo(K0, T0,Block,M:Goal,T):-
  788  remove_samples(M,Block,LS),
  789  copy_term(Goal,Goal1),
  790  (M:Goal1->
  791    Succ=1
  792  ;
  793    Succ=0
  794  ),
  795  T1 is T0 + Succ,
  796  K1 is K0-1,
  797  check_sampled(M,LS),
  798  gibbs_montecarlo(K1,T1,Block,M:Goal,T).
  799
  800remove_samples(M,Block,Samp):-
  801  remove_samp(M,Block,Samp).
  802
  803remove_samp(_M,0,[]):-!.
  804
  805remove_samp(M,Block0,[(R,S)|Samp]):-
  806  retract(M:sampled(R,S,_)),!,
  807  Block is Block0-1,
  808  remove_samp(M,Block,Samp).
  809
  810remove_samp(_M,_Block,[]).
  811
  812
  813% check_sampled(M,R,S):-
  814%   M:sampled(R,S,_),!.
  815
  816check_sampled(M,S):-
  817  maplist(check_sam(M),S).
  818
  819check_sam(M,(R,S)):-
  820  M:samp(R,S,_).
 mc_gibbs_sample(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, -Probability:float, +Options:list) is det
The predicate samples Query a number of Mix+Samples (Mix is set with the options, default value 0) times given that Evidence is true and returns the number of Successes, of Failures and the Probability (Successes/Samples). The first Mix (that is set with the options, default value 0) samples are discarded (mixing time). It performs Gibbs sampling: each sample is obtained from the previous one by resampling a variable given the values of the variables in its Markov blanket. If Query/Evidence are not ground, it considers them as existential queries.

Options is a list of options, the following are recognised by mc_gibbs_sample/5:

block +Block:int
Perform blocked Gibbs: Block variables are sampled together, default value 1
mix(+Mix:int)
The first Mix samples are discarded (mixing time), default value 0
successes(-Successes:int)
Number of succeses
failures(-Failures:int)
Number of failueres /
  844mc_gibbs_sample(M:Goal,M:Evidence,S,P,Options):-
  845  must_be(nonvar,Goal),
  846  must_be(nonvar,Evidence),
  847  must_be(nonneg,S),
  848  must_be(var,P),
  849  must_be(list,Options),
  850  test_prism(M),
  851  option(mix(Mix),Options,0),
  852  option(block(Block),Options,1),
  853  option(successes(T),Options,_T),
  854  option(failures(F),Options,_F),
  855  mc_gibbs_sample(M:Goal,M:Evidence,S,Mix,Block,T,F,P).
  856
  857
  858mc_gibbs_sample(M:Goal,M:Evidence0,S,Mix,Block,T,F,P):-
  859  deal_with_ev(Evidence0,M,Evidence,UpdatedClausesRefs,ClausesToReAdd),
  860  gibbs_sample_cycle(M:Evidence),
  861  copy_term(Goal,Goal1),
  862  (M:Goal1->
  863    Succ=1
  864  ;
  865    Succ=0
  866  ),
  867  (Mix=0->
  868    T1=Succ,
  869    S1 is S-1
  870  ;
  871    T1=0,
  872    S1=S,
  873    Mix1 is Mix-1,
  874    gibbs_montecarlo(Mix1,0,Block,M:Goal,M:Evidence,_TMix)
  875  ),
  876  gibbs_montecarlo(S1,T1,Block,M:Goal,M:Evidence,T),
  877  P is T / S,
  878  F is S - T,
  879  erase_samples(M),
  880  maplist(erase,UpdatedClausesRefs),
  881  maplist(M:assertz,ClausesToReAdd).
  882
  883gibbs_montecarlo(K,T,_Block,_G,_Ev,T):-
  884  K=<0,!.
  885
  886gibbs_montecarlo(K0, T0,Block, M:Goal, M:Evidence,  T):-
  887  remove_samples(M,Block,LS),
  888  save_samples_copy(M,Evidence),
  889  gibbs_sample_cycle(M:Evidence),
  890  delete_samples_copy(M,Evidence),
  891  copy_term(Goal,Goal1),
  892  (M:Goal1->
  893    Succ=1
  894  ;
  895    Succ=0
  896  ),
  897  T1 is T0 + Succ,
  898  K1 is K0-1,
  899  check_sampled(M,LS),
  900  gibbs_montecarlo(K1, T1,Block,M:Goal,M:Evidence, T).
 mc_gibbs_sample_arg(:Query:conjunction_of_literals, +Samples:int, ?Arg:var, -Values:list, +Options:list) is det
The predicate samples Query a number of Samples times. Arg should be a variable in Query. The predicate returns in Values a list of couples L-N where L is the list of values of Arg for which Query succeeds in a world sampled at random and N is the number of samples returning that list of values. The first Mix (that is set with the options, default value 0) samples are discarded (mixing time). It performs Gibbs sampling: each sample is obtained from the previous one by resampling a variable given the values of the variables in its Markov blanket.

Options is a list of options, the following are recognised by mc_gibbs_sample_arg/5:

block +Block:int
Perform blocked Gibbs: Block variables are sampled together, default value 1
mix(+Mix:int)
The first Mix samples are discarded (mixing time), default value 0
bar(-BarChar:dict)
BarChart is a dict for rendering with c3 as a bar chart with a bar for each possible value of L, the list of value of Arg for which Query succeeds in a world sampled at random. /
  927mc_gibbs_sample_arg(M:Goal,S,Arg,ValList,Options):-
  928  must_be(nonvar,Goal),
  929  must_be(nonneg,S),
  930  must_be(var,Arg),
  931  must_be(var,ValList),
  932  must_be(list,Options),
  933  test_prism(M),
  934  option(mix(Mix),Options,0),
  935  option(block(Block),Options,1),
  936  option(bar(Chart),Options,no),
  937  mc_gibbs_sample_arg0(M:Goal,S,Mix,Block,Arg,ValList),
  938  (nonvar(Chart)->
  939    true
  940  ;
  941    argbar(ValList,Chart)
  942  ).
 mc_gibbs_sample_arg(:Query:conjunction_of_literals, +Samples:int, ?Arg:var, -Values:list) is det
Equivalent to mc_gibbs_sample_arg/5 with an empty option list. /
  948mc_gibbs_sample_arg(M:Goal,S,Arg,ValList):-
  949  must_be(nonvar,Goal),
  950  must_be(nonneg,S),
  951  must_be(var,Arg),
  952  must_be(var,ValList),
  953  mc_gibbs_sample_arg(M:Goal,S,Arg,ValList,[]).
 mc_gibbs_sample_arg(:Query:conjunction_of_literals, :Evidence:conjunction_of_groundliterals, +Samples:int, ?Arg:var, -Values:list, +Options:list) is det
The predicate samples Query a number of Samples times given that Evidence is true. Arg should be a variable in Query. The predicate returns in Values a list of couples L-N where L is the list of values of Arg for which Query succeeds in a world sampled at random and N is the number of samples returning that list of values. The first Mix (that is set with the options, default value 0) samples are discarded (mixing time). It performs Gibbs sampling: each sample is obtained from the previous one by resampling a variable given the values of the variables in its Markov blanket.

Options is a list of options, the following are recognised by mc_gibbs_sample_arg/6:

block +Block:int
Perform blocked Gibbs: Block variables are sampled together, default value 1
mix(+Mix:int)
The first Mix samples are discarded (mixing time), default value 0
bar(-BarChar:dict)
BarChart is a dict for rendering with c3 as a bar chart with a bar for each possible value of L, the list of value of Arg for which Query succeeds in a world sampled at random. /
  979mc_gibbs_sample_arg(M:Goal,M:Evidence,S,Arg,ValList,Options):-
  980  must_be(nonvar,Goal),
  981  must_be(nonvar,Evidence),
  982  must_be(nonneg,S),
  983  must_be(var,Arg),
  984  must_be(var,ValList),
  985  must_be(list,Options),
  986  test_prism(M),
  987  option(mix(Mix),Options,0),
  988  option(block(Block),Options,1),
  989  option(bar(Chart),Options,no),
  990  mc_gibbs_sample_arg0(M:Goal,M:Evidence,S,Mix,Block,Arg,ValList),
  991  (nonvar(Chart)->
  992    true
  993  ;
  994    argbar(ValList,Chart)
  995  ).
  996
  997
  998
  999
 1000gibbs_sample_arg(K,_Goals,_Ev,_Block,_Arg,V,V):-
 1001  K=<0,!.
 1002
 1003gibbs_sample_arg(K0,M:Goal, M:Evidence,Block, Arg,V0,V):-
 1004  remove_samples(M,Block,LS),
 1005  save_samples_copy(M,Evidence),
 1006  gibbs_sample_cycle(M:Evidence),
 1007  delete_samples_copy(M,Evidence),
 1008  findall(Arg,M:Goal,La),
 1009  numbervars(La),
 1010  (get_assoc(La, V0, N)->
 1011    N1 is N+1,
 1012    put_assoc(La,V0,N1,V1)
 1013  ;
 1014    put_assoc(La,V0,1,V1)
 1015  ),
 1016  K1 is K0-1,
 1017  check_sampled(M,LS),
 1018  gibbs_sample_arg(K1,M:Goal,M:Evidence,Block,Arg,V1,V).
 mc_gibbs_sample_arg0(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, +Mix:int, +Block:int, +Lag:int, ?Arg:var, -Values:list) is det
The predicate samples Query a number of Samples times given that Evidence is true. Arg should be a variable in Query. The predicate returns in Values a list of couples L-N where L is the list of values of Arg for which Query succeeds in a world sampled at random and N is the number of samples returning that list of values. It performs blocked Gibbs sampling: each sample is obtained from the previous one by resampling Block variables given the values of the variables in its Markov blanket. /
 1035mc_gibbs_sample_arg0(M:Goal,M:Evidence0,S,Mix,Block,Arg,ValList):-
 1036  deal_with_ev(Evidence0,M,Evidence,UpdatedClausesRefs,ClausesToReAdd),
 1037  gibbs_sample_cycle(M:Evidence),
 1038  empty_assoc(Values0),
 1039  findall(Arg,M:Goal,La),
 1040  numbervars(La),
 1041  put_assoc(La,Values0,1,Values1),
 1042  (Mix=0->
 1043    Values2=Values1,
 1044    S1 is S-1
 1045  ;
 1046    Mix1 is Mix-1,
 1047    gibbs_sample_arg(Mix1,M:Goal,M:Evidence,Block,Arg, Values1,_Values),
 1048    S1=S,
 1049    Values2=Values0
 1050  ),
 1051  gibbs_sample_arg(S1,M:Goal,M:Evidence,Block,Arg, Values2,Values),
 1052  erase_samples(M),
 1053  assoc_to_list(Values,ValList0),
 1054  sort(2, @>=,ValList0,ValList),
 1055  maplist(erase,UpdatedClausesRefs),
 1056  maplist(M:assertz,ClausesToReAdd).
 mc_gibbs_sample_arg0(:Query:conjunction_of_literals, +Samples:int, +Mix:int, +Block:int, ?Arg:var, -Values:list) is det
The predicate samples Query a number of Samples times. Arg should be a variable in Query. The predicate returns in Values a list of couples L-N where L is the list of values of Arg for which Query succeeds in a world sampled at random and N is the number of samples returning that list of values. It performs blocked Gibbs sampling: each sample is obtained from the previous one by resampling Block variables given the values of the variables in its Markov blanket. /
 1071mc_gibbs_sample_arg0(M:Goal,S,Mix,Block,Arg,ValList):-
 1072  empty_assoc(Values0),
 1073  findall(Arg,M:Goal,La),
 1074  numbervars(La),
 1075  put_assoc(La,Values0,1,Values1),
 1076  (Mix=0->
 1077    Values2=Values1,
 1078    S1 is S-1
 1079  ;
 1080    Mix1 is Mix-1,
 1081    gibbs_sample_arg(Mix1,M:Goal,Block,Arg, Values1,_Values),
 1082    S1=S,
 1083    Values2=Values0
 1084  ),
 1085  gibbs_sample_arg(S1,M:Goal,Block,Arg, Values2,Values),
 1086  erase_samples(M),
 1087  assoc_to_list(Values,ValList0),
 1088  sort(2, @>=,ValList0,ValList).
 1089
 1090gibbs_sample_arg(K,_Goals,_Block,_Arg,V,V):-
 1091  K=<0,!.
 1092
 1093gibbs_sample_arg(K0,M:Goal,Block, Arg,V0,V):-
 1094  remove_samples(M,Block,LS),
 1095  findall(Arg,M:Goal,La),
 1096  numbervars(La),
 1097  (get_assoc(La, V0, N)->
 1098    N1 is N+1,
 1099    put_assoc(La,V0,N1,V1)
 1100  ;
 1101    put_assoc(La,V0,1,V1)
 1102  ),
 1103  check_sampled(M,LS),
 1104  K1 is K0-1,
 1105  gibbs_sample_arg(K1,M:Goal,Block,Arg,V1,V).
 mc_mh_sample(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, -Probability:float, +Options:list) is det
The predicate samples Query a number of Mix+Samples (Mix is set with the options, default value 0) times given that Evidence is true and returns the number of Successes, of Failures and the Probability (Successes/Samples). The first Mix (that is set with the options, default value 0) samples are discarded (mixing time). It performs Metropolis/Hastings sampling: between each sample, Lag (that is set with the options, default value 1) sampled choices are forgotten and each sample is accepted with a certain probability. If Query/Evidence are not ground, it considers them as existential queries.

Options is a list of options, the following are recognised by mc_mh_sample/5:

mix(+Mix:int)
The first Mix samples are discarded (mixing time), default value 0
lag(+Lag:int)
lag between each sample, Lag sampled choices are forgotten, default value 1
successes(-Successes:int)
Number of succeses
failures(-Failures:int)
Number of failueres /
 1130mc_mh_sample(M:Goal,M:Evidence,S,P,Options):-
 1131  must_be(nonvar,Goal),
 1132  must_be(nonneg,S),
 1133  must_be(var,P),
 1134  must_be(list,Options),
 1135  test_prism(M),
 1136  option(lag(L),Options,1),
 1137  option(mix(Mix),Options,0),
 1138  option(successes(T),Options,_T),
 1139  option(failures(F),Options,_F),
 1140  mc_mh_sample(M:Goal,M:Evidence,S,Mix,L,T,F,P).
 mc_mh_sample(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, -Probability:float) is det
Equivalent to mc_mh_sample/5 with an empty option list. /
 1147mc_mh_sample(M:Goal,M:Evidence,S,P):-
 1148  must_be(nonvar,Goal),
 1149  must_be(nonvar,Evidence),
 1150  must_be(nonneg,S),
 1151  must_be(var,P),
 1152  mc_mh_sample(M:Goal,M:Evidence,S,P,[]).
 mc_mh_sample(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, +Mix:int, +Lag:int, -Successes:int, -Failures:int, -Probability:float) is det
The predicate samples Query a number of Mix+Samples times given that Evidence is true and returns the number of Successes, of Failures and the Probability (Successes/Samples). The first Mix samples are discarded (mixing time). It performs Metropolis/Hastings sampling: between each sample, Lag sampled choices are forgotten and each sample is accepted with a certain probability. If Query/Evidence are not ground, it considers them as existential queries. /
 1167mc_mh_sample(M:Goal,M:Evidence0,S,Mix,L,T,F,P):-
 1168  must_be(nonvar,Goal),
 1169  must_be(nonvar,Evidence0),
 1170  must_be(nonneg,S),
 1171  must_be(nonneg,Mix),
 1172  must_be(nonneg,L),
 1173  must_be(var,T),
 1174  must_be(var,F),
 1175  must_be(var,P),
 1176  deal_with_ev(Evidence0,M,Evidence,UpdatedClausesRefs,ClausesToReAdd),
 1177  initial_sample_cycle(M:Evidence),!,
 1178  copy_term(Goal,Goal1),
 1179  (M:Goal1->
 1180    Succ=1
 1181  ;
 1182    Succ=0
 1183  ),
 1184  count_samples(M,NC),
 1185  Mix1 is Mix-1,
 1186  save_samples_copy(M,Goal),
 1187  mh_montecarlo(L,Mix1,NC,0, Succ,Succ,Succ1,M:Goal, M:Evidence, _NMix, _TMix),
 1188  count_samples(M,NC1),
 1189  mh_montecarlo(L,S,NC1,0, 0,Succ1,_Succ1, M:Goal, M:Evidence, _N, T),
 1190  P is T / S,
 1191  F is S - T,
 1192  erase_samples(M),
 1193  delete_samples_copy(M,Goal),
 1194  maplist(erase,UpdatedClausesRefs),
 1195  maplist(M:assertz,ClausesToReAdd).
 1196
 1197
 1198gibbs_sample_cycle(M:G):-
 1199  copy_term(G,G1),
 1200  (M:G1->
 1201    true
 1202  ;
 1203    erase_samples(M),
 1204    restore_samples(M,G),
 1205    gibbs_sample_cycle(M:G)
 1206  ).
 1207
 1208
 1209initial_sample_cycle(M:G):-
 1210  copy_term(G,G1),
 1211  (initial_sample(M:G1)->
 1212    true
 1213  ;
 1214    erase_samples(M),
 1215    initial_sample_cycle(M:G)
 1216  ).
 1217
 1218initial_sample(_M:true):-!.
 1219
 1220initial_sample(M:(A~= B)):-!,
 1221  add_arg(A,B,A1),
 1222  initial_sample(M:A1).
 1223
 1224initial_sample(M:msw(A,B)):-!,
 1225  msw(M:A,B).
 1226
 1227initial_sample(_M:(sample_head(R,VC,M,HL,NH))):-!,
 1228  sample_head(R,VC,M,HL,NH).
 1229
 1230initial_sample(_M:take_a_sample(R,VC,M,Distr,S)):-!,
 1231  take_a_sample(R,VC,M,Distr,S).
 1232
 1233initial_sample(M:(G1,G2)):-!,
 1234  initial_sample(M:G1),
 1235  initial_sample(M:G2).
 1236
 1237initial_sample(M:(G1;G2)):-!,
 1238  initial_sample(M:G1);
 1239  initial_sample(M:G2).
 1240
 1241initial_sample(M:(\+ G)):-!,
 1242  \+ initial_sample(M:G).
 1243%  initial_sample_neg(M:G).
 1244
 1245initial_sample(M:findall(A,G,L)):-!,
 1246  findall(A,initial_sample(M:G),L).
 1247
 1248initial_sample(M:G):-
 1249  builtin(G),!,
 1250  M:call(G).
 1251
 1252initial_sample(M:G):-
 1253  findall((G,B),M:clause(G,B),L),
 1254  sample_one_back(L,(G,B)),
 1255  initial_sample(M:B).
 1256
 1257initial_sample_neg(_M:true):-!,
 1258  fail.
 1259
 1260initial_sample_neg(_M:(sample_head(R,VC,M,HL,N))):-!,
 1261  sample_head(R,VC,M,HL,NH),
 1262  NH\=N.
 1263
 1264initial_sample_neg(_M:take_a_sample(R,VC,M,Distr,S)):-!,
 1265  take_a_sample(R,VC,M,Distr,S).
 1266
 1267
 1268initial_sample_neg(M:(G1,G2)):-!,
 1269  (initial_sample_neg(M:G1),!;
 1270  initial_sample(M:G1),
 1271  initial_sample_neg(M:G2)).
 1272
 1273initial_sample_neg(M:(G1;G2)):-!,
 1274  initial_sample_neg(M:G1),
 1275  initial_sample_neg(M:G2).
 1276
 1277initial_sample_neg(M:(\+ G)):-!,
 1278  initial_sample(M:G).
 1279
 1280initial_sample_neg(M:G):-
 1281  builtin(G),!,
 1282  \+ M:call(G).
 1283
 1284initial_sample_neg(M:G):-
 1285  findall(B,M:clause(G,B),L),
 1286  initial_sample_neg_all(L,M).
 1287
 1288initial_sample_neg_all([],_M).
 1289
 1290initial_sample_neg_all([H|T],M):-
 1291  initial_sample_neg(M:H),
 1292  initial_sample_neg_all(T,M).
 1293
 1294check(R,VC,M,N):-
 1295  M:sampled(R,VC,N),!.
 1296
 1297check(R,VC,M,N):-
 1298  \+ M:sampled(R,VC,_N),
 1299  assertz(M:sampled(R,VC,N)).
 1300
 1301check_neg(R,VC,M,_LN,N):-
 1302  M:sampled(R,VC,N1),!,
 1303  N\=N1.
 1304
 1305check_neg(R,VC,M,LN,N):-
 1306  \+ M:sampled(R,VC,_N),
 1307  member(N1,LN),
 1308  N1\= N,
 1309  assertz(M:sampled(R,VC,N1)).
 1310
 1311listN(0,[]):-!.
 1312
 1313listN(N,[N1|T]):-
 1314  N1 is N-1,
 1315  listN(N1,T).
 mc_sample_arg(:Query:conjunction_of_literals, +Samples:int, ?Arg:var, -Values:list, +Options:list) is det
The predicate samples Query a number of Samples times. Arg should be a variable in Query. The predicate returns in Values a list of couples L-N where L is the list of values of Arg for which Query succeeds in a world sampled at random and N is the number of samples returning that list of values.

Options is a list of options, the following are recognised by mc_sample_arg/5:

bar(-BarChar:dict)
BarChart is a dict for rendering with c3 as a bar chart with a bar for each possible value of L, the list of value of Arg for which Query succeeds in a world sampled at random. /
 1334mc_sample_arg(M:Goal,S,Arg,ValList,Options):-
 1335  must_be(nonvar,Goal),
 1336  must_be(nonneg,S),
 1337  must_be(var,Arg),
 1338  must_be(var,ValList),
 1339  must_be(list,Options),
 1340  empty_assoc(Values0),
 1341  sample_arg(S,M:Goal,Arg, Values0,Values),
 1342  erase_samples(M),
 1343  assoc_to_list(Values,ValList0),
 1344  sort(2, @>=,ValList0,ValList),
 1345  option(bar(Chart),Options,no),
 1346  (nonvar(Chart)->
 1347    true
 1348  ;
 1349    argbar(ValList,Chart)
 1350  ).
 mc_sample_arg(:Query:conjunction_of_literals, +Samples:int, ?Arg:var, -Values:list) is det
Equivalent to mc_sample_arg/5 with an empty option list. /
 1356mc_sample_arg(M:Goal,S,Arg,ValList):-
 1357  must_be(nonvar,Goal),
 1358  must_be(nonneg,S),
 1359  must_be(var,Arg),
 1360  must_be(var,ValList),
 1361  mc_sample_arg(M:Goal,S,Arg,ValList,[]).
 mc_rejection_sample_arg(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, ?Arg:var, -Values:list, +Options:list) is det
The predicate samples Query a number of Samples times given that Evidence is true. Arg should be a variable in Query. The predicate returns in Values a list of couples L-N where L is the list of values of Arg for which Query succeeds in a world sampled at random and N is the number of samples returning that list of values. Rejection sampling is performed.

Options is a list of options, the following are recognised by mc_rejection_sample_arg/6:

bar(-BarChar:dict)
BarChart is a dict for rendering with c3 as a bar chart with a bar for each possible value of L, the list of value of Arg for which Query succeeds in a world sampled at random. /
 1382mc_rejection_sample_arg(M:Goal,M:Evidence0,S,Arg,ValList,Options):-
 1383  must_be(nonvar,Goal),
 1384  must_be(nonvar,Evidence0),
 1385  must_be(nonneg,S),
 1386  must_be(var,Arg),
 1387  must_be(var,ValList),
 1388  must_be(list,Options),
 1389  test_prism(M),
 1390  deal_with_ev(Evidence0,M,Evidence,UpdatedClausesRefs,ClausesToReAdd),
 1391  empty_assoc(Values0),
 1392  rejection_sample_arg(S,M:Goal,M:Evidence,Arg, Values0,Values),
 1393  erase_samples(M),
 1394  assoc_to_list(Values,ValList0),
 1395  sort(2, @>=,ValList0,ValList),
 1396  maplist(erase,UpdatedClausesRefs),
 1397  maplist(M:assertz,ClausesToReAdd),
 1398  option(bar(Chart),Options,no),
 1399  (nonvar(Chart)->
 1400    true
 1401  ;
 1402    argbar(ValList,Chart)
 1403  ).
 mc_rejection_sample_arg(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, ?Arg:var, -Values:list) is det
Equivalent to mc_rejection_sample_arg/6 with an empty option list. /
 1410mc_rejection_sample_arg(M:Goal,M:Evidence0,S,Arg,ValList):-
 1411  must_be(nonvar,Goal),
 1412  must_be(nonvar,Evidence0),
 1413  must_be(nonneg,S),
 1414  must_be(var,Arg),
 1415  must_be(var,ValList),
 1416  mc_rejection_sample_arg(M:Goal,M:Evidence0,S,Arg,ValList,[]).
 mc_mh_sample_arg(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, ?Arg:var, -Values:list, +Options:list) is det
The predicate samples Query a number of Samples times given that Evidence is true. Arg should be a variable in Query. The predicate returns in Values a list of couples L-N where L is the list of values of Arg for which Query succeeds in a world sampled at random and N is the number of samples returning that list of values. The first Mix (that is set with the options, default value 0) samples are discarded (mixing time). It performs Metropolis/Hastings sampling: between each sample, Lag (that is set with the options, default value 1) sampled choices are forgotten and each sample is accepted with a certain probability.

Options is a list of options, the following are recognised by mc_mh_sample_arg/6:

mix(+Mix:int)
The first Mix samples are discarded (mixing time), default value 0
lag(+Lag:int)
lag between each sample, Lag sampled choices are forgotten, default value 1
bar(-BarChar:dict)
BarChart is a dict for rendering with c3 as a bar chart with a bar for each possible value of L, the list of value of Arg for which Query succeeds in a world sampled at random. /
 1443mc_mh_sample_arg(M:Goal,M:Evidence,S,Arg,ValList,Options):-
 1444  must_be(nonvar,Goal),
 1445  must_be(nonvar,Evidence),
 1446  must_be(nonneg,S),
 1447  must_be(var,Arg),
 1448  must_be(var,ValList),
 1449  must_be(list,Options),
 1450  test_prism(M),
 1451  option(mix(Mix),Options,0),
 1452  option(lag(L),Options,1),
 1453  option(bar(Chart),Options,no),
 1454  mc_mh_sample_arg0(M:Goal,M:Evidence,S,Mix,L,Arg,ValList),
 1455  (nonvar(Chart)->
 1456    true
 1457  ;
 1458    argbar(ValList,Chart)
 1459  ).
 mc_mh_sample_arg(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, ?Arg:var, -Values:list) is det
Equivalent to mc_mh_sample_arg/6 with an empty option list. /
 1466mc_mh_sample_arg(M:Goal,M:Evidence,S,Arg,ValList):-
 1467  must_be(nonvar,Goal),
 1468  must_be(nonvar,Evidence),
 1469  must_be(nonneg,S),
 1470  must_be(var,Arg),
 1471  must_be(var,ValList),
 1472  mc_mh_sample_arg(M:Goal,M:Evidence,S,Arg,ValList,[]).
 mc_mh_sample_arg0(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, +Mix:int, +Lag:int, ?Arg:var, -Values:list) is det
The predicate samples Query a number of Samples times given that Evidence is true. Arg should be a variable in Query. The predicate returns in Values a list of couples L-N where L is the list of values of Arg for which Query succeeds in a world sampled at random and N is the number of samples returning that list of values. It performs Metropolis/Hastings sampling: between each sample, Lag sampled choices are forgotten and each sample is accepted with a certain probability. /
 1487mc_mh_sample_arg0(M:Goal,M:Evidence0,S,Mix,L,Arg,ValList):-
 1488  deal_with_ev(Evidence0,M,Evidence,UpdatedClausesRefs,ClausesToReAdd),
 1489  initial_sample_cycle(M:Evidence),!,
 1490  empty_assoc(Values0),
 1491  findall(Arg,M:Goal,La),
 1492  numbervars(La),
 1493  put_assoc(La,Values0,1,Values1),
 1494  count_samples(M,NC),
 1495  Mix1 is Mix-1,
 1496  save_samples_copy(M,Goal),
 1497  mh_sample_arg(L,Mix1,NC,M:Goal,M:Evidence,Arg, La,La1,Values1,_Values),
 1498  count_samples(M,NC1),
 1499  mh_sample_arg(L,S,NC1,M:Goal,M:Evidence,Arg, La1,_La,Values0,Values),
 1500  erase_samples(M),
 1501  assoc_to_list(Values,ValList0),
 1502  sort(2, @>=,ValList0,ValList),
 1503  delete_samples_copy(M,Goal),
 1504  maplist(erase,UpdatedClausesRefs),
 1505  maplist(M:assertz,ClausesToReAdd).
 1506
 1507
 1508
 1509mh_sample_arg(_L,K,_NC0,_Goals,_Ev,_Arg,AP,AP,V,V):-
 1510  K=<0,!.
 1511
 1512mh_sample_arg(L,K0,NC0,M:Goal, M:Evidence, Arg,AP0,AP,V0,V):-
 1513  resample(M,L),
 1514  copy_term(Evidence,Ev1),
 1515  (M:Ev1->
 1516    findall(Arg,M:Goal,La),
 1517    numbervars(La),
 1518    count_samples(M,NC1),
 1519    (accept(NC0,NC1)->
 1520     (get_assoc(La, V0, N)->
 1521        N1 is N+1,
 1522        put_assoc(La,V0,N1,V1)
 1523      ;
 1524        put_assoc(La,V0,1,V1)
 1525      ),
 1526      delete_samples_copy(M,Goal),
 1527      save_samples_copy(M,Goal),
 1528      K1 is K0-1,
 1529      AP1 = La
 1530    ;
 1531      (get_assoc(AP0, V0, N)->
 1532        N1 is N+1,
 1533        put_assoc(AP0,V0,N1,V1)
 1534      ;
 1535        put_assoc(AP0,V0,1,V1)
 1536      ),
 1537      K1 is K0-1,
 1538      AP1=AP0,
 1539      erase_samples(M),
 1540      restore_samples(M,Goal)
 1541    )
 1542  ;
 1543    (get_assoc(AP0, V0, N)->
 1544      N1 is N+1,
 1545      put_assoc(AP0,V0,N1,V1)
 1546    ;
 1547      put_assoc(AP0,V0,1,V1)
 1548    ),
 1549    K1 is K0-1,
 1550    NC1 = NC0,
 1551    AP1=AP0,
 1552    erase_samples(M),
 1553    restore_samples(M,Goal)
 1554  ),
 1555  mh_sample_arg(L,K1,NC1,M:Goal,M:Evidence,Arg,AP1,AP,V1,V).
 1556
 1557
 1558rejection_sample_arg(0,_Goals,_Ev,_Arg,V,V):-!.
 1559
 1560rejection_sample_arg(K1, M:Goals,M:Ev,Arg,V0,V):-
 1561  erase_samples(M),
 1562  copy_term(Ev,Ev1),
 1563  (M:Ev1->
 1564    copy_term((Goals,Arg),(Goals1,Arg1)),
 1565    findall(Arg1,M:Goals1,L),
 1566    numbervars(L),
 1567    (get_assoc(L, V0, N)->
 1568      N1 is N+1,
 1569      put_assoc(L,V0,N1,V1)
 1570    ;
 1571      put_assoc(L,V0,1,V1)
 1572    ),
 1573    K2 is K1-1
 1574  ;
 1575    V1=V0,
 1576    K2=K1
 1577  ),
 1578  rejection_sample_arg(K2,M:Goals,M:Ev,Arg,V1,V).
 1579
 1580sample_arg(0,_Goals,_Arg,V,V):-!.
 1581
 1582sample_arg(K1, M:Goals,Arg,V0,V):-
 1583  erase_samples(M),
 1584  copy_term((Goals,Arg),(Goals1,Arg1)),
 1585  findall(Arg1,M:Goals1,L),
 1586  numbervars(L),
 1587  (get_assoc(L, V0, N)->
 1588    N1 is N+1,
 1589    put_assoc(L,V0,N1,V1)
 1590  ;
 1591    put_assoc(L,V0,1,V1)
 1592  ),
 1593  K2 is K1-1,
 1594  sample_arg(K2,M:Goals,Arg,V1,V).
 mc_particle_sample(:Query:conjunction_of_literals, :Evidence:list, +Samples:int, -Prob:float) is det
The predicate samples Query a number of Samples times given that Evidence is true. Evidence is a list of goals. The predicate returns in Prob the probability that the query is true. It performs particle filtering with likelihood weighting: each sample is weighted by the likelihood of an element of the Evidence list and constitutes a particle. After weighting, particles are resampled and the next element of Evidence is considered. /
 1608mc_particle_sample(M:Goal,M:Evidence,S,P):-
 1609  must_be(nonvar,Goal),
 1610  must_be(nonvar,Evidence),
 1611  must_be(nonneg,S),
 1612  must_be(var,P),
 1613  M:asserta(('$goal'(1):-Goal,!),Ref1),
 1614  M:asserta('$goal'(0),Ref0),
 1615  mc_particle_sample_arg(M:'$goal'(A),M:Evidence,S,A,ValList),
 1616  foldl(agg_val,ValList,0,Sum),
 1617  foldl(value_cont_single,ValList,0,SumTrue),
 1618  P is SumTrue/Sum,
 1619  erase(Ref1),
 1620  erase(Ref0).
 mc_particle_sample_arg(:Query:conjunction_of_literals, +Evidence:list, +Samples:int, ?Arg:term, -Values:list) is det
The predicate samples Query a number of Samples times given that Evidence is true. It performs particle filtering with likelihood weighting: each sample is weighted by the likelihood of an element of the Evidence list and constitutes a particle. After weighting, particles are resampled and the next element of Evidence is considered. Arg should be a variable in Query. Evidence is a list of goals. Query can be either a single goal or a list of goals. When Query is a single goal, the predicate returns in Values a list of couples V-W where V is a value of Arg for which Query succeeds in a particle in the last set of particles and W is the weight of the particle. For each element of Evidence, the particles are obtained by sampling Query in each current particle and weighting the particle by the likelihood of the evidence element. When Query is a list of goals, Arg is a list of variables, one for each query of Query and Arg and Query must have the same length of Evidence. Values is then list of the same length of Evidence and each of its elements is a list of couples V-W where V is a value of the corresponding element of Arg for which the corresponding element of Query succeeds in a particle and W is the weight of the particle. For each element of Evidence, the particles are obtained by sampling the corresponding element of Query in each current particle and weighting the particle by the likelihood of the evidence element. /
 1650mc_particle_sample_arg(M:Goal,M:Evidence,S,Arg,ValList):-
 1651  must_be(nonvar,Goal),
 1652  must_be(nonvar,Evidence),
 1653  must_be(nonneg,S),
 1654  must_be(var,Arg),
 1655  must_be(var,ValList),
 1656  ValList=[V0|ValList0],
 1657  test_prism(M),
 1658  Goal=[G1|GR],!,
 1659  Evidence=[Ev1|EvR],
 1660  Arg=[A1|AR],
 1661  particle_sample_first_gl(0,S,M:G1,M:Ev1,A1,V0),
 1662  particle_sample_arg_gl(M:GR,M:EvR,AR,1,S,ValList0),
 1663  retractall(M:mem(_,_,_,_)),
 1664  retractall(M:mem(_,_,_,_,_)),
 1665  retractall(M:value_particle(_,_,_)),
 1666  retractall(M:weight_particle(_,_,_)).
 1667
 1668mc_particle_sample_arg(M:Goal,M:Evidence,S,Arg,ValList):-
 1669  Evidence=[Ev1|EvR],
 1670  particle_sample_first(0,S,M:Goal,M:Ev1,Arg),
 1671  particle_sample_arg(M:EvR,M:Goal,1,S,ValList0),
 1672  foldl(agg_val,ValList0,0,Sum),
 1673  Norm is S/Sum,
 1674  retractall(M:mem(_,_,_,_)),
 1675  retractall(M:mem(_,_,_,_,_)),
 1676  retractall(M:value_particle(_,_,_)),
 1677  retractall(M:weight_particle(_,_,_)),
 1678  maplist(norm(Norm),ValList0,ValList).
 mc_particle_expectation(:Query:conjunction_of_literals, :Evidence:list, +N:int, ?Arg:var, -Exp:float) is det
The predicate computes the expected value of Arg in Query given Evidence by particle filtering. It uses N particle and sums up the weighted value of Arg for each particle. The overall sum is divided by the sum of weights to give Exp. Arg should be a variable in Query. /
 1689mc_particle_expectation(M:Goal,M:Evidence,S,Arg,E):-
 1690  must_be(nonvar,Goal),
 1691  must_be(nonvar,Evidence),
 1692  must_be(nonneg,S),
 1693  must_be(var,Arg),
 1694  must_be(var,E),
 1695  mc_particle_sample_arg(M:Goal,M:Evidence,S,Arg,ValList),
 1696  average(ValList,E).
 1697
 1698particle_sample_arg_gl(M:[],M:[],[],I,_S,[]):- !,
 1699  retractall(M:mem(I,_,_,_,_)).
 1700
 1701particle_sample_arg_gl(M:[HG|TG],M:[HE|TE],[HA|TA],I,S,[HV|TV]):-
 1702  I1 is I+1,
 1703  resample_gl(M,I,I1,S),
 1704  retractall(M:value_particle(I,_,_)),
 1705  retractall(M:weight_particle(I,_,_)),
 1706  particle_sample_gl(0,S,M:HG,M:HE,HA,I1,HV),
 1707  particle_sample_arg_gl(M:TG,M:TE,TA,I1,S,TV).
 1708
 1709resample_gl(M,I,I1,S):-
 1710  get_values(M,I,V0),
 1711  foldl(agg_val,V0,0,Sum),
 1712  Norm is 1.0/Sum,
 1713  maplist(norm(Norm),V0,V1),
 1714  numlist(1,S,L),
 1715  maplist(weight_to_prob,L,V1,V2),
 1716  W is 1.0/S,
 1717  take_samples_gl(M,0,S,I,I1,W,V2),
 1718  retractall(M:mem(I,_,_,_,_)).
 1719
 1720weight_to_prob(I,_V-W,I:W).
 1721
 1722take_samples_gl(_M,S,S,_I,_I1,_W,_V):-!.
 1723
 1724take_samples_gl(M,S0,S,I,I1,W,V):-
 1725  S1 is S0+1,
 1726  discrete(V,_M,SInd),
 1727  restore_samples(M,I,SInd),
 1728  save_samples_tab(M,I1,S1),
 1729  take_samples_gl(M,S1,S,I,I1,W,V).
 1730
 1731particle_sample_gl(K,K,M:_G,_Ev,_A,I,L):-!,
 1732  get_values(M,I,L0),
 1733  foldl(agg_val,L0,0,Sum),
 1734  Norm is K/Sum,
 1735  maplist(norm(Norm),L0,L).
 1736
 1737
 1738particle_sample_gl(K0,S,M:Goal,M:Evidence,Arg,I,L):-
 1739  K1 is K0+1,
 1740  restore_samples(M,I,K1),
 1741  copy_term((Goal,Arg),(Goal1,Arg1)),
 1742  lw_sample_cycle(M:Goal1),
 1743  copy_term(Evidence,Ev1),
 1744  lw_sample_weight_cycle(M:Ev1,W),
 1745  save_samples_tab(M,I,K1),
 1746  assert(M:weight_particle(I,K1,W)),
 1747  assert(M:value_particle(I,K1,Arg1)),
 1748  particle_sample_gl(K1,S,M:Goal,M:Evidence,Arg,I,L).
 1749
 1750particle_sample_first_gl(K,K,M:_Goals,_Ev,_Arg,L):-!,
 1751  get_values(M,1,L0),
 1752  foldl(agg_val,L0,0,Sum),
 1753  Norm is K/Sum,
 1754  maplist(norm(Norm),L0,L).
 1755
 1756
 1757particle_sample_first_gl(K0,S,M:Goal, M:Evidence, Arg,V):-
 1758  K1 is K0+1,
 1759  copy_term((Goal,Arg),(Goal1,Arg1)),
 1760  lw_sample_cycle(M:Goal1),
 1761  copy_term(Evidence,Ev1),
 1762  lw_sample_weight_cycle(M:Ev1,W),
 1763  save_samples_tab(M,1,K1),
 1764  assert(M:weight_particle(1,K1,W)),
 1765  assert(M:value_particle(1,K1,Arg1)),
 1766  particle_sample_first_gl(K1,S,M:Goal,M:Evidence,Arg,V).
 1767
 1768
 1769particle_sample_arg(M:[],_Goal,I,_S,L):-!,
 1770  get_values(M,I,L).
 1771
 1772particle_sample_arg(M:[HE|TE],M:Goal,I,S,L):-
 1773  I1 is I+1,
 1774  resample(M,I,I1,S),
 1775  retractall(M:value_particle(I,_,_)),
 1776  retractall(M:weight_particle(I,_,_)),
 1777  particle_sample(0,S, M:HE, I1),
 1778  retractall(M:mem(I,_,_,_,_)),
 1779  particle_sample_arg(M:TE,M:Goal,I1,S,L).
 1780
 1781resample(M,I,I1,S):-
 1782  get_values(M,I,V0),
 1783  foldl(agg_val,V0,0,Sum),
 1784  Norm is 1.0/Sum,
 1785  maplist(norm(Norm),V0,V1),
 1786  numlist(1,S,L),
 1787  maplist(weight_to_prob,L,V1,V2),
 1788  W is 1.0/S,
 1789  take_samples(M,0,S,I,I1,W,V2).
 1790
 1791
 1792take_samples(_M,S,S,_I,_I1,_W,_V):-!.
 1793
 1794take_samples(M,S0,S,I,I1,W,V):-
 1795  S1 is S0+1,
 1796  discrete(V,_M,SInd),
 1797  restore_samples(M,I,SInd),
 1798  save_samples_tab(M,I1,S1),
 1799  M:value_particle(I,SInd,Arg),!,
 1800  assert(M:value_particle(I1,S1,Arg)),
 1801  take_samples(M,S1,S,I,I1,W,V).
 1802
 1803
 1804particle_sample(K,K,_Ev,_I):-!.
 1805
 1806particle_sample(K0,S,M:Evidence,I):-
 1807  K1 is K0+1,
 1808  restore_samples(M,I,K1),
 1809  copy_term(Evidence,Ev1),
 1810  lw_sample_weight_cycle(M:Ev1,W),
 1811  save_samples_tab(M,I,K1),
 1812  assert(M:weight_particle(I,K1,W)),
 1813  particle_sample(K1,S,M:Evidence,I).
 1814
 1815particle_sample_first(K,K,_Goals,_Ev,_Arg):-!.
 1816
 1817particle_sample_first(K0,S,M:Goal, M:Evidence, Arg):-
 1818  K1 is K0+1,
 1819  copy_term((Goal,Arg),(Goal1,Arg1)),
 1820  lw_sample_cycle(M:Goal1),
 1821  copy_term(Evidence,Ev1),
 1822  lw_sample_weight_cycle(M:Ev1,W),
 1823  save_samples_tab(M,1,K1),
 1824  assert(M:weight_particle(1,K1,W)),
 1825  assert(M:value_particle(1,K1,Arg1)),
 1826  particle_sample_first(K1,S,M:Goal,M:Evidence,Arg).
 1827
 1828get_values(M,I,V):-
 1829  findall(A-W,(M:value_particle(I,S,A),M:weight_particle(I,S,W)),V).
 mc_lw_sample(:Query:conjunction_of_literals, :Evidence:conjunction_of_literals, +Samples:int, -Prob:float) is det
The predicate samples Query a number of Samples times given that Evidence is true. The predicate returns in Prob the probability that the query is true. It performs likelihood weighting: each sample is weighted by the likelihood of evidence in the sample. /
 1840mc_lw_sample(M:Goal,M:Evidence0,S,P):-
 1841  must_be(nonvar,Goal),
 1842  must_be(nonvar,Evidence0),
 1843  must_be(nonneg,S),
 1844  must_be(var,P),
 1845  test_prism(M),
 1846  deal_with_ev(Evidence0,M,Evidence,UpdatedClausesRefs,ClausesToReAdd),
 1847  erase_samples(M),
 1848  lw_sample_bool(S,M:Goal,M:Evidence,ValList),
 1849  foldl(agg_val,ValList,0,Sum),
 1850  foldl(value_cont_single,ValList,0,SumTrue),
 1851  P is SumTrue/Sum,
 1852  maplist(erase,UpdatedClausesRefs),
 1853  maplist(M:assertz,ClausesToReAdd).
 1854
 1855
 1856value_cont_single(H-W,S,S+H*W).
 mc_lw_sample_arg(:Query:atom, :Evidence:atom, +Samples:int, ?Arg:var, -Values:list) is det
The predicate samples Query a number of Samples times given that Evidence is true. Arg should be a variable in Query. The predicate returns in Values a list of couples V-W where V is a value of Arg for which Query succeeds in a world sampled at random and W is the weight of the sample. It performs likelihood weighting: each sample is weighted by the likelihood of evidence in the sample. /
 1871mc_lw_sample_arg(M:Goal,M:Evidence0,S,Arg,ValList):-
 1872  must_be(nonvar,Goal),
 1873  must_be(nonvar,Evidence0),
 1874  must_be(nonneg,S),
 1875  must_be(var,Arg),
 1876  must_be(var,ValList),
 1877  test_prism(M),
 1878  deal_with_ev(Evidence0,M,Evidence,UpdatedClausesRefs,ClausesToReAdd),
 1879  lw_sample_arg(S,M:Goal,M:Evidence,Arg,ValList0),
 1880  foldl(agg_val,ValList0,0,Sum),
 1881  Norm is S/Sum,
 1882  maplist(norm(Norm),ValList0,ValList),
 1883  maplist(erase,UpdatedClausesRefs),
 1884  maplist(M:assertz,ClausesToReAdd).
 mc_lw_sample_arg_log(:Query:atom, :Evidence:atom, +Samples:int, ?Arg:var, -Values:list) is det
The predicate samples Query a number of Samples times given that Evidence is true. Arg should be a variable in Query. The predicate returns in Values a list of couples V-W where V is a value of Arg for which Query succeeds in a world sampled at random and W is the natural logarithm of the weight of \ the sample. It performs likelihood weighting: each sample is weighted by the likelihood of evidence in the sample. It differs from mc_lw_sample_arg/5 because the natural logarithm of the weight is returned, useful when the evidence is very unlikely. /
 1901mc_lw_sample_arg_log(M:Goal,M:Evidence0,S,Arg,ValList):-
 1902  must_be(nonvar,Goal),
 1903  must_be(nonvar,Evidence0),
 1904  must_be(nonneg,S),
 1905  must_be(var,Arg),
 1906  must_be(var,ValList),
 1907  test_prism(M),
 1908  deal_with_ev(Evidence0,M,Evidence,UpdatedClausesRefs,ClausesToReAdd),
 1909  lw_sample_arg_log(S,M:Goal,M:Evidence,Arg,ValList),
 1910  maplist(erase,UpdatedClausesRefs),
 1911  maplist(M:assertz,ClausesToReAdd).
 mc_lw_expectation(:Query:atom, :Evidence:atom, +N:int, ?Arg:var, -Exp:float) is det
The predicate computes the expected value of Arg in Query given Evidence by likelihood weighting. It takes N samples of Query and sums up the weighted value of Arg for each sample. The overall sum is divided by the sum of weights to give Exp. Arg should be a variable in Query. /
 1922mc_lw_expectation(M:Goal,M:Evidence,S,Arg,E):-
 1923  must_be(nonvar,Goal),
 1924  must_be(nonvar,Evidence),
 1925  must_be(nonneg,S),
 1926  must_be(var,Arg),
 1927  must_be(var,E),  
 1928  mc_lw_sample_arg(M:Goal,M:Evidence,S,Arg,ValList),
 1929  average(ValList,E).
 1930
 1931
 1932
 1933norm(NF,V-W,V-W1):-
 1934  W1 is W*NF.
 1935
 1936lw_sample_bool(0,_Goals,_Ev,[]):-!.
 1937
 1938lw_sample_bool(K0,M:Goal, M:Evidence, [S-W|V]):-
 1939  copy_term(Goal,Goal1),
 1940  (lw_sample(M:Goal1)->
 1941    S=1
 1942  ;
 1943    S=0
 1944  ),
 1945  copy_term(Evidence,Ev1),
 1946  (lw_sample_weight(M:Ev1,1,W0)->
 1947    W=W0
 1948  ;
 1949    W=0
 1950  ),
 1951  K1 is K0-1,
 1952  erase_samples(M),
 1953  lw_sample_bool(K1,M:Goal,M:Evidence,V).
 1954
 1955lw_sample_arg(0,_Goals,_Ev,_Arg,[]):-!.
 1956
 1957lw_sample_arg(K0,M:Goal, M:Evidence, Arg,[Arg1-W|V]):-
 1958  copy_term((Goal,Arg),(Goal1,Arg1)),
 1959  lw_sample_cycle(M:Goal1),
 1960  copy_term(Evidence,Ev1),
 1961  lw_sample_weight_cycle(M:Ev1,W),
 1962  K1 is K0-1,
 1963  erase_samples(M),
 1964  lw_sample_arg(K1,M:Goal,M:Evidence,Arg,V).
 1965
 1966lw_sample_cycle(M:G):-
 1967  (lw_sample(M:G)->
 1968    true
 1969  ;
 1970    erase_samples(M),
 1971    lw_sample_cycle(M:G)
 1972  ).
 1973
 1974lw_sample_arg_log(0,_Goals,_Ev,_Arg,[]):-!.
 1975
 1976lw_sample_arg_log(K0,M:Goal, M:Evidence, Arg,[Arg1-W|V]):-
 1977  copy_term((Goal,Arg),(Goal1,Arg1)),
 1978  lw_sample_cycle(M:Goal1),
 1979  copy_term(Evidence,Ev1),
 1980  lw_sample_logweight_cycle(M:Ev1,W),
 1981  K1 is K0-1,
 1982  erase_samples(M),
 1983  lw_sample_arg_log(K1,M:Goal,M:Evidence,Arg,V).
 1984
 1985
 1986lw_sample(_M:true):-!.
 1987
 1988lw_sample(M:A~=B):-!,
 1989  add_arg(A,B,A1),
 1990  lw_sample(M:A1).
 1991
 1992lw_sample(M:msw(A,B)):-!,
 1993  msw(M:A,B).
 1994
 1995lw_sample(M:G):-
 1996  G=..[call,P|A],!,
 1997  G1=..[P|A],
 1998  lw_sample(M:G1).
 1999
 2000lw_sample(_M:(sample_head(R,VC,M,HL,N))):-!,
 2001  sample_head(R,VC,M,HL,N).
 2002
 2003lw_sample(_M:take_a_sample(R,VC,M,Distr,S)):-!,
 2004  take_a_sample(R,VC,M,Distr,S).
 2005
 2006lw_sample(M:(G1,G2)):-!,
 2007  lw_sample(M:G1),
 2008  lw_sample(M:G2).
 2009
 2010lw_sample(M:(G1;G2)):-!,
 2011  lw_sample(M:G1);
 2012  lw_sample(M:G2).
 2013
 2014lw_sample(M:(\+ G)):-!,
 2015  \+ lw_sample(M:G).
 2016
 2017lw_sample(M:G):-
 2018  builtin(G),!,
 2019  M:call(G).
 2020
 2021lw_sample(M:G):-
 2022  \+ \+ M:sampled_g(G),!,
 2023  M:sampled_g(G).
 2024
 2025lw_sample(M:G):-
 2026  findall((G,B),M:clause(G,B),L),
 2027  sample_one_back(L,(G,B)),
 2028  lw_sample(M:B),
 2029  assert(M:sampled(r,G,1)).
 2030
 2031
 2032lw_sample_weight_cycle(M:G,W):-
 2033  (lw_sample_weight(M:G,1,W)->
 2034    true
 2035  ;
 2036    erase_samples(M),
 2037    lw_sample_weight_cycle(M:G,W)
 2038  ).
 2039
 2040lw_sample_weight(_M:true,W,W):-!.
 2041
 2042lw_sample_weight(M:A~= B,W0,W):-!,
 2043  add_arg(A,B,A1),
 2044  lw_sample_weight(M:A1,W0,W).
 2045
 2046lw_sample_weight(M:G,W0,W):-
 2047  G=..[call,P|A],!,
 2048  G1=..[P|A],
 2049  lw_sample_weight(M:G1,W0,W).
 2050
 2051lw_sample_weight(M:msw(A,B),W0,W):-!,
 2052  (var(B)->
 2053    msw(M:A,B),
 2054    W=W0
 2055  ;
 2056    msw_weight(M:A,B,W1),
 2057    W is W0*W1
 2058  ).
 2059
 2060lw_sample_weight(_M:(sample_head(R,VC,M,HL,N)),W0,W):-!,
 2061  % sample_head(R,VC,M,HL,N0),
 2062  % N=N0,
 2063  check(R,VC,M,N),
 2064%  W=W0.
 2065  nth0(N,HL,_:P),
 2066  W is W0*P.
 2067
 2068
 2069lw_sample_weight(_M:take_a_sample(R,VC,M,Distr,S),W0,W):-!,
 2070  (var(S)->
 2071    take_a_sample(R,VC,M,Distr,S),
 2072    W=W0
 2073  ;
 2074    (is_discrete(M,Distr)->
 2075      check(R,VC,M,S)
 2076    ;
 2077      true
 2078    ),
 2079    call(Distr,M,S,D),
 2080    W is W0*D
 2081   ).
 2082
 2083lw_sample_weight(M:(G1,G2),W0,W):-!,
 2084  lw_sample_weight(M:G1,W0,W1),
 2085  lw_sample_weight(M:G2,W1,W).
 2086
 2087lw_sample_weight(M:(G1;G2),W0,W):-!,
 2088  lw_sample_weight(M:G1,W0,W);
 2089  lw_sample_weight(M:G2,W0,W).
 2090
 2091lw_sample_weight(M:(\+ G),W0,W):-!,
 2092  lw_sample_weight(M:G,1,W1),
 2093  W is W0*(1-W1).
 2094
 2095lw_sample_weight(M:G,W,W):-
 2096  builtin(G),!,
 2097  M:call(G).
 2098
 2099lw_sample_weight(M:G,W0,W):-
 2100  \+ \+ M:sampled_g(G,_W1),!,
 2101  M:sampled_g(G,W1),
 2102  W is W0*W1.
 2103
 2104lw_sample_weight(M:G,W0,W):-
 2105  findall((G,B),M:clause(G,B),L),
 2106  sample_one_back(L,(G,B)),
 2107  lw_sample_weight(M:B,1,W1),
 2108  assert(M:sampled(rw,G,W1)),
 2109  W is W0*W1.
 2110
 2111
 2112lw_sample_logweight_cycle(M:G,W):-
 2113  (lw_sample_logweight(M:G,0,W)->
 2114    true
 2115  ;
 2116    erase_samples(M),
 2117    lw_sample_logweight_cycle(M:G,W)
 2118  ).
 2119
 2120
 2121lw_sample_logweight(_M:true,W,W):-!.
 2122
 2123lw_sample_logweight(M:A~= B,W0,W):-!,
 2124  add_arg(A,B,A1),
 2125  lw_sample_logweight(M:A1,W0,W).
 2126
 2127lw_sample_logweight(M:msw(A,B),W0,W):-!,
 2128  (var(B)->
 2129    msw(M:A,B),
 2130    W=W0
 2131  ;
 2132    msw_weight(M:A,B,W1),
 2133    W is W0+log(W1)
 2134  ).
 2135
 2136lw_sample_logweight(_M:(sample_head(R,VC,M,HL,N)),W0,W):-!,
 2137  sample_head(R,VC,M,HL,N0),
 2138  N=N0,
 2139  check(R,VC,M,N),
 2140 % W=W0.
 2141  nth0(N,HL,_:P),
 2142  W is W0+log(P).
 2143
 2144lw_sample_logweight(_M:take_a_sample(R,VC,M,Distr,S),W0,W):-!,
 2145  (var(S)->
 2146    take_a_sample(R,VC,M,Distr,S),
 2147    W=W0
 2148  ;
 2149    (is_discrete(M,Distr)->
 2150      check(R,VC,M,S)
 2151    ;
 2152      true
 2153    ),
 2154    call(Distr,M,S,D),
 2155    W is W0+log(D)
 2156   ).
 2157
 2158lw_sample_logweight(M:(G1,G2),W0,W):-!,
 2159  lw_sample_logweight(M:G1,W0,W1),
 2160  lw_sample_logweight(M:G2,W1,W).
 2161
 2162lw_sample_logweight(M:(G1;G2),W0,W):-!,
 2163  lw_sample_logweight(M:G1,W0,W);
 2164  lw_sample_logweight(M:G2,W0,W).
 2165
 2166lw_sample_logweight(M:(\+ G),W0,W):-!,
 2167  lw_sample_logweight(M:G,0,W1),
 2168  W is W0-W1.
 2169
 2170
 2171lw_sample_logweight(M:G,W,W):-
 2172  builtin(G),!,
 2173  M:call(G).
 2174
 2175lw_sample_logweight(M:G,W0,W):-
 2176  findall((G,B),M:clause(G,B),L),
 2177  sample_one_back(L,(G,B)),
 2178  lw_sample_logweight(M:B,W0,W).
 2179
 2180is_var(S):-
 2181  var(S),!.
 2182
 2183is_var(S):-
 2184  maplist(var,S).
 mc_sample_arg_first(:Query:atom, +Samples:int, ?Arg:var, -Values:list, +Options:list) is det
The predicate samples Query a number of Samples times. Arg should be a variable in Query. The predicate returns in Values a list of couples V-N where V is the value of Arg returned as the first answer by Query in a world sampled at random and N is the number of samples returning that value. V is failure if the query fails.

Options is a list of options, the following are recognised by mc_sample_arg_first/5:

bar(-BarChar:dict)
BarChart is a dict for rendering with c3 as a bar chart with with a bar for each value of Arg returned as a first answer by Query in a world sampled at random. The size of the bar is the number of samples that returned that value. /
 2204mc_sample_arg_first(M:Goal,S,Arg,ValList,Options):-
 2205  must_be(nonvar,Goal),
 2206  must_be(nonneg,S),
 2207  must_be(var,Arg),
 2208  must_be(var,ValList),
 2209  must_be(list,Options),
 2210  empty_assoc(Values0),
 2211  sample_arg_first(S,M:Goal,Arg, Values0,Values),
 2212  erase_samples(M),
 2213  assoc_to_list(Values,ValList0),
 2214  sort(2, @>=,ValList0,ValList),
 2215  option(bar(Chart),Options,no),
 2216  (nonvar(Chart)->
 2217    true
 2218  ;
 2219    argbar(ValList,Chart)
 2220  ).
 mc_sample_arg_first(:Query:atom, +Samples:int, ?Arg:var, -Values:list) is det
Equivalent to mc_sample_arg_first/5 with an empty option list. /
 2227mc_sample_arg_first(M:Goal,S,Arg,ValList):-
 2228  must_be(nonvar,Goal),
 2229  must_be(nonneg,S),
 2230  must_be(var,Arg),
 2231  must_be(var,ValList),
 2232  mc_sample_arg_first(M:Goal,S,Arg,ValList,[]).
 2233
 2234sample_arg_first(0,_Goals,_Arg,V,V):-!.
 2235
 2236sample_arg_first(K1, M:Goals,Arg,V0,V):-
 2237  erase_samples(M),
 2238  copy_term((Goals,Arg),(Goals1,Arg1)),
 2239  (M:Goals1->
 2240    numbervars(Arg1),
 2241    Val=Arg1
 2242  ;
 2243    Val=failure
 2244  ),
 2245  (get_assoc(Val, V0, N)->
 2246    N1 is N+1,
 2247    put_assoc(Val,V0,N1,V1)
 2248  ;
 2249    put_assoc(Val,V0,1,V1)
 2250  ),
 2251  K2 is K1-1,
 2252  sample_arg_first(K2,M:Goals,Arg,V1,V).
 mc_sample_arg_one(:Query:atom, +Samples:int, ?Arg:var, -Values:list, +Options:list) is det
The predicate samples Query a number of Samples times. Arg should be a variable in Query. The predicate returns in Values a list of couples V-N where V is a value of Arg sampled with uniform probability from those returned by Query in a world sampled at random and N is the number of samples returning that value. V is failure if the query fails.

Options is a list of options, the following are recognised by mc_sample_arg_one/5:

bar(-BarChar:dict)
BarChart is a dict for rendering with c3 as a bar chart with a bar for each value of Arg returned by sampling with uniform probability one answer from those returned by Query in a world sampled at random. The size of the bar is the number of samples. /
 2273mc_sample_arg_one(M:Goal,S,Arg,ValList,Options):-
 2274  must_be(nonvar,Goal),
 2275  must_be(nonneg,S),
 2276  must_be(var,Arg),
 2277  must_be(var,ValList),
 2278  must_be(list,Options),
 2279  empty_assoc(Values0),
 2280  sample_arg_one(S,M:Goal,Arg, Values0,Values),
 2281  erase_samples(M),
 2282  assoc_to_list(Values,ValList0),
 2283  sort(2, @>=,ValList0,ValList),
 2284  option(bar(Chart),Options,no),
 2285  (nonvar(Chart)->
 2286    true
 2287  ;
 2288    argbar(ValList,Chart)
 2289  ).
 mc_sample_arg_one(:Query:atom, +Samples:int, ?Arg:var, -Values:list) is det
Equivalent to mc_sample_arg_one/5 with an empty option list. /
 2296mc_sample_arg_one(M:Goal,S,Arg,ValList):-
 2297  must_be(nonvar,Goal),
 2298  must_be(nonneg,S),
 2299  must_be(var,Arg),
 2300  must_be(var,ValList),
 2301  mc_sample_arg_one(M:Goal,S,Arg,ValList,[]).
 2302
 2303sample_arg_one(0,_Goals,_Arg,V,V):-!.
 2304
 2305sample_arg_one(K1, M:Goals,Arg,V0,V):-
 2306  erase_samples(M),
 2307  copy_term((Goals,Arg),(Goals1,Arg1)),
 2308  findall(Arg1,M:Goals1,L),
 2309  numbervars(L),
 2310  sample_one(L,Val),
 2311  (get_assoc(Val, V0, N)->
 2312    N1 is N+1,
 2313    put_assoc(Val,V0,N1,V1)
 2314  ;
 2315    put_assoc(Val,V0,1,V1)
 2316  ),
 2317  K2 is K1-1,
 2318  sample_arg_one(K2,M:Goals,Arg,V1,V).
 2319
 2320sample_one([],failure):-!.
 2321
 2322sample_one(List,El):-
 2323  length(List,L),
 2324  random(0,L,Pos),
 2325  nth0(Pos,List,El).
 2326
 2327sample_one_back([],_):-!,
 2328  fail.
 2329
 2330sample_one_back(List,El):-
 2331  length(List,L),
 2332  random(0,L,Pos),
 2333  nth0(Pos,List,El0,Rest),
 2334  sample_backtracking(Rest,El0,El).
 2335
 2336sample_backtracking([],El,El):-!.
 2337
 2338sample_backtracking(_,El,El).
 2339
 2340sample_backtracking(L,_El,El):-
 2341  sample_one_back(L,El).
 mc_sample_arg_raw(:Query:atom, +Samples:int, ?Arg:var, -Values:list) is det
The predicate samples Query a number of Samples times. Arg should be a variable in Query. The predicate returns in Values a list of values of Arg returned as the first answer by Query in a world sampled at random. The value is failure if the query fails. /
 2352mc_sample_arg_raw(M:Goal,S,Arg,Values):-
 2353  must_be(nonvar,Goal),
 2354  must_be(nonneg,S),
 2355  must_be(var,Arg),
 2356  must_be(var,Values),
 2357  sample_arg_raw(S,M:Goal,Arg,Values),
 2358  erase_samples(M).
 2359
 2360sample_arg_raw(0,_Goals,_Arg,[]):-!.
 2361
 2362sample_arg_raw(K1, M:Goals,Arg,[Val|V]):-
 2363  erase_samples(M),
 2364  copy_term((Goals,Arg),(Goals1,Arg1)),
 2365  (M:Goals1->
 2366    numbervars(Arg1),
 2367    Val=Arg1
 2368  ;
 2369    Val=failure
 2370  ),
 2371  K2 is K1-1,
 2372  sample_arg_raw(K2,M:Goals,Arg,V).
 mc_expectation(:Query:atom, +N:int, ?Arg:var, -Exp:float) is det
The predicate computes the expected value of Arg in Query by sampling. It takes N samples of Query and sums up the value of Arg for each sample. The overall sum is divided by N to give Exp. Arg should be a variable in Query. /
 2384mc_expectation(M:Goal,S,Arg,E):-
 2385  must_be(nonvar,Goal),
 2386  must_be(nonneg,S),
 2387  must_be(var,Arg),
 2388  must_be(var,E),
 2389  sample_val(S,M:Goal,Arg, 0,Sum),
 2390  erase_samples(M),
 2391  E is Sum/S.
 mc_gibbs_expectation(:Query:atom, +N:int, ?Arg:var, -Exp:float, +Options:list) is det
The predicate computes the expected value of Arg in Query by sampling. It takes N samples of Query and sums up the value of Arg for each sample. The overall sum is divided by N to give Exp. Arg should be a variable in Query. Options is a list of options, the following are recognised by mc_mh_sample_arg/6:
block +Block:int
Perform blocked Gibbs: Block variables are sampled together, default value 1
mix(+Mix:int)
The first Mix samples are discarded (mixing time), default value 0 /
 2407mc_gibbs_expectation(M:Goal,S,Arg,E,Options):-
 2408  must_be(nonvar,Goal),
 2409  must_be(nonneg,S),
 2410  must_be(var,Arg),
 2411  must_be(var,E),
 2412  must_be(list,Options),
 2413  mc_gibbs_sample_arg(M:Goal,S,Arg,ValList,Options),
 2414  average(ValList,[E]),
 2415  erase_samples(M).
 mc_gibbs_expectation(:Query:atom, +N:int, ?Arg:var, -Exp:float) is det
Equivalent to mc_gibbs_expectation/5 with an empty option list. /
 2422mc_gibbs_expectation(M:Goal,S,Arg,E):-
 2423  must_be(nonvar,Goal),
 2424  must_be(nonneg,S),
 2425  must_be(var,Arg),
 2426  must_be(var,E),
 2427  mc_gibbs_expectation(M:Goal,S,Arg,E,[]).
 mc_rejection_expectation(:Query:atom, :Evidence:atom, +N:int, ?Arg:var, -Exp:float) is det
The predicate computes the expected value of Arg in Query by sampling. It takes N samples of Query and sums up the value of Arg for each sample. The overall sum is divided by N to give Exp. Arg should be a variable in Query. /
 2439mc_rejection_expectation(M:Goal,M:Evidence,S,Arg,E):-
 2440  must_be(nonvar,Goal),
 2441  must_be(nonvar,Evidence),
 2442  must_be(nonneg,S),
 2443  must_be(var,Arg),
 2444  must_be(var,E),
 2445  mc_rejection_sample_arg(M:Goal,M:Evidence,S,Arg,ValList,[]),
 2446  average(ValList,[E]),
 2447  erase_samples(M).
 mc_gibbs_expectation(:Query:atom, :Evidence:atom, +N:int, ?Arg:var, -Exp:float, +Options:list) is det
The predicate computes the expected value of Arg in Query by Gibbs sampling. It takes N samples of Query and sums up the value of Arg for each sample. The overall sum is divided by N to give Exp. Arg should be a variable in Query.

Options is a list of options, the following are recognised by mc_mh_expectation/6:

block +Block:int
Perform blocked Gibbs: Block variables are sampled together, default value 1
mix(+Mix:int)
The first Mix samples are discarded (mixing time), default value 0 /
 2464mc_gibbs_expectation(M:Goal,M:Evidence,S,Arg,E,Options):-
 2465  must_be(nonvar,Goal),
 2466  must_be(nonvar,Evidence),
 2467  must_be(nonneg,S),
 2468  must_be(var,Arg),
 2469  must_be(var,E),
 2470  must_be(list,Options),
 2471  mc_gibbs_sample_arg(M:Goal,M:Evidence,S,Arg,ValList,Options),
 2472  average(ValList,[E]),
 2473  erase_samples(M).
 mc_mh_expectation(:Query:atom, :Evidence:atom, +N:int, ?Arg:var, -Exp:float, +Options:list) is det
The predicate computes the expected value of Arg in Query by Metropolis Hastings sampling. It takes N samples of Query and sums up the value of Arg for each sample. The overall sum is divided by N to give Exp. Arg should be a variable in Query.

Options is a list of options, the following are recognised by mc_mh_expectation/6:

mix(+Mix:int)
The first Mix samples are discarded (mixing time), default value 0
lag(+Lag:int)
lag between each sample, Lag sampled choices are forgotten, default value 1 /
 2490mc_mh_expectation(M:Goal,M:Evidence,S,Arg,E,Options):-
 2491  must_be(nonvar,Goal),
 2492  must_be(nonvar,Evidence),
 2493  must_be(nonneg,S),
 2494  must_be(var,Arg),
 2495  must_be(var,E),
 2496  must_be(list,Options),
 2497  mc_mh_sample_arg(M:Goal,M:Evidence,S,Arg,ValList,Options),
 2498  average(ValList,[E]),
 2499  erase_samples(M).
 mc_mh_expectation(:Query:atom, :Evidence:atom, +N:int, ?Arg:var, -Exp:float) is det
Equivalent to mc_mh_expectation/6 with an empty option list. /
 2506mc_mh_expectation(M:Goal,M:Evidence,S,Arg,E):-
 2507  must_be(nonvar,Goal),
 2508  must_be(nonvar,Evidence),
 2509  must_be(nonneg,S),
 2510  must_be(var,Arg),
 2511  must_be(var,E),
 2512  mc_mh_expectation(M:Goal,M:Evidence,S,Arg,E,[]).
 2513
 2514value_cont([]-_,0):-!.
 2515
 2516value_cont([H|_T]-N,S,S+N*H).
 2517
 2518sample_val(0,_Goals,_Arg,Sum,Sum):-!.
 2519
 2520sample_val(K1, M:Goals,Arg,Sum0,Sum):-
 2521  erase_samples(M),
 2522  copy_term((Goals,Arg),(Goals1,Arg1)),
 2523  (M:Goals1->
 2524    Sum1 is Sum0+Arg1
 2525  ;
 2526    Sum1=Sum
 2527  ),
 2528  K2 is K1-1,
 2529  sample_val(K2,M:Goals,Arg,Sum1,Sum).
 2530
 2531
 2532get_next_rule_number(PName,R):-
 2533  retract(PName:rule_n(R)),
 2534  R1 is R+1,
 2535  assert(PName:rule_n(R1)).
 2536
 2537
 2538assert_all([],_M,[]).
 2539
 2540assert_all([H|T],M,[HRef|TRef]):-
 2541  assertz(M:H,HRef),
 2542  assert_all(T,M,TRef).
 2543
 2544
 2545retract_all([]):-!.
 2546
 2547retract_all([H|T]):-
 2548  erase(H),
 2549  retract_all(T).
 2550
 2551
 2552
 2553add_mod_arg(A,_Module,A1):-
 2554  A=..[P|Args],
 2555  A1=..[P|Args].
 sample_head(+R:int, +Variables:list, +M:module, +HeadList:list, -HeadNumber:int) is det
samples a head from rule R instantiated as indicated by Variables (list of constants, one per variable. HeadList contains the head as a list. HeadNumber is the number of the sample head. Internal predicates used by the transformed input program /
 2564sample_head(R,VC,M,_HeadList,N):-
 2565  M:sampled(R,VC,NH),!,
 2566  N=NH.
 2567
 2568sample_head(R,VC,M,HeadList,N):-
 2569  sample(HeadList,NH),
 2570  assertz(M:sampled(R,VC,NH)),
 2571  N=NH.
 2572
 2573sample(HeadList, HeadId) :-
 2574  random(Prob),
 2575  sample(HeadList, 0, 0, Prob, HeadId), !.
 2576
 2577sample([_HeadTerm:HeadProb|Tail], Index, Prev, Prob, HeadId) :-
 2578	Succ is Index + 1,
 2579	Next is Prev + HeadProb,
 2580	(Prob =< Next ->
 2581		HeadId = Index;
 2582		sample(Tail, Succ, Next, Prob, HeadId)).
 uniform_dens(+Lower:float, +Upper:float, +M:module, -S:float) is det
Returns in S a sample from a distribution uniform in (Lower,Upper) /
 2589uniform_dens(L,U,_M,S):-
 2590  number(L),
 2591  random(V),
 2592  S is L+V*(U-L).
 2593
 2594uniform_dens(L,U,_M,_S,D):-
 2595  D is 1/(U-L).
 uniform(+Values:list, +M:module, -S:float) is det
Returns in S a value from Values chosen uniformly /
 2602uniform(Values,M,S):-
 2603  length(Values,Len),Prob is 1.0/Len,
 2604  maplist(add_prob(Prob),Values,Dist),
 2605  discrete(Dist,M,S).
 2606
 2607
 2608uniform(Values,_M,_S,D):-
 2609  length(Values,L),
 2610  D is 1.0/L.
 take_a_sample(+R:int, +VC:list, +M:module, +Distr:term, -S:term) is det
Returns in S a sample for a random variable with distribution Distr associated to rule R with substitution VC. If the variable has already been sampled, it retrieves the sampled value, otherwise it takes a new sample and records it for rule R with substitution VC. /
 2620take_a_sample(R,VC,M,_Distr,S):-
 2621  M:sampled(R,VC,S0),!,
 2622  S=S0.
 2623
 2624take_a_sample(R,VC,M,Distr,S):-
 2625  call(Distr,M,S0),
 2626  assertz(M:sampled(R,VC,S0)),
 2627  S=S0.
 gaussian(+Mean:float, +Variance:float, +M:module, -S:float) is det
samples a value from a Gaussian with mean Mean and variance Variance and returns it in S /
 2635gaussian(Mean,Variance,_M,S):-
 2636  number(Mean),!,
 2637  random(U1),
 2638  random(U2),
 2639  R is sqrt(-2*log(U1)),
 2640  Theta is 2*pi*U2,
 2641  S0 is R*cos(Theta),
 2642  StdDev is sqrt(Variance),
 2643  S is StdDev*S0+Mean.
 2644
 2645gaussian([Mean1,Mean2],Covariance,_M,[X1,X2]):-!,
 2646  random(U1),
 2647  random(U2),
 2648  R is sqrt(-2*log(U1)),
 2649  Theta is 2*pi*U2,
 2650  S0 is R*cos(Theta),
 2651  S1 is R*sin(Theta),
 2652  cholesky_decomposition(Covariance,A),
 2653  matrix_multiply(A,[[S0],[S1]],Az),
 2654  matrix_sum([[Mean1],[Mean2]],Az,[[X1],[X2]]).
 2655
 2656gaussian(Mean,Covariance,_M,X):-
 2657  length(Mean,N),
 2658  n_gaussian_var(0,N,Z),
 2659  cholesky_decomposition(Covariance,A),
 2660  transpose([Z],ZT),
 2661  matrix_multiply(A,ZT,Az),
 2662  transpose([Mean],MT),
 2663  matrix_sum(MT,Az,XT),
 2664  transpose(XT,[X]).
 2665
 2666n_gaussian_var(N,N,[]):-!.
 2667
 2668n_gaussian_var(N1,N,[Z]):-
 2669  N1 is N-1,!,
 2670  random(U1),
 2671  random(U2),
 2672  R is sqrt(-2*log(U1)),
 2673  Theta is 2*pi*U2,
 2674  Z is R*cos(Theta).
 2675
 2676n_gaussian_var(N1,N,[Z1,Z2|T]):-
 2677  N2 is N1+2,
 2678  random(U1),
 2679  random(U2),
 2680  R is sqrt(-2*log(U1)),
 2681  Theta is 2*pi*U2,
 2682  Z1 is R*cos(Theta),
 2683  Z2 is R*sin(Theta),
 2684  n_gaussian_var(N2,N,T).
 gaussian(+Mean:float, +Variance:float, +M:module, +S:float, -Density:float) is det
Computes the probability density of value S according to a Gaussian with mean Mean and variance Variance and returns it in Density. /
 2693gaussian(Mean,Variance,_M,S,D):-
 2694  number(Mean),!,
 2695  StdDev is sqrt(Variance),
 2696  D is 1/(StdDev*sqrt(2*pi))*exp(-(S-Mean)*(S-Mean)/(2*Variance)).
 2697
 2698gaussian(Mean,Covariance,_M,S,D):-
 2699  determinant(Covariance,Det),
 2700  matrix_diff([Mean],[S],S_M),
 2701  matrix_inversion(Covariance,Cov_1),
 2702  transpose(S_M,S_MT),
 2703  matrix_multiply(S_M,Cov_1,Aux),
 2704  matrix_multiply(Aux,S_MT,[[V]]),
 2705  length(Mean,K),
 2706  D is 1/sqrt((2*pi)^K*Det)*exp(-V/2).
 gamma(+Shape:float, +Scale:float, +M:module, -S:float) is det
samples a value from a Gamma density function with shape Shape and scale Scale returns it in S /
 2715gamma(A,Scale,_M,S):-
 2716  (A>=1->
 2717    gamma_gt1(A,S1)
 2718  ;
 2719    random(U),
 2720    A1 is A+1,
 2721    gamma_gt1(A1,S0),
 2722    S1 is S0*U^(1/A)
 2723  ),
 2724  S is Scale*S1.
 2725
 2726gamma_gt1(A,S):-
 2727  D is A-1.0/3.0,
 2728  C is 1.0/sqrt(9.0*D),
 2729  cycle_gamma(D,C,S).
 2730
 2731cycle_gamma(D,C,S):-
 2732  getv(C,X,V),
 2733  random(U),
 2734  S0 is D*V,
 2735  (U<1-0.0331*X^4->
 2736    S=S0
 2737  ;
 2738    LogU is log(U),
 2739    LogV is log(V),
 2740    (LogU<0.5*X^2+D*(1-V+LogV)->
 2741      S=S0
 2742    ;
 2743      cycle_gamma(D,C,S)
 2744    )
 2745  ).
 2746
 2747getv(C,X,V):-
 2748  gaussian(0.0,1.0,_M,X0),
 2749  V0 is (1+C*X0)^3,
 2750  (V0=<0->
 2751    getv(C,X,V)
 2752  ;
 2753    V=V0,
 2754    X=X0
 2755  ).
 2756
 2757gamma(K,Scale,_M,S,D):-
 2758  D is exp(-lgamma(K))/(Scale^K)*S^(K-1)*exp(-S/Scale).
 beta(+Alpha:float, +Beta:float, +M:module, -S:float) is det
samples a value from a beta probability distribution with parameters Alpha and Beta and returns it in S. Uses the algorithm of van der Waerden, B. L., "Mathematical Statistics", Springer see also https://en.wikipedia.org/wiki/Beta_distribution#Generating_beta-distributed_random_variates /
 2770beta(Alpha,Beta,M,S):-
 2771  gamma(Alpha,1,M,X),
 2772  gamma(Beta,1,M,Y),
 2773  S is X/(X+Y).
 2774
 2775beta(Alpha,Beta,_M,X,D):-
 2776  B is exp(lgamma(Alpha)+lgamma(Beta)-lgamma(Alpha+Beta)),
 2777  D is X^(Alpha-1)*(1-X)^(Beta-1)/B.
 poisson(+Lambda:float, +M:module, -S:int) is det
samples a value from a Poisson probability distribution with parameter Lambda and returns it in S. Uses the inversion by sequential search Devroye, Luc (1986). "Discrete Univariate Distributions" Non-Uniform Random Variate Generation. New York: Springer-Verlag. p. 505. /
 2789poisson(Lambda,_M,X):-
 2790  P is exp(-Lambda),
 2791  random(U),
 2792  poisson_cycle(0,X,Lambda,P,P,U).
 2793
 2794poisson_cycle(X,X,_L,_P,S,U):-
 2795  U=<S,!.
 2796
 2797poisson_cycle(X0,X,L,P0,S0,U):-
 2798  X1 is X0+1,
 2799  P is P0*L/X1,
 2800  S is S0+P,
 2801  poisson_cycle(X1,X,L,P,S,U).
 2802
 2803poisson(Lambda,_M,X,P):-
 2804  fact(X,1,FX),
 2805  P is (Lambda^X)*exp(-Lambda)/FX.
 2806
 2807fact(N,F,F):- N =< 0, !.
 2808
 2809fact(N,F0,F):-
 2810  F1 is F0*N,
 2811  N1 is N-1,
 2812  fact(N1,F1,F).
 binomial(+N:int, +P:float, +M:module, -S:int) is det
samples a value from a binomial probability distribution with parameters N and P and returns it in S. /
 2820binomial(N,1.0,_M,N):-!.
 2821
 2822binomial(N,P,_M,X):-
 2823  Pr0 is (1-P)^N,
 2824  random(U),
 2825  binomial_cycle(0,X,N,P,Pr0,Pr0,U).
 2826
 2827binomial_cycle(X,X,_N,_P,_Pr,CPr,U):-
 2828  U=<CPr,!.
 2829
 2830binomial_cycle(X0,X,N,P,Pr0,CPr0,U):-
 2831  X1 is X0+1,
 2832  Pr is Pr0*P*(N-X0)/(X1*(1-P)),
 2833  CPr is CPr0+Pr,
 2834  binomial_cycle(X1,X,N,P,Pr,CPr,U).
 2835
 2836binomial(N,P,_M,X,Pr):-
 2837  fact(N,1,FN),
 2838  fact(X,1,FX),
 2839  N_X is N-X,
 2840  fact(N_X,1,FN_X),
 2841  Pr is P^X*(1-P)^N_X*FN/(FX*FN_X).
 multinomial(+N:int, +P:list, +M:module, -S:list) is det
samples a value from a multinomial probability distribution with parameters N and P and returns it in S. /
 2850multinomial(N,P,_M,S):-
 2851  length(P,K),
 2852  numlist(1,K,Outcomes),
 2853  length(Distribution,K),
 2854  maplist(pair,Outcomes,P,Distribution),
 2855  findall(0,between(1,K,_),SampleState),
 2856  Sample =..[sample|SampleState],
 2857  multinomial_cycle(0,N,Distribution,Sample),
 2858  Sample=..[_|S].
 2859
 2860pair(A,B,A:B).
 2861
 2862multinomial_cycle(N,N,_Distribution,_Sample):-!.
 2863
 2864multinomial_cycle(N0,N,Distribution,Sample):-
 2865  discrete(Distribution,_,Index),
 2866  arg(Index,Sample,Count),
 2867  Count1 is Count+1,
 2868  setarg(Index,Sample,Count1),
 2869  N1 is N0+1,
 2870  multinomial_cycle(N1,N,Distribution,Sample).
 2871
 2872multinomial(N,P,_M,S,Pr):-
 2873  fact(N,1,FN),
 2874  foldl(factor_multi,P,S,FN,Pr).
 2875
 2876factor_multi(Pi,Xi,Pr0,Pr):-
 2877  fact(Xi,1,FXi),
 2878  Pr is Pr0*Pi^Xi/FXi.
 dirichlet(+Par:list, +M:module, -S:float) is det
samples a value from a Dirichlet probability density with parameters Par and returns it in S /
 2886dirichlet(Par,_M,S):-
 2887  maplist(get_gamma,Par,Gammas),
 2888  sum_list(Gammas,Sum),
 2889  maplist(divide(Sum),Gammas,S).
 2890
 2891divide(S0,A,S):-
 2892  S is A/S0.
 2893
 2894get_gamma(A,G):-
 2895  gamma(A,1.0,_M,G).
 2896
 2897dirichlet(Par,_M,S,D):-
 2898  beta(Par,B),
 2899  foldl(prod,S,Par,1,D0),
 2900  D is D0*B.
 2901
 2902prod(X,A,P0,P0*X^(A-1)).
 geometric(+P:float, +M:module, -I:int) is det
samples a value from a geometric probability distribution with parameters P and returns it in I (I belongs to [1,infinity] /
 2911geometric(P,_M,I):-
 2912  geometric_val(1,P,I).
 2913
 2914geometric_val(N0,P,N):-
 2915  random(R),
 2916  (R=<P->
 2917    N=N0
 2918  ;
 2919    N1 is N0+1,
 2920    geometric_val(N1,P,N)
 2921  ).
 2922
 2923geometric(P,_M,I,D):-
 2924  D is (1-P)^(I-1)*P.
 finite(+Distribution:list, +M:module, -S:float) is det
samples a value from a discrete distribution Distribution (a list of couples Val:Prob) and returns it in S /
 2931finite(D,M,S):-
 2932  discrete(D,M,S).
 2933
 2934finite(D,M,S,Dens):-
 2935  discrete(D,M,S,Dens).
 discrete(+Distribution:list, +M:module, -S:float) is det
samples a value from a discrete distribution Distribution (a list of couples Val:Prob) and returns it in S /
 2943discrete(D,_M,S):-
 2944  random(U),
 2945  discrete_int(D,0,U,S).
 2946
 2947
 2948discrete_int([S:_],_,_,S):-!.
 2949
 2950discrete_int([S0:W|T],W0,U,S):-
 2951  W1 is W0+W,
 2952  (U=<W1->
 2953    S=S0
 2954  ;
 2955    discrete_int(T,W1,U,S)
 2956  ).
 2957
 2958discrete(D,_M,S,Dens):-
 2959  member(S:Dens,D).
 2960% discrete(_D,_S,1.0).
 exponential(+Lambda:float, +M:module, -V:int) is det
Samples a value from exponential distribution with parameter Lambda

**/

 2967exponential(Lambda,_M,V):-
 2968  random(U),
 2969  V is - log(U)/Lambda.
 2970
 2971exponential(Lambda,_M,X,V):-
 2972  V is Lambda*exp(-Lambda*X).
 negative_binomial(+N:int, +P:float, +M:module, -Value:int) is det
samples a value from a negative binomial probability distribution with parameters N and P and returns it in Value. N is the number of successes P is the success probability /
 2984% N number of successes
 2985% P probability of success
 2986negative_binomial(N,P,_M,Value):-
 2987  random(RandomVal),
 2988  negative_binomial_prob_(0,N,P,0,RandomVal,Value).
 2989
 2990negative_binomial_prob_(I,_,_,CurrentProb,RandomVal,I1):-
 2991  RandomVal =< CurrentProb, !,
 2992  I1 is I - 1.
 2993negative_binomial_prob_(I,R,P,CurrentProb,RandomVal,V):-
 2994  negative_binomial_int(I,R,P,V0),
 2995  I1 is I+1,
 2996  CurrentProb1 is V0 + CurrentProb,
 2997  negative_binomial_prob_(I1,R,P,CurrentProb1,RandomVal,V).
 2998
 2999/*
 3000* N number of successes
 3001* X number of failures
 3002* P probability of success
 3003*/
 3004negative_binomial_int(X,N,P,Value):-
 3005  XN1 is X+N-1,
 3006  binomial_coeff(XN1,X,Bin),
 3007  Value is Bin*(P**N)*(1-P)**X.
 3008
 3009binomial_coeff(N,K,0):- N < K, !.
 3010binomial_coeff(N,K,Val):-
 3011  fact(N,1,NF),
 3012  fact(K,1,KF),
 3013  NK is N-K,
 3014  fact(NK,1,NKF),
 3015  Val is NF/(KF*NKF).
 user(+Distr:atom, +M:module, -Value:int) is det
samples a value from a user defined distribution /
 3022user(Distr,M,S):-
 3023  call(M:Distr,S).
 3024
 3025user(Distr,M,S,D):-
 3026  call(M:Distr,S,D).
 3027
 3028generate_rules_fact([],HeadList,VC,M,R,_N,[Rule]):-
 3029  Rule=(samp(R,VC,N):-(sample_head(R,VC,M,HeadList,N))).
 3030
 3031generate_rules_fact([],_HL,_VC,_M,_R,_N,[]).
 3032
 3033generate_rules_fact([Head:_P1,'':_P2],HeadList,VC,M,R,N,[Clause]):-!,
 3034  Clause=(Head:-(sample_head(R,VC,M,HeadList,N))).
 3035
 3036generate_rules_fact([Head:_P|T],HeadList,VC,M,R,N,[Clause|Clauses]):-
 3037  Clause=(Head:-(sample_head(R,VC,M,HeadList,N))),
 3038  N1 is N+1,
 3039  generate_rules_fact(T,HeadList,VC,M,R,N1,Clauses).
 3040
 3041
 3042generate_clause_distr(Head,Body,VC,M,R,Var,Distr,Clause):-
 3043  Clause=[(Head:-(Body,take_a_sample(R,VC,M,Distr,Var))),
 3044  (samp(R,VC,Var):-take_a_sample(R,VC,M,Distr,Var))].
 3045
 3046
 3047generate_clause_samp(Head,Body,HeadList,VC,M,R,N,[Clause,Rule]):-
 3048  generate_clause(Head,Body,HeadList,VC,M,R,N,Clause),
 3049  Rule=(samp(R,VC,Val):-sample_head(R,VC,M,HeadList,Val)).
 3050
 3051generate_clause(Head,Body,HeadList,VC,M,R,N,Clause):-
 3052  Clause=[(Head:-(Body,sample_head(R,VC,M,HeadList,N)))].
 3053
 3054
 3055generate_rules([],_Body,HeadList,VC,M,R,_N,[Rule]):-
 3056  Rule=(samp(R,VC,N):-sample_head(R,VC,M,HeadList,N)).
 3057
 3058generate_rules([Head:_P1,'':_P2],Body,HeadList,VC,M,R,N,[Clause]):-!,
 3059  generate_clause(Head,Body,HeadList,VC,M,R,N,Clause).
 3060
 3061generate_rules([Head:_P|T],Body,HeadList,VC,M,R,N,[Clause|Clauses]):-
 3062  generate_clause(Head,Body,HeadList,VC,M,R,N,Clause),
 3063  N1 is N+1,
 3064  generate_rules(T,Body,HeadList,VC,M,R,N1,Clauses).
 3065
 3066
 3067
 3068
 3069
 3070
 3071process_head(HeadList, GroundHeadList) :-
 3072  ground_prob(HeadList), !,
 3073  process_head_ground(HeadList, 0, GroundHeadList).
 3074
 3075process_head(HeadList0, HeadList):-
 3076  get_probs(HeadList0,PL),
 3077  foldl(minus,PL,1,PNull),
 3078  append(HeadList0,['':PNull],HeadList).
 3079
 3080minus(A,B,B-A).
 3081
 3082prob_ann(_:P,P):-!.
 3083prob_ann(P::_,P).
 add_prob(?Prob:float, :Goal:atom, ?AnnGoal:atom) is det
From Prob and Goal builds the annotated atom AnnGoal=Goal:Prob. /
 3090add_prob(P,A,A:P).
 3091
 3092/* process_head_ground([Head:ProbHead], Prob, [Head:ProbHead|Null])
 3093 * ----------------------------------------------------------------
 3094 */
 3095process_head_ground([H], Prob, [Head:ProbHead1|Null]) :-
 3096  (H=Head:ProbHead;H=ProbHead::Head),!,
 3097  ProbHead1 is ProbHead,
 3098  ProbLast is 1 - Prob - ProbHead1,
 3099  prolog_load_context(module, M),mc_input_mod(M),
 3100  M:local_mc_setting(epsilon_parsing, Eps),
 3101  EpsNeg is - Eps,
 3102  ProbLast > EpsNeg,
 3103  (ProbLast > Eps ->
 3104    Null = ['':ProbLast]
 3105  ;
 3106    Null = []
 3107  ).
 3108
 3109process_head_ground([H|Tail], Prob, [Head:ProbHead1|Next]) :-
 3110  (H=Head:ProbHead;H=ProbHead::Head),
 3111  ProbHead1 is ProbHead,
 3112  ProbNext is Prob + ProbHead1,
 3113  process_head_ground(Tail, ProbNext, Next).
 3114
 3115
 3116ground_prob([]).
 3117
 3118ground_prob([_Head:ProbHead|Tail]) :-!,
 3119  ground(ProbHead), % Succeeds if there are no free variables in the term ProbHead.
 3120  ground_prob(Tail).
 3121
 3122ground_prob([ProbHead::_Head|Tail]) :-
 3123  ground(ProbHead), % Succeeds if there are no free variables in the term ProbHead.
 3124  ground_prob(Tail).
 3125
 3126get_probs(Head, PL):-
 3127  maplist(prob_ann,Head,PL).
 set_mc(:Parameter:atom, +Value:term) is det
The predicate sets the value of a parameter For a list of parameters see https://friguzzi.github.io/cplint/ /
 3138set_mc(M:Parameter,Value):-
 3139  must_be(atom,Parameter),
 3140  must_be(nonvar,Value),
 3141  retract(M:local_mc_setting(Parameter,_)),
 3142  assert(M:local_mc_setting(Parameter,Value)).
 setting_mc(:Parameter:atom, ?Value:term) is det
The predicate returns the value of a parameter For a list of parameters see https://friguzzi.github.io/cplint/ /
 3151setting_mc(M:P,V):-
 3152  must_be(atom,P),
 3153  M:local_mc_setting(P,V).
 3154
 3155
 3156delete_equal([],_,[]).
 3157
 3158delete_equal([H|T],E,T):-
 3159  H == E,!.
 3160
 3161delete_equal([H|T],E,[H|T1]):-
 3162  delete_equal(T,E,T1).
 3163
 3164add_arg(A,Arg,A1):-
 3165  A=..L,
 3166  append(L,[Arg],L1),
 3167  A1=..L1.
 set_sw(:Var:term, +List:lit) is det
Sets the domain of the random variable Var to List. This is a predicate for programs in the PRISM syntax /
 3175set_sw(M:A,B):-
 3176  M:local_mc_setting(prism_memoization,false),!,
 3177  assert(M:sw(A,B)).
 3178
 3179set_sw(M:A,B):-
 3180  M:local_mc_setting(prism_memoization,true),
 3181  get_next_rule_number(M,R),
 3182  assert(M:sw(R,A,B)).
 msw(:Var:term, ?Value:term) is det
Gets or tests the Value of the random variable Var. This is a predicate for programs in the PRISM syntax /
 3190msw(M:A,B):-
 3191  M:values(A,V),
 3192  M:sw(A,D),!,
 3193  sample_msw(V,D,B).
 3194
 3195msw(M:A,B):-
 3196  M:values(A,V),
 3197  M:sw(R,A,D),
 3198  sample_msw(V,M,R,A,D,B).
 3199
 3200sample_msw(real,norm(Mean,Variance),V):-!,
 3201  gaussian(Mean,Variance,_M,S),
 3202  S=V.
 3203
 3204sample_msw(Values,Dist,V):-
 3205  maplist(combine,Values,Dist,VD),
 3206  sample(VD,N),
 3207  nth0(N,Values,V).
 3208
 3209sample_msw(real,M,R,A,norm(Mean,Variance),V):-!,
 3210  take_a_sample(R,A,M,gaussian(Mean,Variance),V).
 3211
 3212sample_msw(Values,M,R,A,Dist,V):-
 3213  maplist(combine,Values,Dist,VD),
 3214  take_a_sample(R,A,M,discrete(VD),V).
 3215
 3216combine(V,P,V:P).
 3217
 3218msw_weight(M:A,B,W):-
 3219  M:values(A,V),
 3220  M:sw(A,D),!,
 3221  msw_weight(V,D,B,W).
 3222
 3223msw_weight(M:A,B,W):-
 3224  M:values(A,V),
 3225  M:sw(_R,A,D),
 3226  msw_weight(V,D,B,W).
 3227
 3228msw_weight(real,norm(Mean,Variance),V,W):-!,
 3229  gaussian(Mean,Variance,V,W).
 3230
 3231msw_weight(Values,Dist,V,W):-
 3232  maplist(combine,Values,Dist,VD),
 3233  member(V:W,VD).
 3234
 3235
 3236test_prism(M):-
 3237  (M:local_mc_setting(prism_memoization,false),M:values(_,_)->
 3238    throw(error("This predicate doesn't support PRISM programs without memoization."))
 3239  ;
 3240    true
 3241  ).
 3242
 3243act(M,A/B):-
 3244  M:(dynamic A/B).
 3245
 3246tab(A/B,A/B1):-
 3247  B1 is B + 2.
 3248
 3249
 3250mcintyre_expansion((:- mcaction Conj), []) :-!,
 3251  prolog_load_context(module, M),
 3252  mc_input_mod(M),!,
 3253  list2and(L,Conj),
 3254  maplist(act(M),L).
 3255
 3256mcintyre_expansion((:- table(Conj)), [:- table(Conj1)]) :-!,
 3257  prolog_load_context(module, M),
 3258  mc_input_mod(M),!,
 3259  list2and(L,Conj),
 3260  maplist(tab,L,L1),
 3261  list2and(L1,Conj1).
 3262
 3263mcintyre_expansion((:- begin_lpad), []) :-
 3264  prolog_load_context(module, M),
 3265  mc_input_mod(M),!,
 3266  assert(M:mc_on).
 3267
 3268mcintyre_expansion((:- end_lpad), []) :-
 3269  prolog_load_context(module, M),
 3270  mc_input_mod(M),!,
 3271  retractall(M:mc_on).
 3272
 3273mcintyre_expansion((Head:=Body),(H1:-Body)) :-
 3274  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3275% disjunctive fact with guassia distr
 3276  (Head \= ((mcintyre_expansion(_,_)) :- _ )),
 3277  Head=(H~val(Var)), !,
 3278  add_arg(H,Var,H1).
 3279
 3280
 3281mcintyre_expansion((Head:=Body),Clause) :-
 3282  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3283% fact with distr
 3284  (Head \= ((mcintyre_expansion(_,_)) :- _ )),
 3285  Head=(H~Distr0), !,
 3286  add_arg(H,Var,H1),
 3287  switch_finite(Distr0,Distr),
 3288  term_variables([Head,Body],VC),
 3289  get_next_rule_number(M,R),
 3290  (M:local_mc_setting(single_var,true)->
 3291    generate_clause_distr(H1,Body,[],M,R,Var,Distr,Clause)
 3292  ;
 3293    generate_clause_distr(H1,Body,VC,M,R,Var,Distr,Clause)
 3294  ).
 3295
 3296mcintyre_expansion((Head:=Body),(Head:- Body)) :-
 3297  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,!.
 3298
 3299mcintyre_expansion((Head :- Body), Clauses):-
 3300  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3301% disjunctive clause with more than one head atom senza depth_bound
 3302  Head = (_;_), !,
 3303  list2or(HeadListOr, Head),
 3304  process_head(HeadListOr, HeadList),
 3305  term_variables((Head :- Body),VC),
 3306  get_next_rule_number(M,R),
 3307  (M:local_mc_setting(single_var,true)->
 3308    generate_rules(HeadList,Body,HeadList,[],M,R,0,Clauses)
 3309  ;
 3310    generate_rules(HeadList,Body,HeadList,VC,M,R,0,Clauses)
 3311  ).
 3312
 3313mcintyre_expansion((Head:-Body),Clause) :-
 3314  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3315% rule with distr
 3316  (Head \= ((mcintyre_expansion(_,_)) :- _ )),
 3317  Head=(H:Distr0),
 3318  nonvar(Distr0),
 3319  \+ number(Distr0),
 3320  Distr0=..[D,Var|Pars],
 3321  is_dist(M,D),!,
 3322  Distr=..[D|Pars],
 3323  term_variables([Head,Body],VC0),
 3324  delete_equal(VC0,Var,VC),
 3325  get_next_rule_number(M,R),
 3326  (M:local_mc_setting(single_var,true)->
 3327    generate_clause_distr(H,Body,[],M,R,Var,Distr,Clause)
 3328  ;
 3329    generate_clause_distr(H,Body,VC,M,R,Var,Distr,Clause)
 3330  ).
 3331
 3332mcintyre_expansion((Head :- Body), []) :-
 3333% disjunctive clause with a single head atom con prob. 0 senza depth_bound --> la regola e' non  caricata nella teoria e non e' conteggiata in NR
 3334  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3335  ((Head:-Body) \= ((mcintyre_expansion(_,_) ):- _ )),
 3336  (Head = (_:P);Head=(P::_)),
 3337  ground(P),
 3338  P=:=0.0, !.
 3339
 3340
 3341mcintyre_expansion((Head :- Body), Clauses) :-
 3342% disjunctive clause with a single head atom senza DB, con prob. diversa da 1
 3343  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3344  ((Head:-Body) \= ((mcintyre_expansion(_,_) ):- _ )),
 3345  (Head = (H:_);Head = (_::H)), !,
 3346  list2or(HeadListOr, Head),
 3347  process_head(HeadListOr, HeadList),
 3348  term_variables((Head :- Body),VC),
 3349  get_next_rule_number(M,R),
 3350  (M:local_mc_setting(single_var,true)->
 3351    generate_clause_samp(H,Body,HeadList,[],M,R,0,Clauses)
 3352  ;
 3353    generate_clause_samp(H,Body,HeadList,VC,M,R,0,Clauses)
 3354  ).
 3355
 3356
 3357mcintyre_expansion(Head,Clauses) :-
 3358  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3359% disjunctive fact with more than one head atom senza db
 3360  Head=(_;_), !,
 3361  list2or(HeadListOr, Head),
 3362  process_head(HeadListOr, HeadList),
 3363  term_variables(Head,VC),
 3364  get_next_rule_number(M,R),
 3365  (M:local_mc_setting(single_var,true)->
 3366    generate_rules_fact(HeadList,HeadList,[],M,R,0,Clauses)
 3367  ;
 3368    generate_rules_fact(HeadList,HeadList,VC,M,R,0,Clauses)
 3369  ).
 3370
 3371mcintyre_expansion(Head,Clause) :-
 3372  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3373% fact with distr
 3374  (Head \= ((mcintyre_expansion(_,_)) :- _ )),
 3375  Head=(H~Distr0),
 3376  nonvar(Distr0),
 3377  !,
 3378  switch_finite(Distr0,Distr),
 3379  add_arg(H,Var,H1),
 3380  term_variables(Head,VC),
 3381  get_next_rule_number(M,R),
 3382  (M:local_mc_setting(single_var,true)->
 3383    generate_clause_distr(H1,true,[],M,R,Var,Distr,Clause)
 3384  ;
 3385    generate_clause_distr(H1,true,VC,M,R,Var,Distr,Clause)
 3386  ).
 3387
 3388mcintyre_expansion(Head,Clause) :-
 3389  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3390% disjunctive fact with dirichlet distr
 3391  (Head \= ((mcintyre_expansion(_,_)) :- _ )),
 3392  Head=(H:Distr0),
 3393  nonvar(Distr0),
 3394  \+ number(Distr0),
 3395  Distr0=..[D,Var|Pars],
 3396  is_dist(M,D),!,
 3397  Distr=..[D|Pars],
 3398  term_variables([Head],VC0),
 3399  delete_equal(VC0,Var,VC),
 3400  get_next_rule_number(M,R),
 3401  (M:local_mc_setting(single_var,true)->
 3402    generate_clause_distr(H,true,[],M,R,Var,Distr,Clause)
 3403  ;
 3404    generate_clause_distr(H,true,VC,M,R,Var,Distr,Clause)
 3405  ).
 3406
 3407
 3408mcintyre_expansion(Head,[]) :-
 3409  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3410% disjunctive fact with a single head atom con prob. 0
 3411  (Head \= ((mcintyre_expansion(_,_)) :- _ )),
 3412  (Head = (_:P); Head = (P::_)),
 3413  ground(P),
 3414  P=:=0.0, !.
 3415
 3416
 3417mcintyre_expansion(Head,H) :-
 3418  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3419% disjunctive fact with a single head atom con prob. 1, senza db
 3420  (Head \= ((mcintyre_expansion(_,_)) :- _ )),
 3421  (Head = (H:P);Head =(P::H)),
 3422  ground(P),
 3423  P=:=1.0, !.
 3424
 3425
 3426mcintyre_expansion(Head,Clause) :-
 3427  prolog_load_context(module, M),mc_input_mod(M),M:mc_on,
 3428% disjunctive fact with a single head atom e prob. generiche, senza db
 3429  (Head \= ((mcintyre_expansion(_,_)) :- _ )),
 3430  (Head=(H:_);Head=(_::H)), !,
 3431  list2or(HeadListOr, Head),
 3432  process_head(HeadListOr, HeadList),
 3433  term_variables(HeadList,VC),
 3434  get_next_rule_number(M,R),
 3435  (M:local_mc_setting(single_var,true)->
 3436    generate_clause_samp(H,true,HeadList,[],M,R,0,Clause)
 3437  ;
 3438    generate_clause_samp(H,true,HeadList,VC,M,R,0,Clause)
 3439  ).
 begin_lpad_pred is det
Initializes LPAD loading. /
 3448begin_lpad_pred:-
 3449  assert(mc_input_mod(user)),
 3450  assert(user:mc_on).
 end_lpad_pred is det
Terminates the cplint inference module. /
 3457end_lpad_pred:-
 3458  retractall(mc_input_mod(_)),
 3459  retractall(user:mc_on).
 3460
 3461list2or([],true):-!.
 3462
 3463list2or([X],X):-
 3464    X\=;(_,_),!.
 3465
 3466list2or([H|T],(H ; Ta)):-!,
 3467    list2or(T,Ta).
 3468
 3469
 3470list2and([],true):-!.
 3471
 3472list2and([X],X):-
 3473    X\=(_,_),!.
 3474
 3475list2and([H|T],(H,Ta)):-!,
 3476    list2and(T,Ta).
 3477
 3478
 3479builtin(average(_L,_Av)) :- !.
 3480builtin(mc_prob(_,_,_)) :- !.
 3481builtin(mc_prob(_,_)) :- !.
 3482builtin(mc_sample(_,_,_,_)) :- !.
 3483builtin(db(_)) :- !.
 3484builtin(G) :-
 3485  swi_builtin(G).
 3486
 3487
 3488is_dist(_M,D):-
 3489  member(D,[
 3490    finite,
 3491    discrete,
 3492    uniform,
 3493    uniform_dens,
 3494    gaussian,
 3495    dirichlet,
 3496    gamma,
 3497    beta,
 3498    poisson,
 3499    binomial,
 3500    multinomial,
 3501    geometric,
 3502    exponential,
 3503    negative_binomial,
 3504    user
 3505    ]),!.
 3506
 3507
 3508switch_finite(D0,D):-
 3509  D0=..[F,Arg0],
 3510  (F=finite;F=discrete),!,
 3511  maplist(swap,Arg0,Arg),
 3512  D=..[F,Arg].
 3513
 3514switch_finite(D,D).
 3515
 3516is_discrete(_M,D):-
 3517  functor(D,F,A),
 3518  member(F/A,[
 3519    finite/1,
 3520    discrete/1,
 3521    uniform/1
 3522    ]),!.
 3523
 3524is_discrete(M,D):-
 3525  functor(D,F,_A),
 3526  M:disc(F).
 swap(?Term1:term, ?Term2:term) is det
If Term1 is of the form A:B, then Term2 is of the form B:A. /
 3532swap(A:B,B:A).
 :Term:term ~= +B:term is det
equality predicate for distributional clauses /
 3539(M:A) ~= B :-
 3540  A=..L,
 3541  append(L,[B],L1),
 3542  A1=..L1,
 3543  M:A1.
 3544
 3545:- multifile sandbox:safe_meta/2. 3546
 3547sandbox:safe_meta(mcintyre:s(_,_), []).
 3548sandbox:safe_meta(mcintyre:mc_prob(_,_,_), []).
 3549sandbox:safe_meta(mcintyre:mc_prob(_,_), []).
 3550sandbox:safe_meta(mcintyre:mc_sample(_,_,_,_,_), []).
 3551sandbox:safe_meta(mcintyre:mc_rejection_sample(_,_,_,_,_), []).
 3552sandbox:safe_meta(mcintyre:mc_rejection_sample(_,_,_,_), []).
 3553sandbox:safe_meta(mcintyre:mc_mh_sample(_,_,_,_,_), []).
 3554sandbox:safe_meta(mcintyre:mc_mh_sample(_,_,_,_,_,_,_,_), []).
 3555sandbox:safe_meta(mcintyre:mc_gibbs_sample(_,_,_), []).
 3556sandbox:safe_meta(mcintyre:mc_gibbs_sample(_,_,_,_), []).
 3557sandbox:safe_meta(mcintyre:mc_gibbs_sample(_,_,_,_,_), []).
 3558sandbox:safe_meta(mcintyre:mc_sample(_,_,_), []).
 3559sandbox:safe_meta(mcintyre:mc_sample(_,_,_,_), []).
 3560sandbox:safe_meta(mcintyre:mc_sample_arg(_,_,_,_,_), []).
 3561sandbox:safe_meta(mcintyre:mc_sample_arg_first(_,_,_,_), []).
 3562sandbox:safe_meta(mcintyre:mc_sample_arg_first(_,_,_,_,_), []).
 3563sandbox:safe_meta(mcintyre:mc_sample_arg_one(_,_,_,_), []).
 3564sandbox:safe_meta(mcintyre:mc_sample_arg_one(_,_,_,_,_), []).
 3565sandbox:safe_meta(mcintyre:mc_sample_arg_raw(_,_,_,_), []).
 3566sandbox:safe_meta(mcintyre:mc_rejection_sample_arg(_,_,_,_,_,_), []).
 3567sandbox:safe_meta(mcintyre:mc_rejection_sample_arg(_,_,_,_,_), []).
 3568sandbox:safe_meta(mcintyre:mc_mh_sample_arg(_,_,_,_,_,_), []).
 3569sandbox:safe_meta(mcintyre:mc_mh_sample_arg(_,_,_,_,_), []).
 3570sandbox:safe_meta(mcintyre:mc_gibbs_sample_arg(_,_,_,_,_,_), []).
 3571sandbox:safe_meta(mcintyre:mc_gibbs_sample_arg(_,_,_,_,_), []).
 3572sandbox:safe_meta(mcintyre:mc_gibbs_sample_arg(_,_,_,_), []).
 3573sandbox:safe_meta(mcintyre:mc_expectation(_,_,_,_), []).
 3574sandbox:safe_meta(mcintyre:mc_mh_expectation(_,_,_,_,_), []).
 3575sandbox:safe_meta(mcintyre:mc_mh_expectation(_,_,_,_,_,_), []).
 3576sandbox:safe_meta(mcintyre:mc_rejection_expectation(_,_,_,_,_), []).
 3577sandbox:safe_meta(mcintyre:mc_gibbs_expectation(_,_,_,_), []).
 3578sandbox:safe_meta(mcintyre:mc_gibbs_expectation(_,_,_,_,_), []).
 3579sandbox:safe_meta(mcintyre:mc_gibbs_expectation(_,_,_,_,_,_), []).
 3580
 3581sandbox:safe_meta(mcintyre:mc_lw_sample(_,_,_,_), []).
 3582sandbox:safe_meta(mcintyre:mc_lw_sample_arg(_,_,_,_,_), []).
 3583sandbox:safe_meta(mcintyre:mc_lw_sample_arg_log(_,_,_,_,_), []).
 3584sandbox:safe_meta(mcintyre:mc_lw_expectation(_,_,_,_,_), []).
 3585sandbox:safe_meta(mcintyre:mc_particle_expectation(_,_,_,_,_), []).
 3586sandbox:safe_meta(mcintyre:mc_particle_sample_arg(_,_,_,_,_), []).
 3587sandbox:safe_meta(mcintyre:mc_particle_sample(_,_,_,_), []).
 3588sandbox:safe_meta(mcintyre:msw(_,_), []).
 3589
 3590sandbox:safe_meta(mcintyre:set_mc(_,_), []).
 3591sandbox:safe_meta(mcintyre:setting_mc(_,_), []).
 3592sandbox:safe_meta(mcintyre:set_sw(_,_) ,[]).
 3593
 3594:- license(artisticv2). 3595
 3596:- thread_local mcintyre_file/1. 3597
 3598user:term_expansion((:- mc), []) :-!,
 3599  prolog_load_context(source, Source),
 3600  asserta(mcintyre_file(Source)),
 3601  prolog_load_context(module, M),
 3602  retractall(local_mc_setting(_,_)),
 3603  findall(local_mc_setting(P,V),default_setting_mc(P,V),L),
 3604  assert_all(L,M,_),
 3605  assert(mc_input_mod(M)),
 3606  retractall(M:rule_n(_)),
 3607  assert(M:rule_n(0)),
 3608  dynamic((M:samp/3,M:mem/4,M:mc_on/0,M:sw/2,M:sw/3,M:sampled/3, M:sampled_g/2, M:sampled_g/1, M:disc/1,M:values/2)),
 3609  retractall(M:samp(_,_,_)),
 3610  style_check(-discontiguous).
 3611
 3612user:term_expansion(end_of_file, end_of_file) :-
 3613  mcintyre_file(Source),
 3614  prolog_load_context(source, Source),
 3615  retractall(mcintyre_file(Source)),
 3616  prolog_load_context(module, M),
 3617  mc_input_mod(M),!,
 3618  retractall(mc_input_mod(M)),
 3619  style_check(+discontiguous).
 3620
 3621user:term_expansion(In, Out) :-
 3622  \+ current_prolog_flag(xref, true),
 3623  mcintyre_file(Source),
 3624  prolog_load_context(source, Source),
 3625  mcintyre_expansion(In, Out)