function [prev_state_index, input_index, trellis_out, A] = siva_init ( next_state, output, N_IN , N_OUT )

% Vectorized Soft-Input Viterbi Algorithm Initialization
% [prev_state_index, input_index, trellis_out, A] = siva_init ( next_state, output, N_IN , N_OUT )
%
% prev_state_index : N_STATE by num_trans_max matrix. prev_state_index( i , : ) indicates the indexes of states whose next
%                    state transition is i. If num_trans(i) < num_trans_max, prev_state_index( i, num_trans(i)+ 1 : num_trans_max ) = 1;
%                    
% input_index : N_STATE by num_trans_max matrix which indicates the associated inputs with the state transition
%                       
%
% trellis_out : N_STATE*num_trans_max by N_OUT matrix which indicates the asociated trellis outputs in bipolar forms.
%
%                 [ O(1,1) O(1,2) .. O(1,N_OUT)   ----------------
%                 O(1,1) O(1,2) .. O(1,N_OUT)   |               |
%                                 ..            |  num_trans(1) |
%                 O(1,1) O(1,2) .. O(1,N_OUT)   ---             | num_trans_max
%                   0      0         0                          | 
%                                 ..                            |
%                   0      0         0          ------------------
%                 O(2,1) O(2,2) .. O(2,N_OUT)   ---
%                 O(2,1) O(2,2) .. O(2,N_OUT)   |
%                               ..              |  num_trans(2)
%                 O(2,1) O(2,2) .. O(2,N_OUT)   ---
%                               ..
%                 O(N_STATE,1)  .. O(N_STATE,N_OUT) ]
%
%               Output associated with the state tansition
%               
%  A : Additional matrix to be added to make irrelevant cost infinite
%     A returns 0 for the simple case of  num_trans = constant
%
% Sep 3, 2002
%

N_STATE=size(next_state,1);     %Number of states

if prod(double(size(next_state)==size(output)))==0
    error('Dimensions of next_state and output should be the same');
    return;
end

for i = 1 : N_STATE
    
    [ row_index , col_index ] = find( next_state == i - 1 );
    num_trans( i ) = length(col_index ) ;
    
end

num_trans_max = max(num_trans);

trellis_out = zeros( N_STATE * num_trans_max , N_OUT );
prev_state_index = zeros( N_STATE , num_trans_max );

for i=1:N_STATE

    [ row_index ,col_index] = find( next_state == i-1 );
    
    K = num_trans(i);
    
    for j=1:N_OUT
        trellis_out( (i-1)*K +1 : i*K , j ) = 2*bitget( diag(output( row_index, col_index)), j)-1;
    end
    
    prev_state_index ( i , 1 : K ) = row_index';
    input_index ( i , 1 : K)  = col_index';
    
end

[ I , J ] = find(prev_state_index == 0 );

if isempty(I)
    %Simple case : the number of transition to each state is the same
    
    A = 0;
    
else
    
    BRANCH_MAX = 1e5;
    
    A = zeros( N_STATE , num_trans_max ); % Additional matrix to be added to make cost infinite
    
    for i = 1 : length(I)
    
        prev_state_index(I(i),J(i)) = 1;
        B( I(i) , J(i) ) = BRANCH_MAX;
        
    end
    
end

num_trans_max = max(num_trans);
total_trans = sum (num_trans);
row_trellis_out = size ( trellis_out , 1);
cum_num_trans = cumsum (num_trans) - num_trans;