function [fees, regionInfees] = testFeasibility(m, p, dp, yData, yEps)

m = m.update(p, m.ic);

a1 = m.a1;
a2 = m.a2;
a1P = m.a1P;
a2P = m.a2P;


%% Formulate the primal problem
for i = 1:m.nX

    v   = 1/2*a1(i,:);
    q   = full(reshape(a2(i,:), m.nX, m.nX));
    k   = m.k(i);
    Q{i} = [ k, v;
            v', q];
             
end

for i = 1:m.nY
    L{2*(i-1)+1}  =  [-yData(i)* (1-yEps(i)), m.c(i,:)];
    L{2*(i-1)+2}  = [ yData(i)* (1+yEps(i)),  -m.c(i,:)];
end
    

%% Now Get Dual
e1 =  double([(1:(m.nX+1))==1]'); %HAXOR

for  i  = 1:length(L)
    R{i} = e1*L{i} + L{i}'*e1'; 
end

for i = 1:length(L);
    for j = (i+1):length(L);
        P{i,j} = L{i}'*L{j}+ L{j}'*L{i};  
    end
end



%% Make the Sudumi Matricies

% Build up At
At = vec(e1*e1');

for i = 1:length(Q);
    At = [At, vec(Q{i})];
end

for i = 1:length(R);
    At = [At, vec(R{i})];
end

for i = 1:length(P);
    for j = (i+1):length(P);
        At = [At vec(P{i,j})];  
    end
end

nR = length(R);
nP = (length(L)-1)*length(L)/2;
nQ = length(Q);

Qobj = sparse(m.nX+1, m.nX+1);
A = [sparse(nR+nP, nQ+1) speye(nR+nP)];
b = sparse(nR+nP, 1);
c = [-1;  sparse(nR+nP+nQ, 1)];

Att = [-A; At];
btt = -c;
ct   = vec(Qobj);
ctt = [-b; ct];

K.l = size(A,1);
K.s = size(Qobj,1);
par.fid = 0;

[x, y, info] = sedumi(Att, btt, ctt, K, par);
% figure(99)
% yEQ = m.c*x;
% loglog(yEQ(1), yEQ, 'x', yEQ(1),x,'o');
% hold on;
% drawnow;

if(info.numerr == 1)
    fprintf(1, 'Feasibility test failed due to numerical errors! %d %d\n', p(1), p(2));
    fees = -1;
else
    fees = double(info.pinf == 0);
end


%% If primal-infeasible, test the whole region for primal-infeasiblilty
if fees == 0
    nP    = length(p);
    gamma = y(1);
    lambda = y(2:(m.nX+1));
    r = y((1:(m.nY*2))+(m.nX+1));
    s = y((m.nY*2+m.nX+1):end);
    Lmat = cell2mat(L');
    S = sparse(size(Lmat,1),size(Lmat,1));
    for i = 2:size(Lmat,1)
        for j = 1:(i-1)
            S(i,j) = s(1);
            s = s(2:end);
        end
    end

    for delta = 0:((2^nP)-1)
        M = gamma*e1*e1' + e1*r'*Lmat + Lmat'*r*e1' + Lmat'*S*Lmat;
        Qmod = Q;
        for i = 1:m.nX
            for j = 1:m.nP

                a1Pj = a1P(((j-1)*m.nX+1):(j*m.nX),:);
                a2Pj = a2P(((j-1)*m.nX+1):(j*m.nX),:);
                kPj  = m.kP(((j-1)*m.nX+1):(j*m.nX),:);
                v   = 1/2*a1Pj(i,:);
                q   = full(reshape(a2Pj(i,:), m.nX, m.nX));
                k   = kPj(i);
                dQ{i,j} = [k, v;
                           v', q];

                Qmod{i} = Qmod{i} + dQ{i,j}*dp(j)*(double(bitget(delta,j))-1/2);

            end
            M = M + lambda(i)*Qmod{i};
        end
        if negativeSemiDef(M)
            regionInfees = 1;
        else
            regionInfees = 0;
            return
        end
    end
else
    regionInfees = 0;   %% center-point is primal-feasible
end

