function[] = DST_init_rec(cur,maxent,infer_only)

global D;

%chain
if( D(cur).type==1)   

    %use slightly perturbed max entropy initialization
    if(maxent)
        
        if(~infer_only)
            %no parent
            if(cur==1)
		
                D(cur).Phi = mk_stochastic(rand(D(cur).S,D(cur).S)/20 + .5)';
                D(cur).phi = normalise(rand(D(cur).S,1)/20 + .5);
     
            %has parent    
            else
      
                for i=1:D(D(cur).parent).S
                    D(cur).Phi(:,:,i) = mk_stochastic(rand(D(cur).S,D(cur).S)/20 + .5)';
                    D(cur).phi(:,i) = normalise(rand(D(cur).S,1)/20 + .5,1);
		        end
            end
        end
        
        for n=1:D(cur).N
        
            for t=1:D(cur).T(n)
	
                D(cur).s{n}(:,t) = normalise(rand(D(cur).S,1)/20 + .5,1);

            end
	    end
        
    %use random initialization    
    else
        if(~infer_only)
            %no parent
            if(cur==1)   
                D(cur).Phi = mk_stochastic(rand(D(cur).S,D(cur).S))';
                D(cur).phi = normalise(rand(D(cur).S,1));
            
            %has parent
            else
	
                for i=1:D(D(cur).parent).S
                    D(cur).Phi(:,:,i) = mk_stochastic(rand(D(cur).S,D(cur).S))';
                    D(cur).phi(:,i) = normalise(rand(D(cur).S,1));
                end
            end
        end
            
        for n=1:D(cur).N
        
            for t=1:D(cur).T(n)
	
                D(cur).s{n}(:,t) = normalise(rand(D(cur).S,1));

            end
	    end
        
    end
    
    for c=D(cur).child
        
        DST_init_rec(c,maxent,infer_only);
        
    end
    
elseif( D(cur).type==2)
    
    %need to write sometime
    
elseif( D(cur).type==3)
    
    if(~infer_only)
        xdim = D(cur).xdim;
        
        ydim = D(cur).ydim;
	
		D(cur).C = rand(ydim,xdim);
		
		D(cur).R = rand(ydim);
		D(cur).R = D(cur).R * D(cur).R' + 1e-8 * eye(ydim);
        
        A = rand(xdim);
        Q = rand(xdim);
        Q = Q * Q' + 1e-8 * eye(xdim);
        Qinv = safeinv(Q);
		Q1 = rand(xdim);
        Q1 = Q1 * Q1' + 1e-8 * eye(xdim);
        Q1inv = safeinv(Q1);
        mu= rand(xdim,1);
        
        if(cur==1)
            D(cur).A = A;
		
            D(cur).Q = Q;
		
            D(cur).Qinv = Qinv;
            
            D(cur).Q1 = Q1;
	
            D(cur).Q1inv = Q1inv;
		
            D(cur).mu = mu;
		
        else
            
            Dtemp(1).type=1;
            Dtemp(1).S=1;
            Dtemp(1).phi=1;
            Dtemp(1).Phi=1;
            Dtemp(1).T = D(cur).T;
            Dtemp(1).N = D(cur).N;
            Dtemp(1).child = 2;
            for n=1:Dtemp(1).N            
                Dtemp(1).s{n} = ones(1, Dtemp(1).T(n));
	%             Dtemp(1).xi{n} = ones(1, Dtemp(1).T(n)-1);
            end
            
            Dtemp(2).parent=1;
            Dtemp(2).type=3;
            Dtemp(2).A = A;
            Dtemp(2).Q = Q;   
            Dtemp(2).Qinv = Qinv;
            Dtemp(2).Q1 = Q1;
            Dtemp(2).Q1inv = Q1inv;
            Dtemp(2).mu = mu;
            Dtemp(2).y = D(cur).y; 
            Dtemp(2).xdim = D(cur).xdim;
            Dtemp(2).ydim = D(cur).ydim;
		    Dtemp(2).C = D(cur).C;
		    Dtemp(2).R = D(cur).R;
            Dtemp(2).T = D(cur).T;
            Dtemp(2).N = D(cur).N;
            Dtemp(2).missing = D(cur).missing;
                    
            [Dtemp, bound] = DST_em2(Dtemp, 10, 1, 0);
	
            D(cur).C = Dtemp(2).C;
            D(cur).R = Dtemp(2).R;
	
            %D(cur).R = symmat(D(cur).R);
            
            [spldata, time] = splitdata(D(cur).y,D(D(cur).parent).S);
            missing = splitdata((D(cur).missing),D(D(cur).parent).S,1);
      
	
			for i=1:D(D(cur).parent).S
	
                D(cur).mu(:,i) = Dtemp(2).mu;
                D(cur).Q1(:,:,i) = Dtemp(2).Q1;
                D(cur).Q1inv(:,:,i) = Dtemp(2).Q1inv;
            
                Dtemp(2).y = {spldata{i,:}};
                Dtemp(2).missing = {missing{i,:}};
     
                Dtemp(2).T = time{i};
                Dtemp(1).T = time{i};
                
                [Dtemp2, bound] = DST_em2(Dtemp, 10, 1, 0);
                
                D(cur).A(:,:,i) = Dtemp2(2).A;
                D(cur).Q(:,:,i) = Dtemp2(2).Q;
                D(cur).Qinv(:,:,i) = Dtemp2(2).Qinv;
	%             D(cur).Q1(:,:,i) = Dtemp2(2).Q1;
	%             D(cur).Q1inv(:,:,i) = Dtemp2(2).Q1inv;
	%             D(cur).mu(:,:,i) = Dtemp2(2).mu;
	
                
	%             [D(cur).A(:,:,i), C, D(cur).Q(:,:,i), R, D(cur).mu(:,i), D(cur).Q1(:,:,i), LL] = ...
	%     learn_kalman(trans_cell({spldata{i,:}}), A, D(cur).C, Q, D(cur).R, mu, Q1, 10);
	
	%             D(cur).Q(:,:,i) = symmat(D(cur).Q(:,:,i));
	% 
	%             D(cur).Qinv(:,:,i) = symmat(safeinv(D(cur).Q(:,:,i)));
	%    
	%             D(cur).Q1(:,:,i) = symmat(D(cur).Q1(:,:,i));
	%             
	%             D(cur).Q1inv(:,:,i) = symmat(safeinv(D(cur).Q1(:,:,i)));
                
			end
        end
    end
    
end