%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.datamatrix=xlsread('braindata');
    
    %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.matrix=xlsread('conditie'); %function to read the time covariate of Rijmen, vansteelandt, De Boeck
        %covariates.names={'treatment_A','treatment_C'};
        

        
%%%%%       
%2.define bayes net (model structure)
%%%%%

    %specify bayes net using one of the construct_bnet functions or use your own
        %here: rijmen, ip et al paper on brain cancer patients, with random effect within states 
    %specifications of options for the construct_bnet function that is used. 
    S=4;% number of latent states 
    D=5; %number of measurements
    M=13; %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_hmm(D,S,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
    nr_nodes=length(onames);
    equiv_class=zeros(nr_nodes,1);
    equiv_class(strmatch('timepoint_1',onames))=1;
    equiv_class(mysetdiff(strmatch('timepoint',onames),strmatch('timepoint_1',onames)))=2;
    equiv_class(strmatch('Y',onames))=3;
     
    param_equiv_class=equiv_class;

    nr_equiv=max(equiv_class);
    for i=1:nr_equiv
        link{i}='multinomial';
    end
    
    %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 
            %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 which is stored in file
            design=xlsread('domaindesign');
            pred_mat(eclass)={design};
        
        %second, define the structure of the linear predictor for each equiv class (no restrictions, only main effects for parents etc)    
        %saturated: no restrictions
        lin_pred_struct=define_lin_pred_struct_sat(equiv_class,parents);
        
        
   
        %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);
        %further changes to design matrices can be done manually or
        %with code provided by user
        %here item main effects have to be included, so we need to add
        %eye(M) for design{3} and design{4}
        %domain effects of first class have to be put to zero for
        %identifiability (skip columns 2:4)
        ss=size(design{3},1)/M;
        design{3}(:,2:4)=[];
        design{3}=[design{3} repmat(eye(M),ss,1)];
        
        %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
            %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   
            %in rijmen, ip et al. only main efffects of covariates
            %lin_pred_struct_cov = define_lin_pred_struct_cov_main(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
            
  %define restrictions on parameters (gauss kwad nodes are
  %restricted automatically in gen_rand_start function below)
  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>.00000001 && logl_diff/logl_avg>.00000001 )
            %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