function [s_dec, f_sdr] = sdr_detector(Hi, yi, opts) 
% =========================================================================
% Implementation of SDR Detector based on DSDP convex optimization package. 
% Purpose: approximation algorithm for solving maximum-likelihood detection 
% problem given by:
%
%                              OR      min   Trace(Q*X)
%    min   || y - H * s ||    <==>     s.t.  X(i,i) = 1,
%    s.t.  s(i) = -1 or +1.                  X >= 0 (positive semidefinite)
%                                            X = [s; 1]*[s, 1] (X is rank-1)
%                                            Q = [H'*H, -H'*y; -y'*H, y'*y]
%
% =========================================================================
%
%               [s_dec, f_sdr] = sdr_detector(H, y, opts)  
%
% -------------------------------------------------------------------------
% Channel model:                 
%                   y = H * s + v,    H = sqrt(snr/n) * A,
% where 
%   n is the number of transmitters,
%   m is the number of receivers,
%   snr is signal-to-noise ratio [abs. units],
%   H is an [m by n] matrix, can be real or complex,
%   v is an [m by 1] vector of Gaussian noise, can be real or complex,
%     channel model is normalized such that v ~ N(0,1) or CN(0,1), 
%   s is the vector of transmitted signals, BPSK or QPSK; 
%     use an equivalent channel representation with BPSK or QPSK signals 
%     to model other PAM or QAM constellation.
%
% -------------------------------------------------------------------------
% INPUT:
%   REQUIRED PARAMETERS:
%           H               -   real or complex channel matrix,
%           y               -   real or complex vector of received signals,
%   OPTIONAL PARAMETERS: 
%           opts.snr        -   signal-to-noise ratio in the channel [abs. units],
%           opts.con        -   signal constellation, valid options:
%                               opts.con = 'bpsk' 
%                                  [DEFAULT] BPSK modulation: 
%                                  for any i, s_i = -1 or +1,
%                               opts.con = 'qpsk'
%                                  QPSK modulation:
%                                  for any i, s_i = ((+/-)1 + (+/-)1*j)/sqrt(2),
%           opts.alg        -   quasi-maximum-likelihood algorithm to use:
%                               opts.alg = 0 
%                                  [DEFAULT for TALL or SQUARE channel matrices] 
%                                  dual-scaling algorithm with dimension reduction 
%                                  and a truncated Sphere Decoder as a warm start
%                               opts.alg = 1
%                                   dual-scaling algorithm with the eigenvector-based 
%                                   randomized rounding procedure and a truncated 
%                                   Sphere Decoder as a warm start
%                               opts.alg = 2
%                                   [DEFAULT for FAT channel matrices]
%                                   dual-scaling algorithm with dimension reduction
%                               opts.alg = 3
%                                   dual-scaling algorithm with the eigenvector-based 
%                                   randomized rounding procedure
%   OUTPUT:   
%           s_dec           -   quasi-maximum-likelihood estimate of
%                               transmitted signals, s
%           f_sdr           -   objective value achieved by s_dec:
%                                   f_sdr := || y - H * s_dec ||^2, 
% =========================================================================
% EXAMPLE:
% Given the channel with
%          H = [2, -1; -4, 7];  y = [5; -14];  SNR = 35; (15.4 dB)
%
% Run SDR Detector:
% >> opts.snr = 35;
% >> [s_est] = sdr_detector(H, y, opts);
%
% The quasi-ML estimate at the output of SDR Detector is s_est = [1; -1];
% =========================================================================
% Developed by  :   Nick (Mikalai) Kisialiou,
%                   Prof. Zhi-Quan Luo,
% Contacts      :   kisi0004@umn.edu, luoxx032@umn.edu
% Affiliations  :   Prof. Z.-Q. Luo's research laboratory,
%                   Department of Electrical and Computer Engineering,
%                   University of Minnesota, MN, United States
% Release date  :   October 15, 2006
% Version       :   0.9.2
% License       :   Please, see the website
%                   http://www.ece.umn.edu/users/luozq/software/sw_about.html
% Compatibility :   Windows(R) 2000, XP+ platform, Matlab(R) 6.0+ versions.
% Uses          :   DSDP package for convex optimization problems developed
%                   by S.J. Benson, http://www-unix.mcs.anl.gov/DSDP/
% =========================================================================

% -------------------------------------------------------------------------
% Constants
% -------------------------------------------------------------------------
% Generic constants and default parameters
YES = 0;
NO = 1;
% Quasi-ML with dimension reduction
ALG_ML_DIM_RED = 0; 
% Quasi-ML with SDP solver followed by randomized rounding based on the
% eigenvector corresponding to the largest eigenvalue of Xopt.
ALG_ML_EIG_RND = 1; 
% Run pure DSDP with dimension reduction, no initialization phase
ALG_ML_DSDP_DIM_RED = 2;
% Run pure DSDP followed by the eigenvector rounding
ALG_ML_DSDP_EIG_RND = 3;
% Default threshold for dimension reduction
DEF_THRESHOLD = 0.8;
% Default fraction of the total number of lanczos iterations in dimension reduction
DEF_LANC_ITER = 1;
% Default number of Dual-Scaling algorithm iterations before the dimension
% reduction occurs
DEF_DS_DIM_RED_MAXIT = 4;
DEF_DS_EIG_RND_MAXIT = 15;
% Default constellation
CON_BPSK = 'bpsk';
CON_QPSK = 'qpsk';
DEF_CONSTELLATION = CON_BPSK;

% -------------------------------------------------------------------------
% Switch it off if input verification takes too long
% -------------------------------------------------------------------------
CHECK_INPUT = YES;

% -------------------------------------------------------------------------
% Check that the required input is present
% -------------------------------------------------------------------------
% Verify number of input/output parameters 
if (nargin < 2),
    error('INPUT ERROR -> Not enough input parameters.');
end
if (nargin > 3),
    error('INPUT ERROR -> Too many input parameters.');
end

% -------------------------------------------------------------------------
% Set the constellation either provided by the user or default
% -------------------------------------------------------------------------
if (nargin > 2),
    if ~isfield(opts,'con'),
        % No constellation is specified
        con = DEF_CONSTELLATION;
    else
        con = opts.con;
        % Verify the validity of provided constellation
        if ~ischar(con) | isempty(con),
            error('INPUT ERROR -> Constellation is specified by the string ''bpsk'' or ''qpsk''.');
        end
        if ~(strcmp(con,CON_BPSK) | strcmp(con,CON_QPSK)),
            error('INPUT ERROR -> Currently supported constellations are ''bpsk'' and ''qpsk''.');
        end
    end
else % options structure is not provided by the user, use default
    con = DEF_CONSTELLATION;
end

% -------------------------------------------------------------------------
% Check whether the input channel is real or complex
% -------------------------------------------------------------------------
% If we deal with BPSK constellation
if strcmp(con,CON_BPSK),
    if isreal(Hi) & isreal(yi),
        % We are given the real channel
        H = Hi;
        y = yi;
        isChReal = YES;
    else 
        % There are complex entries in H and/or y
        H = [real(Hi); imag(Hi)];
        y = [real(yi); imag(yi)];
        isChReal = NO;
    end
end
% If we deal with QPSK constellation
if strcmp(con,CON_QPSK),
    % There are complex entries in H and/or y
    H = [real(Hi), -imag(Hi); imag(Hi), real(Hi)]/sqrt(2);
    y = [real(yi); imag(yi)];
    isChReal = NO;
end

% -------------------------------------------------------------------------
% Verify dimension consistency
% -------------------------------------------------------------------------
% Fat channel matrix at the input -> needs to skip truncated Sphere Decoder
if (size(H,1) < size(H,2)),
    % If no algorithm has been specified, switch from default to DSDP + dimension
    % reduction strategy
    if ((nargin <= 2) | (~isfield(opts,'alg'))),
        opts.alg = ALG_ML_DSDP_DIM_RED;
    end
    % If the algorithm has been specified to use the warm start with Sphere
    % Decoder - warn the user and switch to the one without Sphere Decoder
    if (isfield(opts,'alg')),
        if (ALG_ML_DIM_RED == opts.alg),
            fprintf('Warning : warm start with the truncated Sphere Decoder can not be used for fat channel matrices...\n');
            fprintf('Action : switching to (opts.alg = 2) with the dual-scaling algorithm and the dimension reduction strategy...\n');
            opts.alg = ALG_ML_DSDP_DIM_RED;
        elseif (ALG_ML_EIG_RND == opts.alg),
            fprintf('Warning : warm start with the truncated Sphere Decoder can not be used for fat channel matrices...\n');
            fprintf('Action : switching to (opts.alg = 3) with the dual-scaling algorithm and the randomized rounding procedure...\n');
            opts.alg = ALG_ML_DSDP_EIG_RND;
        end
    end
end
if (size(H,1) ~= size(y,1)),
    error('INPUT ERROR -> Input matrix/vector dimensions do not agree.');
end
% Number of transmitters
n = size(H,2);
% Number of receivers
m = size(H,1);

% -------------------------------------------------------------------------
% Check whether SNR is provided by the user or compute it explicitly
% -------------------------------------------------------------------------
if (nargin > 2),
    if ~isfield(opts,'snr'),
        % No SNR is provided, estimate instantaneous SNR
        snr = trace(H'*H) / m;
    else
        % Pick the provided SNR
        snr = opts.snr;
    end
else
    snr = trace(H'*H) / m;
end

% -------------------------------------------------------------------------
% Set the default values
% -------------------------------------------------------------------------
% Check if opts structure is provided
if (nargin > 2),
    if ~isfield(opts,'alg'),
        % No algorithm is selected, pick quasi-ML
        alg = ALG_ML_DIM_RED;
    else
        % Pick the provided algorithm
        alg = opts.alg;
    end
    if ~isfield(opts,'th'),
        % No threshold for dimension reduction is selected, pick default
        th = DEF_THRESHOLD;
    else
        % Pick the provided threshold
        th = opts.th;
    end
    if ~isfield(opts,'lanc'),
        % Default fraction of the total number of lanczos iterations in
        % dimension reduction
        lanc = DEF_LANC_ITER;
    else
        % Pick the provided number of lanczos iterations
        lanc = opts.lanc;
    end
    if ~isfield(opts,'maxit'),
        % Default number of Dual-Scaling algorithm iterations before the dimension
        % reduction occurs
        if ((alg == ALG_ML_DIM_RED) | (alg == ALG_ML_DSDP_DIM_RED)),
            maxit = DEF_DS_DIM_RED_MAXIT;
        else
            maxit = DEF_DS_EIG_RND_MAXIT;
        end
    else
        % Pick the provided number of Dual-Scaling iterations
        maxit = opts.maxit;
    end
    if ~isfield(opts,'ninit'), 
        % Default number of initialization vectors to check
        ninit = max(2, 10*round(2^(-2*n/snr) * 120/snr));
    else
        % Pick the provided number of initialization vectors
        if isnumeric(snr) & (snr > 0),
            ninit = opts.ninit;
        else
            error('INPUT ERROR -> SNR (Signal-to-Noise Ratio) must be a positive scalar');
        end
    end
else
    % Set all parameters to defaults
    alg = ALG_ML_DIM_RED;
    th = DEF_THRESHOLD;
    lanc = DEF_LANC_ITER;
    maxit = DEF_DS_DIM_RED_MAXIT;
    ninit = max(2, 10*round(2^(-2*n/snr) * 120/snr));
end

% -------------------------------------------------------------------------
% Check input validity
% -------------------------------------------------------------------------
if (CHECK_INPUT == YES),
    % Verify the validity of provided SNR
    if ~isnumeric(snr) | (snr <= 0) | (length(snr) > 1) | isempty(snr),
        error('INPUT ERROR -> Signal-to-Noise Ratio must be a positive scalar.');
    end
    % Verify that the algorithm is recognized
    AlgRecognized = (alg == ALG_ML_DIM_RED) ...
                  | (alg == ALG_ML_EIG_RND) ...
                  | (alg == ALG_ML_DSDP_DIM_RED) ...
                  | (alg == ALG_ML_DSDP_EIG_RND);
    if (~AlgRecognized) | (~isnumeric(alg)) | isempty(alg),
        error('INPUT ERROR -> Unrecognized algorithm requested.');
    end
    % Verify the value of dimension reduction threshold is valid
    if ~isnumeric(th) | (th > 1) | (th <= 0) | isempty(th),
        error('INPUT ERROR -> Threshold for dimension reduction must be in (0,1] range.');
    end
    % Verify the fraction of the total number of lanczos iterations in
    % dimension reduction
    if ~isnumeric(lanc) | (lanc <= 0) | (lanc > 1) | isempty(lanc),
        error('INPUT ERROR -> Fraction of lanczos iterations must be in (0,1] range.');
    end
    % Verify the number of Dual-Scaling iterations before dimension
    % reduction takes place
    if ~isnumeric(maxit) | (maxit <= 0) | isempty(maxit),
        error('INPUT ERROR -> Number of Dual-Scaling iterations must be a positive scalar.');
    end
    % Verify the number of initialization vectors
    if ~isnumeric(ninit) | (ninit <= 0) | isempty(ninit),
        error('INPUT ERROR -> Number of initialization vectors must be a positive scalar.');
    end
end

% -------------------------------------------------------------------------
% Take the QR decomposition of the channel
% -------------------------------------------------------------------------
[Q,R] = qr(H,0);

% -------------------------------------------------------------------------
% Find the zero-forcing estimate
% -------------------------------------------------------------------------
ZF = R\(Q'*y);

% -------------------------------------------------------------------------
% Make all the diagonal entries nonnegative
% -------------------------------------------------------------------------
R = sign(diag(diag(R)))*R;

% -------------------------------------------------------------------------
% Set the detector optimization parameters
% -------------------------------------------------------------------------
RadIni = 0.1;
dsdp_t = (0.0017 * snr + 0.52) * n^3 * 1.0e-6;
rad_exp = dsdp_t * snr * 2^(1.2 * 2 * n / snr) / 0.3;
RadFin = rad_exp^(snr * 0.6 / n) * 2;

% -------------------------------------------------------------------------
% Call the decoder
% -------------------------------------------------------------------------
[Xout, f_sdr, Xopt, Fsdp] = sdr(H, y, R, ZF, alg, ninit, snr, th, lanc, maxit, RadIni, RadFin);

% -------------------------------------------------------------------------
% Constellation specific output
% -------------------------------------------------------------------------
if strcmp(con,CON_QPSK),
    % If we deal with QPSK constellation, conversion to complex is needed 
    s_dec = ( Xout(1:n/2) + sqrt(-1)*Xout(n/2+1:n) ) / sqrt(2);
else
    % If we deal with BPSK constellation, no need to do anything, just output
    s_dec = Xout;
end

% =========================================================================
% End of the program
% =========================================================================
