function design=...
    construct_design_mats(parents,node_sizes,equiv_class,lin_pred_struct,pred_mat,gausskwadnodes)
%design=construct_design constructs design matrix, includes gauss quad
%nodes
%columns (and parms) are ordered according to order of nodenr: last column
%corresponds to highest nodenumber (children). for interactions: lowest nodenr of
%interactions counts
%rows of design are ordered as in sas proc mixed: highest nodenrs (children) change
%first
%this is all quite complex and might not work bugfree all of the time
%checks are thus needed

K=max(equiv_class);
design=cell(K,1);
for k=1:K
    node_nr=find(equiv_class==k,1);
    par=parents{node_nr};
    siz=node_sizes(par);
    des=pred_mat{k};    
%incorporate effects for parents
        if ~isempty(par) 
          term=lin_pred_struct{k};%term is cell or vector, see define_lin_pred_... functions
          if ~iscell(term) %term is vector, one type of effects
             if term(1)==0 %interactions with variables in pred_mat
                repfactor=prod(siz(term(2:end)));%size of state space of the set of parents 
                dc=eye(repfactor);                  %starts from second element, first element of term is 0
                des= kron(dc,des);   
             else %no interactions with covariates
                repfactor=prod(siz(term));%starts from first element, first element of term is nonzero
                dc=dummycode(0:repfactor-1);
                des2= kron(dc,ones(size(des,1),1));
                des= [des2 repmat(des,repfactor,1)];
             end
          else %term is cell, several sets of effects, tricky
            
             repfactor2=1;
             check_order(term);%check wheter lin_pred_struc is correctly defined
             first=1;
             for j=1:length(term)
                 term2=term{j};%first cell element lowest nodes ('oldest'parents)
                 
                 if term2(1)==0 
                    %repfactor1=repfactor1+prod(siz(term2(2:end)));
                    repfactor=prod(siz(term2(2:end)));
                    dc=eye(repfactor);
                    des= kron(dc,des); 
                    dd=des;
                 else
                    repfactor=prod(siz(term2));
                    repfactor2=repfactor2*repfactor;
                    if ~mysubset(par(term2),gausskwadnodes.nodenrs) 
                        dc=dummycode(0:repfactor-1);
                    else dc=herzo(siz(term2));
                    end
                    if ~first
                        dd=[kron(dd,ones(repfactor,1)) repmat(dc,size(dd,1),1)];
                    else
                        dd=dc;
                        first=0;
                    end
                    %des2= kron(dc,ones(size(des,1),1));
                    %des= [des2 repmat(des,repfactor,1)];
                 end
              end
             %dc=eye(repfactor1);
             %des= kron(dc,des);
             if repfactor2>0 
                 %dc=dummycode(0:repfactor2-1);
                 des2= kron(dd,ones(size(des,1),1));
                 des= [des2 repmat(des,repfactor2,1)];
             
             end
             
          end
        end
  design{k}=des;    
        
end