% Copyright 2005 William S. Krasker

% This program requires the array S from the MATLAB program
% http://www.footballcommentary.com/dataandcode20040517/dyn2min.m.txt
% for all three values of scorefactor (0, 1, and 0.5). These must be in the MATLAB output files
% May1down4cjustS, May1down2cjustS, and May1down3cjustS respectively.
% See http://www.footballcommentary.com/twomindrill.htm
% for a discussion of that model.
%----------------------------------------------------------------------------------------------
tmax=181;

maxd=50;
dtop=2*maxd+1;
yardline=15;   % units of yardage in the two-minute drill model, dyn2min.m
p2pt=0.4;
p1pt=0.985;
p3out=0.5;     % probability of a three-and-out when the offense will definitely run the ball
tiehaircut=1;

ptdmarginal=0.2;   % regular offense
pfgmarginal=0.12;  % regular offense

% set up the vector ptime of probabilities for the length of the possession, regular offense
% first if no score, then if fg, then if td
maxdrivetime=42;
mindrivetime=3;

mean=12.5;
stddev=11;
[mu, sigma]=lognorparams(mean,stddev);
ptimenoscore=lognordensity(mu,sigma,mindrivetime,maxdrivetime);

mean=36.65;
stddev=25;
[mu, sigma]=lognorparams(mean,stddev);
ptimefg=lognordensity(mu,sigma,mindrivetime,maxdrivetime);

mean=30.2;
stddev=25;
[mu, sigma]=lognorparams(mean,stddev);
ptimetd=lognordensity(mu,sigma,mindrivetime,maxdrivetime);

ptime=(1-ptdmarginal-pfgmarginal)*ptimenoscore + ptdmarginal*ptimetd + pfgmarginal*ptimefg;
ptd=[zeros(1,mindrivetime-1) ...
       (ptimetd(mindrivetime:maxdrivetime)*ptdmarginal)./ptime(mindrivetime:maxdrivetime)];
pfg=[zeros(1,mindrivetime-1) ...
       (ptimefg(mindrivetime:maxdrivetime)*pfgmarginal)./ptime(mindrivetime:maxdrivetime)];

a=zeros(maxdrivetime,1);

PTDMARGINAL=0.07;   % slowdown offense
PFGMARGINAL=0.042;  % slowdown offense

% set up the vector PTIME of probabilities for the length of the possession, slowdown offense
% first if no score, then if fg, then if td
MAXDRIVETIME=42;
MINDRIVETIME=13;

mean=12.5;
stddev=11;
[mu, sigma]=lognorparams(mean,stddev);
ptimenoscore=lognordensity(mu,sigma,MINDRIVETIME,MAXDRIVETIME);

mean=36.65;
stddev=25;
[mu, sigma]=lognorparams(mean,stddev);
ptimefg=lognordensity(mu,sigma,MINDRIVETIME,MAXDRIVETIME);

mean=30.2;
stddev=25;
[mu, sigma]=lognorparams(mean,stddev);
ptimetd=lognordensity(mu,sigma,MINDRIVETIME,MAXDRIVETIME);

PTIME=(1-PTDMARGINAL-PFGMARGINAL)*ptimenoscore + PTDMARGINAL*ptimetd + PFGMARGINAL*ptimefg;
PTD=[zeros(1,MINDRIVETIME-1) ...
       (ptimetd(MINDRIVETIME:MAXDRIVETIME)*PTDMARGINAL)./PTIME(MINDRIVETIME:MAXDRIVETIME)];
PFG=[zeros(1,MINDRIVETIME-1) ...
       (ptimefg(MINDRIVETIME:MAXDRIVETIME)*PFGMARGINAL)./PTIME(MINDRIVETIME:MAXDRIVETIME)];

A=zeros(MAXDRIVETIME,1);

% parameters related to hurry-up offense
intercept=0.98;
slope=0.0;
intercept7=1.05;
slope7=0.0;
hurryperiods=18;

% boundary conditions
R=zeros(181,dtop,4,4,2);

Strat=ones(181,dtop,4,4);

for T1=1:4
   for T2=1:4
      R(1,maxd+1,T1,T2,1)=0.5;
      for d=maxd+2:dtop
         R(1,d,T1,T2,1)=1;
      end
   end
end

for T1=1:4
   for T2=1:4
      for d=1:dtop
         R(1,d,T1,T2,2)=1-R(1,dtop-d+1,T2,T1,1);
      end
   end
end

% last two minutes

for t=2:13

%                                    lead by two scores. (if trail by two scores, win prob is 0)
   for d=maxd+1+9:dtop
      for T1=1:4
         for T2=1:4
            R(t,d,T1,T2,1)=1;
         end
      end
   end
 
   load May1down4cjustS

%                                    lead by 8
   for T2=1:4
      R(t,maxd+1+8,1:4,T2,1)=1-0.5*p2pt*p3out*S(yardline,max(10*(t-1-14+(T2-1)*4)+1,1),1,2,1);
   end

%                                    lead by 7
   for T2=1:4
      R(t,maxd+1+7,1:4,T2,1)=1-0.5*p1pt*p3out*S(yardline,max(10*(t-1-14+(T2-1)*4)+1,1),1,2,1);
   end

%                                    lead by 6
   for T2=1:4
      R(t,maxd+1+6,1:4,T2,1)=1-(p1pt+(1-p1pt)*0.5)*p3out*S(yardline,max(10*(t-1-14+(T2-1)*4)+1,1),1,2,1);
   end
 
%                                    lead by 5
   for T2=1:4
      R(t,maxd+1+5,1:4,T2,1)=1-p3out*S(yardline,max(10*(t-1-14+(T2-1)*4)+1,1),1,2,1);
   end
 
%                                    lead by 4
   for T2=1:4
      R(t,maxd+1+4,1:4,T2,1)=1-p3out*S(yardline,max(10*(t-1-14+(T2-1)*4)+1,1),1,2,1);
   end
  
   load May1down3cjustS
  
%                                    lead by 3
   for T2=1:4
      R(t,maxd+1+3,1:4,T2,1)=1-p3out*S(yardline,max(10*(t-1-14+(T2-1)*4)+1,1),1,2,1);
   end
 
   load May1down2cjustS
  
%                                    lead by 2
   for T2=1:4
      R(t,maxd+1+2,1:4,T2,1)=1-p3out*S(yardline,max(10*(t-1-14+(T2-1)*4)+1,1),1,2,1);
   end
  
%                                    lead by 1
   for T2=1:4
      R(t,maxd+1+1,1:4,T2,1)=1-p3out*S(yardline,max(10*(t-1-14+(T2-1)*4)+1,1),1,2,1);
   end
  
%                                    tied
   for T1=1:4
      R(t,maxd+1,T1,1:4,1)=S(yardline, round(tiehaircut*10*(t-1))+1 ,1,2,T1) ...
                             +(1-S(yardline, round(tiehaircut*10*(t-1))+1 ,1,2,T1))*0.5;
   end
  
%                                    trail by 1
   for T1=1:4
      R(t,maxd+1-1,T1,1:4,1)=S(yardline,10*(t-1)+1,1,2,T1);
   end
 
%                                    trail by 2
   for T1=1:4
      R(t,maxd+1-2,T1,1:4,1)=S(yardline,10*(t-1)+1,1,2,T1);
   end
 
   load May1down3cjustS
 
%                                    trail by 3
   for T1=1:4
      R(t,maxd+1-3,T1,1:4,1)=S(yardline,10*(t-1)+1,1,2,T1);
   end
 
   load May1down4cjustS
 
%                                    trail by 4
   for T1=1:4
      R(t,maxd+1-4,T1,1:4,1)=S(yardline,10*(t-1)+1,1,2,T1);
   end
 
%                                    trail by 5
   for T1=1:4
      R(t,maxd+1-5,T1,1:4,1)=S(yardline,10*(t-1)+1,1,2,T1);
   end
 
%                                    trail by 6
   for T1=1:4
      R(t,maxd+1-6,T1,1:4,1)=(p1pt+(1-p1pt)*0.5)*S(yardline,10*(t-1)+1,1,2,T1);
   end
 
%                                    trail by 7
   for T1=1:4
      R(t,maxd+1-7,T1,1:4,1)=0.5*p1pt*S(yardline,10*(t-1)+1,1,2,T1);
   end
 
%                                    trail by 8
   for T1=1:4
      R(t,maxd+1-8,T1,1:4,1)=0.5*p2pt*S(yardline,10*(t-1)+1,1,2,T1);
   end
 
% fill in values for when opponents have ball
   for T1=1:4
      for T2=1:4
         for d=1:dtop
            R(t,d,T1,T2,2)=1-R(t,dtop-d+1,T2,T1,1);
         end
      end
   end
  
end


S=zeros(181,dtop,4,4,2);           % re-initialize S, this time as value before attempting try

for T1=1:4
   for T2=1:4
      for d=1:maxd-2
         S(1,d,T1,T2,1)=0;
      end
      S(1,maxd-1,T1,T2,1)=p2pt*0.5;
      S(1,maxd,T1,T2,1)=p1pt*0.5;
      S(1,maxd+1,T1,T2,1)=p1pt+(1-p1pt)*0.5;
      for d=maxd+2:dtop
         S(1,d,T1,T2,1)=1;
      end
   end
end

for T1=1:4
   for T2=1:4
      for d=1:dtop
         S(1,d,T1,T2,2)=1-S(1,dtop-d+1,T2,T1,1);
      end
   end
end

for t=2:13
   for T1=1:4
      for T2=1:4
         for d=1:dtop
            S(t,d,T1,T2,1)= max( p1pt*R(t-1,min(d+1,dtop),T1,T2,2)+(1-p1pt)*R(t-1,d,T1,T2,2) , ...
                                 p2pt*R(t-1,min(d+2,dtop),T1,T2,2)+(1-p2pt)*R(t-1,d,T1,T2,2) );
         end
      end
   end
   for T1=1:4
      for T2=1:4
         for d=1:dtop
            S(t,d,T1,T2,2)=1-S(t,dtop-d+1,T2,T1,1);
         end
      end
   end
end

% before the two-minute warning

for t=14:tmax
   
%                                           variables for hurry-up   
   mindrivetimeh=2;
   maxdrivetimeh=min(hurryperiods,t-1);
   ptdmarginalh=ptdmarginal*(intercept+slope*min(t,19));
   pfgmarginalh=pfgmarginal*(intercept+slope*min(t,19));
   ptimetdh=[0 (1:maxdrivetimeh-1)];
   ptimetdh=ptimetdh/sum(ptimetdh);
   ptimefgh=ptimetdh;
   ptimenoscoreh=[0 ones(1,maxdrivetimeh-1)];
   ptimenoscoreh=ptimenoscoreh/sum(ptimenoscoreh);
   ptimeh=(1-ptdmarginalh-pfgmarginalh)*ptimenoscoreh + ptdmarginalh*ptimetdh + pfgmarginalh*ptimefgh;
   ptdh=[zeros(1,mindrivetimeh-1) ...
       (ptimetdh(mindrivetimeh:maxdrivetimeh)*ptdmarginalh)./ptimeh(mindrivetimeh:maxdrivetimeh)];
   pfgh=[zeros(1,mindrivetimeh-1) ...
       (ptimefgh(mindrivetimeh:maxdrivetimeh)*pfgmarginalh)./ptimeh(mindrivetimeh:maxdrivetimeh)];

   ah=zeros(maxdrivetimeh,1);

%                                           variables for hurry-up, td only
   ptdmarginal7=ptdmarginal*(intercept7+slope7*min(t,19));
   ptime7=(1-ptdmarginal7)*ptimenoscoreh + ptdmarginal7*ptimetdh;
   ptd7=[zeros(1,mindrivetimeh-1) ...
       (ptimetdh(mindrivetimeh:maxdrivetimeh)*ptdmarginal7)./ptime7(mindrivetimeh:maxdrivetimeh)];

   a7=zeros(maxdrivetimeh,1);
   
   for d=1:dtop
      for T1=1:4
         for T2=1:4
            
            S(t,d,T1,T2,1)= max( p1pt*R(t-1,min(d+1,dtop),T1,T2,2)+(1-p1pt)*R(t-1,d,T1,T2,2) , ...
                                 p2pt*R(t-1,min(d+2,dtop),T1,T2,2)+(1-p2pt)*R(t-1,d,T1,T2,2) );
            
            for k=mindrivetime:maxdrivetime               % periods used in regular drive        
               for j=1:T2                                 % end-of-drive remaining defensive timeouts
                  used=T2-j;                              % defensive timeouts used
                  finishtime=min(t-k+used*4,t-3);
                  if finishtime >= 1
                     v(j)=ptd(k)*S( finishtime , min(d+6,dtop) , T1, j, 1) ...
                         +pfg(k)*R( finishtime , min(d+3,dtop) , T1, j, 2) ...
                          +(1-ptd(k)-pfg(k))*R( finishtime , d , T1, j, 2);
                  else
                     v(j)=R(1,d,T1,T2,1);
                  end
               end
               a(k)=min(v(1:T2));             
            end           
            regularvalue=ptime*a;
            
            for k=MINDRIVETIME:MAXDRIVETIME               % periods used in slowdown drive        
               for j=1:T2                                 % end-of-drive remaining defensive timeouts
                  used=T2-j;                              % defensive timeouts used
                  finishtime=min(t-k+used*4,t-3);
                  if finishtime >= 1
                     v(j)=PTD(k)*S( finishtime , min(d+6,dtop) , T1, j, 1) ...
                         +PFG(k)*R( finishtime , min(d+3,dtop) , T1, j, 2) ...
                          +(1-PTD(k)-PFG(k))*R( finishtime , d , T1, j, 2);
                  else
                     v(j)=R(1,d,T1,T2,1);
                  end
               end
               A(k)=min(v(1:T2));             
            end           
            SLOWDOWNVALUE=PTIME*A;
            
         if d<=maxd+1
            for k=mindrivetimeh:maxdrivetimeh             % periods used in hurry drive        
               for j=1:T1                                 % end-of-drive remaining offensive timeouts
                  used=T1-j;                              % offensive timeouts used
                  finishtime=min(t-k+used,t-2);             
                  v(j)=ptdh(k)*S( finishtime , min(d+6,dtop) , j, T2, 1) ...
                      +pfgh(k)*R( finishtime , min(d+3,dtop) , j, T2, 2) ...
                      +(1-ptdh(k)-pfgh(k))*R( finishtime , d , j, T2, 2);                         
               end
               ah(k)=max(v(1:T1));             
            end           
            hurryvalue=ptimeh*ah;

            for k=mindrivetimeh:maxdrivetimeh             % periods used in hurry drive for td only       
               for j=1:T1                                 % end-of-drive remaining offensive timeouts
                  used=T1-j;                              % offensive timeouts used
                  finishtime=min(t-k+used,t-2);             
                  v(j)=ptd7(k)*S( finishtime , min(d+6,dtop) , j, T2, 1) ...
                              +(1-ptd7(k))*R( finishtime , d , j, T2, 2);                         
               end
               a7(k)=max(v(1:T1));             
            end           
            hurryvalue7=ptime7*a7;
         else
            hurryvalue=0;
            hurryvalue7=0;
         end
            
            R(t,d,T1,T2,1)=max([regularvalue SLOWDOWNVALUE hurryvalue hurryvalue7]);
            
            if R(t,d,T1,T2,1)==SLOWDOWNVALUE
               Strat(t,d,T1,T2)=2;
            elseif R(t,d,T1,T2,1)==hurryvalue
               Strat(t,d,T1,T2)=3;
            elseif R(t,d,T1,T2,1)==hurryvalue7
               Strat(t,d,T1,T2)=4;
            end
               
         end  
      end     
   end
   
                                                     % fill in values for when opponents have ball
   for T1=1:4
      for T2=1:4
         for d=1:dtop
            R(t,d,T1,T2,2)=1-R(t,dtop-d+1,T2,T1,1);
            S(t,d,T1,T2,2)=1-S(t,dtop-d+1,T2,T1,1);
         end
      end
   end
  t 
end