%% prSAMP package demo
% by Boshra Rajaei 29 July 2016

%% Print some initial info
clear;
clc;
close all;

fprintf('............prSAMP package demo...............\n')
fprintf('prSAMP: phase retrieval swept approximate message passing\n')
fprintf('Using this demo, you can run prSAMP for both calibration and recovery phases.\n')
fprintf('In phase retrieval, y=|Hx|, calibration is about estimating transmission matrix, H, using some train input-output (x,y) samples.\n')
fprintf('Recovery is solving y=|Hx| for known or estimated H and unknown complex or binary x\n');
fprintf('NOTE: this demo in calibration step uses mpiexec command which is provided by MPICH2 \n')
fprintf('..............................................\n\n\n')

fprintf('Choose one scenario:\n')
fprintf('1. Calibration\n')
fprintf('2. Recovery (using known TM) \n')
sc = input('');

%% Calibration
if(sc==1)
    fprintf('\n............prSAMP calibration...............\n')
    fprintf('Here we estimate $H in C^{MxN}$ using $X in {0,1}^{NxP}$ and $Y in R^{MxP}$ known input-output samples from Y=|HX|.\n')
    n = input('Input (x) dimension: N=');
    m = input('Measurement (y) dimension: M=');
    p = input('Number of calibration patterns: P=');
    snr = input('Noise snr: ');
    %num_blk = 1;  %number of blocks
    %if(n>=512)
    %    fprintf('In high dimensions it is recommended to use block-based prSAMP.\n')
    %    blkopt = input('Do you want to try block-based prSAMP? (y/n)   ', 's');
    %    if blkopt=='y'
    %        blk_siz = input('Block size:');
    %        num_blk = n/blk_siz;
    %        if(num_blk~=floor(num_blk))
    %            fprintf('Error: N should be divisible to block size\n');
    %            return;
    %        end
    %    end
    %end
    
    %if(num_blk==1)
    Ho = randn(m,n) + randn(m,n)*1i; Ho = Ho/sqrt(n); %TM
    X = double(rand(n,p)>0.5); %input binary patterns
    Y = addNoise(abs(Ho*X), snr);
    prblmid = save2file_calib(X, Y, snr, 0); %write the problem data to text file
    
    fprintf('\nRunning prSAMP (it takes few seconds to multiple minutes according to problem size and processor power) ..........\n');
    np=feature('numCores');
    %c = parcluster('local'); % build the 'local' cluster object
    %np = c.NumWorkers        % get the number of workers
    st = system(['mpiexec -n ' num2str(np) ' ./prSAMP_Calib ' num2str(prblmid) '0']);
    
    if(st)
        fprintf('Error in running prSAMP_Calib.exe\n');
        fprintf('Run <<mpirun -np  <number_of_processors> calibration/prSAMP_calib %d0>> in an unix shell\n', prblmid);
        fprintf('Copy-paste resulting output%d0.txt in the same directory as demo.m\n', prblmid)
        fprintf('Press any key when done...');
        pause;
    end
    
    H_hat = loadfile_calib([num2str(prblmid) '0'], m, n); %load output.txt file which contains matrix H estimation
    
    % calculate performance over measurements
    for i=1:size(X,2)
        y1 = Y(:,i);
        y2 = abs(H_hat*X(:,i));
        msey(i) = mean((y2-y1).^2)/mean(y1.^2);
    end
    fprintf('Calibration performance (in terms of average MSE): %0.4f \n',mean(msey));
    %else
    %    fprintf('Block-based prSAMP is not yet implemented in this demo.\n')
    %    return;
    %end
    
    
end
%% Recovery
if(sc==2)
    fprintf('\n............prSAMP recovery...............\n')
    fprintf('Here we estimate $x in C^N$ using $H in {0,1}^{MxN}$ and $y in R^M$ from y=|Hx|.\n')
    n = input('Input (x) dimension: N=');
    m = input('Measurement (y) dimension: M=');
    snr = input('Noise snr: ');
    
    xo = randn(n,1) + randn(n,1)*1i; xo = xo/sqrt(n); %unknown input vector
    H = double(rand(m,n)>0.5); %TM
    y = addNoise(abs(H*xo), snr);
    prblmid = save2file_recovery(H, y, snr, 0); %write the problem data to text file
    
    fprintf('\nRunning prSAMP (it takes few seconds to multiple minutes according to problem size and processor power) ..........\n');
    
    st = system(['prSAMP_Rec_BinaryTM ' num2str(prblmid) '0']);
    
    if(st)
        fprintf('Error in running prSAMP_Rec_BinaryTM.exe\n');
        return
    end
    
    x_hat = loadfile_recovery([num2str(prblmid) '0'], n); %load output.txt file which contains x estimation
    
    % calculate performance 
    msex = mean(abs(x_hat-xo).^2);
    fprintf('Recovery performance (in terms of MSE): %0.4f \n',msex);
end