%main file for estimation prob graphical models
%here, the hierarchical hidden markov model with time interval as a covariate for transitions between signals is taken as an example
%see the examples in the example map for other specific models
%or define your own model
%to be done for each model
    %1.load data 
    %2.define bayes net (model structure)
    %3.define equivalence sets of nodes, and link function and design
        %matrix for each equivalent set
    %4. estimation, automatic after speciofication of 1-3
    %5. postestimation: compute standard errors, save results 

%%%%%
%1.load data
%%%%%
    %obs_variables: the observed variables that are nodes in the
    %graphical model
    %obs_variables is structure :
        %obs_variables.names: cell array of names (does not have to be
        %defined)
        %obs_variables.datamatrix: 'wide' format datamatrix, there is a
            %separate row for each independent case (e.g. two repeated
            %measurements defines two separate variables)
            %missings are coded by the value -1
        
    %matlab import wizard can be used to import data from e.g. excel
    obs_variables=loadtimesamplingdata; %function to read the time sampling data of Rijmen, vansteelandt, De Boeck
    
    %covariates: the covariates, variables that vary over cases
    %but that are themselves not included as separate nodes in the network
    %covariates.names: cell array of names
    %covariates.matrix: again 'wide' format.
    %covariates=struct([]); %default
    covariates=loadtime; %function to read the time covariate of Rijmen, vansteelandt, De Boeck
        
%%%%%       
%2.define bayes net (model structure)
%%%%%

    %specify bayes net using one of the construct_bnet functions or use your own
        %here: the hierarchical latent
        %markov model with multiple indicators of rijmen et al, 2005
    %specifications of options for the construct_bnet function that is used. 
    S_signal=2;% number of latent states at signal level
    S_day=2; %number of latent states at day-level
    B=9; %number of signals in day
    D=7; %number of days
    M=12; %number of items at each measurement occasion
    S_item=2;%number of response categories for items
    N=size(obs_variables.datamatrix,1); %number of cases
    %bnet construction
    [bnet,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,...
    names,onames,order,inv_o]=...
    construct_bnet_hier_hmm(B,S_day,S_signal,D,M,S_item,obs_variables.datamatrix);

    %convert bnet object from BNT (Murphy) into what we need 
    [parents,child,node_sizes,postorder,preorder,cliquetable,septable,clq_ass_to_node,pot_to_CPT]= franks_from_BNT(bnet,size(obs_variables.datamatrix,1));
    

%%%%
%3.define equivalence sets of nodes, link function,  design
    %matrix and restrictions on parameters for each equivalent (parameter) set
%%%%
    %equivalent nodes are nodes with the same design matrix, the same link and governed by the same set of parameters 
    %param_equivalent nodes are nodes with a differnetn design matrix, but
    %the same link and governed by the same set of parameters
    [equiv_class,param_equiv_class,link]=equiv_classes_hier_hmm(onames,B,D);
    
    %design matrix for all equiv classes
        %first, define pred_mat which is a cell array. each cell contains
        %the design matrix of the
        %variable(s) in a node without modelling the parents. e.g. 
        %for a four-category variable, a 3*3 design matrix
        %case covariates are not yet included here 
        %default is an identity matrix
        pred_mat=construct_predmat(equiv_class,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes);
        % item design matrices can be incorporated here as follows (here i use an
        % identity matrix as design matrix because i do not want restrictions)
            %first find the class of equivalent nodes that contain the observed variables
            %n=strmatch('Y',onames);
            %n=n(1);
            %eclass=equiv_class(n);
            %then input design matrix
            %design=eye(M);
            %pred_mat{eclass}=design;
        
        %second, define the structure of the linear predictor for each node (no restrictions, only main effects for parents etc)    
        lin_pred_struct=define_lin_pred_struct_hier_hmm_main(equiv_class,parents,terminal_merged_nodes.nodenrs);
        %third, construct design matrix for each node based on its pred_mat and  linear
        %predictor
        % no covariates yet
        design=construct_design_mats(parents,node_sizes,equiv_class,lin_pred_struct,pred_mat,gausskwadnodes);
        
        %fourth, inclusion of case covariates, makes design three dimensional, last
        %dimension represents cases

            %specify which covariates belong to which equivalent sets of
            %nodes
            %cov_nodes is cell array that contains as rows a pair of
            %equiv_class number,column numbers of the covariate matrix
            cov_nodes={};%default
            %cov_nodes=link_covariates_to_nodes_hier_hmm_time(equiv_class,covariates.names,onames,B,D);
            %specify main and or interaction effects of the
            %covariates
            %default: interactions at highest level (covariates have different effect 
            %for each category of the dependent variable and combination of the parents    
            %lin_pred_struct_cov=define_lin_pred_struct_cov_default(cov_nodes);
            %include covariates in design matrices
            %design=cov_into_design(covariates.matrix,cov_nodes,lin_pred_struct_cov,design,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,equiv_class);
            
       %fifth, further changes to design matrices can be done manually or
       %with code provided by user
       %example: the model in rijmen et al, where there are main
       %effects of items, and interactions between positive vs
       %negative emotions and latent states, design{4} is the appropriate design matrix
       %first four pos emotions, then 8 neg emotions
        aa=[ [ones(4,1);zeros(8,1)] [zeros(4,1); ones(8,1)]];
        aa=[kron([0 0 1 1]',aa) kron([0 1 0 1]',aa)];
        design{4}=[aa design{4}(:,3:end)];
   
   %define restrictions on parameters
   restparms=cell(max(param_equiv_class),1);
   %default =0, no restrictions (1 =restricted at fixed value)
   for i=1:max(param_equiv_class);
       node_nr=find(param_equiv_class==i,1);
       restparms{i}=zeros(size(design{equiv_class(node_nr)},2),1);
   end
   
%%%%
%4. estimation, automatic after specification of 1-3
%%%%
    %estimation with EM
    %can be altered: number of runs, starting values, convergence criteria

    %number of runs
    nr=1; 
    for run=1:nr
        loglikelihood=-1e10;
        v=Inf;
        it=1;
        %starting values   (here: random start)
        [parms,restparms]=gen_random_start(design,parents,equiv_class,param_equiv_class,gausskwadnodes,link);
        loglikelihood=-1/eps;
        logl_diff=1;
        logl_avg=1;
        %EM loop
        while (v>.000000001 && logl_diff/logl_avg>.000000001 )
            %1 iteration
            parms_old=parms;
            loglikelihood_old=loglikelihood;
            [parms, restparms, loglikelihood ]=EM_iteration(link,parms,restparms,design,parents,node_sizes...
            ,equiv_class,param_equiv_class,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,...
            gausskwadnodes,preorder, postorder ,cliquetable, septable,pot_to_CPT,N);    
            %check convergence
            v=0;
            for i=1:max(param_equiv_class)
                vvec=abs(parms_old{i}-parms{i});
                v=[v max(vvec)];
                v=max(v);
            end
            v
            it=it+1;
            logl_diff=loglikelihood-loglikelihood_old;
            if logl_diff<-.0001
                disp('decrease in loglikelihood. press any key to continue');
                pause;
            end
            logl_avg=(abs(loglikelihood)+abs(loglikelihood_old))/2
        end
%%%%
%5. postestimation: compute standard errors, save results ?
%%%%

        %compute final CPTs?
    	lin_pred=construct_lin_pred(parms,design,parents,node_sizes,equiv_class,param_equiv_class,terminal_merged_nodes,N);
        equiv_class_CPTs=construct_equiv_class_CPT(link,lin_pred,N);
    
        %save?
        %eval(['save D:\f\mijnmatlabfuncties\bayesnetfullgraph\kristof\test\beepitem_main_posneg2_tijd_juist_run',int2str(run),'_personstate',int2str(S_person),'_beepstate',int2str(S_beep),' parms loglikelihood equiv_class_CPTs ']);

        %compute standard errors?
        %info=num_infomatrix_anal_score(link,parms,restparms,design,parents,node_sizes,equiv_class,...
        %postorder,septable, cliquetable,preorder,param_equiv_class,clq_ass_to_node,evidence_nodes,...
        %partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,N,.000001);
    end