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

commiting public release

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cpack/dml/lib/grammars.pl	Tue Feb 09 21:05:06 2016 +0100
@@ -0,0 +1,425 @@
+/* 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).
+