annotate cpack/dml/lib/grammars.pl @ 0:718306e29690 tip

commiting public release

/* Part of DML (Digital Music Laboratory)
	Copyright 2014-2015 Samer Abdallah, University of London
	 
	This program is free software; you can redistribute it and/or
	modify it under the terms of the GNU General Public License
	as published by the Free Software Foundation; either version 2
	of the License, or (at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public
	License along with this library; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/

:- module(grammars,
      [  model_name/2
      ,  model_module_prep/3
      ,  build_options/2
      ,  build_subset_options/2
      ,  learn/6
      ,  learn_bpe/6
      ,  learn_model/7
      ,  model_sequence_parses/4
      ,  nnums_ivals/2
      ,  restart_prism/0
      ,  best_first/6
      ,  method_model_dataset_results/4
      ,  dataset_num_events/2
      ]).

:- multifile dataset_sequences/2.

:- use_module(library(memo)).
:- use_module(library(typedef)).
:- use_module(library(lambda)).
:- use_module(library(plml)).
:- use_module(library(prism/prism)).
:- use_module(library(argutils)).
:- use_module(library(snobol)).

:- type pmodule == ground.
:- type natural == nonneg.
:- type prep    == callable.
:- type options == ordset.
:- type model  ---> model(atom, list, prep, list(pair(ground,list(number)))).
:- type method ---> vb(ground) ; map(ground).
:- type matlab_method ---> vb ; map.

% % hmmmn.. maybe module issues here...
% error:has_type(\Checker,Term) :- call(Checker,Term).
error:has_type(ordset,Term) :- is_ordset(Term).

:- persistent_memo 
      learn( +method, +pmodule, +prep, +dataset:ground, +options, -scores:list),
      learn_model( +method, +pmodule, +prep, +dataset:ground, +options, -scores:list, -model),
      dataset_num_events( +dataset:ground, -num_events:nonneg).

dataset_num_events(Dataset,NumEvents) :-
   dataset_sequences(Dataset,Seqs),
   maplist(length,Seqs,Lens),
   sumlist(Lens,NumEvents).
   % aggregate_all(sum(L),(member(S,Seqs),length(S,L)),NumEvents).

:- initialization memo_attach(memo(learned),[]).
:- setting(timeout, number, 900, 'Time limit for learning in seconds').

user:file_search_path(prism,'psm').
user:matlab_path(grammars,['stats/hmm']).

partition_options([],[],[],[]).
partition_options([O|OX],[O|MO],IO,LO) :- option_class(O,model), partition_options(OX,MO,IO,LO).
partition_options([O|OX],MO,[O|IO],LO) :- option_class(O,init), partition_options(OX,MO,IO,LO).
partition_options([O|OX],MO,IO,[O|LO]) :- option_class(O,learn), partition_options(OX,MO,IO,LO).

option_class(switch_mode(_),model).
option_class(prior_weight(_),model).
option_class(gamut(_),model).
option_class(leap_range(_),model).
option_class(log_scale(_),learn).
option_class(O,init) :- \+option_class(O,model), \+option_class(O,learn).

option_mlopt(gamut(X-Y),gamut:[X,Y]) :- !.
option_mlopt(init(none),perturb:0) :- !.
option_mlopt(init(perturb(U)),perturb:U) :- !.
option_mlopt(F,N:V) :- F=..[N,V].

learn(vb(InitMeth),matlab(Function),Prepare,DataSet,Opts,[free_energy(FE)]) :- !,
   dataset_sequences(DataSet, D1),
	maplist(Prepare,D1,D),
   maplist(option_mlopt,[init(InitMeth)|Opts],Opts1),
   compileoptions(Opts1,Opts2),
   [float(FE)]===feval(@Function,cell(D),Opts2).

   % method dependent options

learn(Method,Module,Prepare,DataSet,Opts,Scores) :-
   member(Method,[vb(_),map(_)]),
   % prepare data
   dataset_sequences(DataSet, D1),
	maplist(Prepare,D1,D),
   % method dependent options
   method_switch_mode(Method,Mode),
   option(log_scale(LS), Opts, on),
   % load up prism and initialise
   restart_prism, load_prism(prism(Module)), 
   #init_model([switch_mode(Mode)|Opts]),
   #init_switches(Opts),
   set_prism_flag(log_scale,LS),
   % allow 15 minutes for learning
   setting(timeout, TimeLimit),
   get_time(Time),
   debug(learn,'Calling PRISM with time limit ~w at ~@...',
         [TimeLimit, format_time(current_output,'%T',Time)]),
   call_with_time_limit(TimeLimit, prism_learn(Method, D, [], Scores)).

method_switch_mode(vb(_),a).
method_switch_mode(map(_),d).

learn_model(Method,Module,Prepare,DataSet,Opts,Scores,model(Module,ModelOpts,Prepare,Counts)) :-
   % must re-compute to get final state
   compute(learn(Method,Module,Prepare,DataSet,Opts,Scores)),
   partition_options(Opts,ModelOpts,_,_),
   get_prism_state(ps(_,_,Counts1,_,_,_)),
   map_filter(unfixed_count,Counts1,Counts).

unfixed_count(sw(SW,a,set(unfixed,Counts)), SW-Counts).

model_initial_state(model(Module,MO,_,_),State) :-
   restart_prism,
   load_prism(prism(Module)),
   #init_model([switch_mode(a)|MO]),
   get_prism_state(State).

with_model_data(model(Module,MO,Prepare,Counts),Data,Goal) :-
   restart_prism, load_prism(prism(Module)),
   #init_model([switch_mode(a)|MO]),
   #maplist(SW-C,set_sw_a(SW,C),Counts),
   call(Prepare,Data,D),
   call(Goal,D).

learn_bpe(Method,Module,Prepare,DataSet,Opts,BPE) :-
   browse(learn(Method,Module,Prepare,DataSet,Opts,Scores)),
   catch( bits_per_event(Method,DataSet,Scores,BPE), _, fail).

bits_per_event(Method,DS,Scores,BPE) :-
   dataset_num_events(DS,NumEvents),
   score(Method,NumEvents,Scores,BPE).

score(map(_),NumEvents,Scores,BitsPerEvent) :-
	member(log_lik(LL), Scores),
	BitsPerEvent is -(LL/NumEvents)/log(2).

score(vb(_),NumEvents,Scores,BitsPerEvent) :-
	member(free_energy(LL), Scores),
	BitsPerEvent is -(LL/NumEvents)/log(2).

score(vb_pm(_),NumEvents,Scores,BitsPerEvent) :-
	member(free_energy(LL), Scores),
	BitsPerEvent is -(LL/NumEvents)/log(2).

%% tree_syntax(+Mod:module,+Tree:prism_tree,-Syntax:tree) is det.
%
%  Create a parse tree from a PRISM Viterbi tree.
%  Works for models gilbert1, gilbert2, gilbert2a, gilbert3 and gilbert2m.
tree_syntax(Mod,[s(_),TT],T2) :- tree_parse_tree(Mod,TT,T2).

tree_parse_tree(_,msw(i(I),terminal),node(t(I),[])).
tree_parse_tree(Mod,[pp(s,_,_)|Children],Term) :- member(Mod,[gilbert2,gilbert2a,gilbert3,gilbert2m]), !, 
   member(msw(s,Rule),Children),
   map_filter(tree_parse_tree(Mod),Children,CN),
   (  Rule=grow -> CN=[Child1,T1], T1=node(s,Tail), Term=node(s,[Child1|Tail])
   ;  Rule=last -> CN=[Child], Term=node(s,[Child])
   ).
tree_parse_tree(Mod,[pp(s,_,_)|Children],Term) :- Mod=gilbert1, !,
   member(msw(s,Rule),Children),
   map_filter(tree_parse_tree(Mod),Children,CN),
   (  Rule=grow -> CN=[Child1,T1], T1=node(s,Tail), Term=node(s,[Child1|Tail])
   ;  Rule=first -> CN=[], Term=node(s,[])
   ).
tree_parse_tree(Mod,[pp(H,_,_)|Children],Term) :- !, 
   map_filter(tree_parse_tree(Mod),Children,CN),
   member(msw(H,Rule1),Children),
   ( Rule1=terminal -> Rule=t; Rule=Rule1),
   Term = node(H-Rule,CN).

:- volatile_memo model_sequence_parses(+ground,+list(ground),+natural,-ground).

model_sequence_parses(Model,Seq,N,Parses) :- 
   Model=model(Mod,_,_,_), 
   with_model_data(Model,Seq,parses(Mod,N,Parses)).

parses(Mod,N,Parses,Goal) :- 
   succ(N,M),
   findall(P-T,viterbi_tree(M,Goal,P,[],T),ProbsTrees),
   append(NProbsTrees,[P0-_],ProbsTrees),
   maplist(tree_parse(Mod,P0),NProbsTrees,Parses).

tree_parse(Mod,P0,P-T,RP-S) :- tree_syntax(Mod,T,S), RP is P/P0.


% model declarations
decl_model(markov(nnum,0), p1gram,   markovp, with_nnums(s0)).
decl_model(markov(nnum,1), p2gram,   markovp, with_nnums(s1)).
decl_model(markov(ival,0), i1gram,   markovi, with_ivals(s0)).
decl_model(markov(ival,1), i2gram,   markovi, with_ivals(s1)).
decl_model(gilbert(1),     gilbert1, gilbert1, with_pre_ivals(s)).
decl_model(gilbert(2),     gilbert2, gilbert2, with_ivals(s)).
decl_model(gilbert(3),     gilbert3, gilbert3, with_ivals(s)).
decl_model(gilbert(2-a),   gilbert2a, gilbert2a, with_ivals(s)).
decl_model(gilbert(2-m),   gilbert2m, gilbert2m, with_ivals(s)).
% decl_model(gilbert(4),     gilbert1, gilbert1a, with_pre_ivals(s)).
% decl_model(gilbert(5),     gilbert2, gilbert2a, with_ivals(s)).
% decl_model(gilbert(6),     gilbert3, gilbert3a, with_ivals(s)).
decl_model(hmm(nnum),      phmm,     phmm,     with_nnums(s)).
decl_model(hmm(nnum,NS),   Name,     phmm,     with_nnums(s)) :- atom_concat(phmm,NS,Name).
%decl_model(hmm(ival),      ihmm,     ihmm,     with_ivals(s)).
decl_model(matlab(p1gram), 'ml-p1gram', matlab(p1gram), (=)). 
decl_model(matlab(p2gram), 'ml-p2gram', matlab(p2gram), (=)). 
decl_model(matlab(phmm),   'ml-phmm', matlab(phmm), (=)). 

model_name(Model,Name) :- decl_model(Model,Name,_,_).
model_module_prep(Model,Module,Prepare) :- decl_model(Model,_,Module,Prepare).

with_nnums(Head, Seq, Head1) :-
	addargs(Head,[Seq],Head1).

with_ivals(Head, Seq, Head1) :-
	nnums_ivals(Seq,Seq1),
	addargs(Head,[Seq1],Head1).

with_pre_ivals(Head, Seq, Head1) :-
	nnums_pre_ivals(Seq,Seq1),
	addargs(Head,[Seq1],Head1).

nnums_ivals(NNums,Ivals) :- nnums_post_ivals(NNums,Ivals).
nnums_post_ivals(NNums,Ivals) :- phrase((ivals(NNums),[end]),Ivals,[]).
nnums_pre_ivals(NNums,Ivals) :- phrase(([start],ivals(NNums)),Ivals,[]).

ivals([X0,X1|Xs]) --> {I1 is X1-X0}, [I1], ivals([X1|Xs]).
ivals([_]) --> [].

% NB option defaults here must match those in PRISM source files.
model_options(markov(nnum,_)) -->
   optopt(prior_shape, [binomial+0.1*uniform, binomial+uniform, uniform]),
   optopt(gamut,[40-100]).

model_options(markov(ival,_)) -->
   optopt(prior_shape, [uniform, binomial+uniform, binomial+0.1*uniform]).

model_options(gilbert(_)) -->
   [leap_range((-20)-(20))],
   optopt(leap_shape, [uniform, binomial+uniform, binomial+0.1*uniform]),
   optopt(pass_shape, [binomial, binomial+uniform, binomial+0.1*uniform]).

model_options(hmm(nnum)) -->
   optopt(prior_shape, [binomial+0.1*uniform, binomial+uniform, uniform]),
   anyopt(num_states, [1,2,3,5,7,12,18]),
   optopt(trans_self, [1]),
   % optopt(gamut,[40-100]).
   % anyopt(trans_self, [1]),
   anyopt(gamut,[40-100]).

model_options(hmm(nnum,NS)) -->
   [num_states(NS)],
   optopt(prior_shape, [binomial+0.1*uniform, binomial+uniform, uniform]),
   optopt(trans_self, [1]),
   % optopt(gamut,[40-100]).
   % anyopt(trans_self, [1]),
   anyopt(gamut,[40-100]).

model_options(hmm(ival)) -->
   optopt(prior_shape, [binomial+0.1*uniform, binomial+uniform, uniform]),
   anyopt(num_states, [5,7,12,24]),
   [leap_range((-20)-(20))].

model_options(matlab(p1gram)) -->
   optopt(prior_shape, [binomial+0.1*uniform, binomial+uniform, uniform]),
   optopt(gamut,[40-100]).

model_options(matlab(p2gram)) -->
   optopt(prior_shape, [binomial+0.1*uniform, binomial+uniform, uniform]),
   optopt(gamut,[40-100]).

model_options(matlab(phmm)) -->
   optopt(prior_shape, [binomial+0.1*uniform, binomial+uniform, uniform]),
   anyopt(num_states, [1,2,3,5,7,12,18]),
   optopt(gamut,[40-100]).

model_subset_options(markov(nnum,_)) --> [prior_shape(binomial+0.1*uniform)].
model_subset_options(markov(ival,_)) --> [prior_shape(binomial+0.1*uniform)].
model_subset_options(gilbert(_)) --> [].
model_subset_options(hmm(nnum)) --> 
   [prior_shape(binomial+0.1*uniform)], 
   anyopt(num_states,[2,5,7,12,18]).
model_subset_options(hmm(ival)) --> 
   [prior_shape(binomial+0.1*uniform), leap_range((-20)-(20))],
   anyopt(num_states,[2,5,7,12]).

build_options(Model,Opts) :-
   build_options(Model,Opts1,[]),
   sort(Opts1,Opts).

build_options(Model) -->
   optopt(prior_weight,[0.3,3,10]), % NB have removed 0.1 and 30
   model_options(Model).

build_subset_options(Model,Opts) :- 
   build_subset_options(Model,Opts1,[]),
   sort(Opts1,Opts).

build_subset_options(Model) -->
   optopt(prior_weight,[0.1,0.3,3,10,30]), 
   model_subset_options(Model).

anyopt(Name,Vals) --> {maplist(\X^Y^(Y=..[Name,X]),Vals,Opts)}, any(Opts).
optopt(Name,Vals) --> []; anyopt(Name,Vals).

best_first(Meth,Mod,Prepare,DS,Opts,learn(Meth,Mod,Prepare,DS,Opts,F)) :-
   order_by([asc(F)], learn_bpe(Meth,Mod,Prepare,DS,Opts,F)).

/* TODO: 
  Sort out HMM init options. They are all wrong:
      init_shape : shape of prior over initial state {uniform}
      trans_shape : shape of prior over transition distribution {uniform}
      trans_persistence : add self transtion counts {0,1,3}
      init_noise : perturbation of initial obs counts {0,0.1}
      restarts: {3}

   

   compare p2gram with hmm(1)
      */


%% map_filter(+P:pred(A,B),+Xs:list(A),-Ys:list(B)) is det.
%
%  map_filter(P,Xs,Ys) is similar to maplist(P,Xs,Ys), except that P is allowed to
%  fail and the resulting list Ys contains only those elements Y for which call(P,X,Y).
%  P is used as if it was semi-deterministic: only the first solution is accepted.
map_filter(_,[],[]).
map_filter(P,[X|XX],[Y|YY]) :- call(P,X,Y), !, map_filter(P,XX,YY).
map_filter(P,[_|XX],YY) :- map_filter(P,XX,YY).

user:file_search_path(home,X) :- expand_file_name('~',[X]).
user:file_search_path(prism,home('src/probprog/prism/psm')).
restart_prism :- prism_start('prism.log').

retry_on_error(G,M) :-
   catch((restart_prism,G,Status=ok), Ex, Status=error(Ex)),
   (  Status=ok -> true
   ;  format('*** terminated with ~w.\n',[Ex]),
      shell('say ALERT: THERE WAS AN ERROR!'),
      (  succ(N,M) -> retry_on_error(G,N)
      ;  writeln('failing after too many retries.'), 
         shell('say ALERT: I AM GIVING UP!'),
         shell('say ALERT: I AM GIVING UP!'),
         fail
      )
   ).

%% summary_array(+Meth:method,+Models:list(model),+Datasets:list(dataset),-Summary) is det.
summary_array(Meth,Models,Datasets,arr(Summary)) :-
   maplist( \Model^RR^maplist(\DS^{arr(R)}^method_model_dataset_results(Meth,Model,DS,R),Datasets,RR), Models,Summary).

%% method_model_dataset_results(+Method,+Model,+Dataset,Results:list(list(number))) is det.
method_model_dataset_results(Meth,Model,DS,Results) :- bagof( Row, datum(Meth,Model,DS,Row), Results).

%% datum(+Method:method,+Model:model,+Dataset,-Datum:list(number)) is nondet.
%
%  Maps trained models to a numerical tuple consisting of
%  the free energy (in bits per event), the prior weight parameter, and one or two 
%  shape parameters, depending on the model. The shape parameter is 0 for a uniform
%  prior, 1 for a binomial prior, and between 0 and 1 for linear interpolation between
%  the two.
%
%  Grammar models have leap_shape and pass_shape parameters, while the others
%  have a single prior_shape parameter.
datum(Meth,gilbert(I),DS,[FE,W,K1,K2]) :-
   model_module_prep(gilbert(I),Mod,Prep),
   learn_bpe(Meth,Mod,Prep,DS,Opts,FE),
   option(prior_weight(W),Opts,1),
   option(leap_shape(Sh1),Opts,binomial),
   option(pass_shape(Sh2),Opts,uniform),
   shape_param(Sh1,K1),
   shape_param(Sh2,K2).

datum(Meth,markov(Rep,Ord),DS,[FE,W,K]) :-
   model_module_prep(markov(Rep,Ord),Mod,Prep),
   learn_bpe(Meth,Mod,Prep,DS,Opts,FE),
   option(prior_weight(W),Opts,1),
   option(prior_shape(Sh),Opts,binomial),
   shape_param(Sh,K).

datum(Meth,hmm(Rep),DS,[FE,W,K]) :-
   model_module_prep(hmm(Rep),Mod,Prep),
   learn_bpe(Meth,Mod,Prep,DS,Opts,FE),
   option(prior_weight(W),Opts,1),
   option(prior_shape(Sh),Opts,binomial),
   shape_param(Sh,K).

datum(Meth,hmm(Rep,NS),DS,[FE,W,K]) :-
   model_module_prep(hmm(Rep),Mod,Prep),
   learn_bpe(Meth,Mod,Prep,DS,Opts,FE),
   option(num_states(NS),Opts),
   option(prior_weight(W),Opts,1),
   option(prior_shape(Sh),Opts,binomial),
   shape_param(Sh,K).

datum(Meth,matlab(F),DS,[FE,W,K]) :-
   model_module_prep(matlab(F),Mod,Prep),
   learn_bpe(Meth,Mod,Prep,DS,Opts,FE),
   option(prior_weight(W),Opts,1),
   option(prior_shape(Sh),Opts,binomial),
   shape_param(Sh,K).

shape_param(uniform,0).
shape_param(binomial,1).
shape_param(binomial+uniform,0.5).
shape_param(uniform+binomial,0.5).
shape_param(binomial+K*uniform,Lam) :- Lam = 1/(1+K).