%figure;
hold on;

%primka
a = rand(1);
b = rand(1);
xLine = (-5:5);
yLine = a*xLine + b;

plot(xLine, yLine, 'y--')

%body v normalnim rozdeleni
epsilon = 0.05;

v = randn(2, 100);

sx = v(1, 1:length(v));
x = sx(:);

y = v(2, 1:length(v));
y = a * x + b+ y(:);

plot(x, y, 'r+')

%robustni prokladani primarni ulohy
heigth = length(x);

y1 = y+epsilon; 
y2 = -y+epsilon;
y3 = 0;

A = [ x ones(heigth,1) -1*eye(heigth); -1*x -1*ones(heigth,1) -1*eye(heigth); zeros(heigth,1) zeros(heigth,1) -1*eye(heigth) ];

b = [y1;y2;zeros(heigth,1)];

f = [0; 0; ones(heigth,1)];

v = linprog(f, A, b );

PrimaryOptimum = f'*v


plot(x,y,'r+');

lineX = [-5:0.1:5];
lineX = lineX(:);
lineY = v(1)*lineX + v(2);

plot(lineX,lineY, 'b');
plot(lineX,lineY+epsilon, 'g--');
plot(lineX,lineY-epsilon, 'g--');

hold off;
%robustni prokladani dualni ulohy
fD = -b';
AD = A';
bD = f';

AeqD = (A'); %tady skromnej trik...nema to bejt omezeny na nerovnost, ale na rovnost....
AeqD = AeqD(1:2,:); %protoze jsem linej psat znova ty nerovnice (aD, bD), tak je tam necham, jak jsou.
beqD = f'; %Protoze je rovnost tvrdsi omezeni, tak se jejim pridanim nic neporusi a docilim
beqD = beqD(1:2); %stejnych vlasnosti, jako kdyby aD, bD ony dve rovnice nebyly a byly jen v AeqD a beqD

v2 = linprog(fD, AD, bD,AeqD,beqD,[],zeros(length(b), 1));
SecondaryOptimum = -fD*v2

%vykresleni divergenci a odpovidajicich dualnich promennych
figure
hold on
xAxis = [1:1:length(v2)];
plot(xAxis, v2, 'r+')
hold off;