# cs_decoding.run
## compressed-sensing based decoding
## application in LP book by Matousek and Gaertner

option randseed 0;
param eps default 1e-6;

## original message to be sent
include "jll/dat/hello.dat";

## encoding matrix Q and encoded vector z
param n integer, > 0, default 2*k;
set N := 1..n;
param Q{N,K} default Uniform(-1,1);
param z{j in N} default sum{i in K} Q[j,i]*w[i];

## transmit (=tamper with) z
param e default 0.1;
param xerr{j in N} default if (Uniform01() < e) then Uniform(-1,1) else 0;
param zbar{j in N} default z[j] + xerr[j];


# formulate and solve the problem of finding A orthogonal to Q
var a{K,N} >= -1, <= 1;
maximize nontrivial: sum{i in K, j in K} Uniform(-1,1)*a[i,j];
subject to zero{i in K, h in K}: sum{j in N} a[i,j]*Q[j,h] = 0;
problem findA: a, nontrivial, zero;
option solver cplex;
option cplex_options "baropt bardisplay=1";
#option cplex_options "display=1";
solve findA;
param A{i in K, j in N};
let {i in K, j in N} A[i,j] := a[i,j];

# compute rhs vector b
param b{i in K} default sum{j in N} A[i,j]*zbar[j];



# formulate and solve the ell_1 norm minimization LP
var x{N};
var s{N};
minimize norm1: sum{j in N} s[j];
subject to n1lin1{j in N}: x[j] >= -s[j];
subject to n1lin2{j in N}: x[j] <=  s[j];
subject to enc{i in K}: sum{j in N} A[i,j]*x[j] = b[i];
problem norm1LP: x,s,norm1,n1lin1,n1lin2,enc;
option solver cplex;
option cplex_options "baropt bardisplay=1";
#option cplex_options "display=1";
solve norm1LP;




# compute error-free solution
param y{j in N};
let {j in N} y[j] := zbar[j] - x[j];

# verify it
param count integer, default 0;
for {j in N : abs(z[j] - y[j]) > eps} {
  let count := count + 1;
}
printf "[CS] sent/received: %d different components over %d\n", count, n;
param yzerr := sum{j in N} abs(z[j] - y[j]);
printf "[CS] sent/received ell-1 error = %g\n", yzerr;
printf "[CS] sent/rec mean ell-1 error = %g\n", yzerr / n;



# recover original message
param QTQ{i in K, h in K} default sum{j in N} Q[j,i]*Q[j,h];
param QTy{i in K} default sum{j in N} Q[j,i]*y[j];
# invert QTQ
param M default 1e6;
var q{K,K} <= M, >= -M;
subject to inverse{i in K, h in K}:
  sum{l in K} q[i,l]*QTQ[l,h] = if (i==h) then 1 else 0;
problem invertQTQ: q, inverse;
option solver cplex;
option cplex_options "baropt bardisplay=1";
#option cplex_options "display=1";
solve invertQTQ;



# w_recovered = round(Q^-1 Q^T y)
param wrec{K};
## (reason this works so well: rounding operator - poorer on floats)
let {i in K} wrec[i] := round(sum{h in K} q[i,h]*QTy[h]);
# cap to [0,1]
let {i in K} wrec[i] := if wrec[i] < 0 then 0 else if wrec[i] > 1 then 1 else wrec[i];

# display results
printf "[CS] w_org = ";
for {i in K} {
   printf "%d", w[i];
}
printf "\n";
printf "[CS] w_rec = ";
for {i in K} {
   printf "%d", wrec[i];
}
printf "\n";
param werr >= 0;
let werr := sum{i in K} abs(w[i] - wrec[i]);
printf "[CS] ||w_org - w_rec||_1 = %g\n", werr;

include "jll/bin2str.run";