% Gauss-Markov full-trellis vs. RSSE
%
%
%%%
%%
% global variables for Viterbi
%
global HUGE_VAL; % 1e99
global DEBUG;    % for debugging
global V;        % main state array
global Vl;       % array for L-best indices
global Vw;       % cost array
global Vparam;   % parameters
global Vprev;    % previous state indices
global Vlprev;   % previous state el indices
global VW;       % final path costs

HUGE_VAL = 1e99;

DEBUG = 0;
%
%%
%%%

%%%
%%
% global variables for the cost functions
%
% use big GG_tmp and qq_tmp arrays or not
%
global USE_GGqq_tmp; USE_GGqq_tmp = 1
global G;
global q;

global z;
global trans;
global states;

global FIR_length;
global Memory;
global FIR_0;
global FIR_save;

global Nchips;
global Nusers;
global Nsyms;

global TYPE;     % problem type

global alpha;
global alpha2;
global var_v;
global stdev_v;
global C_inv;

global var_noise;
global SNR;

global GG_tmp;
global qq_tmp;
global GG;
global qq;

global ex_known;
global ex_gm;
global ex_p;
%
%%
%%%

format short g;

% initialize rand() and randn() with time-of-day
%
  seed=floor(sum(100*clock))
  rand( 'state', seed);
  randn( 'state', seed);

%%%
%%
% input data setup
%
Nusers = 2
Nchips = 4
FIR_length = 3

Nsyms = 5 % number of received data samples (symbols) per user

alpha = 0.995
alpha2 = alpha*alpha;

var_v = 0.01
stdev_v = sqrt( var_v);

Memory = FIR_length - 1;     % channel memory size

% full inv(C) matrix would have FIR_length*Nusers*Nchips rows and columns
% but here we are assuming that it is block diagonal with these blocks:
%
C_inv = eye( FIR_length*Nusers, FIR_length*Nusers) / var_v;

% states, trans, and M depend on full vs. RSSE trellis
%
  states_full = FSM( [-1 1], Nusers*Memory);
  trans_full = get_trans_multi( states_full); % state transition feasibility matrix
  M_full = size(states_full,2);               % number of states

  states_rsse = FSM( [-1 1], Memory);
  trans_rsse = ~get_trans( states_rsse, zeros(FIR_length,1));
  M_rsse = size(states_rsse,2);

% user#1 has higher SNR
%
% s = [10 15 20 25 30];
s = [5]
SNR_array = [ s; ones(Nusers-1,1)*(s-10)]

trials_array = 1000 * ones(size(s))

L_full = 3

L_array = [3]

Piter = 3

for SNR_index = 1:size(SNR_array,2)

  SNR = SNR_array(:,SNR_index)

  BER_full = zeros( Nusers, 1);
  BER_rsse = zeros( Nusers, length(L_array));

  trials = trials_array( SNR_index)

  for trial_count = 1:trials

    % generate users' sequences, FIRs, and received data
    %
    get_data;

    % full trellis
    %
    Nv = 1;

    TYPE = 31;       % multi-user, Gauss-Markov cost

    states = states_full;
    trans = trans_full;
    M = M_full;

    L = L_full;

    vit_init( L, M, Nsyms, Nv);

    for t=1:Nsyms
      vit_update;
    end

    paths = vit_paths;

    ex_gm = zeros( Nsyms, L, Nusers);

    for u=1:Nusers
      ex_gm(:,:,u) = get_ex_multi( paths, u);
    end

    tmp = (x ~= reshape(ex_gm(:,1,:),size(x)));
    BER_full = BER_full + sum(tmp)';

    % RSSE - separate trellises
    %
    Nv = Nusers;

    TYPE = 41;       % multi-user, Gauss-Markov cost, separate trellises

    states = states_rsse;
    trans = trans_rsse;
    M = M_rsse;

    for L_index = 1:length(L_array)

      L = L_array(L_index);

      LM = L*M; % L or L*M

      ex_gm = zeros( Nsyms, LM, Nusers);

      ex_p = zeros( LM, Nusers);

      vit_init( L, M, Nsyms, Nv);

      cost_init( 0); % special cost_init() call for this case

      for t = 1:Nsyms
        Vw_save = Vw;
        GG_save = GG;
        qq_save = qq;
        for P = 1:Piter % iterate
          Vparam( 4 + [1:Nusers]) = t - 1;     % reset Viterbi t
          Vw = Vw_save;                        %  and cumulative costs
          GG = GG_save;
          qq = qq_save;
          for u = 1:Nusers
            vit_update( u);
            Vparam(1) = LM; % LM = L*M to get all paths
              paths = vit_paths( u);
            Vparam(1) = L;
            ex_gm(1:t,:,u) = get_ex( paths);
            p = exp( -(VW - min(VW)));   % change costs to probabilities
            ex_p(:,u) = p/sum(p);        %  and normalize so sum is 1.
          end
        end
      end

      tmp = (x ~= reshape(ex_gm(:,1,:),size(x)));
      BER_rsse(:,L_index) = BER_rsse(:,L_index) + sum(tmp)';

    end % L_index

  end % trials

  BER_full = BER_full / (trials*Nsyms)

  BER_rsse = BER_rsse / (trials*Nsyms)

end % SNR
