Daniel@0: 
Daniel@0: function [eg] = dag_to_essential_graph(dag)
Daniel@0: cpdag = dag_to_cpdag(dag);
Daniel@0: eg = dag + dag .* (cpdag + cpdag');
Daniel@0: 
Daniel@0: return;
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: % Coverts a DAG into Essential Graph where edges are coded by 2 and 3, 2 is
Daniel@0: % directed edge and 3 is bidirected edge and is at one (the same as the original DAG) of the two
Daniel@0: % symetrical places. 
Daniel@0: 
Daniel@0: % Is implemented by the algorithm of Max Chickering in D.M.Chickering (1995). 
Daniel@0: % A transformational characterization of equivalent Bayesian network structures. 
Daniel@0: % In Proceedings of Eleventh Conference on Uncertainty in Artificial Intelligence, Montreal, QU,
Daniel@0: % pages 87-98. Morgan Kaufmann 
Daniel@0: % http://research.microsoft.com/~dmax/publications/uai95.pdf 
Daniel@0: 
Daniel@0: % Implemented by Tomas Kocka, AAU.
Daniel@0: 
Daniel@0: function [eg] = dag_to_essential_graph(dagx)
Daniel@0: 
Daniel@0: %print_dag(dagx); % Just checking input
Daniel@0: 
Daniel@0: order = topological_sort(dagx); % get the topological order of nodes and their number
Daniel@0: 
Daniel@0: % fprintf('the topological order is: %d',order);
Daniel@0: % fprintf('\n');
Daniel@0: 
Daniel@0: [nx,ny] = size(dagx); % gets the number of nodes, note that nx == ny
Daniel@0: [I,J] = find(dagx); % finds all nonzero elements in the adjacency matrix, i.e. arcs in the DAG - however we will overwrite it in a special order
Daniel@0: % we will sort the arcs from lowest possible y and highest possible x, arcs are x->y
Daniel@0: e = 1;
Daniel@0: for y = 1:ny
Daniel@0:     for x = nx:-1:1
Daniel@0:         %fprintf('x %d ',order(x)); fprintf('y %d ',order(y));
Daniel@0:         if dagx(order(x),order(y)) == 1 
Daniel@0:             I(e) = order(x);
Daniel@0:             J(e) = order(y);
Daniel@0:             e = e + 1;
Daniel@0:             %fprintf('x order %d',x);
Daniel@0:             %fprintf('y order %d',y);
Daniel@0:             %fprintf('\n');
Daniel@0:         end
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: % fprintf('the arcs are: %d',I);
Daniel@0: % fprintf('\n');
Daniel@0: % fprintf('the arcs are: %d',J);
Daniel@0: % fprintf('\n');
Daniel@0: 
Daniel@0: 
Daniel@0: % Now we have to decide which arcs are part of the essential graph and
Daniel@0: % which are undirected edges in the essential graph.
Daniel@0: % Undecided arc in the DAG are 1, directed in EG are 2 and undirected in EG
Daniel@0: % are 3.
Daniel@0: 
Daniel@0: 
Daniel@0: for e = 1:length(I)
Daniel@0:     if dagx(I(e),J(e)) == 1
Daniel@0:         cont = true;
Daniel@0:         for w = 1:nx 
Daniel@0:             if dagx(w,I(e)) == 2
Daniel@0:                 if dagx(w,J(e)) ~= 0
Daniel@0:                     dagx(w,J(e)) = 2;
Daniel@0:                 else
Daniel@0:                     for ww = 1:nx
Daniel@0:                         if dagx(ww,J(e)) ~= 0
Daniel@0:                            dagx(ww,J(e)) = 2;
Daniel@0:                         end
Daniel@0:                     end % and now skip the rest and start with another arc from the list
Daniel@0:                     w = nx;
Daniel@0:                     cont = false;
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if cont
Daniel@0:            exists = false;
Daniel@0:            for z = 1:nx
Daniel@0:                %fprintf('test %d',dagx(z,J(e)));
Daniel@0:                if dagx(z,J(e)) ~= 0 & z ~= I(e) & dagx(z,I(e)) == 0
Daniel@0:                   exists = true; 
Daniel@0:                   for ww = 1:nx
Daniel@0:                         if dagx(ww,J(e)) == 1
Daniel@0:                            dagx(ww,J(e)) = 2;
Daniel@0:                         end 
Daniel@0:                   end
Daniel@0:                end
Daniel@0:            end
Daniel@0:            if ~ exists
Daniel@0:                for ww = 1:nx
Daniel@0:                    if dagx(ww,J(e)) == 1
Daniel@0:                       dagx(ww,J(e)) = 3;
Daniel@0:                    end 
Daniel@0:                end  
Daniel@0:            end
Daniel@0:         end
Daniel@0:     end            
Daniel@0: end
Daniel@0: 
Daniel@0: %print_dag(dagx); % Just checking output
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: