% This code is for inv control problem pg 26 prob 3 in Denardo ONLY
% obtain cij and X matrix from cost for inventory_Denardo_pg26_prob_3.m
% The program can detect upto 12 paths in any size network. 
% change min to max for longest path
% 0 means nodes not connected
% 0.00001 means that the arc exsits but the cost value is 0. 
% numpath= maximum number of paths. Some paths may be redundant.
clc
%set up the model parameters
%**********************
prodc = [3 5 3]; % monthly producton cost
%**********************
f=0;
F=0;
P=0;
Path=0;
Length=0;
c = [
   0     3     8    13     0     0     0     0
     0     0     0     0   100   100   100     0
     0     0     0     0    17   100   100     0
     0     0     0     0    12    15   100     0
     0     0     0     0     0     0     0     5
     0     0     0     0     0     0     0     6
     0     0     0     0     0     0     0   100
     0     0     0     0     0     0     0     0
];

X = [
 
     0     1     2     3     0     0     0     0
     0     0     0     0     4     5     6     0
     0     0     0     0     3     4     5     0
     0     0     0     0     2     3     4     0
     0     0     0     0     0     0     0     1
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0    -1
     0     0     0     0     0     0     0     0
];

count=0;
temp=0;
d=0;
f=0;
D=size(c);
f(D(1,1),1)=0;
for i= D(1,1)-1:-1:1
    count=0;
    for j=i:1: D(1,1)
         if c(i,j) ~= 0
             count=count+1;
             temp(count,1)= c(i,j)+f(j,1);
             temp(count,2)= j;
         end
    end
    f(i,1)=min(temp (:,1)); % change min to max for longest path
    cz=2;
    G=size(temp);
    for k=1:1:G(1,1)
        if f(i,1)==temp(k,1)
            f(i,cz)=temp(k,2);
            cz=cz+1;
        end
    end
    temp=0;
end
f;

% the code below tells you the number of optimal paths in the solution
cv=0;
     H=size(f);
     numpath=1;
for y=1:1:H(1,1)-1
       cv(y,1)=0;  
    for x=2:1:H(1,2)
        if f(y,x)~= 0
            cv(y,1)=cv(y,1)+1;
        end
              
    end 
    numpath=numpath*cv(y,1)  ;
end        
jj=max(cv);  
numpath

F=0; 
    o=1;
    oo=1;
for y=1:1:H(1,1)-1
       cv(y,1)=0;  
    for x=2:1:H(1,2)
        if f(y,x)~= 0
            cv(y,1)=cv(y,1)+1;
        end
    end
    mm=numpath/cv(y,1);
    m=min(mm,jj);
    if mm==numpath
        F(y,1)=f(y,1); 
        for u=2:1:numpath+1
        F(y,u)=f(y,2);
        end
    elseif m==1;
         for u=1:1:numpath+1
         F(y,u)=f(y,u);
         end
    elseif m==2
         if o==1
          F(y,1)=f(y,1);
              for u=2:cv(y,1):numpath+1
              F(y,u)=f(y,2);
              F(y,u+1)=f(y,3);
              end
         end
         if o==2
          F(y,1)=f(y,1);
              for u=2:cv(y,1)+2:numpath+1
              F(y,u)=f(y,2);
              F(y,u+1)=f(y,2);
              F(y,u+2)=f(y,3);
              F(y,u+3)=f(y,3);
              end
          end
         if o==3
          F(y,1)=f(y,1);
              for u=2:cv(y,1)+6:numpath+1
                F(y,u)=f(y,2);
                F(y,u+1)=f(y,2);
                F(y,u+2)=f(y,2);
                F(y,u+3)=f(y,2);
                F(y,u+4)=f(y,3);
                F(y,u+5)=f(y,3);
                F(y,u+6)=f(y,3);
                F(y,u+7)=f(y,3);
              end  
         end
         o=o+1;
    elseif m==3
         if oo==1;
          F(y,1)=f(y,1);
            for u=2:cv(y,1):numpath+1
            F(y,u)=f(y,2);
            F(y,u+1)=f(y,3);
            F(y,u+2)=f(y,4);
            end 
         end
         if oo==2;
            F(y,1)=f(y,1);
            for u=2:cv(y,1)+9:numpath+1
            F(y,u)=f(y,2);
            F(y,u+1)=f(y,2);
            F(y,u+2)=f(y,2);
            F(y,u+3)=f(y,3);
            F(y,u+4)=f(y,3);
            F(y,u+5)=f(y,3);
            F(y,u+6)=f(y,4);
            F(y,u+7)=f(y,4);
            F(y,u+8)=f(y,4);      
            end 
         end   
          oo=oo+1;
      elseif m==4
          F(y,1)=f(y,1);
          for u=2:cv(y,1):numpath+1
          F(y,u)=f(y,2);
          F(y,u+1)=f(y,3);
          F(y,u+2)=f(y,4);
          F(y,u+3)=f(y,5);
          end 
      elseif m==5
          F(y,1)=f(y,1);
          for u=2:cv(y,1):numpath+1
          F(y,u)=f(y,2);
          F(y,u+1)=f(y,3);
          F(y,u+2)=f(y,4);
          F(y,u+3)=f(y,5);
          F(y,u+4)=f(y,6);
          end 
          elseif m==6
    F(y,1)=f(y,1);
          for u=2:cv(y,1):numpath+1
          F(y,u)=f(y,2);
          F(y,u+1)=f(y,3);
          F(y,u+2)=f(y,4);
          F(y,u+3)=f(y,5);
          F(y,u+4)=f(y,6);
          F(y,u+5)=f(y,7);
          end        
     end
 end   
HH=size(F);
Z=HH(1,2);
if HH(1,2)>13
   Z=13;
end
   for L=2:1:Z
        rr=3;
        P(1,L-1)=1;
        P(2,L-1)=F(1,L);  
        while P(rr-1,L-1)<D(1,1)
        P(rr, L-1)=F(P(rr-1,L-1),L);
        rr=rr+1;
        end
    end
   f;
   F;
  Length=F(1,1)
  Path=P 

  szp=size(Path);
  action=0;
  for j=1:1:szp(1,2)
    
for i=1:1:szp(1,1)-1
  rown(i,1)=Path(i,j);
  rown(i,2)=Path(i+1,j);
  if rown(i,2)~=0
  action(i,j)=X(rown(i,1),rown(i,2));
  end
end
  end
action

%verification for the correct answer.
% there are several paths that yield the correctlength but not all paths
% are correct if both cost and path are traced together.
%the cost matrix has taken into account the cost due to change in
%production qty.
% The following steps will select the path the gives the correct answer.
% this problem is different from the other inv control problems beacuse the
%cij for stage t depends  on the production at stage t and the
%change in production from stage t-1. While tracing paths the above
%conditions must be considered.
szaction=size(action);
costf=0;
startinv=1;

finalaction=0;
for j=1:1:szaction(1,2)
for i=1:1:szaction(1,1)
 costf = costf + (prodc(1,i)*action (i,j))+ 2*(abs(startinv-action(i,j)));
 startinv=action(i,j);
end
if costf==Length
    for k=1:1:szaction(1,1)
    finalaction(k,1)= action(k,j);
    end
    costf
end
costf=0;
end
finalaction

% with max hold 4
% c = [
%      0     3     8    13   100   100     0     0     0     0     0     0
%      0     0     0     0     0     0   100   100   100   100   100     0
%      0     0     0     0     0     0    17   100   100   100   100     0
%      0     0     0     0     0     0    12    15   100   100   100     0
%      0     0     0     0     0     0     9    12    15   100   100     0
%      0     0     0     0     0     0     6     9    12    15   100     0
%      0     0     0     0     0     0     0     0     0     0     0     5
%      0     0     0     0     0     0     0     0     0     0     0     2
%      0     0     0     0     0     0     0     0     0     0     0   100
%      0     0     0     0     0     0     0     0     0     0     0   100
%      0     0     0     0     0     0     0     0     0     0     0   100
%      0     0     0     0     0     0     0     0     0     0     0     0
% ];
% 
% X = [
%      0     1     2     3     4     5     0     0     0     0     0     0
%      0     0     0     0     0     0     4     5     6     7     8     0
%      0     0     0     0     0     0     3     4     5     6     7     0
%      0     0     0     0     0     0     2     3     4     5     6     0
%      0     0     0     0     0     0     1     2     3     4     5     0
%      0     0     0     0     0     0     0     1     2     3     4     0
%      0     0     0     0     0     0     0     0     0     0     0     1
%      0     0     0     0     0     0     0     0     0     0     0     0
%      0     0     0     0     0     0     0     0     0     0     0    -1
%      0     0     0     0     0     0     0     0     0     0     0    -2
%      0     0     0     0     0     0     0     0     0     0     0    -3
%      0     0     0     0     0     0     0     0     0     0     0     0
% ];
%  

% 
% with max hold 2
% c = [
%    0     3     8    13     0     0     0     0
%      0     0     0     0   100   100   100     0
%      0     0     0     0    17   100   100     0
%      0     0     0     0    12    15   100     0
%      0     0     0     0     0     0     0     5
%      0     0     0     0     0     0     0     6
%      0     0     0     0     0     0     0   100
%      0     0     0     0     0     0     0     0
% ];
% 
% X = [
%  
%      0     1     2     3     0     0     0     0
%      0     0     0     0     4     5     6     0
%      0     0     0     0     3     4     5     0
%      0     0     0     0     2     3     4     0
%      0     0     0     0     0     0     0     1
%      0     0     0     0     0     0     0     0
%      0     0     0     0     0     0     0    -1
%      0     0     0     0     0     0     0     0
% ];
