list={'alarm50','alarm100','alarm200','alarm500','alarm50mis','alarm100mis','alarm200mis','alarm500mis','alarm25000'};
%main file for estimation prob graphical models for artificial intelligence
%paper

for s=1%:length(list)-1
    set=list{s}
%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
%%%%%
    
    eval(['load c:\frank\alarm\twee\',set,'.txt -ascii;'])
    
    %load alarm dag
    load c:\frank\alarm\dag;
    dag=data;
    order = topological_sort(dag);
    adj=dag(order,order);
    inv_o=inv_order(order);
    [y,i]=sort(labels);
    restalf=i(1:2:end);
    rest=inv_o(restalf);
    nn=[50 100 200 500];
    restnodes=[29 37 12 16 1 3  35 33 9 19 10 23 5 25 4 21 30 22  17];
    names=labels;
    eval(['data= ',set,''';']);
    N=size(data,2);
    %number of added prior cases to each freq table
    prior_count=N/1000;
          
%%%%%       
%2.define bayes net (model structure)
%%%%%


    [bnet,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,...
    names,onames,order,inv_o]=...
    construct_alarm(dag,data,names);

    %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,N);
    

%%%%
%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=1:length(names);
    param_equiv_class=equiv_class;
   for i=1:length(design)
    link{i}='multinomial';
   end
    for i=1:length(restnodes)
    link{restnodes(i)}='cumulative';
    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
        %default is an identity matrix
        pred_mat=construct_predmat(equiv_class,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes);

        %second, define the structure of the linear predictor for each node (no restrictions, only main effects for parents etc)    
        %sa: no restrictions (never use default with gaussian quadrature nodes,they only have main effects)    
        lin_pred_struct=define_lin_pred_struct_sat(equiv_class,parents);
        %only main effects for rest nodes
        for i=1:length(restnodes)
            j=restnodes(i);
            a=cell(1,1+length(parents{j}));
            for k=0:length(parents{j})
                a(k+1)={k};
            end
            lin_pred_struct(j)={a};
        end
        design=construct_design_mats(parents,node_sizes,equiv_class,lin_pred_struct,pred_mat,gausskwadnodes);
            
        %no 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   
        [parms,restparms]=gen_alarm_start(design,parents,equiv_class,param_equiv_class,gausskwadnodes,link);
        loglikelihood=-1/eps;
        logl_diff=1;
        logl_avg=1;
        %EM loop
        it=1
        %while (v>.000000001 && logl_diff/logl_avg>.000000001 )
        while (v>.00001 || logl_diff/logl_avg>.00001 )
            %1 iteration
            parms_old=parms;
            loglikelihood_old=loglikelihood;
            [parms, restparms, loglikelihood ]=EM_iteration_prior(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,prior_count,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
            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
            it=it+1
        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 c:\frank\alarm\twee\main_cum50_50_',set,' N design parms loglikelihood equiv_class_CPTs parents;']);

        %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
end