% Part (b)
% introduce a slack var x3; the problem is now:
% 
%      min      -x1-x2
%    x1,x2,x3
%    
%    s.t. x1+x2+x3 = 1
%        x1,x2,x3 >= 0
% 
A = [1 1 1];
b = 1;
c = [-1; -3; 0];

x0 = [0.3; 0.3; 0.4];

%%%%%% Part 1(b) %%%%%%
if (ismember(input('Part 1(b): ','s'),['y' 'Y']))
  format compact;
  [xstar,ubflag,fconv,xconv] = affscale(A,b,c,x0,eps,0.5);
  disp('  The solution is:');
  disp(xstar(1:2));
  disp('  The optimal objective value is:');
  disp(fconv(length(fconv)));
end;

%%%%%% Part 1(c) %%%%%%

if (ismember(input('Part 1(c): ','s'),['y' 'Y']))
  beta_v = [0.01 0.5 0.99];

  for i = 1:length(beta_v)
    [xstar,ubflag,fconv,xconv] = affscale(A,b,c,x0,eps,beta_v(i));
    fconv_v{i}=fconv;
  end;

  figure;
  hold on;
  plot(fconv_v{1},'r-+');
  plot(fconv_v{2},'b-o');
  plot(fconv_v{3},'g-x');
  set(gca,'Xscale','log');
  legend({'\beta = 0.01','\beta = 0.5','\beta = 0.99'});
  xlabel('Number of Iterations');
  ylabel('Objective Value');
  title('Convergence plot for 1(c)');
  hold off;
  print -depsc -tiff conv-plot-1c
end;

%%%%%% Part 1(d) %%%%%%

if (ismember(input('Part 1(d): ','s'),['y' 'Y']))
  x0_v = [1/3*[1;1;1] [0.99;0.005;0.005]];

  for i = 1:2
    [xstar,ubflag,fconv,xconv] = affscale(A,b,c,x0_v(:,i),eps,0.9);
    xconv_v{i} = xconv;
  end;

  figure;
  hold on;
  plot(xconv_v{1}(1,:),xconv_v{1}(2,:),'r-+','LineWidth',1);
  plot(xconv_v{2}(1,:),xconv_v{2}(2,:),'b-o','LineWidth',1);
  axis square;
  axis([-0.2 1.2 -0.2 1.2]);
  % feasible region border
  line([1 0],[0 1],'LineWidth',1,'LineStyle','--','Color','black');
  line([0 0],[1 0],'LineWidth',1,'LineStyle','--','Color','black');
  line([0 1],[0 0],'LineWidth',1,'LineStyle','--','Color','black');
  legend({'x_0 = [1/3;1/3;1/3]','x_0 = [0.99;0.005;0.005]'});
  title('Solver path for 1(d); \beta = 0.9');
  hold off;
  print -depsc -tiff path-plot-1d
end;

%%%%%% Part 2(a,b) %%%%%%

% variables
% [x1..x4 y1..y4 z2..z4 s1..s4]

prod_cost = [35*ones(1,4)];
purch_cost = [50*ones(1,4)];
inv_cost = [5*ones(1,3)];
slack_cost = zeros(1,4);

costs = [prod_cost purch_cost inv_cost slack_cost]';
%    x1 x2 x3 x4 y1 y2 y3 y4 z2 z3 z4 s1 s2 s3 s4
Ademand = [eye(4) eye(4) [-1 0 0;   
                           1 -1 0;
                           0  1 -1;
                           0  0  1] zeros(4,4)];
Bdemand = [150 160 225 180]';

Aprod   = [eye(4) zeros(4,4) zeros(4,3) eye(4)];
Bprod   = [160 160 160 160]';

Atot = [Ademand;Aprod];
Btot = [Bdemand;Bprod];

Atot_base = Atot;
Btot_base = Btot;
costs_base = costs;

x0 = [100 100 100 100 70  60 125  60  20  20  20 60  60  60  60]';
   
[xstar,ubflag,fconv,xconv] = affscale(Atot,Btot,costs,x0,0.0001,0.5);

disp(' Optimal Production:');
disp([xstar(1:4)]');
disp(' Optimal Purchasing:');
disp([xstar(5:8)]');
disp(' Optimal Inventory:');
disp([0;xstar(9:11)]');
disp(' Optimal Cost');
disp(fconv(length(fconv)));

fval_base = fconv(length(fconv));

%%%%%% Part 2(c) %%%%%%

month = {'Jan' 'Feb' 'Mar'};
max_level = [151 153 155];

for i=1:3
    disp('If maintainence is in:');
    disp(month{i});
    Bprod(i) = max_level(i);
    Btot_man = [Bdemand;Bprod];

    [xstar,fval] = linprog(costs,[],[],Atot,Btot_man,zeros(15,1),[]);
    
    disp(' Optimal Production:');
    disp([xstar(1:4)]');
    disp(' Optimal Purchasing:');
    disp([xstar(5:8)]');
    disp(' Optimal Inventory:');
    disp([0;xstar(9:11)]');
    disp(' Optimal Cost');
    disp(fval);
end;

%%%%%% Part 2(d) %%%%%%

disp('\n\n%%%%%% Part 2(d) %%%%%%\n\n');

Lcols= [eye(3);zeros(5,3);ones(1,3)];
sLcol = [zeros(8,1);1];

for i=1:3
    disp('If purchase from D in month:');
    disp(month{i});
    Atot_D = [[Atot_base;zeros(1,15)] Lcols(:,i) sLcol];
    Btot_D = [Btot_base;50];
    costs_D = [costs_base;45;0];

    [xstar,fval] = linprog(costs_D,[],[],Atot_D,Btot_D,zeros(17,1),[]);
    
    disp(' Optimal Production:');
    disp([xstar(1:4)]');
    disp(' Optimal Purchasing (C):');
    disp([xstar(5:8)]');
    disp(' Optimal Purchasing (D):');
    disp([zeros(1,i-1) xstar(16)]);
    disp(' Optimal Inventory:');
    disp([0;xstar(9:11)]');
    disp(' Optimal Cost');
    disp(fval);
end;

%%%%%% Part 2(e) %%%%%%

fval_Coff = [];

for delta = [0:1:10]
    costs_Coff = costs_base;
    costs_Coff(6) = 50-delta;
    
    [xstar,fval] = linprog(costs_Coff,[],[],Atot_base,Btot_base,zeros(15,1),[]);
    
    fval_Coff = [fval_Coff;fval(length(fval))];
end;

figure;
plot([0:1:10],fval_Coff-fval_base);
xlabel('Delta cost in Feb [$]');
ylabel('Change in opt cost [$]');
title('Part 2(e) Break point for new price from C in Feb');
print -depsc -tiff break-plot-2e;

%%%%%% Part 2(f) %%%%%%

disp(sprintf('\n\n%%%%%% Part 2(f) %%%%%%\n\n'));

costs_newInt = costs_base;
costs_newInt(9:11) = 8;

[xstar,fval] = linprog(costs_newInt,[],[],Atot_base,Btot_base,zeros(15,1),[]);

disp(' Optimal Production:');
disp([xstar(1:4)]');
disp(' Optimal Purchasing:');
disp([xstar(5:8)]');
disp(' Optimal Inventory:');
disp([0;xstar(9:11)]');
disp(' Optimal Cost');
disp(fval);

%%%%%% Part 2(g) %%%%%%

disp(sprintf('\n\n%%%%%% Part 2(g) %%%%%%\n\n'));

B_g = Btot_base;
B_g(1) = 90;

[xstar,fval] = linprog(costs_base,[],[],Atot_base,B_g,zeros(15,1),[]);

disp(' Optimal Production:');
disp([xstar(1:4)]');
disp(' Optimal Purchasing:');
disp([xstar(5:8)]');
disp(' Optimal Inventory:');
disp([0;xstar(9:11)]');
disp(' Optimal Cost');
disp(fval);