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 ]).
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
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 ).
253mc_load_file(File):-
254 must_be(atom,File),
255 begin_lpad_pred,
256 user:consult(File),
257 end_lpad_pred.
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:), 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.
Options is a list of options, the following are recognised by mc_prob/3:
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 ).
524mc_prob(M:Goal,P):-
525 must_be(nonvar,Goal),
526 must_be(var,P),
527 mc_prob(M:Goal,P,[]).
Options is a list of options, the following are recognised by mc_sample/4:
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 ).
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,[]).
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).
Options is a list of options, the following are recognised by mc_rejection_sample/5:
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).
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,[]).
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(\+ _)).
Options is a list of options, the following are recognised by mc_gibbs_sample/4:
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).
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,_).
Options is a list of options, the following are recognised by mc_gibbs_sample/5:
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).
Options is a list of options, the following are recognised by mc_gibbs_sample_arg/5:
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 ).
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,[]).
Options is a list of options, the following are recognised by mc_gibbs_sample_arg/6:
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).
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).
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).
Options is a list of options, the following are recognised by mc_mh_sample/5:
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).
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,[]).
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).
Options is a list of options, the following are recognised by mc_sample_arg/5:
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 ).
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,[]).
Options is a list of options, the following are recognised by mc_rejection_sample_arg/6:
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 ).
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,[]).
Options is a list of options, the following are recognised by mc_mh_sample_arg/6:
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 ).
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,[]).
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).
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).
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).
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).
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).
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).
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).
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).
Options is a list of options, the following are recognised by mc_sample_arg_first/5:
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 ).
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).
Options is a list of options, the following are recognised by mc_sample_arg_one/5:
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 ).
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).
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).
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.
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).
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,[]).
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).
Options is a list of options, the following are recognised by mc_mh_expectation/6:
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).
Options is a list of options, the following are recognised by mc_mh_expectation/6:
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).
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:,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].
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)).
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).
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.
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.
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).
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).
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).
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.
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).
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).
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.
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)).
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.
2931finite(D,M,S):- 2932 discrete(D,M,S). 2933 2934finite(D,M,S,Dens):- 2935 discrete(D,M,S,Dens).
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).
**/
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).
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).
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).
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).
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)).
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.
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)).
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 ).
3448begin_lpad_pred:-
3449 assert(mc_input_mod(user)),
3450 assert(user:mc_on).
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).
3532swap(A:B,B:A).
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 3598userterm_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 3612userterm_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 3621userterm_expansion(In, Out) :- 3622 \+ current_prolog_flag(xref, true), 3623 mcintyre_file(Source), 3624 prolog_load_context(source, Source), 3625 mcintyre_expansion(In, Out)
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.