function [bnet,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,...
    names,onames,order,inv_o]=...
   construct_bnet_hmm(B,S,D,M,S_item,data)
%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 columns 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 rows 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: onames=names(order)
%inv_o: vector: gives for the unordered node 1...end the corresponding nr of
    %the ordered nodes names=onames(inv_o)

%

%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=[];


disp('it is assumed that the data have a ');
disp(' separate row for each case, and');
disp('the columns are variables organised as follows:');
disp('first all items of first measurement occasion on first day');
disp('then all items of second measurement occasion of first day, etc');
disp('pres a key to continue');
pause;
T=B*D; %total number of timepoints

names=cell(2*T,1);

for i=1:T
    names{i}=eval(['''Xsignal_',int2str(i),'''']);
end

for i=1:T
    
        names{T+i}=eval(['''Y_',int2str(i),'''']);
        
end

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(M*T+(D-1)+D*(B-1),2);
connections=cell(T+D*(B-1),2);
for j=1:D
    for i=1:B-1
        connections{(j-1)*(B-1)+i,1}=eval(['''Xsignal_',int2str((j-1)*(B)+i),'''']);
        connections{(j-1)*(B-1)+i,2}=eval(['''Xsignal_',int2str((j-1)*(B)+i+1),'''']);
        
    end
end

for i=1:T
        connections{D*(B-1)+i,1}=eval(['''Xsignal_',int2str(i),'''']);
        connections{D*(B-1)+i,2}=eval(['''Y_',int2str(i),'''']);
end
    % 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, each of the T Ys contains M observables
    unord_merged_nodes=cell(T,1);
    unord_merged_nodes_respcat=cell(T,1);
    unord_nr_merged_nodes=M*ones(T,1);
    for i=1:T
        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=mysetdiff(1:nr_nodes,terminal_merged_nodes.nodenrs);
hid_nodes.nodesizes=S*ones(1,length(hid_nodes.nodenrs));
[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:T 
    terminal_merged_nodes.data{i}=(data(:,merged_nodes{i}))';
end

unord_nr_merged_nodes=M*ones(T,1);
nr_merged_nodes(ind)= unord_nr_merged_nodes; 
terminal_merged_nodes.nrvars=nr_merged_nodes;
terminal_merged_nodes.nodesizes=ones(T,1);






%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(hid_nodes.nodenrs)=S;

%murphy's
bnet = mk_bnet(adj, node_sizes, 'discrete', discrete_nodes, 'names','onames');