1/*
2
3CEM
4
5Copyright (c) 2011, Fabrizio Riguzzi
6
7*/
8%:- set_prolog_flag(unknown,error).
9%:- set_prolog_flag(profiling,on).
10%:- set_prolog_flag(debug,on).
11:- set_prolog_flag(discontiguous_warnings,on).
12:- set_prolog_flag(single_var_warnings,on).
13:- set_prolog_flag(unknown,fail).
14%:-source.
15%:-yap_flag(gc_trace,very_verbose).
16:- use_module(inference,
17[find_deriv_inf1/3]).
18%:-consult(inference).
19:-use_module(library(rbtrees)).
20:-use_module(library(random)).
21:-use_module(library(avl)).
22:-use_module(library(lists)).
23
24%:-use_module(library(lpadsld)).
25:-load_foreign_files(['cplint'],[],init_my_predicates).
26
27:-dynamic setting/2,rule/5.
28
29
30setting(depth,3).
31setting(single_var,false).
32
33setting(sample_size,1000).
34/* Total number of examples in case in which the models in the kb file contain
35a prob(P). fact. In that case, one model corresponds to sample_size*P examples
36*/
37setting(equivalent_sample_size,100).
38/* equivalent samaple size for computing the BD score of the network refinements
39It is indicated with NPrime in the formulas on Heckerman, Geiger & Chickering
40paper */
41setting(epsilon_em,0.1).
42setting(epsilon_em_fraction,0.01).
43/* if the difference in log likelihood in two successive em iteration is smaller
44than epsilon_em, then em stops */
45setting(epsilon_sem,2).
46setting(random_restarts_number,1).
47/* number of random restarts of em */
48setting(verbosity,1).
49
50
51
52em(File):-
53  generate_file_names(File,FileKB,FileOut,FileL,FileLPAD),
54  reconsult(FileL),
55  load_models(FileKB,DB),
56  load_initial_model(FileLPAD,Model0),!,
57  set(verbosity,3),
58  statistics(cputime,[_,_]),
59  gen_ex(DB,[],DBE),
60  compute_parameters_EM(Model0,Model,SuffStats,CLL,DBE),
61  statistics(cputime,[_,CT]),
62  CTS is CT/1000,
63  format("Final CLL ~f~n",[CLL]),
64  format("Execution time ~f~n",[CTS]),
65  write_stats(user_output,SuffStats),
66  listing(setting/2),
67  format("Model:~n",[]),
68  write_model(Model,user_output),
69  open(FileOut,write,Stream),
70  format(Stream,"/* Final CLL ~f~n",[CLL]),
71  format(Stream,"Execution time ~f~n",[CTS]),
72  tell(Stream),
73  listing(setting/2),
74  write_stats(Stream,SuffStats),
75  format(Stream,"*/~n",[]),
76  write_model(Model,Stream),
77  told.
78
79gen_ex([],DBE,DBE).
80
81gen_ex([H|T],DB0,DB1):-
82        get_ouptut_atoms(O),
83  generate_goal(O,H,[],GL),
84  append(DB0,GL,DB2),
85  gen_ex(T,DB2,DB1).
86
87
88cycle_head([],[],_NR,_S,_NH,_PG,_CSetList,_N):-!.
89
90cycle_head([SSH0|T],[SSH1|T1],NR,S,NH,PG,CSetList,N):-
91  extract_relevant_C_sets(NR,S,NH,CSetList,CSL1),
92  (CSL1=[]->
93    SSH1 is SSH0
94  ;
95    build_formula(CSL1,Formula,[],Var),
96    var2numbers(Var,0,NewVar),
97    compute_prob(NewVar,Formula,Prob,0),
98    SSH1 is SSH0 +Prob/PG*N
99  ),
100  NH1 is NH+1,
101  cycle_head(T,T1,NR,S,NH1,PG,CSetList,N).
102
103cycle_head_neg([],[],_NR,_S,_NH,_NA,_PG,_CSetList,_N):-!.
104
105cycle_head_neg([SSH0|T],[SSH1|T1],NR,S,NH,NA,PG,CSetList,N):-
106        extract_relevant_C_sets_neg(NR,S,NH,NA,CSetList,CSL1),
107        (CSL1=[]->
108                SSH1 is SSH0%+0.001
109        ;
110                build_formula(CSL1,Formula,[],Var),
111                var2numbers(Var,0,NewVar),
112                compute_prob(NewVar,Formula,Prob,0),
113                (Prob>1 ->write(cyc),write(Prob),write(NewVar),nl;true),
114                SSH1 is SSH0 +(1-Prob)/PG*N
115        ),
116        NH1 is NH+1,
117        cycle_head_neg(T,T1,NR,S,NH1,NA,PG,CSetList,N).
118
119extract_relevant_C_sets_neg(NR,S,NH,NA,CS,CS1):-
120        neg_choice(0,NA,NH,NR,S,C),
121        append(CS,C,CS1).
122
123neg_choice(N,N,_NH,_NR,_S,[]):-!.
124
125neg_choice(NH,NA,NH,NR,S,L):-!,
126        N1 is NH+1,
127        neg_choice(N1,NA,NH,NR,S,L).
128
129neg_choice(N,NA,NH,NR,S,[[(N,NR,S)]|L]):-
130        N1 is N+1,
131        neg_choice(N1,NA,NH,NR,S,L).
132
133extract_relevant_C_sets(_NR,_S,_NH,[],[]):-!.
134
135extract_relevant_C_sets(NR,S,NH,[H|T],CS):-
136  member((NH1,NR,S),H),!,
137  extract_relevant_C_sets1(NR,S,NH,NH1,H,T,CS).
138
139extract_relevant_C_sets(NR,S,NH,[H|T],[H1|CS]):-
140  append(H,[(NH,NR,S)],H1),
141  extract_relevant_C_sets(NR,S,NH,T,CS).
142
143extract_relevant_C_sets1(NR,S,NH,NH1,_H,T,CS):-
144  NH1\=NH,!,
145  extract_relevant_C_sets(NR,S,NH,T,CS).
146
147extract_relevant_C_sets1(NR,S,NH,_NH1,H,T,[H|CS]):-
148  extract_relevant_C_sets(NR,S,NH,T,CS).
149
150
151
152/* EM start */
153compute_parameters_EM([],[],SuffStats,-1e200,_DB):-!,
154  rb_new(SuffStats).
155
156compute_parameters_EM(Model0,Model1,SuffStats1,CLL1,DB):-
157  setting(verbosity,Ver),
158  (Ver>0->
159    format("EM computation ~nInitial model:~n",[]),
160    write_model(Model0,user_output),
161    flush_output
162  ;
163    true
164  ),
165  (Ver>2->
166    format("Initial EM Iteration ~n",[]),
167    flush_output
168  ;
169    true
170  ),
171  randomize(Model0,ModelR),
172  em_iteration(ModelR,Model,SuffStats,CLL,DB),
173  (Ver>2->
174    format("CLL ~f~n",[CLL])
175  ;
176    true
177  ),
178  flush_output,
179  setting(random_restarts_number,N),
180  random_restarts(N,Model,SuffStats,CLL,Model1,SuffStats1,CLL1,DB),
181  (Ver>0->
182    format("Final CLL ~f~n",[CLL1]),
183    flush_output
184  ;
185    true
186  ).
187
188random_restarts(1,Model,SS,CLL,Model,SS,CLL,_DB):-!.
189
190random_restarts(N,Model0,SS0,CLL0,Model1,SS1,CLL1,DB):-
191  setting(verbosity,Ver),
192  (Ver>2->
193    setting(random_restarts_number,NMax),
194    Num is NMax-N+1,
195    format("Restart number ~d~n",[Num]),
196    flush_output
197  ;
198    true
199  ),
200  randomize(Model0,ModelR),
201  em_iteration(ModelR,ModelR1,SSR,CLLR,DB),
202  setting(verbosity,Ver),
203  (Ver>2->
204    format("CLL ~f~n",[CLLR])
205  ;
206    true
207  ),
208  N1 is N-1,
209  (CLLR>CLL0->
210    random_restarts(N1,ModelR1,SSR,CLLR,Model1,SS1,CLL1,DB)
211  ;
212    random_restarts(N1,Model0,SS0,CLL0,Model1,SS1,CLL1,DB)
213  ).
214
215randomize([],[]):-!.
216
217randomize([rule(N,V,NH,HL,BL,LogF)|T],[rule(N,V,NH,HL1,BL,LogF)|T1]):-
218  length(HL,L),
219  Int is 1.0/L,
220  randomize_head(Int,HL,0,HL1),
221  randomize(T,T1).
222
223randomize_head(_Int,['':_],P,['':PNull1]):-!,
224  PNull is 1.0-P,
225  (PNull>=0.0->
226    PNull1 =PNull
227  ;
228    PNull1=0.0
229  ).
230
231randomize_head(Int,[H:_|T],P,[H:PH1|NT]):-
232  PMax is 1.0-P,
233  random(0,PMax,PH1),
234  P1 is P+PH1,
235  randomize_head(Int,T,P1,NT).
236
237
238
239em_iteration(Model0,ModelPar,SuffStats1,CLL1,DB):-
240  compute_CLL_stats(Model0,DB,CLL0,SuffStats0),
241/*  setting(verbosity,Ver),
242  (Ver>2->
243    format("EM Iteration numer ~d~nCLL ~f~n",[N,CLL0]),
244    write_stats(user_output,SuffStats0)
245  ;
246    true
247  ),*/
248  cycle_EM(Model0,SuffStats0,CLL0,ModelPar,SuffStats1,CLL1,DB,1).
249
250cycle_EM(Model0,SuffStats0,CLL0,ModelPar,SuffStats,CLL,DB,N):-
251  m_step(Model0,SuffStats0,Model1),
252  compute_CLL_stats(Model1,DB,CLL1,SuffStats1),
253  setting(verbosity,Ver),
254  (Ver>2->
255    format("Iteration: ~d CLL ~f~n",[N,CLL1])
256  ;
257    true
258  ),
259  flush_output,
260%  write_stats(user_output,SuffStats1),
261%  statistics,
262  setting(epsilon_em,Epsilon_EM),
263  setting(epsilon_em_fraction,Epsilon_EM_Frac),
264  ((CLL1-CLL0<Epsilon_EM;(CLL1-CLL0)< - CLL0*Epsilon_EM_Frac)->
265    ModelPar=Model1,
266    SuffStats=SuffStats1,
267    CLL=CLL1,!
268  ;
269    N1 is N+1,!,
270    cycle_EM(Model1,SuffStats1,CLL1,ModelPar,SuffStats,CLL,DB,N1)
271  ).
272
273write_stats(S,SS):-
274  rb_visit(SS,Pairs),
275  format(S,"Suff stats~n",[]),
276  write_stats_list(S,Pairs).
277
278write_stats_list(S,[]):-nl(S),nl(S),!.
279
280write_stats_list(S,[R-d(D,N,I)|T]):-
281  format(S,"~d,~p,~f,~d~n",[R,D,N,I]),
282  write_stats_list(S,T).
283
284m_step([],_SS,[]):-!.
285
286m_step([rule(N,V,NH,HL,BL,LogF)|T],SS,[rule(N,V,NH,HL1,BL,LogF)|T1]):-
287  (rb_lookup(N,d(Distr,_NBT,_NI),SS)->
288    sum_list(Distr,NBT),
289    update_head(HL,Distr,NBT,HL1)
290  ;
291    HL1=HL
292  ),
293  m_step(T,SS,T1).
294
295update_head([],[],_N,[]).
296
297update_head([H:_P|T],[PU|TP],N,[H:P|T1]):-
298  P is PU/N,
299  update_head(T,TP,N,T1).
300
301
302/* EM end */
303
304
305/* Start of computation of log likelihood and sufficient stats */
306compute_CLL_stats(Model,DB,CLL,SuffStats1):-
307  assert_model(Model),
308  compute_CLL_stats_examples(DB,CLL,SuffStats1),
309  retract_model.
310
311assert_model([]):-!.
312
313assert_model([rule(N,V,NH,HL,BL,_LogF)|T]):-
314  assert_rules(HL,0,HL,BL,NH,N,V),
315  assertz(rule_by_num(N,V,NH,HL,BL)),
316  assert_model(T).
317
318retract_model:-
319  retractall(rule_by_num(_,_,_,_,_)),
320  retractall(rule(_,_,_,_,_,_,_,_)).
321
322assert_rules([],_Pos,_HL,_BL,_Nh,_N,_V1):-!.
323
324assert_rules(['':_P],_Pos,_HL,_BL,_Nh,_N,_V1):-!.
325
326assert_rules([H:P|T],Pos,HL,BL,NH,N,V1):-
327        assertz(rule(H,P,Pos,N,V1,NH,HL,BL)),
328        Pos1 is Pos+1,
329        assert_rules(T,Pos1,HL,BL,NH,N,V1).
330
331compute_CLL_stats_examples(DB,CLL,SuffStats1):-
332  rb_new(SuffStats0),
333  compute_CLL_stats_cplint(DB,0,CLL,SuffStats0,SuffStats1).
334
335get_ouptut_atoms(O):-
336  findall((A/Ar),output((A/Ar)),O).
337
338generate_goal([],_H,G,G):-!.
339
340generate_goal([P/A|T],H,G0,G1):-
341  functor(Pred,P,A),
342  Pred=..[P|Rest],
343  Pred1=..[P,H|Rest],
344  findall(Pred1,call(Pred1),L),
345  findall(\+ Pred1,call(neg(Pred1)),LN),
346  append(G0,L,G2),
347  append(G2,LN,G3),
348  generate_goal(T,H,G3,G1).
349
350compute_CLL_stats_cplint([],CLL,CLL,S,S):-!.
351
352compute_CLL_stats_cplint([\+ H|T],CLL0,CLL1,Stats0,Stats1):-!,
353        setting(verbosity,V),
354        (V>3->
355                write(user_error,(\+ H)),nl(user_error),flush_output
356        ;
357                true
358        ),
359        s([H],CL,CSetList,PG),!,
360        (PG=:=1.0->
361                CLL2=CLL0,
362                Stats2=Stats0
363        ;
364                (prob(H,P)->
365                        setting(sample_size,NTot),
366                        N is P*NTot
367                ;
368                        N=1
369                ),
370                PG1 is 1-PG,
371                CLL2 is CLL0+log(PG1)*N,
372                collect_stats_cplint_neg(CL,PG1,CSetList,N,Stats0,Stats2)
373        ),
374        compute_CLL_stats_cplint(T,CLL2,CLL1,Stats2,Stats1).
375
376compute_CLL_stats_cplint([H|T],CLL0,CLL1,Stats0,Stats1):-
377        setting(verbosity,V),
378        (V>3->
379                write(user_error,H),nl(user_error),flush_output
380        ;
381                true
382        ),
383  s([H],CL,CSetList,PG),!,
384  (PG=0.0->
385    CLL2=CLL0,
386    Stats2=Stats0
387  ;
388    (prob(H,P)->
389      setting(sample_size,NTot),
390      N is P*NTot
391    ;
392      N=1
393    ),
394    CLL2 is CLL0+log(PG)*N,
395    collect_stats_cplint(CL,PG,CSetList,N,Stats0,Stats2)
396  ),
397  compute_CLL_stats_cplint(T,CLL2,CLL1,Stats2,Stats1).
398
399
400
401s(GoalsList,GroundLpad,CSets,Prob):-
402  solve(GoalsList,GroundLpad,CSets,Prob).
403
404solve(GoalsList,GroundLpad,LDup,Prob):-
405        setting(depth,D),
406        findall(Deriv,inference:find_deriv_inf1(GoalsList,D,Deriv),LDup),
407        (LDup=[]->
408                Prob=0.0,
409                GroundLpad=[]
410        ;
411                append(LDup,L0),
412                remove_head(L0,L1),
413                remove_duplicates(L1,L2),
414                build_ground_lpad(L2,GroundLpad),
415                build_formula(LDup,Formula,[],Var),
416                var2numbers(Var,0,NewVar),
417                compute_prob(NewVar,Formula,Prob,0),
418                true
419        ).
420
421collect_stats_cplint([],_PG,_CSetList,_N,Stats,Stats):-!.
422
423collect_stats_cplint([(R,S,Head,_Body)|T],PG,CSetList,N,Stats0,Stats1):-
424  (rb_lookup(R,d(Distr0,N1,NInst1),Stats0)->
425    cycle_head(Distr0,Distr,R,S,0,PG,CSetList,N),
426    N2 is N+N1,
427    rb_update(Stats0,R,d(Distr,N2,NInst1),Stats2)
428  ;
429    length(Head,LH),
430    list0(0,LH,Distr0),
431    cycle_head(Distr0,Distr,R,S,0,PG,CSetList,N),
432    rb_insert(Stats0,R,d(Distr,N,1),Stats2)
433  ),
434  collect_stats_cplint(T,PG,CSetList,N,Stats2,Stats1).
435
436collect_stats_cplint_neg([],_PG,_CSetList,_N,Stats,Stats):-!.
437
438collect_stats_cplint_neg([(R,S,Head,_Body)|T],PG,CSetList,N,Stats0,Stats1):-
439        length(Head,NA),
440        (rb_lookup(R,d(Distr0,N1,NInst1),Stats0)->
441                cycle_head_neg(Distr0,Distr,R,S,0,NA,PG,CSetList,N),
442                N2 is N+N1,
443                rb_update(Stats0,R,d(Distr,N2,NInst1),Stats2)
444        ;
445                length(Head,LH),
446                list0(0,LH,Distr0),
447                cycle_head_neg(Distr0,Distr,R,S,0,NA,PG,CSetList,N),
448                rb_insert(Stats0,R,d(Distr,N,1),Stats2)
449        ),
450        collect_stats_cplint_neg(T,PG,CSetList,N,Stats2,Stats1).
451
452/* build_formula(LC,Formula,VarIn,VarOut) takes as input a set of C sets
453LC and a list of Variables VarIn and returns the formula and a new list
454of variables VarOut
455Formula is of the form [Term1,...,Termn]
456Termi is of the form [Factor1,...,Factorm]
457Factorj is of the form (Var,Value) where Var is the index of
458the multivalued variable Var and Value is the index of the value
459*/
460build_formula([],[],Var,Var,C,C).
461
462build_formula([D|TD],[F|TF],VarIn,VarOut,C0,C1):-
463        length(D,NC),
464        C2 is C0+NC,
465        build_term(D,F,VarIn,Var1),
466        build_formula(TD,TF,Var1,VarOut,C2,C1).
467
468build_formula([],[],Var,Var).
469
470build_formula([D|TD],[F|TF],VarIn,VarOut):-
471        build_term(D,F,VarIn,Var1),
472        build_formula(TD,TF,Var1,VarOut).
473
474build_term([],[],Var,Var).
475
476build_term([(_,pruned,_)|TC],TF,VarIn,VarOut):-!,
477        build_term(TC,TF,VarIn,VarOut).
478
479build_term([(N,R,S)|TC],[[NVar,N]|TF],VarIn,VarOut):-
480        (nth0_eq(0,NVar,VarIn,(R,S))->
481                Var1=VarIn
482        ;
483                append(VarIn,[(R,S)],Var1),
484                length(VarIn,NVar)
485        ),
486        build_term(TC,TF,Var1,VarOut).
487
488/* nth0_eq(PosIn,PosOut,List,El) takes as input a List,
489an element El and an initial position PosIn and returns in PosOut
490the position in the List that contains an element exactly equal to El
491*/
492nth0_eq(N,N,[H|_T],El):-
493        H==El,!.
494
495nth0_eq(NIn,NOut,[_H|T],El):-
496        N1 is NIn+1,
497        nth0_eq(N1,NOut,T,El).
498
499/* var2numbers converts a list of couples (Rule,Substitution) into a list
500of triples (N,NumberOfHeadsAtoms,ListOfProbabilities), where N is an integer
501starting from 0 */
502var2numbers([],_N,[]).
503
504var2numbers([(R,S)|T],N,[[N,ValNumber,Probs]|TNV]):-
505        find_probs(R,S,Probs),
506        length(Probs,ValNumber),
507        N1 is N+1,
508        var2numbers(T,N1,TNV).
509
510find_probs(R,S,Probs):-
511        rule_by_num(R,S,_N,Head,_Body),
512        get_probs(Head,Probs).
513
514get_probs(uniform(_A:1/Num,_P,_Number),ListP):-
515        Prob is 1/Num,
516        list_el(Num,Prob,ListP).
517
518get_probs([],[]).
519
520get_probs([_H:P|T],[P1|T1]):-
521        P1 is P,
522        get_probs(T,T1).
523
524list_el(0,_P,[]):-!.
525
526list_el(N,P,[P|T]):-
527        N1 is N-1,
528        list_el(N1,P,T).
529
530sum(_NS,[],[],[]):-!.
531
532sum(NS,[H0|T0],[H1|T1],[H2|T2]):-
533  H2 is H0+H1*NS,
534  sum(NS,T0,T1,T2).
535
536times(_NS,[],[]):-!.
537
538times(NS,[H0|T0],[H1|T1]):-
539  H1 is H0*NS,
540  times(NS,T0,T1).
541
542/* End of computation of log likelihood and sufficient stats */
543
544/* Utility predicates */
545generate_file_names(File,FileKB,FileOut,FileL,FileLPAD):-
546    generate_file_name(File,".kb",FileKB),
547    generate_file_name(File,".rules",FileOut),
548    generate_file_name(File,".cpl",FileLPAD),
549    generate_file_name(File,".l",FileL).
550
551generate_file_name(File,Ext,FileExt):-
552    name(File,FileString),
553    append(FileString,Ext,FileStringExt),
554    name(FileExt,FileStringExt).
555
556
557set(Parameter,Value):-
558  retract(setting(Parameter,_)),
559  assert(setting(Parameter,Value)).
560
561load_initial_model(File,Model):-
562  open(File,read,S),
563  read_clauses(S,C),
564  close(S),
565  process_clauses(C,1,_N,[],Model).
566
567process_clauses([(end_of_file,[])],N,N,Model,Model).
568
569process_clauses([((H:-B),_V)|T],N,N2,Model0,Model1):-
570        H=(db(A)),!,
571  assert((A:-B)),
572  process_clauses(T,N,N2,Model0,Model1).
573
574process_clauses([((H:-B),V)|T],N,N2,Model0,[rule(N,V1,NH,HL,BL,0)|Model1]):-
575  H=(_;_),!,
576  list2or(HL1,H),
577  process_head(HL1,HL,VI),
578  list2and(BL0,B),
579  add_int_atom(BL0,BL,VI),
580  length(HL,LH),
581  listN(0,LH,NH),
582  N1 is N+1,
583  (setting(single_var,true)->
584    V1=[]
585  ;
586    V1=V
587  ),
588%  assertz(rule(N,V,NH,HL,BL)),
589  process_clauses(T,N1,N2,Model0,Model1).
590
591process_clauses([((H:-B),V)|T],N,N2,Model0,[rule(N,V1,NH,HL,BL,0)|Model1]):-
592  H=(_:_),!,
593  list2or(HL1,H),
594  process_head(HL1,HL,VI),
595  list2and(BL0,B),
596  add_int_atom(BL0,BL,VI),
597  length(HL,LH),
598  listN(0,LH,NH),
599  (setting(single_var,true)->
600    V1=[]
601  ;
602    V1=V
603  ),
604  N1 is N+1,
605%  assertz(rule(N,V1,NH,HL,BL)),
606  process_clauses(T,N1,N2,Model0,Model1).
607
608process_clauses([((H:-B),V)|T],N,N2,Model0,[rule(N,V1,NH,HL,BL,0)|Model1]):-!,
609  process_head([H:1.0],HL,VI),
610  list2and(BL0,B),
611  add_int_atom(BL0,BL,VI),
612  length(HL,LH),
613  listN(0,LH,NH),
614  (setting(single_var,true)->
615    V1=[]
616  ;
617    V1=V
618  ),
619  N1 is N+1,
620%  assertz(rule(N,V1,NH,HL,BL)),
621  process_clauses(T,N1,N2,Model0,Model1).
622
623process_clauses([(H,V)|T],N,N2,Model0,[rule(N,V1,NH,HL,[],0)|Model1]):-
624  H=(_;_),!,
625  list2or(HL1,H),
626  process_head(HL1,HL,_VI),
627  length(HL,LH),
628  listN(0,LH,NH),
629  (setting(single_var,true)->
630    V1=[]
631  ;
632    V1=V
633  ),
634  N1 is N+1,
635%  assertz(rule(N,V,NH,HL,[])),
636  process_clauses(T,N1,N2,Model0,Model1).
637
638process_clauses([(H,V)|T],N,N2,Model0,[rule(N,V1,NH,HL,[],0)|Model1]):-
639  H=(_:_),!,
640  list2or(HL1,H),
641  process_head(HL1,HL,_VI),
642  length(HL,LH),
643  listN(0,LH,NH),
644  (setting(single_var,true)->
645    V1=[]
646  ;
647    V1=V
648  ),
649  N1 is N+1,
650%  assertz(rule(N,V,NH,HL,[])),
651  process_clauses(T,N1,N2,Model0,Model1).
652
653process_clauses([(H,V)|T],N,N2,Model0,[rule(N,V1,NH,HL,[],0)|Model1]):-
654  process_head([H:1.0],HL,_VI),
655  length(HL,LH),
656  listN(0,LH,NH),
657  (setting(single_var,true)->
658    V1=[]
659  ;
660    V1=V
661  ),
662  N1 is N+1,
663%  assertz(rule(N,V,NH,HL,[])),
664  process_clauses(T,N1,N2,Model0,Model1).
665
666
667/* if the annotation in the head are not ground, the null atom is not added
668and the eventual formulas are not evaluated */
669
670process_head([H:P|T],NHL,VI):-!,
671    process_head_prob([H:P|T],0.0,NHL,VI).
672
673process_head(HL,NHL,VI):-
674    process_head_random(HL,0.0,NHL,VI).
675
676process_head_random([],P,['':PNull1],_VI):-
677  PNull is 1.0-P,
678  (PNull>=0.0->
679    PNull1 =PNull
680  ;
681    PNull1=0.0
682  ).
683
684process_head_random([H|T],P,[H1:PH1|NT],VI):-
685  add_int_atom([H],[H1],VI),
686  PMax is 1.0-P,
687  random(0,PMax,PH1),
688  P1 is P+PH1,
689  process_head_random(T,P1,NT,VI).
690
691
692process_head_prob([H:PH],P,[H1:PH1,'':PNull1],VI):-
693  add_int_atom([H],[H1],VI),
694  PH1 is PH,
695  PNull is 1.0-P-PH1,
696  (PNull>=0.0->
697    PNull1 =PNull
698  ;
699    PNull1=0.0
700  ).
701
702process_head_prob([H:PH|T],P,[H1:PH1|NT],VI):-
703  add_int_atom([H],[H1],VI),
704  PH1 is PH,
705  P1 is P+PH1,
706  process_head_prob(T,P1,NT,VI).
707
708
709add_int_atom([],[],_VI).
710
711add_int_atom([\+ H|T],[\+ H|T1],VI):-
712    inference:builtin(H),!,
713  add_int_atom(T,T1,VI).
714
715add_int_atom([\+ H|T],[\+ H1|T1],VI):-!,
716  H=..[F|Args],
717  H1=..[F,VI|Args],
718  add_int_atom(T,T1,VI).
719
720add_int_atom([H|T],[H|T1],VI):-
721  inference:builtin(H),!,
722  add_int_atom(T,T1,VI).
723
724add_int_atom([H|T],[H1|T1],VI):-
725  H=..[F|Args],
726  H1=..[F,VI|Args],
727  add_int_atom(T,T1,VI).
728
729/* predicates for reading in the program clauses */
730read_clauses(S,Clauses):-
731    read_clauses_ground_body(S,Clauses).
732
733
734read_clauses_ground_body(S,[(Cl,V)|Out]):-
735  read_term(S,Cl,[variable_names(V)]),
736  (Cl=end_of_file->
737    Out=[]
738  ;
739    read_clauses_ground_body(S,Out)
740  ).
741
742
743
744
745listN(N,N,[]):-!.
746
747listN(NIn,N,[NIn|T]):-
748  N1 is NIn+1,
749  listN(N1,N,T).
750
751list0(N,N,[]):-!.
752
753list0(NIn,N,[0|T]):-
754  N1 is NIn+1,
755  list0(N1,N,T).
756
757/* end of predicates for parsing an input file containing a program */
758
759
760load_models(File,ModulesList):-
761    open(File,read,Stream),
762    read_models(Stream,ModulesList),
763    close(Stream).
764
765read_models(Stream,[Name1|Names]):-
766    read(Stream,begin(model(Name))),!,
767    (number(Name)->
768        name(Name,NameStr),
769        append("i",NameStr,Name1Str),
770        name(Name1,Name1Str)
771    ;
772        Name1=Name
773    ),
774    read_all_atoms(Stream,Name1),
775    read_models(Stream,Names).
776
777read_models(_S,[]).
778
779read_all_atoms(Stream,Name):-
780    read(Stream,At),
781    At \=end(model(_Name)),!,
782    (At=neg(Atom)->
783      Atom=..[Pred|Args],
784      Atom1=..[Pred,Name|Args],
785      assertz(neg(Atom1))
786    ;
787      At=..[Pred|Args],
788      Atom1=..[Pred,Name|Args],
789      assertz(Atom1)
790    ),
791    read_all_atoms(Stream,Name).
792
793read_all_atoms(_S,_N).
794
795
796list2or([],true):-!.
797
798list2or([X],X):-
799    X\=;(_,_),!.
800
801list2or([H|T],(H ; Ta)):-!,
802    list2or(T,Ta).
803
804list2and([],true):-!.
805
806list2and([X],X):-
807    X\=(_,_),!.
808
809list2and([H|T],(H,Ta)):-!,
810    list2and(T,Ta).
811
812
813write_model([],_Stream):-!.
814
815write_model([rule(_N,_V,_NH,HL,BL,_LogF)|Rest],Stream):-
816  copy_term((HL,BL),(HL1,BL1)),
817    numbervars((HL1,BL1),0,_M),
818    write_disj_clause(Stream,(HL1:-BL1)),
819    format(Stream,".~n~n",[]),
820    write_model(Rest,Stream).
821
822
823write_disj_clause(S,(H:-[])):-!,
824  write_head(S,H).
825
826write_disj_clause(S,(H:-B)):-
827  write_head(S,H),
828  write(S,' :-'),
829  nl(S),
830  write_body(S,B).
831
832write_head(S,[A:1.0|_Rest]):-!,
833  remove_int_atom(A,A1),
834    format(S,"~p",[A1]).
835
836write_head(S,[A:P,'':_P]):-!,
837  remove_int_atom(A,A1),
838    format(S,"~p:~f",[A1,P]).
839
840write_head(S,[A:P|Rest]):-
841  remove_int_atom(A,A1),
842    format(S,"~p:~f ; ",[A1,P]),
843    write_head(S,Rest).
844
845write_body(S,[\+ A]):-!,
846  remove_int_atom(A,A1),
847    format(S,"\t\\+ ~p",[A1]).
848
849write_body(S,[A]):-!,
850  remove_int_atom(A,A1),
851    format(S,"\t~p",[A1]).
852
853write_body(S,[\+ A|T]):-!,
854  remove_int_atom(A,A1),
855    format(S,"\t\\+ ~p,~n",[A1]),
856    write_body(S,T).
857
858write_body(S,[A|T]):-
859  remove_int_atom(A,A1),
860    format(S,"\t~p,~n",[A1]),
861    write_body(S,T).
862
863
864remove_int_atom(A,A1):-
865  A=..[F,_|T],
866  A1=..[F|T].
867
868build_ground_lpad([],[]):-!.
869
870build_ground_lpad([(R,S)|T],[(R,S,Head,Body)|T1]):-
871  user:rule_by_num(R,S,_,Head,Body),
872  build_ground_lpad(T,T1).
873
874
875remove_head([],[]).
876
877remove_head([(_N,R,S)|T],[(R,S)|T1]):-
878  remove_head(T,T1).
879
880
881append_all([],L,L):-!.
882
883append_all([LIntH|IntT],IntIn,IntOut):-
884  append(IntIn,LIntH,Int1),
885  append_all(IntT,Int1,IntOut).
886
887