slipcover

   16/*
   17
   18SLIPCOVER
   19
   20Copyright (c) 2016, Fabrizio Riguzzi and Elena Bellodi
   21
   22*/
   23:-module(liftcover,[set_lift/2,setting_lift/2,
   24  induce_lift/2,induce_par_lift/2,test_lift/7,
   25  filter_rules/2,filter_rules/3,sort_rules/2,
   26  remove_zero/2,
   27  op(500,fx,#),op(500,fx,'-#'),
   28  test_prob_lift/6,
   29  prob_lift/2,prob_lift/3,
   30  explain_lift/2,explain_lift/3,
   31  ranked_answers/3,
   32  ranked_answers/4,
   33  rank/3,
   34  rank_answer/3,
   35  rank_answer/4,
   36  hits_at_k/6,
   37  hits_at_k/7,
   38  rank_ex/3,
   39  rank_exs/4,
   40  inst_exs/4
   41  ]).   42:-use_module(library(auc)).   43:-use_module(library(lists)).   44:-use_module(library(random)).   45:-use_module(library(system)).   46:-use_module(library(terms)).   47:-use_module(library(ordsets)).   48:-use_module(library(apply)).   49:-use_module(library(settings)).   50:-use_module(library(clpfd), [transpose/2]).   51:-absolute_file_name(library(lbfgs),F,[solutions(all)]),atomic_concat(F,'.pl',Fpl),exists_file(Fpl),use_module(library(lbfgs));true.   52%:-use_foreign_library(foreign(bddem),install).
   53:-set_prolog_flag(unknown,warning).   54
   55%:-multifile setting_lift/2.
   56%:-use_module(library(sandbox)).
   57
   58
   59:- dynamic lift_input_mod/1.   60
   61:- thread_local v/3, lift_input_mod/1, local_setting/2, rule_lift_n/1.   62
   63:- meta_predicate induce_lift(:,-).   64:- meta_predicate induce_rules(:,-).   65:- meta_predicate induce_par_lift(:,-).   66:- meta_predicate induce_parameters(:,-).   67
   68
   69
   70:- meta_predicate test_lift(:,+,-,-,-,-,-).   71:- meta_predicate test_prob_lift(:,+,-,-,-,-).   72:- meta_predicate prob_lift(:,-).   73:- meta_predicate prob_lift(:,+,-).   74:- meta_predicate explain_lift(:,-).   75:- meta_predicate explain_lift(:,+,-).   76:- meta_predicate ranked_answers(:,-,-).   77:- meta_predicate ranked_answers(:,-,+,-).   78:- meta_predicate rank_answer(:,+,-).   79:- meta_predicate rank_answer(:,+,+,-).   80:- meta_predicate hits_at_k(:,+,+,+,-,-).   81:- meta_predicate hits_at_k(:,+,+,+,+,-,-).   82:- meta_predicate rank_ex(:,+,+).   83:- meta_predicate rank_exs(:,+,+,+).   84:- meta_predicate inst_exs(:,+,+,+).   85:- meta_predicate set_lift(:,+).   86:- meta_predicate setting_lift(:,-).   87:- meta_predicate filter_rules(:,-).   88
   89
   90
   91
   92:- setting(max_threads, integer, 256,
   93    'Maximum number of threads to use in scoring clause refinements and parameter learning').   94:- setting(allow_gpu, boolean, true,
   95    'Whether the use of GPU is allowed').   96
   97
   98
   99default_setting_lift(eps,0.0001).
  100default_setting_lift(eps_f,0.00001).
  101
  102
  103default_setting_lift(random_restarts_number,1). % number of random restarts of parameter learning
  104default_setting_lift(random_restarts_number_str_learn,1). % number of random restarts during structure learning of parameter learning for single clauses
  105default_setting_lift(iter,100).
  106default_setting_lift(d,1).
  107default_setting_lift(verbosity,1).
  108default_setting_lift(logzero,log(0.000001)).
  109default_setting_lift(megaex_bottom,1).
  110default_setting_lift(initial_clauses_per_megaex,1).
  111default_setting_lift(max_iter,10).
  112default_setting_lift(max_var,4).
  113default_setting_lift(maxdepth_var,2).
  114default_setting_lift(beamsize,100).
  115default_setting_lift(max_clauses,1000).
  116default_setting_lift(max_body_length,100).
  117default_setting_lift(neg_literals,false).
  118
  119default_setting_lift(specialization,bottom).
  120/* allowed values: mode,bottom */
  121
  122default_setting_lift(seed,rand(10,1231,3032)).
  123default_setting_lift(neg_ex,cw).
  124
  125
  126default_setting_lift(epsilon_parsing, 1e-5).
  127
  128
  129
  130
  131default_setting_lift(zero,0.000001).
  132default_setting_lift(minus_infinity,-1.0e20).
  133
  134default_setting_lift(regularization,l1). % regularization: no, l1, l2, bayesian 
  135default_setting_lift(gamma,10). % set the value of gamma for regularization l1 and l2
  136default_setting_lift(ab,[0,10]). % set the value of a and b for regularization baysian
  137default_setting_lift(min_probability,1e-5).  % Threshold of the probability under which the clause is dropped
  138default_setting_lift(parameter_learning,em). % parameter learning algorithm: em_python, em, lbfgs, gd_python, gd 
  139default_setting_lift(max_initial_weight,0.5). % initial weights of dphil in [-0.5 0.5]
  140
  141default_setting_lift(parameter_update,fixed_learning_rate). % values: fixed_learning_rate,adam
  142default_setting_lift(eta,0.01). % fixed learning rate
  143default_setting_lift(adam_params,[0.001,0.9,0.999,1e-8]). % default Adam hyper-pameters
  144default_setting_lift(processor,cpu). % where to run em_python and gd_python: cpu or gpu
  145default_setting_lift(threads,1). % number of threads to use in scoring clause refinements and parameter learning
  146default_setting_lift(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.

 induce_lift(:TrainFolds:list_of_atoms, -P:probabilistic_program) is det

The predicate performs structure learning using the folds indicated in TrainFolds for training. It returns in P the learned probabilistic program. /

  155induce_lift(TrainFolds,P):-
  156  induce_rules(TrainFolds,P0),
  157  rules2terms(P0,P).

 test_lift(:P:probabilistic_program, +TestFolds:list_of_atoms, -LL:float, -AUCROC:float, -ROC:dict, -AUCPR:float, -PR:dict) is det

The predicate takes as input in P a probabilistic program, tests P on the folds indicated in TestFolds and returns the log likelihood of the test examples in LL, the area under the Receiver Operating Characteristic curve in AUCROC, a dict containing the points of the ROC curve in ROC, the area under the Precision Recall curve in AUCPR and a dict containing the points of the PR curve in PR /

  168test_lift(P,TestFolds,LL,AUCROC,ROC,AUCPR,PR):-
  169  test_prob_lift(P,TestFolds,_NPos,_NNeg,LL,LG),
  170  compute_areas_diagrams(LG,AUCROC,ROC,AUCPR,PR).

 test_prob_lift(:P:probabilistic_program, +TestFolds:list_of_atoms, -NPos:int, -NNeg:int, -LL:float, -Results:list) is det

The predicate takes as input in P a probabilistic program, tests P on the folds indicated in TestFolds and returns the number of positive examples in NPos, the number of negative examples in NNeg, the log likelihood in LL and in Results a list containing the probabilistic result for each query contained in TestFolds. /

  181test_prob_lift(M:P,TestFolds,NPos,NNeg,CLL,Results) :-
  182  write2(M,'Testing\n'),
  183  make_dynamic(M),
  184  process_clauses(P,M,PRules),
  185  generate_clauses(PRules,M,0,Prog),
  186  (M:bg(RBG0)->
  187    process_clauses(RBG0,M,RBG),
  188    generate_clauses_bg(RBG,ClBG),
  189    assert_all(ClBG,M,ClBGRef)
  190  ;
  191    true
  192  ),
  193  findall(Exs,(member(F,TestFolds),M:fold(F,Exs)),L),
  194  append(L,DB),
  195  test_no_area(DB,M,Prog,NPos,NNeg,CLL,Results),
  196  (M:bg(RBG0)->
  197    retract_all(ClBGRef)
  198  ;
  199    true
  200  ).
  201
  202induce_rules(M:Folds,R):-
  203  load_python_module(M),
  204  make_dynamic(M),
  205  M:local_setting(seed,Seed),
  206  setrand(Seed),
  207  findall(Exs,(member(F,Folds),M:fold(F,Exs)),L),
  208  append(L,DB),
  209  (M:bg(RBG0)->
  210    process_clauses(RBG0,M,RBG),
  211    generate_clauses_bg(RBG,ClBG),
  212    assert_all(ClBG,M,ClBGRef)
  213  ;
  214    true
  215  ),
  216  find_ex(DB,M,Pos,Neg,NPos,_Neg),
  217  M:local_setting(megaex_bottom, NumMB),
  218  (NPos >= NumMB ->
  219      true
  220    ;
  221      format2(M,"~nWARN: Number of required bottom clauses is greater than the number of training examples!~n. The number of required bottom clauses will be equal to the number of training examples", []),
  222      M:set_lift(megaex_bottom, NPos)
  223  ),
  224  statistics(walltime,[_,_]),
  225  (M:local_setting(specialization,bottom)->
  226    M:local_setting(megaex_bottom,MB),
  227    deduct(MB,M,DB,[],InitialTheory),
  228    remove_duplicates(InitialTheory,R1)
  229  ;
  230    get_head_atoms(O,M),
  231    generate_top_cl(O,M,R1)
  232  ),
  233  learn_struct(Pos,Neg,M,R1,R2,Score),
  234  sort_rules_int(R2,R),
  235  statistics(walltime,[_,WT]),
  236  WTS is WT/1000,
  237  write2(M,'\n\n'),
  238  format2(M,'/* Final score ~f~n',[Score]),
  239  format2(M,'Wall time ~f */~n',[WTS]),
  240  write_rules2(M,R,user_output),
  241  (M:bg(RBG0)->
  242    retract_all(ClBGRef)
  243  ;
  244    true
  245  ),
  246  retractall(M:ref_clause(_)),
  247  retractall(M:ref(_)).

 sort_rules(+RulesIn:list_of_rules, -RulesOut:list_of_rules) is det

The predicate sorts RulesIn according to the probability of the rules /

  254sort_rules(R0,R):-
  255  rules2terms(P0,R0),
  256  sort_rules_int(P0,P),
  257  rules2terms(P,R).
  258
  259sort_rules_int(P0,P):-
  260  maplist(to_pair,P0,P1),
  261  sort(1,@>=,P1,P2),
  262  maplist(to_pair,P,P2).
  263
  264to_pair(rule(N,[H:P|R],BL,Lit),P-rule(N,[H:P|R],BL,Lit)).
  265
  266
  267make_dynamic(M):-
  268  M:(dynamic int/1),
  269  findall(O,M:output(O),LO),
  270  findall(I,M:input(I),LI),
  271  findall(I,M:input_cw(I),LIC),
  272  findall(D,M:determination(D,_DD),LDH),
  273  findall(DD,M:determination(_D,DD),LDD),
  274  findall(DH,(M:modeh(_,_,_,LD),member(DH,LD)),LDDH),
  275  append([LO,LI,LIC,LDH,LDD,LDDH],L0),
  276  remove_duplicates(L0,L),
  277  maplist(to_dyn(M),L).
  278
  279to_dyn(M,P/A):-
  280  A1 is A+1,
  281  M:(dynamic P/A1),
  282  A2 is A1+1,
  283  M:(dynamic P/A2).
  284
  285
  286
  287learn_struct(Pos,Neg,Mod,Beam,R,Score):-  %+Beam:initial theory of the form [rule(NR,[h],[b]],...], -R:final theory of the same form, -CLL
  288  format2(Mod,"Clause search~n~n",[]),
  289  Mod:local_setting(max_iter,M),
  290  Mod:local_setting(beamsize,BS),
  291  Mod:local_setting(max_clauses,MC),
  292  cycle_beam(Beam,Mod,Pos,Neg,[],CL,0,M,BS,MC),
  293  maplist(get_cl,CL,LC,MIC,MINC),
  294  maplist(append,MIC,MIC1),
  295  transpose(MIC1,MI),
  296  append(MINC,MIN),
  297  length(LC,NumCL),
  298  write2(Mod,"Final parameter learning"),nl2(Mod),
  299  Mod:local_setting(random_restarts_number,RR),
  300  learn_param_int(MI,MIN,NumCL,Mod,RR,Par,Score),
  301  update_theory(LC,Par,Program1),
  302  remove_zero(Program1,R),
  303  format2(Mod,"Best target theory~n~n",[]),
  304  write_rules2(Mod,R,user_output).
  305
  306get_cl(c(_,C,[MI,MIN]),C,MI,MIN).
  307
  308pick_first(0,_,[]):-!.
  309
  310pick_first(_,[],[]):-!.
  311
  312pick_first(N,[(H,_S)|T],[H|T1]):-
  313  N1 is N-1,
  314  pick_first(N1,T,T1).
  315
  316remove_score([],[]).
  317
  318remove_score([(H,_S)|T],[H|T1]):-
  319  remove_score(T,T1).
  320
  321init_gd_par(0,_Max,[]):-!.
  322
  323init_gd_par(I,Max,[W|TW]):-
  324  I1 is I-1,
  325  W is -Max+random_float*2*Max,
  326  init_gd_par(I1,Max,TW).
  327
  328init_par(_Env,0):-!.
  329
  330init_par(Env,I):-
  331  I1 is I-1,
  332  optimizer_set_x(Env,I1,0.5),
  333  init_par(Env,I1).
  334
  335evaluate_L(Env,M,MIP,MI,L):-
  336  compute_likelihood_pos(MIP,Env,M,0,0,LP),
  337  compute_likelihood_neg(MI,Env,M,LN),
  338  compute_likelihood(LN,M,LP,L).
  339  
  340gen_initial_counts(0,[]):-!.
  341
  342gen_initial_counts(N0,[0|MIP0]):-
  343  N1 is N0-1,
  344  gen_initial_counts(N1,MIP0).
  345
  346evaluate(Env,L,_N,_Step,[M,MIP,MI]):-
  347%  M:mip(MIP),
  348%  M:mi(MI),
  349  compute_likelihood_pos(MIP,Env,M,0,0,LP),
  350  compute_likelihood_neg(MI,Env,M,LN),
  351  compute_likelihood(LN,M,LP,L),
  352  compute_grad(MIP,Env,M,0,MI,LN).
  353
  354compute_grad([],_Env,_M,_N,_MI,_LN):-!.
  355
  356compute_grad([HMIP|TMIP],Env,M,N0,MI,LN):-
  357  compute_sum_neg(MI,M,LN,N0,0,S),
  358  optimizer_get_x(Env,N0,P0),
  359  M:local_setting(zero,Zero),
  360  (P0=<0 ->
  361    PI=Zero
  362  ;
  363    (P0>=1.0->
  364       PI is 1.0-Zero
  365     ;
  366       PI=P0
  367     )
  368  ),
  369
  370 (PI=:= 1.0->
  371    G is 1.0/Zero
  372  ;
  373    G is (HMIP-S)/(1.0-PI)
  374  ),
  375  optimizer_set_g(Env,N0,G),
  376  N1 is N0+1,
  377  compute_grad(TMIP,Env,M,N1,MI,LN).
  378
  379compute_sum_neg([],_M,_LN,_I,S,S).
  380
  381compute_sum_neg([HMI|TMI],M,[HLN|TLN],I,S0,S):-
  382  nth0(I,HMI,MIR),
  383  Den is 1.0-exp(-HLN),
  384  M:local_setting(zero,Zero),
  385  (Den=<0.0->
  386    Den1 is Zero
  387  ;
  388    Den1 = Den
  389  ),
  390  S1 is S0+MIR*exp(-HLN)/Den1,
  391  compute_sum_neg(TMI,M,TLN,I,S1,S).
  392
  393compute_likelihood([],_M,L,L).
  394
  395compute_likelihood([HP|TP],M,L0,L):-
  396  A is 1.0-exp(-HP),
  397  M:local_setting(zero,Zero),
  398  (A=<0.0->
  399    A1 is Zero
  400  ;
  401    A1=A
  402  ),
  403  L1 is L0-log(A1),
  404  compute_likelihood(TP,M,L1,L).
  405
  406compute_likelihood_neg([],_Env,_M,[]).
  407
  408compute_likelihood_neg([HMI|TMI],Env,M,[HLN|TLN]):-
  409  compute_likelihood_pos(HMI,Env,M,0,0,HLN),
  410  compute_likelihood_neg(TMI,Env,M,TLN).
  411
  412compute_likelihood_pos([],_Env,_M,_,LP,LP).
  413
  414compute_likelihood_pos([HMIP|TMIP],Env,M,I,LP0,LP):-
  415  optimizer_get_x(Env,I,P0),
  416  M:local_setting(zero,Zero),
  417  (P0=<0.0 ->
  418    P=Zero
  419  ;
  420    (P0>=1.0->
  421       P is 1-Zero
  422     ;
  423       P=P0
  424     )
  425  ),
  426  LP1 is LP0-log(1-P)*HMIP,
  427  I1 is I+1,
  428  compute_likelihood_pos(TMIP,Env,M,I1,LP1,LP).
  429
  430progress(_Env,FX,X_Norm,G_Norm,Step,_N,Iteration,Ls,0,[M|_]) :-
  431  format4(M,'~d. Iteration :  f(X)=~4f  |X|=~4f
  432                |g(X)|=~4f  Step=~4f  Ls=~4f~n',
  433                [Iteration,FX,X_Norm,G_Norm,Step,Ls]).
  434
  435gen_par(NC,NC,[]):-!.
  436
  437gen_par(N0,NC,[[N0,[0.5,0.5]]|T]):-
  438  N1 is N0+1,
  439  gen_par(N1,NC,T).
  440
  441
  442logistic(X,Sigma_X):-
  443  Sigma_X is 1/(1+exp(-X)).
  444
  445
  446load_python_module(M):-
  447  M:local_setting(parameter_learning,PL),
  448  (PL=em_python;PL=gd_python),!,
  449  absolute_file_name(library('liftcover.pl'), F),
  450  file_directory_name(F,Dir),
  451  py_add_lib_dir(Dir),
  452  processor(M,Proc),
  453  M:local_setting(verbosity,Verb),
  454  py_call(liftcover:init(PL,Proc,Verb)).
  455
  456load_python_module(_).
  457
  458processor(M,Proc):-
  459  M:local_setting(processor,Proc0),
  460  (setting(allow_gpu,true)->
  461    Proc=Proc0
  462  ;
  463    Proc = cpu
  464  ).

 induce_par_lift(:TrainFolds:list_of_atoms, -P:probabilistic_program) is det

The predicate learns the parameters of the program stored in the in/1 fact of the input file using the folds indicated in TrainFolds for training. It returns in P the input program with the updated parameters. /

  473induce_par_lift(Folds,ROut):-
  474  induce_parameters(Folds,R),
  475  rules2terms(R,ROut).
  476
  477induce_parameters(M:Folds,R):-
  478  load_python_module(M),
  479  make_dynamic(M),
  480  M:local_setting(seed,Seed),
  481  setrand(Seed),
  482  findall(Exs,(member(F,Folds),M:fold(F,Exs)),L),
  483  append(L,DB),
  484  statistics(walltime,[_,_]),
  485  (M:bg(RBG0)->
  486    process_clauses(RBG0,M,RBG),
  487    generate_clauses_bg(RBG,ClBG),
  488    assert_all(ClBG,M,ClBGRef)
  489  ;
  490    true
  491  ),
  492  M:in(R00),
  493  process_clauses(R00,M,R0),
  494  statistics(walltime,[_,_]),
  495  find_ex(DB,M,Pos,Neg,_NPos,_NNeg),
  496  M:local_setting(random_restarts_number_str_learn,RR),
  497  number_of_threads(M,Th),
  498  learn_param(R0,M,Pos,Neg,RR,Th,R1,Score,_MI,_MIN),
  499  sort_rules_int(R1,R),
  500  statistics(walltime,[_,CT]),
  501  CTS is CT/1000,
  502  format2(M,'/* Final score ~f~n',[Score]),
  503  format2(M,'Wall time ~f */~n',[CTS]),
  504  write_rules2(M,R,user_output),
  505  (M:bg(RBG0)->
  506    retract_all(ClBGRef)
  507  ;
  508    true
  509  ).
  510
  511number_of_threads(M,Th):-
  512  M:local_setting(threads,Th0),
  513  current_prolog_flag(cpu_count,Cores),
  514  ((Th0=cpu;Th0>Cores)->
  515    Th1 = Cores
  516  ;
  517    Th1 = Th0
  518  ),
  519  setting(max_threads,ThMax),
  520  Th is min(Th1,ThMax).

 filter_rules(:RulesIn:list_of_rules, -RulesOut:list_of_rules) is det

The predicate removes the rules with a probability below or equal to the min_prob parmeter. /

  527filter_rules(M:R0,R):-
  528  M:local_setting(min_probability,Min_prob),
  529  filter_rules(R0,R,Min_prob).

 filter_rules(+RulesIn:list_of_rules, -RulesOut:list_of_rules, +Min_prob:float) is det

The predicate removes from the rules with a probability below or equal to Min_prob. /

  537filter_rules(R0,R,Min_prob):-
  538  (R0=[(_ :- _)|_]->
  539    rules2terms(R0At,R0),
  540    remove_clauses(R0At,Min_prob,RAt,_Num),
  541    rules2terms(RAt,R)
  542  ;
  543    remove_clauses(R0,Min_prob,R,_Num)  
  544  ).

 remove_zero(+RulesIn:list_of_rules, -RulesOut:list_of_rules) is det

The predicate removes the rules with a probability of 0.0. /

  551remove_zero(R0,R1):-
  552  filter_rules(R0,R1,0.0).
  553
  554
  555remove_clauses(Rules,Prob,RulesOut,Num):-
  556  remove_clauses_loop(Rules,Prob,0,Num,[],RulesOut).
  557
  558remove_clauses_loop([],_,Num,Num,Rules,Rules).
  559remove_clauses_loop([Rule|Rest],Prob,NumCur,Num,RulesCur,RulesOut):-
  560  Rule=rule(_N,[_Head:Par|_],_,_),
  561  Par =< Prob,!,
  562  NumCur1 is NumCur+1,
  563  remove_clauses_loop(Rest, Prob, NumCur1,Num,RulesCur, RulesOut).
  564
  565remove_clauses_loop([Rule|Rest],Prob,NumCur,Num,RulesCur,[Rule|RulesOut]):-
  566  remove_clauses_loop(Rest, Prob, NumCur,Num,RulesCur, RulesOut).
  567
  568
  569test_theory_neg_prob(Ex,M,Theory,MIP0,MIP):-
  570  (M:local_setting(single_var,false)->
  571    test_clause_prob(Theory,M,Ex,MIP0,MIP)
  572  ;
  573    test_clause_prob_sv(Theory,M,Ex,MIP0,MIP)
  574  ).
  575
  576test_clause_prob([],_M,_Exs,MIP,MIP).
  577
  578test_clause_prob([(H,B,V,_P)|Rest],M,Exs,[MIPH0|MIPT0],[MIPH|MIPT]):-
  579  maplist(test_ex(V,H,B,M),Exs,L),
  580  sum_list(L,MIP),
  581  MIPH is MIPH0+MIP,
  582  test_clause_prob(Rest,M,Exs,MIPT0,MIPT).
  583
  584test_ex(_V,H,B,M,E,N):-
  585  findall(1,(H=E,M:B),L),
  586  length(L,N).
  587
  588test_clause_prob_sv([],_M,_Exs,MIP,MIP).
  589
  590test_clause_prob_sv([(H,B,V,_P)|Rest],M,Exs,[MIPH0|MIPT0],[MIPH|MIPT]):-
  591  maplist(test_ex_sv(V,H,B,M),Exs,L),
  592  sum_list(L,MIP),
  593  MIPH is MIPH0+MIP,
  594  test_clause_prob_sv(Rest,M,Exs,MIPT0,MIPT).
  595
  596
  597test_ex_sv(_V,H,B,M,E,N):-
  598  (\+ (H=E,M:B)->
  599    N=0
  600  ;
  601    N=1
  602  ).  
  603
  604test_theory_pos_prob(Ex,M,Th,N,LMI):-
  605  (M:local_setting(single_var,false)->
  606    test_theory_pos_prob_mv(Ex,M,Th,N,LMI)
  607  ;
  608    test_theory_pos_prob_sv(Ex,M,Th,N,LMI)
  609  ).
  610
  611test_theory_pos_prob_mv([],_M,_Theory,_N,[]).
  612
  613test_theory_pos_prob_mv([Ex|Rest],M,Th,N,[MI|LMI]):-
  614  gen_initial_counts(N,MI0),
  615  test_clause_prob(Th,M,[Ex],MI0,MI),
  616  test_theory_pos_prob_mv(Rest,M,Th,N,LMI).
  617
  618test_theory_pos_prob_sv([],_M,_Theory,_N,[]).
  619
  620test_theory_pos_prob_sv([Ex|Rest],M,Th,N,[MI|LMI]):-
  621  gen_initial_counts(N,MI0),
  622  test_clause_prob_sv(Th,M,[Ex],MI0,MI),
  623  test_theory_pos_prob_sv(Rest,M,Th,N,LMI).
  624
  625learn_param([],M,_,_,_,_,[],MInf,[],[]):-!,
  626  M:local_setting(minus_infinity,MInf).
  627
  628learn_param(Program0,M,Pos,Neg,RR,Th,Program,LL,MI,MIN):-
  629  generate_clauses(Program0,M,0,Pr1),
  630  length(Program0,N),
  631  format4(M,'Computing clause statistics~n',[]),
  632  gen_initial_counts(N,MIN0),
  633  clauses_statistics(Pr1,N,M,Pos,Neg,MIN0,MI,MIN,Th),
  634  format4(M,'Updating parameters~n',[]),
  635  learn_param_int(MI,MIN,N,M,RR,Par,LL),
  636  update_theory(Program0,Par,Program1),
  637  remove_zero(Program1,Program).
  638
  639
  640clauses_statistics(Pr,N,M,Pos,Neg,MIN0,MI,MIN,Th):-
  641  (Th=1->
  642    test_theory_neg_prob(Neg,M,Pr,MIN0,MIN),
  643    test_theory_pos_prob(Pos,M,Pr,N,MI)
  644  ;
  645    current_prolog_flag(cpu_count,Cores),
  646    ((Th=cpu;Th>Cores)->
  647      Chunks = Cores
  648    ;
  649      Chunks = Th
  650    ),
  651    chunks(Pos,Chunks,PosC),
  652    chunks(Neg,Chunks,NegC),
  653    concurrent_maplist(test_theory_neg_prob_conc(Pr,M,MIN0),NegC,MINC),
  654    concurrent_maplist(test_theory_pos_prob_conc(Pr,M,N),PosC,MIC),
  655    append(MIC,MI),
  656    transpose(MINC,MINT),
  657    maplist(sum_list,MINT,MIN)
  658  ).
  659
  660test_theory_neg_prob_conc(Pr,M,MIN0,Neg,MIN):-
  661  test_theory_neg_prob(Neg,M,Pr,MIN0,MIN).
  662
  663test_theory_pos_prob_conc(Pr,M,N,Pos,MI):-
  664  test_theory_pos_prob(Pos,M,Pr,N,MI).
  665
  666chunks(L,N,Chunks):-
  667  length(L,Len),
  668  LenChunks is round(Len/N),
  669  split_list(L,N,LenChunks,Chunks).
  670
  671split_list(L,1,_,[L]):-!.
  672
  673split_list(L0,N,NL,[H|L]):-
  674  N>1,
  675  N1 is N-1,
  676  length(H1,NL),
  677  (append(H1,T,L0)->
  678    H=H1,
  679    split_list(T,N1,NL,L)
  680  ;
  681    H=L0,
  682    L=[]
  683  ).
  684
  685learn_param_int(MI,MIN,N,M,NR,Par,LL):-
  686  M:local_setting(parameter_learning,em),!,
  687  random_restarts(0,NR,M,-1e20,LL,N,initial,Par,MI,MIN),
  688  format3(M,"Final LL ~f~n",[LL]).
  689
  690learn_param_int(MI,MIN,_N,M,NR,Par,LL):-
  691  M:local_setting(parameter_learning,em_python),!,
  692  M:local_setting(eps,EA),
  693  M:local_setting(eps_f,ER),
  694  M:local_setting(iter,Iter),
  695  M:local_setting(regularization,Reg),
  696  M:local_setting(gamma,Gamma),
  697  M:local_setting(zero,Zero),
  698  M:local_setting(ab,[A,B]),
  699  M:local_setting(verbosity,Verb),
  700  processor(M,Device),
  701  py_call(liftcover:random_restarts(MI,MIN,Device,NR,Iter,EA,ER,Reg,Zero,Gamma,A,B,Verb),-(Par,LL)),
  702  format3(M,"Final LL ~f~n",[LL]).
  703
  704learn_param_int(MI,MIN,N,M,NR,Par,LL):-
  705  M:local_setting(parameter_learning,gd),!,
  706  random_restarts_gd(0,NR,M,-1e20,PLL,N,initial,ParR,MI,MIN),  %computes new parameters Par
  707  maplist(logistic,ParR,Par),
  708  LL is -PLL,
  709  format3(M,"Final LL ~f~n",[LL]).
  710
  711
  712learn_param_int(MI,MIN,_N,M,NR,Par,LL):-
  713  M:local_setting(parameter_learning,gd_python),!,
  714  M:local_setting(verbosity,Verb),
  715  M:local_setting(parameter_update,UpdateMethod),
  716  M:local_setting(iter,Iter),
  717  M:local_setting(eps,Eps),
  718  M:local_setting(adam_params,[Eta,Beta1,Beta2,Epsilon]),
  719  M:local_setting(eta,LearningRate),
  720  M:local_setting(gamma,Gamma),
  721  M:local_setting(regularization,Reg),
  722  M:local_setting(zero,Zero),
  723  processor(M,Device),
  724  py_call(liftcover:random_restarts_gd(MI,MIN,Device,NR,UpdateMethod,
  725    Iter,Eps,Reg,Gamma,LearningRate,Eta,-(Beta1,Beta2),Epsilon,Zero,Verb),-(Par,LL)),
  726  format3(M,"Final LL ~f~n",[LL]).
  727
  728
  729learn_param_int(MI,MIN,N,M,_,Par,LL):-
  730  M:local_setting(parameter_learning,lbfgs),
  731  optimizer_initialize(N,liftcover,evaluate,progress,[M,MIN,MI],Env),
  732  init_par(Env,N),
  733  evaluate_L(Env,M,MIN,MI,L),
  734% parte da modificare fine
  735  IL is -L,
  736  format4(M,"~nInitial L ~f~n",[IL]),
  737  optimizer_run(Env,_LL,Status),
  738  format4(M,"Status ~p~n",[Status]),
  739  new_pars_lbfgs(Env,M,0,N,Par),
  740  evaluate_L(Env,M,MIN,MI,NewL),
  741  LL is -NewL,
  742  optimizer_finalize(Env).
  743
  744
  745new_pars_lbfgs(_Env,_M,N,N,[]):-!.
  746
  747new_pars_lbfgs(Env,M,N,NMax,[P|Rest1]):-
  748    optimizer_get_x(Env,N,P0),
  749    M:local_setting(zero,Zero),
  750    (P0=<0.0->
  751      P=Zero
  752    ;
  753      (P0>=1.0->
  754        P is 1.0-Zero
  755      ;
  756        P=P0
  757      )
  758    ),
  759    N1 is N+1,
  760    new_pars_lbfgs(Env,M,N1,NMax,Rest1).
  761
  762
  763random_restarts(N,N,_M,Score,Score,_N,Par,Par,_MI,_MIN):-!.
  764
  765random_restarts(N,NMax,M,Score0,Score,NR,Par0,Par,MI,MIN):-
  766  N1 is N+1,
  767  format3(M,"Restart number ~d~n~n",[N1]),
  768  length(Par1,NR),
  769  maplist(random,Par1),
  770  M:local_setting(eps,EA),
  771  M:local_setting(eps_f,ER),
  772  M:local_setting(iter,Iter),
  773  em(EA,ER,0,Iter,M,NR,Par1,-1e20,MI,MIN,ParR,ScoreR),
  774  format3(M,"Random_restart: Score ~f~n",[ScoreR]),
  775  (ScoreR>Score0->
  776    random_restarts(N1,NMax,M,ScoreR,Score,NR,ParR,Par,MI,MIN)
  777  ;
  778    random_restarts(N1,NMax,M,Score0,Score,NR,Par0,Par,MI,MIN)
  779  ).
  780
  781random_restarts_gd(N,N,_M,Score,Score,_N,Par,Par,_MI,_MIN):-!.
  782
  783random_restarts_gd(N,NMax,M,_Score0,Score,NR,_Par0,Par,MI,MIN):-
  784  N1 is N+1,
  785  format3(M,"Restart number ~d~n~n",[N1]),
  786  M:local_setting(max_initial_weight,Max),
  787  init_gd_par(NR,Max,Par1),
  788  evaluate_L_gd(M,MIN,MI,Par1,L),
  789  ScoreIn is -L,
  790  format3(M,"GD Random_restart: initial score ~f~n",[ScoreIn]),
  791  M:local_setting(eps,EA),
  792  M:local_setting(eps_f,ER),
  793  M:local_setting(iter,Iter),
  794  (M:local_setting(parameter_update,adam)->
  795    findall(0,between(1,NR,_),M0),
  796    findall(0,between(1,NR,_),M1),
  797    gd_adam(EA,ER,0,Iter,M,NR,Par1,M0,M1,-1e20,MI,MIN,ParR,ScoreR)
  798  ;
  799    gd(EA,ER,0,Iter,M,NR,Par1,-1e20,MI,MIN,ParR,ScoreR)
  800  ),
  801  format3(M,"GD Random_restart: Score ~f~n",[ScoreR]),
  802  random_restarts_gd(N1,NMax,M,ScoreR,Score,NR,ParR,Par,MI,MIN).
  803
  804
  805gd(_EA,_ER,MaxIter,MaxIter,_M,_NR,Par,Score,_MI,_MIP,Par,Score):-!.
  806
  807gd(EA,ER,Iter0,MaxIter,M,NR,Par0,Score0,MI,MIP,Par,Score):-
  808  compute_gradient_gd(MIP,MI,M,Par0,G,LL),
  809  Score1 is -LL,
  810  Iter is Iter0+1,
  811  format4(M,"GD Iteration ~d LL ~f~n",[Iter,Score1]),
  812  Diff is Score1-Score0,
  813  Fract is -Score1*ER,
  814  (( Diff<EA;Diff<Fract)->
  815    Score=Score1,
  816    Par=Par0
  817  ;
  818    M:local_setting(eta,Eta),
  819    M:local_setting(gamma,Gamma),
  820    (M:local_setting(regularization,l2)->
  821      maplist(l2,Par0,Reg)
  822    ;
  823      (M:local_setting(regularization,l1)->
  824        maplist(l1,Par0,Reg)
  825      ;
  826        findall(0,between(1,NR,_),Reg)
  827      )
  828    ),
  829    maplist(update_par(Eta,Gamma),Par0,Reg,G,Par1),
  830    gd(EA,ER,Iter,MaxIter,M,NR,Par1,Score1,MI,MIP,Par,Score)
  831  ).
  832
  833gd_adam(_EA,_ER,Iter,Iter,_M,_NR,Par,_M0,_M1,Score,_MI,_MIP,Par,Score):-!.
  834
  835gd_adam(EA,ER,Iter0,MaxIter,M,NR,Par0,M00,M10,Score0,MI,MIP,Par,Score):-
  836  compute_gradient_gd(MIP,MI,M,Par0,G,LL),
  837  Score1 is -LL,
  838  Iter is Iter0+1,
  839  format4(M,"Iteration ~d LL ~f~n",[Iter,Score1]),
  840  Diff is Score1-Score0,
  841  Fract is -Score1*ER,
  842  (( Diff<EA;Diff<Fract)->
  843    Score=Score1,
  844    Par=Par0
  845  ;
  846    M:local_setting(gamma,Gamma),
  847    M:local_setting(adam_params,[Eta,Beta1,Beta2,Epsilon]),
  848    (M:local_setting(regularization,l2)->
  849      maplist(l2,Par0,Reg)
  850    ;
  851      (M:local_setting(regularization,l1)->
  852        maplist(l1,Par0,Reg)
  853      ;
  854        findall(0,between(1,NR,_),Reg)
  855      )
  856    ),
  857    maplist(update_grad(Gamma),G,Reg,G1),
  858    maplist(update_M0(Beta1),M00,G1,M0),
  859    maplist(update_M1(Beta2),M10,G1,M1),
  860    EtaIter is Eta*sqrt(1-Beta2^Iter)/(1-Beta1^Iter),
  861    maplist(update_par_adam(EtaIter,Epsilon),Par0,M0,M1,Par1),
  862    gd_adam(EA,ER,Iter,MaxIter,M,NR,Par1,M0,M1,Score1,MI,MIP,Par,Score)
  863  ).
  864
  865update_par_adam(EtaIter,Epsilon,Par0,M0,M1,Par1):-
  866  Par1 is Par0-EtaIter*M0/(sqrt(M1)+Epsilon).
  867
  868
  869update_M0(Beta1,M00,G,M0):-
  870  M0 is Beta1*M00+(1-Beta1)*G.
  871
  872update_M1(Beta2,M10,G,M1):-
  873  M1 is Beta2*M10+(1-Beta2)*G^2.
  874
  875update_grad(Gamma,G,Reg,G1):-
  876  G1 is G+Gamma*Reg.
  877
  878
  879l1(W,R):-
  880  logistic(W,S),
  881  R is S*(1-S).
  882
  883l2(W,R):-
  884  logistic(W,S),
  885  R is 2*S^2*(1-S).
  886
  887update_par(Eta,Gamma,Par0,Reg,G,Par1):-
  888  Par1 is Par0-Eta*(G+Gamma*Reg).
  889
  890
  891
  892evaluate_L_gd(M,MIP,MI,Par,L):-
  893  maplist(logistic,Par,Prob),
  894  compute_likelihood_pos_gd(MIP,Prob,M,0,LP),
  895%  write(lpos),nl,
  896  compute_likelihood_neg_gd(MI,Prob,M,LN),
  897%  write(lneg),nl,
  898  compute_likelihood_gd(LN,M,LP,L).
  899
  900
  901
  902compute_gradient_gd(MIP,MI,M,Par,G,L):-
  903  maplist(logistic,Par,Prob),
  904  compute_likelihood_pos_gd(MIP,Prob,M,0,LP),
  905  compute_likelihood_neg_gd(MI,Prob,M,LN),
  906  compute_likelihood_gd(LN,M,LP,L),
  907  compute_grad_gd(MIP,Prob,M,0,MI,LN,G).
  908%  write(grad),nl.
  909
  910compute_likelihood_neg_gd([],_Prob,_M,[]).
  911
  912compute_likelihood_neg_gd([HMI|TMI],Prob,M,[HLN|TLN]):-
  913  compute_likelihood_pos_gd(HMI,Prob,M,0,HLN),
  914  compute_likelihood_neg_gd(TMI,Prob,M,TLN).
  915
  916compute_likelihood_pos_gd([],[],_M,LP,LP).
  917
  918compute_likelihood_pos_gd([HMIP|TMIP],[P|TP],M,LP0,LP):-
  919  LP1 is LP0-log(1-P)*HMIP,
  920  compute_likelihood_pos_gd(TMIP,TP,M,LP1,LP).
  921
  922
  923compute_likelihood_gd([],_M,L,L).
  924
  925compute_likelihood_gd([HP|TP],M,L0,L):-
  926  A is 1.0-exp(-HP),
  927  (A=:=0.0->
  928    M:local_setting(logzero,LZ),
  929    L1 is L0-LZ
  930  ;
  931    L1 is L0-log(A)
  932  ),
  933  compute_likelihood_gd(TP,M,L1,L).
  934
  935
  936compute_grad_gd([],[],_M,_N,_MI,_LN,[]):-!.
  937
  938compute_grad_gd([HMIP|TMIP],[P|TP],M,N0,MI,LN,[G|TG]):-
  939%  write(prima_comp_grad),nl,
  940  compute_sum_neg(MI,M,LN,N0,0,S),
  941%  write(opt),nl,
  942  G is (HMIP-S)*P,
  943  N1 is N0+1,
  944  compute_grad_gd(TMIP,TP,M,N1,MI,LN,TG).
  945
  946compute_sum_neg_gd([],_M,_LN,_I,S,S).
  947
  948compute_sum_neg_gd([HMI|TMI],M,[HLN|TLN],I,S0,S):-
  949%  write(HMI),write(hmi),nl,
  950%  write(I),write('I'),nl,
  951  nth0(I,HMI,MIR),
  952%  write(MIR),write(mir),nl,
  953%  write(HLN),write(hln),nl,
  954  Den is 1.0-exp(-HLN),
  955  S1 is S0+MIR*exp(-HLN)/Den,
  956  compute_sum_neg_gd(TMI,M,TLN,I,S1,S).
  957
  958
  959em(_EA,_ER,MaxIter,MaxIter,_M,_NR,Par,Score,_MI,_MIN,Par,Score):-!.
  960
  961em(EA,ER,Iter0,MaxIter,M,NR,Par0,Score0,MI,MIN,Par,Score):-
  962  length(Eta0,NR),
  963  expectation_quick(Par0,M,MI,MIN,Eta0,Eta,Score1),
  964  Iter is Iter0+1,
  965  format4(M,"Iteration ~d LL ~f~n",[Iter,Score1]),
  966  maximization_quick(Eta,M,Par1),
  967  Diff is Score1-Score0,
  968  Fract is -Score1*ER,
  969  (( Diff<EA;Diff<Fract)->
  970    Score=Score1,
  971    Par=Par1
  972  ;
  973    em(EA,ER,Iter,MaxIter,M,NR,Par1,Score1,MI,MIN,Par,Score)
  974  ).
  975
  976expectation_quick(Par,M,MI,MIN,Eta0,Eta,Score):-
  977 /* LLO is the negative examples contribution in the LL*/
  978  M:local_setting(logzero,LogZero),
  979  foldl(llm(LogZero),Par,MIN,0,LL0),
  980  maplist(eta0,MIN,Eta0),
  981  /* positive examples contibution in LL*/
  982  scan_pos(MI,M,Par,LL0,Eta0,Score,Eta).
  983
  984maximization_quick(Eta,M,Par):-
  985  (M:local_setting(regularization,l1)->
  986    maplist(maximize_L1(M),Eta,Par)
  987  ;
  988    (M:local_setting(regularization,l2)->
  989      maplist(maximize_L2(M),Eta,Par)
  990  ;
  991      (M:local_setting(regularization,bayesian)->
  992        maplist(maximize_bayesian(M),Eta,Par)
  993      ;
  994        maplist(maximize(M),Eta,Par)
  995  )
  996    )
  997  ).
  998
  999maximize(_M,[Eta0,Eta1],Par):-
 1000  (Eta0+Eta1=:=0.0->
 1001    Par=0.0
 1002  ;
 1003    Par is Eta1/(Eta0+Eta1)  
 1004  ).
 1005
 1006
 1007maximize_L1(M,[Eta0,Eta1],Par):-
 1008  M:local_setting(gamma,Gamma),
 1009/*(((Eta0+Eta1)^2+Gamma^2+2*Gamma*(Eta0-Eta1))<0 ->
 1010writeln((Eta0+Eta1)^2+Gamma^2+2*Gamma*(Eta0-Eta1)),
 1011  Par is 4*Eta1/(2*(Gamma+Eta0+Eta1))
 1012;((2*(Gamma+Eta0+Eta1+sqrt((Eta0+Eta1)^2+Gamma^2+2*Gamma*(Eta0-Eta1)))=:=0.0)->
 1013writeln((2*(Gamma+Eta0+Eta1+sqrt((Eta0+Eta1)^2+Gamma^2+2*Gamma*(Eta0-Eta1)))))
 1014
 1015;*/
 1016  Par is 4*Eta1/(2*(Gamma+Eta0+Eta1+sqrt((Eta0+Eta1)^2+Gamma^2+2*Gamma*(Eta0-Eta1))))
 1017%)
 1018  %)
 1019  .
 1020
 1021maximize_L2(M,[Eta0,Eta1],Par):-
 1022  M:local_setting(gamma,Gamma),
 1023  Sum is 3*Eta0+3*Eta1+Gamma,
 1024  Arccos is acos(sqrt(Gamma/Sum)*(9*Eta0/2-9*Eta1+Gamma)/(3*Eta0+3*Eta1+Gamma)),
 1025  Par is 2*sqrt(Sum/Gamma)*cos(Arccos/3-2*pi/3)/3+1/3.
 1026
 1027maximize_bayesian(M,[Eta0,Eta1],Par):-
 1028  M:local_setting(ab,[A,B]),
 1029  Par is (Eta1+A)/(Eta0+Eta1+A+B).
 1030
 1031/*....scan_pos predicate..... */
 1032
 1033scan_pos([],_M,_Par,LL,Eta,LL,Eta).
 1034
 1035scan_pos([MIH|MIT],M,Par,LL0,Eta0,LL,Eta):-
 1036  M:local_setting(logzero,LogZero),
 1037  foldl(rule_contrib,MIH,Par,1,Prod),
 1038  ProbEx is 1-Prod,
 1039  (ProbEx=:=0.0->
 1040    LLCurrent is LL0+LogZero
 1041   ;
 1042    LLCurrent is LL0+log(ProbEx)
 1043  ),
 1044  maplist(update_eta(ProbEx,M),Eta0,Par,MIH,EtaCurrent),
 1045  scan_pos(MIT,M,Par,LLCurrent,EtaCurrent,LL,Eta).
 1046
 1047
 1048update_eta(ProbEx,M,[Etai00,Etai10],Pi,MIR,[Etai0,Etai1]):-
 1049  M:local_setting(zero,Zero),
 1050  ( ProbEx=:=0.0->
 1051    CondP is Pi/Zero
 1052  ;
 1053    CondP is Pi/ProbEx
 1054  ),
 1055  OCondP0 is 1-CondP,
 1056  ( OCondP0 <0.0->
 1057    OCondP = 0.0
 1058  ;
 1059    OCondP = OCondP0
 1060  ),
 1061  Etai0 is Etai00+MIR*OCondP,
 1062  Etai1 is Etai10+MIR*CondP.
 1063
 1064rule_contrib(MIR,Pi,P0,P):-
 1065  P is P0*(1-Pi)^MIR.
 1066
 1067llm(LogZero,Pi,MI,LL0,LL):-
 1068  ((1-Pi)=:=0.0 ->
 1069
 1070    ( MI==0 ->
 1071      LL is LL0
 1072     ;
 1073       LL is LL0+MI*LogZero
 1074    )
 1075   ;
 1076    ( MI==0 ->
 1077       LL is LL0
 1078     ;
 1079       LL is LL0+MI*log(1-Pi)
 1080    )
 1081
 1082  ).
 1083
 1084eta0(MIN,[MIN,0]).
 1085
 1086
 1087update_theory([],_N,[]):-!.
 1088
 1089update_theory([rule(Name,[H:_,_],B,L)|Rest],[P|T],[rule(Name,[H:P,'':PN],B,L)|Rest1]):-
 1090    PN is 1-P,
 1091    update_theory(Rest,T,Rest1).
 1092
 1093
 1094
 1095cycle_beam([],_Mod,_Pos,_Neg,CL,CL,_M0,_M,_BS,_MC):-!.
 1096
 1097cycle_beam(_Beam,_Mod,_Pos,_Neg,CL,CL,M,M,_BS,_MC):-!.
 1098
 1099cycle_beam(Beam,Mod,Pos,Neg,CL0,CL,M0,M,BS,MC):-
 1100  M1 is M0+1,%decreases the number of max_iter M
 1101  format2(Mod,"Clause iteration ~d~n~n",[M1]),
 1102  /*write('\n\ncurrent beam\n\n'),
 1103  write(Beam),
 1104  write('\n\n'),*/
 1105  cycle_clauses(Beam,Mod,Pos,Neg,[],NB0,CL0,CL1),
 1106  (length(NB,BS),append(NB,_,NB0)->
 1107    true
 1108  ;
 1109   NB=NB0 
 1110  ),
 1111  (MC=:= inf ->
 1112    CL2=CL1
 1113  ;
 1114    (length(CL2,MC),append(CL2,_,CL1)->
 1115      true
 1116    ;
 1117      CL2=CL1
 1118    )
 1119  ),
 1120  cycle_beam(NB,Mod,Pos,Neg,CL2,CL,M1,M,BS,MC).
 1121
 1122cycle_clauses([],_M,_Pos,_Neg,NB,NB,CL,CL):-!.
 1123
 1124cycle_clauses([c(_ScoreH,RH,_)|T],M,Pos,Neg,NB0,NB,CL0,CL):-
 1125%  write3('\n\nRevising clause\n'),
 1126%  write_rules3([RH],user_output),
 1127%  RH=rule(_,H,B,Lits),
 1128%  write3(H),write3(' '),write3(B),
 1129%  write3('\n'),write3(Lits),write3('\n'),
 1130  findall(RS,specialize_rule(RH,M,RS,_L),LR),!,   %-LR:list of lists, each one correponding to a different revised theory; specialize_rule defined in revise.pl
 1131  length(LR,NR),
 1132  write3(M,'Number of revisions '),write3(M,NR),write3(M,'\n'),
 1133  score_clause_refinements(LR,M,Pos,Neg,NB0,NB1,CL0,CL1),
 1134  cycle_clauses(T,M,Pos,Neg,NB1,NB,CL1,CL).
 1135
 1136
 1137score_clause_refinements(LR,M,Pos,Neg,NB0,NB,CL0,CL):-  %scans the list of revised theories
 1138  number_of_threads(M,Th),
 1139  (Th=1->
 1140     score_clause_refinements_int(M,1,Pos,Neg,LR,NB1,CL1),
 1141     list_to_ord_set(NB1,NB1OS),
 1142     ord_union(NB1OS,NB0,NB),
 1143     list_to_ord_set(CL1,CL1OS),
 1144     ord_union(CL1OS,CL0,CL)
 1145  ;
 1146    chunks(LR,Th,LRC),
 1147    concurrent_maplist(score_clause_refinements_int(M,1,Pos,Neg),LRC,NBs,CLs),
 1148    merge_ordsets(NBs,NB0,NB),
 1149    merge_ordsets(CLs,CL0,CL)
 1150  ).
 1151
 1152merge_ordsets([],OS,OS).
 1153
 1154merge_ordsets([H|T],OS0,OS):-
 1155  list_to_ord_set(H,HOS),
 1156  ord_union(HOS,OS0,OS1),
 1157  merge_ordsets(T,OS1,OS).
 1158
 1159
 1160score_clause_refinements_int(_M,_N,_Pos,_Neg,[],[],[]):-!.
 1161
 1162score_clause_refinements_int(M,Nrev,Pos,Neg,[R1|T],[c(Score,R3,_)|NB],CL):-  %scans the list of revised theories
 1163  already_scored_clause(R1,R3,M,Score),!,
 1164  format3(M,'Score ref.  ~d~n',[Nrev]),
 1165  write3(M,'Already scored, updated refinement\n'),
 1166  write_rules3(M,[R3],user_output),
 1167  write3(M,'Score '),write3(M,Score),write3(M,'\n\n\n'),
 1168  Nrev1 is Nrev+1,
 1169  score_clause_refinements_int(M,Nrev1,Pos,Neg,T,NB,CL).
 1170
 1171score_clause_refinements_int(M,Nrev,Pos,Neg,[R1|T],NB,CL):-
 1172  format3(M,'Score ref.  ~d~n',[Nrev]),
 1173  write_rules3(M,[R1],user_output),
 1174  M:local_setting(random_restarts_number_str_learn,NR),
 1175  learn_param([R1],M,Pos,Neg,NR,1,NewR,Score,MI,MIN),
 1176  write3(M,'Updated refinement\n'),
 1177  write_rules3(M,NewR,user_output),
 1178  write3(M,'Score (CLL) '),write3(M,Score),write3(M,'\n\n\n'),
 1179  (NewR=[R3]->
 1180    NB=[c(Score,R3,_)|NB0],
 1181    (range_restricted(R3)->
 1182      CL=[c(Score,R3,[MI,MIN])|CL0]
 1183    ;
 1184      CL=CL0
 1185    ),
 1186    format2(M,"Added a target clause~n",[]),
 1187    store_clause_refinement(R1,R3,M,Score),
 1188    Nrev1 is Nrev+1
 1189  ;
 1190    NB=NB0,
 1191    CL=CL0,
 1192    Nrev1=Nrev+1
 1193  ),
 1194  score_clause_refinements_int(M,Nrev1,Pos,Neg,T,NB0,CL0).
 1195
 1196range_restricted(rule(_N,HL,BL,_Lit)):-
 1197  term_variables(HL,VH),
 1198  term_variables(BL,VB),
 1199  sublisteq(VH,VB).
 1200
 1201sublisteq([],_).
 1202
 1203sublisteq([H|T],L):-
 1204  member_eq(H,L),
 1205  sublisteq(T,L).
 1206
 1207target(R,M):-
 1208  get_output_preds(R,O),
 1209  member(T,O),
 1210  M:output(T),!.
 1211
 1212get_output_preds(rule(_N,HL,_BL,_Lit),O):-
 1213  scan_head(HL,[],O).
 1214
 1215scan_head(['':_],O,O):-!.
 1216scan_head([],O,O):-!.
 1217scan_head([H:_P|T],O0,O):-
 1218  functor(H,F,N),
 1219  (member(F/N,O0)->
 1220    O1=O0
 1221  ;
 1222    O1=[F/N|O0]
 1223  ),
 1224  scan_head(T,O1,O).
 1225
 1226
 1227store_clause_refinement(Ref,RefP,M,Score):-
 1228  elab_clause_ref(Ref,Ref1),
 1229  assert(M:ref_clause(r(Ref1,RefP,Score))).
 1230
 1231store_refinement(Ref,RefP,M,Score):-
 1232  elab_ref(Ref,Ref1),
 1233  assert(M:ref(r(Ref1,RefP,Score))).
 1234
 1235already_scored_clause(R,R1,M,Score):-
 1236  elab_ref([R],[rule(H,B)]),
 1237  M:ref_clause(r(rule(H,B1),R1,Score)),
 1238  permutation(B,B1).
 1239
 1240already_scored(R,R1,M,Score):-
 1241  elab_ref(R,RR),
 1242  M:ref(r(RR,R1,Score)).
 1243
 1244
 1245elab_clause_ref(rule(_NR,H,B,_Lits),rule(H1,B1)):-
 1246  copy_term((H,B),(H1,B1)).
 1247
 1248elab_ref([],[]).
 1249
 1250elab_ref([rule(_NR,H,B,_Lits)|T],[rule(H1,B1)|T1]):-!,
 1251  copy_term((H,B),(H1,B1)),
 1252  numbervars((H1,B1),0,_N),
 1253  elab_ref(T,T1).
 1254
 1255elab_ref([def_rule(H,B,_Lits)|T],[rule(H1,B1)|T1]):-
 1256  copy_term((H,B),(H1,B1)),
 1257  numbervars((H1,B1),0,_N),
 1258  elab_ref(T,T1).
 1259
 1260%insertion in the beam
 1261insert_in_order([],C,BeamSize,[C]):-
 1262  BeamSize>0,!.
 1263
 1264insert_in_order(Beam,_New,0,Beam):-!.
 1265
 1266insert_in_order([[Th1,Heuristic1|Rest1]|RestBeamIn],[Th,Heuristic|Rest],BeamSize,BeamOut):-
 1267  Heuristic>Heuristic1,!,
 1268  % larger heuristic, insert here
 1269  NewBeam=[[Th,Heuristic|Rest],[Th1,Heuristic1|Rest1]|RestBeamIn],
 1270  length(NewBeam,L),
 1271  (L>BeamSize->
 1272    nth1(L,NewBeam,_Last,BeamOut)
 1273  ;
 1274    BeamOut=NewBeam
 1275  ).
 1276
 1277insert_in_order([[Th1,Heuristic1|Rest1]|RestBeamIn],[Th,Heuristic|Rest],BeamSize,
 1278[[Th1,Heuristic1|Rest1]|RestBeamOut]):-
 1279  BeamSize1 is BeamSize -1,
 1280  insert_in_order(RestBeamIn,[Th,Heuristic|Rest],BeamSize1,
 1281  RestBeamOut).
 1282
 1283
 1284
 1285remove_int_atom_list([],[]).
 1286
 1287remove_int_atom_list([\+ A|T],[\+ A1|T1]):-!,
 1288  A=..[F,_|Arg],
 1289  A1=..[F|Arg],
 1290  remove_int_atom_list(T,T1).
 1291
 1292remove_int_atom_list([A|T],[A1|T1]):-
 1293  A=..[F,_|Arg],
 1294  A1=..[F|Arg],
 1295  remove_int_atom_list(T,T1).
 1296
 1297
 1298
 1299remove_int_atom(\+ A,\+ A1):-!,
 1300  A=..[F,_|T],
 1301  A1=..[F|T].
 1302
 1303remove_int_atom(A,A1):-
 1304  A=..[F,_|T],
 1305  A1=..[F|T].
 1306
 1307
 1308get_heads([],[]).
 1309
 1310get_heads([_-H|T],[H|TN]):-
 1311  get_heads(T,TN).
 1312
 1313randomize([],[]):-!.
 1314
 1315randomize([rule(N,V,NH,HL,BL,LogF)|T],[rule(N,V,NH,HL1,BL,LogF)|T1]):-
 1316  length(HL,L),
 1317  Int is 1.0/L,
 1318  randomize_head(Int,HL,0,HL1),
 1319  randomize(T,T1).
 1320
 1321randomize_head(_Int,['':_],P,['':PNull1]):-!,
 1322  PNull is 1.0-P,
 1323  (PNull>=0.0->
 1324    PNull1 =PNull
 1325  ;
 1326    PNull1=0.0
 1327  ).
 1328
 1329randomize_head(Int,[H:_|T],P,[H:PH1|NT]):-
 1330  PMax is 1.0-P,
 1331  random(0,PMax,PH1),
 1332  P1 is P+PH1,
 1333  randomize_head(Int,T,P1,NT).
 1334
 1335
 1336
 1337update_head([],[],_N,[]):-!.
 1338
 1339update_head([H:_P|T],[PU|TP],N,[H:P|T1]):-
 1340  P is PU/N,
 1341  update_head(T,TP,N,T1).
 1342
 1343/* utilities */
 1344
 1345
 1346rules2terms(R,T):-
 1347  maplist(rule2term,R,T).
 1348
 1349rule2term(rule(_N,HL,BL,_Lit),(H:-B)):-!,
 1350  list2or(HL,H),
 1351  list2and(BL,B).
 1352
 1353rule2term(def_rule(H,BL,_Lit),((H:1.0):-B)):-
 1354  list2and(BL,B).
 1355
 1356
 1357write_rules([],_S).
 1358
 1359write_rules([rule(_N,HL,BL,Lit)|T],S):-!,
 1360  copy_term((HL,BL,Lit),(HL1,BL1,Lit1)),
 1361  numbervars((HL1,BL1,Lit1),0,_M),
 1362  write_disj_clause(S,(HL1:-BL1)),
 1363  write_rules(T,S).
 1364
 1365write_rules([def_rule(H,BL,Lit)|T],S):-
 1366  copy_term((H,BL,Lit),(H1,BL1,Lit1)),
 1367  numbervars((H1,BL1,Lit1),0,_M),
 1368  write_disj_clause(S,([H1:1.0]:-BL1)),
 1369  write_rules(T,S).
 1370
 1371
 1372new_par([],[],[]).
 1373
 1374new_par([HP|TP],[Head:_|TO],[Head:HP|TN]):-
 1375  new_par(TP,TO,TN).
 1376
 1377
 1378
 1379write_disj_clause(S,(H:-[])):-!,
 1380  write_head(S,H),
 1381  format(S,".~n~n",[]).
 1382
 1383write_disj_clause(S,(H:-B)):-
 1384  write_head(S,H),
 1385  format(S,' :-',[]),
 1386  nl(S),
 1387  write_body(S,B).
 1388
 1389
 1390write_head(S,[A:1.0|_Rest]):-!,
 1391  format(S,"~q:1.0",[A]).
 1392
 1393write_head(S,[A:P,'':_P]):-!,
 1394  format(S,"~q:~g",[A,P]).
 1395
 1396write_head(S,[A:P]):-!,
 1397  format(S,"~q:~g",[A,P]).
 1398
 1399write_head(S,[A:P|Rest]):-
 1400  format(S,"~q:~g ; ",[A,P]),
 1401  write_head(S,Rest).
 1402
 1403write_body(S,[]):-!,
 1404  format(S,"  true.~n~n",[]).
 1405
 1406write_body(S,[A]):-!,
 1407  format(S,"  ~q.~n~n",[A]).
 1408
 1409write_body(S,[A|T]):-
 1410  format(S,"  ~q,~n",[A]),
 1411  write_body(S,T).