Daniel@0: % dataset      -> (1=>user data) or (2=>toy example)
Daniel@0: % type         -> (1=> Regression model) or (2=>Classification model)
Daniel@0: % num_glevel   -> number of hidden nodes in the net (gating levels)
Daniel@0: % num_exp      -> number of experts in the net
Daniel@0: % branch_fact  -> dimension of the hidden nodes in the net
Daniel@0: % cov_dim      -> root node dimension
Daniel@0: % res_dim      -> output node dimension
Daniel@0: % nodes_info   -> 4 x num_glevel+2 matrix that contain all the info about the nodes:
Daniel@0: %                 nodes_info(1,:) = nodes type: (0=>gaussian)or(1=>softmax)or(2=>mlp)
Daniel@0: %                 nodes_info(2,:) = nodes size: [cov_dim   num_glevel x branch_fact   res_dim]
Daniel@0: %                 nodes_info(3,:) = hidden units number (for mlp nodes)
Daniel@0: %                                  |- optimizer iteration number (for softmax & mlp CPD)
Daniel@0: %                 nodes_info(4,:) =|- covariance type (for gaussian CPD)-> 
Daniel@0: %                                  | (1=>Full)or(2=>Diagonal)or(3=>Full&Tied)or(4=>Diagonal&Tied)
Daniel@0: % fh1 -> Figure: data & decizion boundaries; fh2 -> confusion matrix; fh3 -> LL trace                                                                       
Daniel@0: % test_data    -> test data matrix
Daniel@0: % train_data   -> training data matrix
Daniel@0: % ntrain       -> size(train_data,2)
Daniel@0: % ntest        -> size(test_data,2)
Daniel@0: % cases        -> (cell array) training data formatted for the learning engine
Daniel@0: % bnet         -> bayesian net before learning
Daniel@0: % bnet2        -> bayesian net after learning
Daniel@0: % ll           -> log-likelihood before learning
Daniel@0: % LL2          -> log-likelihood trace
Daniel@0: % onodes       -> obs nodes in bnet & bnet2
Daniel@0: % max_em_iter  -> maximum number of interations of the EM algorithm
Daniel@0: % train_result -> prediction on the training set (as test_result)
Daniel@0: % 
Daniel@0: % IMPORTANT: CHECK the loading path (lines 64 & 364)
Daniel@0: % ----------------------------------------------------------------------------------------------------
Daniel@0: % -> pierpaolo_b@hotmail.com   or   -> pampo@interfree.it
Daniel@0: % ----------------------------------------------------------------------------------------------------
Daniel@0: 
Daniel@0: error('this no longer works with the latest version of BNT')
Daniel@0: 
Daniel@0: clear all;
Daniel@0: clc;
Daniel@0: disp('---------------------------------------------------');
Daniel@0: disp('  Hierarchical Mixtures of Experts models builder ');
Daniel@0: disp('---------------------------------------------------');
Daniel@0: disp(' ')
Daniel@0: disp('   Using this script you can build both an HME model')
Daniel@0: disp('as in [Wat94] and [Jor94] i.e. with ''softmax'' gating')
Daniel@0: disp('nodes and ''gaussian'' ( for regression ) or ''softmax''') 
Daniel@0: disp('( for classification ) expert node, and its variants')
Daniel@0: disp('called ''gated nets'' where we use ''mlp'' models in')
Daniel@0: disp('place of a number of ''softmax'' ones [Mor98], [Wei95].')
Daniel@0: disp('  You can decide to train and test the model on your')
Daniel@0: disp('datasets  or  to evaluate its  performance on  a toy')
Daniel@0: disp('example.')
Daniel@0: disp(' ')
Daniel@0: disp('Reference')
Daniel@0: disp('[Mor98] P. Moerland (1998):')
Daniel@0: disp('        Localized mixtures of experts. (http://www.idiap.ch/~perry/)')
Daniel@0: disp('[Jor94] M.I. Jordan, R.A. Jacobs (1994):')
Daniel@0: disp('        HME and the EM algorithm. (http://www.cs.berkeley.edu/~jordan/)')
Daniel@0: disp('[Wat94] S.R. Waterhouse, A.J. Robinson (1994):') 
Daniel@0: disp('        Classification using HME. (http://www.oigeeza.com/steve/)')
Daniel@0: disp('[Wei95] A.S. Weigend, M. Mangeas (1995):') 
Daniel@0: disp('        Nonlinear gated experts for time series.')
Daniel@0: disp(' ')
Daniel@0: 
Daniel@0: if 0
Daniel@0: disp('(See the figure)')
Daniel@0: pause(5);
Daniel@0: %%%%%WARNING!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: im_path=which('HMEforMatlab.jpg');
Daniel@0: fig=imread(im_path, 'jpg');
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: figure('Units','pixels','MenuBar','none','NumberTitle','off', 'Name', 'HME model');
Daniel@0: image(fig); 
Daniel@0: axis image;
Daniel@0: axis off;
Daniel@0: clear fig;
Daniel@0: set(gca,'Position',[0 0 1 1])
Daniel@0: disp('(Press any key to continue)')
Daniel@0: pause
Daniel@0: end
Daniel@0: 
Daniel@0: clc
Daniel@0: disp('---------------------------------------------------');
Daniel@0: disp('              Specify the Architecture             ');
Daniel@0: disp('---------------------------------------------------');
Daniel@0: disp(' ');
Daniel@0: disp('What kind of model do you need?')
Daniel@0: disp(' ')
Daniel@0: disp('1) Regression ')
Daniel@0: disp('2) Classification')
Daniel@0: disp(' ')
Daniel@0: type=input('1 or 2?: ');
Daniel@0: if (isempty(type)|(~ismember(type,[1 2]))), error('Invalid value'); end
Daniel@0: clc
Daniel@0: disp('----------------------------------------------------');
Daniel@0: disp('               Specify the Architecture             ');
Daniel@0: disp('----------------------------------------------------');
Daniel@0: disp(' ')
Daniel@0: disp('Now you have to set the number of experts and gating')
Daniel@0: disp('levels in the net.  This script builds only balanced')
Daniel@0: disp('hierarchy with the same branching factor (>1)at each')
Daniel@0: disp('(gating) level. So remember that: ')
Daniel@0: disp(' ')
Daniel@0: disp('         num_exp = branch_fact^num_glevel           ')
Daniel@0: disp(' ')
Daniel@0: disp('with branch_fact >=2.')
Daniel@0: disp('You can also set to zeros the number of gating level')
Daniel@0: disp('in order to obtain a classical GLM model.           ')
Daniel@0: disp(' ')
Daniel@0: disp('----------------------------------------------------');
Daniel@0: disp(' ')
Daniel@0: num_glevel=input('Insert the number of gating levels {0,...,20}: ');
Daniel@0: if (isempty(num_glevel)|(~ismember(num_glevel,[0:20]))), error('Invalid value'); end
Daniel@0: nodes_info=zeros(4,num_glevel+2);
Daniel@0: if num_glevel>0, %------------------------------------------------------------------------------------
Daniel@0:     for i=2:num_glevel+1,
Daniel@0:         clc
Daniel@0:         disp('----------------------------------------------------');
Daniel@0:         disp('               Specify the Architecture             ');
Daniel@0:         disp('----------------------------------------------------');
Daniel@0:         disp(' ')   
Daniel@0:         disp(['-> Gating network ', num2str(i-1), ' is a: '])
Daniel@0:         disp(' ')
Daniel@0:         disp('   1) Softmax model');
Daniel@0:         disp('   2) Two layer perceptron model')
Daniel@0:         disp(' ')
Daniel@0:         nodes_info(1,i)=input('1 or 2?: ');
Daniel@0:         if (isempty(nodes_info(1,i))|(~ismember(nodes_info(1,i),[1 2]))), error('Invalid value'); end
Daniel@0:         disp(' ')
Daniel@0:         if nodes_info(1,i)==2,
Daniel@0:            nodes_info(3,i)=input('Insert the number of units in the hidden layer: ');
Daniel@0:            if (isempty(nodes_info(3,i))|(floor(nodes_info(3,i))~=nodes_info(3,i))|(nodes_info(3,i)<=0)), 
Daniel@0:               error(['Invalid value: ', num2str(nodes_info(3,i)), ' is not a positive integer!']);
Daniel@0:            end
Daniel@0:            disp(' ')
Daniel@0:         end
Daniel@0:         nodes_info(4,i)=input('Insert the optimizer iteration number: ');
Daniel@0:         if (isempty(nodes_info(4,i))|(floor(nodes_info(4,i))~=nodes_info(4,i))|(nodes_info(4,i)<=0)), 
Daniel@0:            error(['Invalid value: ', num2str(nodes_info(4,i)), ' is not a positive integer!']);
Daniel@0:         end    
Daniel@0:     end
Daniel@0:     clc
Daniel@0:     disp('---------------------------------------------------------');
Daniel@0:     disp('                 Specify the Architecture                ');
Daniel@0:     disp('---------------------------------------------------------');
Daniel@0:     disp(' ')
Daniel@0:     disp('Now you have to set the number  of experts in the network');
Daniel@0:     disp('The value will be adjusted in order to obtain a hierarchy');
Daniel@0:     disp('as said above.')
Daniel@0:     disp(' ');    
Daniel@0:     num_exp=input(['Insert the approximative number of experts (>=', num2str(2^num_glevel), '): ']);
Daniel@0:     if (isempty(num_exp)|(num_exp<=0)|(num_exp<2^num_glevel)), 
Daniel@0:         error('Invalid value');
Daniel@0:     end
Daniel@0:     app1=0; base=2;
Daniel@0:     while app1<num_exp,
Daniel@0:         app1=base^num_glevel;
Daniel@0:         base=base+1;
Daniel@0:     end
Daniel@0:     app2=(base-2)^num_glevel;
Daniel@0:     branch_fact=base-1;
Daniel@0:     if app2>=(2^num_glevel)&(abs(app2-num_exp)<abs(app1-num_exp)),
Daniel@0:         branch_fact=base-2;
Daniel@0:     end
Daniel@0:     clear app1 app2 base;
Daniel@0:     disp(' ')
Daniel@0:     disp(['The effective number of experts in the net is: ', num2str(branch_fact^num_glevel), '.'])
Daniel@0:     disp(' ');
Daniel@0: else
Daniel@0:     clc
Daniel@0:     disp('---------------------------------------------------------');
Daniel@0:     disp('        Specify the Architecture (GLM model)             ');
Daniel@0:     disp('---------------------------------------------------------');
Daniel@0:     disp(' ')
Daniel@0: end % END of: if num_glevel>0-------------------------------------------------------------------------
Daniel@0: 
Daniel@0: if type==2,
Daniel@0:     disp(['-> Expert node is a: '])
Daniel@0:     disp(' ')
Daniel@0:     disp('   1) Softmax model');
Daniel@0:     disp('   2) Two layer perceptron model')
Daniel@0:     disp(' ')
Daniel@0:     nodes_info(1,end)=input('1 or 2?: ');
Daniel@0:     if (isempty(nodes_info(1,end))|(~ismember(nodes_info(1,end),[1 2]))), 
Daniel@0:         error('Invalid value'); 
Daniel@0:     end
Daniel@0:     disp(' ')
Daniel@0:     if nodes_info(1,end)==2,
Daniel@0:        nodes_info(3,end)=input('Insert the number of units in the hidden layer: ');
Daniel@0:        if (isempty(nodes_info(3,end))|(floor(nodes_info(3,end))~=nodes_info(3,end))|(nodes_info(3,end)<=0)), 
Daniel@0:            error(['Invalid value: ', num2str(nodes_info(3,end)), ' is not a positive integer!']);
Daniel@0:        end
Daniel@0:        disp(' ')
Daniel@0:     end
Daniel@0:     nodes_info(4,end)=input('Insert the optimizer iteration number: ');
Daniel@0:     if (isempty(nodes_info(4,end))|(floor(nodes_info(4,end))~=nodes_info(4,end))|(nodes_info(4,end)<=0)), 
Daniel@0:         error(['Invalid value: ', num2str(nodes_info(4,end)), ' is not a positive integer!']);
Daniel@0:     end
Daniel@0: elseif type==1,
Daniel@0:     disp('What kind of covariance matrix structure do you want?')
Daniel@0:     disp(' ')
Daniel@0:     disp('   1) Full');
Daniel@0:     disp('   2) Diagonal')
Daniel@0:     disp('   3) Full & Tied');
Daniel@0:     disp('   4) Diagonal & Tied')
Daniel@0: 
Daniel@0:     disp(' ')
Daniel@0:     nodes_info(4,end)=input('1, 2, 3 or 4?: ');
Daniel@0:     if (isempty(nodes_info(4,end))|(~ismember(nodes_info(4,end),[1 2 3 4]))), 
Daniel@0:         error('Invalid value'); 
Daniel@0:     end  
Daniel@0: end
Daniel@0: clc
Daniel@0: disp('----------------------------------------------------');
Daniel@0: disp('                    Specify the Input               ');
Daniel@0: disp('----------------------------------------------------');
Daniel@0: disp(' ')
Daniel@0: disp('Do you want to...')
Daniel@0: disp(' ')
Daniel@0: disp('1) ...use your own dataset?')
Daniel@0: disp('2) ...apply the model on a toy example?')
Daniel@0: disp(' ')
Daniel@0: dataset=input('1 or 2?: ');
Daniel@0: if (isempty(dataset)|(~ismember(dataset,[1 2]))), error('Invalid value'); end
Daniel@0: if dataset==1,
Daniel@0:     if type==1,
Daniel@0:         clc
Daniel@0:         disp('-------------------------------------------------------');
Daniel@0:         disp('        Specify the Input - Regression problem         ');
Daniel@0:         disp('-------------------------------------------------------');
Daniel@0:         disp(' ')
Daniel@0:         disp('Be sure that each row of your data matrix is an example');
Daniel@0:         disp('with the covariate values that precede the respond ones')
Daniel@0:         disp(' ')
Daniel@0:         disp('-------------------------------------------------------');
Daniel@0:         disp(' ')
Daniel@0:         cov_dim=input('Insert the covariate space dimension: ');
Daniel@0:         if (isempty(cov_dim)|(floor(cov_dim)~=cov_dim)|(cov_dim<=0)), 
Daniel@0:           error(['Invalid value: ', num2str(cov_dim), ' is not a positive integer!']);
Daniel@0:         end
Daniel@0:         disp(' ')
Daniel@0:         res_dim=input('Insert the dimension of the respond variable: ');
Daniel@0:         if (isempty(res_dim)|(floor(res_dim)~=res_dim)|(res_dim<=0)), 
Daniel@0:             error(['Invalid value: ', num2str(res_dim), ' is not a positive integer!']);
Daniel@0:         end 
Daniel@0:         disp(' ');
Daniel@0:     elseif type==2
Daniel@0:         clc
Daniel@0:         disp('-------------------------------------------------------');
Daniel@0:         disp('      Specify the Input - Classification problem       ');
Daniel@0:         disp('-------------------------------------------------------');
Daniel@0:         disp(' ')
Daniel@0:         disp('Be sure that each row of your data matrix is an example');
Daniel@0:         disp('with the covariate values that precede the class labels');
Daniel@0:         disp('(integer value >=1).                                   ');
Daniel@0:         disp(' ')
Daniel@0:         disp('-------------------------------------------------------');
Daniel@0:         disp(' ')
Daniel@0:         cov_dim=input('Insert the covariate space dimension: ');
Daniel@0:         if (isempty(cov_dim)|(floor(cov_dim)~=cov_dim)|(cov_dim<=0)), 
Daniel@0:           error(['Invalid value: ', num2str(cov_dim), ' is not a positive integer!']);
Daniel@0:         end
Daniel@0:         disp(' ')
Daniel@0:         res_dim=input('Insert the number of classes: ');
Daniel@0:         if (isempty(res_dim)|(floor(res_dim)~=res_dim)|(res_dim<=0)), 
Daniel@0:           error(['Invalid value: ', num2str(res_dim), ' is not a positive integer!']);
Daniel@0:         end        
Daniel@0:         disp(' ')               
Daniel@0:     end    
Daniel@0:     % ------------------------------------------------------------------------------------------------
Daniel@0:     % Loading training data --------------------------------------------------------------------------
Daniel@0:     % ------------------------------------------------------------------------------------------------
Daniel@0:     train_path=input('Insert the complete (with extension) path of the training data file:\n >> ','s');    
Daniel@0:     if isempty(train_path), error('You must specify a data set for training!'); end
Daniel@0:     if ~isempty(findstr('.mat',train_path)),
Daniel@0:         ap=load(train_path); app=fieldnames(ap); train_data=eval(['ap.', app{1,1}]);
Daniel@0:         clear ap app;
Daniel@0:     elseif ~isempty(findstr('.txt',train_path)),
Daniel@0:         train_data=load(train_path, '-ascii');
Daniel@0:     else
Daniel@0:         error('Invalid data format: not a .mat or a .txt file')
Daniel@0:     end
Daniel@0:     if (size(train_data,2)~=cov_dim+res_dim)&(type==1),
Daniel@0:         error(['Invalid data matrix size: ', num2str(size(train_data,2)), ' columns rather than ',...
Daniel@0:             num2str(cov_dim+res_dim),'!']);
Daniel@0:     elseif (size(train_data,2)~=cov_dim+1)&(type==2),
Daniel@0:         error(['Invalid data matrix size: ', num2str(size(train_data,2)), ' columns rather than ',...
Daniel@0:             num2str(cov_dim+1),'!']);    
Daniel@0:     elseif (~isempty(find(ismember(intersect([train_data(:,end)' 1:res_dim],...
Daniel@0:             train_data(:,end)'),[1:res_dim])==0)))&(type==2),
Daniel@0:         error('Invalid class label');
Daniel@0:     end    
Daniel@0:     ntrain=size(train_data,1);
Daniel@0:     train_d=train_data(:,1:cov_dim);
Daniel@0:     if type==2,
Daniel@0:         train_t=zeros(ntrain, res_dim);
Daniel@0:         for m=1:res_dim,
Daniel@0:             train_t((find(train_data(:,end)==m))',m)=1;
Daniel@0:         end
Daniel@0:     else
Daniel@0:         train_t=train_data(:,cov_dim+1:end);
Daniel@0:     end        
Daniel@0:     disp(' ')
Daniel@0:     % ------------------------------------------------------------------------------------------------
Daniel@0:     % Loading test data ------------------------------------------------------------------------------
Daniel@0:     % ------------------------------------------------------------------------------------------------
Daniel@0:     disp('(If you don''t want to specify a test-set press ''return'' only)');
Daniel@0:     test_path=input('Insert the complete (with extension) path of the test data file:\n >> ','s');  
Daniel@0:     if ~isempty(test_path),
Daniel@0:         if ~isempty(findstr('.mat',test_path)),        
Daniel@0:             ap=load(test_path); app=fieldnames(ap); test_data=eval(['ap.', app{1,1}]);
Daniel@0:             clear ap app;
Daniel@0:         elseif ~isempty(findstr('.txt',test_path)),
Daniel@0:             test_data=load(test_path, '-ascii');
Daniel@0:         else
Daniel@0:             error('Invalid data format: not a .mat or a .txt file')
Daniel@0:         end
Daniel@0:         if (size(test_data,2)~=cov_dim)&(size(test_data,2)~=cov_dim+res_dim)&(type==1),
Daniel@0:             error(['Invalid data matrix size: ', num2str(size(test_data,2)), ' columns rather than ',...
Daniel@0:                 num2str(cov_dim+res_dim), ' or ', num2str(cov_dim), '!']);
Daniel@0:         elseif (size(test_data,2)~=cov_dim)&(size(test_data,2)~=cov_dim+1)&(type==2),
Daniel@0:             error(['Invalid data matrix size: ', num2str(size(test_data,2)), ' columns rather than ',...
Daniel@0:                 num2str(cov_dim+1), ' or ', num2str(cov_dim), '!']);
Daniel@0:         elseif (~isempty(find(ismember(intersect([test_data(:,end)' 1:res_dim],...
Daniel@0:                 test_data(:,end)'),[1:res_dim])==0)))&(type==2)&(size(test_data,2)==cov_dim+1),
Daniel@0:             error('Invalid class label');
Daniel@0:         end
Daniel@0:         ntest=size(test_data,1);        
Daniel@0:         test_d=test_data(:,1:cov_dim);
Daniel@0:         if (type==2)&(size(test_data,2)>cov_dim),
Daniel@0:             test_t=zeros(ntest, res_dim);
Daniel@0:             for m=1:res_dim,
Daniel@0:                 test_t((find(test_data(:,end)==m))',m)=1;
Daniel@0:             end
Daniel@0:         elseif (type==1)&(size(test_data,2)>cov_dim),
Daniel@0:             test_t=test_data(:,cov_dim+1:end);
Daniel@0:         end
Daniel@0:         disp(' ');
Daniel@0:     end
Daniel@0: else    
Daniel@0:     clc
Daniel@0:     disp('----------------------------------------------------');
Daniel@0:     disp('                  Specify the Input                 ');
Daniel@0:     disp('----------------------------------------------------');
Daniel@0:     disp(' ')
Daniel@0:     ntrain = input('Insert the number of examples in training (<500): ');
Daniel@0:     if (isempty(ntrain)|(floor(ntrain)~=ntrain)|(ntrain<=0)|(ntrain>500)), 
Daniel@0:           error(['Invalid value: ', num2str(ntrain), ' is not a positive integer <500!']);
Daniel@0:     end        
Daniel@0:     disp(' ')
Daniel@0:     test_path='toy';
Daniel@0:     ntest = input('Insert the number of examples in test (<500): ');
Daniel@0:     if (isempty(ntest)|(floor(ntest)~=ntest)|(ntest<=0)|(ntest>500)), 
Daniel@0:           error(['Invalid value: ', num2str(ntest), ' is not a positive integer <500!']);
Daniel@0:     end        
Daniel@0: 
Daniel@0:     if type==2,
Daniel@0:         cov_dim=2;
Daniel@0:         res_dim=3;
Daniel@0:         seed = 42;
Daniel@0:         [train_d, ntrain1, ntrain2, train_t]=gen_data(ntrain, seed);
Daniel@0:         for m=1:ntrain
Daniel@0:             q=[]; q = find(train_t(m,:)==1);
Daniel@0:             train_data(m,:)=[train_d(m,:) q];
Daniel@0:         end
Daniel@0:         [test_d, ntest1, ntest2, test_t]=gen_data(ntest);
Daniel@0:         for m=1:ntest
Daniel@0:             q=[]; q = find(test_t(m,:)==1);
Daniel@0:             test_data(m,:)=[test_d(m,:) q];
Daniel@0:         end
Daniel@0:     else
Daniel@0:         cov_dim=1;
Daniel@0:         res_dim=1;
Daniel@0:         global HOME
Daniel@0:         %%%%%WARNING!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:         load([HOME '/examples/static/Misc/mixexp_data.txt'], '-ascii');
Daniel@0:         %%%%%WARNING!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:         train_data = mixexp_data(1:ntrain, :);
Daniel@0:         train_d=train_data(:,1:cov_dim); train_t=train_data(:,cov_dim+1:end);
Daniel@0:         test_data = mixexp_data(ntrain+1:ntrain+ntest, :);
Daniel@0:         test_d=test_data(:,1:cov_dim); 
Daniel@0:         if size(test_data,2)>cov_dim,
Daniel@0:             test_t=test_data(:,cov_dim+1:end);
Daniel@0:         end
Daniel@0:     end    
Daniel@0: end
Daniel@0: % Set the nodes dimension-----------------------------------
Daniel@0: if num_glevel>0,
Daniel@0:     nodes_info(2,2:num_glevel+1)=branch_fact;
Daniel@0: end
Daniel@0: nodes_info(2,1)=cov_dim; nodes_info(2,end)=res_dim;
Daniel@0: %-----------------------------------------------------------
Daniel@0: % Prepare the training data for the learning engine---------
Daniel@0: %-----------------------------------------------------------
Daniel@0: cases = cell(size(nodes_info,2), ntrain);
Daniel@0: for m=1:ntrain,
Daniel@0:     cases{1,m}=train_data(m,1:cov_dim)';
Daniel@0:     cases{end,m}=train_data(m,cov_dim+1:end)';
Daniel@0: end
Daniel@0: %-----------------------------------------------------------------------------------------------------
Daniel@0: [bnet onodes]=hme_topobuilder(nodes_info);
Daniel@0: engine = jtree_inf_engine(bnet, onodes);
Daniel@0: clc
Daniel@0: disp('---------------------------------------------------------------------');
Daniel@0: disp('                         L  E  A  R  N  I  N  G                      ');
Daniel@0: disp('---------------------------------------------------------------------');
Daniel@0: disp(' ')
Daniel@0: ll = 0;
Daniel@0: for l=1:ntrain
Daniel@0:   scritta=['example number: ', int2str(l),'---------------------------------------------'];
Daniel@0:   disp(scritta);
Daniel@0:   ev = cases(:,l);
Daniel@0:   [engine, loglik] = enter_evidence(engine, ev);
Daniel@0:   ll = ll + loglik;
Daniel@0: end
Daniel@0: disp(' ')
Daniel@0: disp(['Log-likelihood before learning: ', num2str(ll)]);
Daniel@0: disp(' ')
Daniel@0: disp('(Press any key to continue)');
Daniel@0: pause
Daniel@0: %-----------------------------------------------------------
Daniel@0: clc
Daniel@0: disp('---------------------------------------------------------------------');
Daniel@0: disp('                         L  E  A  R  N  I  N  G                      ');
Daniel@0: disp('---------------------------------------------------------------------');
Daniel@0: disp(' ')
Daniel@0: max_em_iter=input('Insert the maximum number of the EM algorithm iterations: ');
Daniel@0: if (isempty(max_em_iter)|(floor(max_em_iter)~=max_em_iter)|(max_em_iter<=1)), 
Daniel@0:           error(['Invalid value: ', num2str(ntest), ' is not a positive integer >1!']);
Daniel@0: end 
Daniel@0: disp(' ')
Daniel@0: disp(['Log-likelihood before learning: ', num2str(ll)]);
Daniel@0: disp(' ')
Daniel@0: 
Daniel@0: [bnet2, LL2] = learn_params_em(engine, cases, max_em_iter);
Daniel@0: disp(' ')
Daniel@0: fprintf('HME: loglik before learning %f, after %d iters %f\n', ll, length(LL2),  LL2(end));
Daniel@0: disp(' ')
Daniel@0: disp('(Press any key to continue)');
Daniel@0: pause
Daniel@0: %-----------------------------------------------------------------------------------
Daniel@0: % Classification problem: plot data & decision boundaries if the input data size = 2
Daniel@0: % Regression problem: plot data & prediction if the input data size = 1
Daniel@0: %-----------------------------------------------------------------------------------
Daniel@0: if (type==2)&(nodes_info(2,1)==2)&(~isempty(test_path)),
Daniel@0:     fh1=hme_class_plot(bnet2, nodes_info, train_data, test_data);
Daniel@0:     disp(' ');
Daniel@0:     disp('(See the figure)');
Daniel@0: elseif (type==2)&(nodes_info(2,1)==2)&(isempty(test_path)),
Daniel@0:     fh1=hme_class_plot(bnet2, nodes_info, train_data);
Daniel@0:     disp(' ');
Daniel@0:     disp('(See the figure)');
Daniel@0: elseif (type==1)&(nodes_info(2,1)==1)&(~isempty(test_path)),
Daniel@0:     fh1=hme_reg_plot(bnet2, nodes_info, train_data, test_data);
Daniel@0:     disp(' ');
Daniel@0:     disp('(See the figure)');
Daniel@0: elseif (type==1)&(nodes_info(2,1)==1)&(isempty(test_path)),
Daniel@0:     fh1=hme_reg_plot(bnet2, nodes_info, train_data);
Daniel@0:     disp(' ')
Daniel@0:     disp('(See the figure)');
Daniel@0: end
Daniel@0: %-----------------------------------------------------------------------------------
Daniel@0: % Classification problem: plot confusion matrix
Daniel@0: %-----------------------------------------------------------------------------------
Daniel@0: if (type==2)
Daniel@0:     ztrain=fhme(bnet2, nodes_info, train_d, size(train_d,1));  
Daniel@0:     [Htrain, trainRate]=confmat(ztrain, train_t); % CM on the training set
Daniel@0:     fh2=figure('Name','Confusion matrix', 'MenuBar', 'none', 'NumberTitle', 'off');
Daniel@0:     if (~isempty(test_path))&(size(test_data,2)>cov_dim),
Daniel@0:         ztest=fhme(bnet2, nodes_info, test_d, size(test_d,1));
Daniel@0:         [Htest, testRate]=confmat(ztest, test_t);   % CM on the test set
Daniel@0:         subplot(1,2,1);
Daniel@0:     end
Daniel@0:     plotmat(Htrain,'b','k',12)
Daniel@0:     tick=[0.5:1:(0.5+nodes_info(2,end)-1)];
Daniel@0:     set(gca,'XTick',tick)
Daniel@0:     set(gca,'YTick',tick)
Daniel@0:     grid('off')
Daniel@0:     ylabel('True')
Daniel@0:     xlabel('Prediction')
Daniel@0:     title(['Confusion Matrix: training set (' num2str(trainRate(1)) '%)'])
Daniel@0:     if (~isempty(test_path))&(size(test_data,2)>cov_dim),
Daniel@0:         subplot(1,2,2)
Daniel@0:         plotmat(Htest,'b','k',12)
Daniel@0:         set(gca,'XTick',tick)
Daniel@0:         set(gca,'YTick',tick)
Daniel@0:         grid('off')
Daniel@0:         ylabel('True')
Daniel@0:         xlabel('Prediction')
Daniel@0:         title(['Confusion Matrix: test set (' num2str(testRate(1)) '%)'])
Daniel@0:     end
Daniel@0:     disp(' ')
Daniel@0:     disp('(Press any key to continue)');
Daniel@0:     pause
Daniel@0: end
Daniel@0: %-----------------------------------------------------------------------------------
Daniel@0: % Regression & Classification problem: calculate the predictions & plot the LL trace
Daniel@0: %-----------------------------------------------------------------------------------
Daniel@0: train_result=fhme(bnet2,nodes_info,train_d,size(train_d,1));
Daniel@0: if ~isempty(test_path),
Daniel@0:     test_result=fhme(bnet2,nodes_info,test_d,size(test_d,1));
Daniel@0: end
Daniel@0: fh3=figure('Name','Log-likelihood trace', 'MenuBar', 'none', 'NumberTitle', 'off')
Daniel@0: plot(LL2,'-ro',...
Daniel@0:                 'MarkerEdgeColor','k',...
Daniel@0:                 'MarkerFaceColor',[1 1 0],...
Daniel@0:                 'MarkerSize',4)
Daniel@0: title('Log-likelihood trace')
Daniel@0: %-----------------------------------------------------------------------------------
Daniel@0: % Regression & Classification problem: save the predictions
Daniel@0: %-----------------------------------------------------------------------------------
Daniel@0: clc
Daniel@0: disp('------------------------------------------------------------------');
Daniel@0: disp('                           Save the results                       ');
Daniel@0: disp('------------------------------------------------------------------');
Daniel@0: disp(' ')
Daniel@0: %-----------------------------------------------------------------------------------
Daniel@0: save_quest_m=input('Do you want to save the HME model (Y/N)? [Y default]: ', 's');
Daniel@0: if isempty(save_quest_m),
Daniel@0:     save_quest_m='Y';
Daniel@0: end
Daniel@0: if ~findstr(save_quest_m, ['Y', 'N']), error('Invalid input'); end
Daniel@0: if save_quest_m=='Y',
Daniel@0:     disp(' ');
Daniel@0:     m_save=input('Insert the complete path for save the HME model (.mat):\n >> ', 's');
Daniel@0:     if isempty(m_save), error('You must specify a path!'); end
Daniel@0:     save(m_save, 'bnet2');
Daniel@0: end
Daniel@0: %-----------------------------------------------------------------------------------    
Daniel@0: disp(' ')
Daniel@0: save_quest=input('Do you want to save the HME predictions (Y/N)? [Y default]: ', 's');
Daniel@0: disp(' ')
Daniel@0: if isempty(save_quest),
Daniel@0:     save_quest='Y';
Daniel@0: end
Daniel@0: if ~findstr(save_quest, ['Y', 'N']), error('Invalid input'); end
Daniel@0: if save_quest=='Y',
Daniel@0:     tr_save=input('Insert the complete path for save the training data prediction (.mat):\n >> ', 's');    
Daniel@0:     if isempty(tr_save), error('You must specify a path!'); end
Daniel@0:     save(tr_save, 'train_result');  
Daniel@0:     if ~isempty(test_path),
Daniel@0:         disp(' ')
Daniel@0:         te_save=input('Insert the complete path for save the test data prediction (.mat):\n >> ', 's');
Daniel@0:         if isempty(te_save), error('You must specify a path!'); end
Daniel@0:         save(te_save, 'test_result');
Daniel@0:     end
Daniel@0: end
Daniel@0: clc
Daniel@0: disp('----------------------------------------------------');
Daniel@0: disp('                      B  Y  E !                     ');
Daniel@0: disp('----------------------------------------------------');
Daniel@0: pause(2)
Daniel@0: %clear 
Daniel@0: clc