function [bnet,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,names,onames,order,inv_o]=...
    construct_bnet_mixlltm(S,nrquadpoints,M,S_item,data)
%construct bnet for mixture rasch model
%S=nr classes
%bnet: bnet constructed from bnet software of murphy (struct array)
%evidence_nodes: struct(4): NOT YET IMPLEMENTED: ARE TREATED AS PARTIALLY
    %OBSERVED EVIDENCE NODES
    %evidence_nodes.nodenrs vector of (ordered) nodenrs of nodes that are observed for all cases
    %evidence_nodes.respcat nr of resp categories (not equal to nodesizes,
    %%evidence_nodes.nodesizes 
    %evidence_nodes.data cell array of corresponding columns in dataset  
%partial_evidence_nodes: struct(3): 
    %partial_evidence_nodes.nodenrs vector of (ordered) nodenrs of nodes that are partially observed 
    %partial_evidence_nodes.nodesizes number of resp categories (equals
    %nodesizes in the graph
    %partial_evidence_nodes.data cell array of of corresponding rows in dataset  
%terminal_merged_nodes: struct(5): 
    %terminal_merged_nodes.nodenrs vector of (ordered) nodenrs of terminal (partially) observed nodes 
    %terminal_merged_nodes.nodesizes
    %terminal_merged_nodes.data cell array of corresponding sets of columns in dataset 
    %that are merged into single node 
    %terminal_merged_nodes.nrvars vector of numbers of nodes merged
    %into each (ordered) terminal node
    %terminal_merged_nodes.respcat cell array of nr of respcat of merged
    %nodes
%hid_nodes struct(2)
    %hid_nodes.nodenrs vector of (ordered) nodenrs of nodes that are
    %hidden
    %hid_nodes.nodesizes
%gausskwadnodes: struct(2)
    %gausskwadnodes.nodenrs (ordered) nrs continuous normal hidden nodes that are discretised
 
    %gausskwadnodes.nodesizes

%the following output is not needed but can be useful during debugging
%names: names of the original order of the nodes
%onames: names for the reordered nodes (reordering is done during the
    %construction of the graphical model)
%order: vector: gives for the ordered node 1...end the corresponding nr of
    %the unordered nodes: names=onames(order)
%inv_o: vector: gives for the unordered node 1...end the corresponding nr of
    %the ordered nodes onames=names(order)

%

%default values
evidence_nodes.nodenrs =[];
evidence_nodes.nodesizes =[];
evidence_nodes.data ={};  
evidence_nodes.respcat =[];
partial_evidence_nodes.nodenrs =[];
partial_evidence_nodes.nodesizes =[];
partial_evidence_nodes.data ={}; 
terminal_merged_nodes.nodenrs =[];
terminal_merged_nodes.nodesizes =[];
terminal_merged_nodes.data ={};
terminal_merged_nodes.nrvars =[];
terminal_merged_nodes.respcat ={};
hid_nodes.nodenrs=[];
hid_nodes.nodesizes=[];
gausskwadnodes.nodenrs=[];
gausskwadnodes.nodesizes=[];

names=cell(3,1);
names{1}='class';
names{2}='theta';
names{3}='Y';

nr_nodes=length(names);
%We can specify the topology using a cell array as follows, where each row specifies a connection between two named nodes: 

connections=cell(2,2);
connections{1,1}='class';
connections{1,2}='Y';
connections{2,1}='theta';
connections{2,2}='Y';
    % convert this cell array to an adjacency matrix using the following function: 
[adj ,onames,order ] = franks_mk_adj_mat(connections, names,1);

inv_o=inv_order(order);
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%alternatively, specify the adjacency matrix directly 
%%%%%%%%%%%%%%%%%%%%%%%%%%%


%define terminal_merged_nodes
    %specify which terminal observed nodes are merged 
    %for unordered nodes, 
    unord_merged_nodes=cell(1,1);
    unord_merged_nodes_respcat=cell(1,1);
    unord_nr_merged_nodes=M*ones(1,1);
    for i=1:1
        unord_merged_nodes{i}=(i-1)*M+[1:M];
        unord_merged_nodes_respcat{i}=S_item*ones(M,1);
    end
terminal_merged_nodes.nodenrs=strmatch('Y', onames);
hid_nodes.nodenrs=strmatch('class', onames);
[y,ind]=sort(order(terminal_merged_nodes.nodenrs));
terminal_merged_nodes.respcat(ind)=unord_merged_nodes_respcat;
merged_nodes(ind)=unord_merged_nodes;
for i=1:1
    terminal_merged_nodes.data{i}=(data(:,merged_nodes{i})');
end

unord_nr_merged_nodes=M*ones(1,1);
nr_merged_nodes(ind)= unord_nr_merged_nodes; 
terminal_merged_nodes.nrvars=nr_merged_nodes;
terminal_merged_nodes.nodesizes=ones(1,1);
gausskwadnodes.nodenrs=strmatch('theta', onames);
gausskwadnodes.nodesizes=nrquadpoints;



%define bnet
discrete_nodes = 1:nr_nodes;
node_sizes(terminal_merged_nodes.nodenrs) = 1;%for HMM:type models: either they are observed, and actual node size is then one
%or they are unobserved and thus can be left out because terminal hidden nodes have no influence
node_sizes(strmatch('class', onames)) = S;
node_sizes(strmatch('theta', onames)) = nrquadpoints;
hid_nodes.nodesizes=S;
%murphy's
bnet = mk_bnet(adj, node_sizes, 'discrete', discrete_nodes, 'names','onames');