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