function figures
%% parameters
s_tp=1;
K_T=1.5;
gamma=0.1;
beta=0.15;
eta=1e-2;
mu=0.05;
s_fp=0; C_fS=50;
C_t=200;
C_S=15;
C_T=400;
n=10;
N=1e3;

omega=@(p) gamma*(1-p).*(beta*p+eta);
f=@(p,lambda,delta) (1-p).*(beta*p+eta)-mu*p-(1+delta*K_T).*(lambda*s_tp.*p+omega(p));
C=@(lambda,delta,p) omega(p)*C_t+lambda.*(C_S+C_t*s_tp*p+C_fS*s_fp*(1-p))+delta*C_T.*(lambda*s_tp.*p+omega(p));
pvals=linspace(0,1,10000).^2; %square to have more points when functions change more rapidly
pvals=pvals(2:end); %to avoid division by zero
pstar=(C_S+C_fS*s_fp)/(s_tp*C_T/K_T+C_fS*s_fp-C_t*s_tp); %this is P' from the text
P_0=0.5-eta/(2*beta)-mu/(2*beta*(1-gamma))...
	+sqrt((0.5-eta/(2*beta)-mu/(2*beta*(1-gamma)))^2+eta/beta);

pvals1=pvals(pvals<pstar); %p s.t. delta=1 is optimal
pvals0=pvals(pvals>=pstar); %p s.t. delta=0 is optimal

graycolor=[0.7 0.7 0.7];

lambda_L=@(p,delta) 1./(p*s_tp).*((beta*p.*(1-p)+eta*(1-p)-mu*p)./(1+delta*K_T)-omega(p));
C_L=@(p,delta) C(lambda_L(p,delta),delta,p);

%% unchanged prevalence
lambda_U=@(p,delta) 1./(N*p*s_tp).*(n./(1+delta*K_T)-N*omega(p));

figure(1)
clf
hold on
plot(pvals1,lambda_U(pvals1,1),'k')
plot(pvals0,lambda_U(pvals0,0),'k')
plot([pstar pstar],[0,lambda_U(pstar,0)],'k--')
hold off
axis([0 1 0 0.1])
box on
set(gca,'YTick',[],'Fontsize',12)
set(gca,'XTick',[0,pstar,1])
set(gca,'XTickLabel',{'0','P''','1'})
xlabel('p')
ylabel('\lambda^*_U(n;p)','Rotation',0,'HorizontalAlignment','Right')
text(pstar/2,3e-2,'\delta_U^*=1')
text(0.5*(pstar+1),3e-2,'\delta_U^*=0')
set(findall(gcf,'Type','line'),'linewidth',1)
set(findall(gcf,'Type','text'),'fontsize',12)
savefig('schematic_lambdap');

%%
C_U=@(p,delta) C(lambda_U(p,delta),delta,p);
figure(2)
clf
hold on
plot(pvals1,C_U(pvals1,0),'k','color',graycolor)
plot(pvals0,C_U(pvals0,0),'k','linewidth',1.5)
plot(pvals0,C_U(pvals0,1),'k','color',graycolor)
plot(pvals1,C_U(pvals1,1),'k','linewidth',1.5)
plot([pstar pstar],[0,C_U(pstar,0)],'k--');
hold off
axis([0 1 0 5])
box on
set(gca,'XTick',[0,pstar,1])
set(gca,'XTickLabel',{'0','P''','1'})
set(gca,'YTick',[],'Fontsize',12)
text(pstar/2,1.5,'\delta_U^*=1')
text(0.5*(pstar+1),1.5,'\delta_U^*=0')
xlabel('p')
ylabel('C^*_U(n;p)','Rotation',0,'HorizontalAlignment','Right')
set(findall(gcf,'Type','line','linewidth',0.5),'linewidth',1)
set(findall(gcf,'Type','text'),'fontsize',12)
savefig('schematic_cost');

%% long term costs
figure(3)
clf
plot(pvals1,lambda_L(pvals1,1),'k',pvals0,lambda_L(pvals0,0),'k',[pstar pstar],[0,lambda_L(pstar,0)],'k--')
axis([0 1 0 0.2])
set(gca,'YTick',[],'Fontsize',12)
set(gca,'XTick',[0,pstar,P_0,1])
set(gca,'XTickLabel',{'0','P''','','1'})
text(P_0,-0.009,'P_0','HorizontalAlignment','center')
text(pstar/2,0.05,'\delta_L^*=1')
text(0.5*(pstar+1),0.05,'\delta_L^*=0')
xlabel P
ylabel('\lambda^*_L(P)','Rotation',0,'HorizontalAlignment','Right')
set(findall(gcf,'Type','line','linewidth',0.5),'linewidth',1)
set(findall(gcf,'Type','text'),'fontsize',12)
savefig('schematic_lambdaE');

%%
figure(4)
clf
hold on
plot(pvals1,C_L(pvals1,0),'k','color',graycolor)
plot(pvals0(pvals0<=P_0),C_L(pvals0(pvals0<=P_0),0),'k','linewidth',1.5)
plot(pvals0(pvals0<=P_0),C_L(pvals0(pvals0<=P_0),1),'k','color',graycolor)
plot(pvals1,C_L(pvals1,1),'k','linewidth',1.5)
plot([pstar pstar],[0,C_L(pstar,0)],'k--');
hold off
axis([0 1 0 10])
box on
set(gca,'YTick',[],'Fontsize',12)
set(gca,'XTick',[0,pstar,P_0,1])
set(gca,'XTickLabel',{'0','P''','','1'})
text(P_0,-0.48,'P_0','HorizontalAlignment','center')
text(pstar/2,3,'\delta_L^*=1')
text(0.5*(pstar+1),3,'\delta_L^*=0')
xlabel P
ylabel('C^*_L(P)','Rotation',0,'HorizontalAlignment','Right')
set(findall(gcf,'Type','line','linewidth',0.5),'linewidth',1)
set(findall(gcf,'Type','text'),'fontsize',12)
savefig('schematic_costE');

%% minimizing total cost
%the first order condition of the optimization problem eq 27 
%is a cubic equation in p_1
%we write it as a vector [linear term, const, 1/p, 1/p^2]
r=0.05;

% C_L'(p_1) (unfortunately the parameters are kind of hard coded here)
C_Ldiff=@(delta) [-60-120*delta,15.975+36.3375*delta,0,-0.135+0.0225*delta]/(1+1.5*delta);

cubic=@(delta,alpha) C_Ldiff(delta)+[0,-(C_t*s_tp-C_fS*s_fp+delta*s_tp*C_T)/(1+delta*K_T)+alpha/r,...
	-(C_S+C_fS*s_fp)/(1+delta*K_T),0];
alphas=1e3:1e2:1e4;
figure(5)
clf
for i=1:length(alphas)
	temproots1=roots(cubic(1,alphas(i)));
	temproots0=roots(cubic(0,alphas(i)));
	rvals{i}=[temproots1(real(temproots1)<=pstar & real(temproots1)>=0);temproots0(real(temproots0)>=pstar & real(temproots0)<=1)];
end
plot(alphas,cat(1,rvals{:}),'k')
xlabel('\alpha')
ylabel('p_1^*','Rotation',0)
set(findall(gcf,'Type','line','linewidth',0.5),'linewidth',1)
set(findall(gcf,'Type','text','fontsize',10),'fontsize',12)
savefig('optcontrol2');

%% figure requested by referee
%specific case lambdabar=0.8,r=0.1,p0=0.05,alpha=500
%p_0 and p1 found by running above and looking at plot(pvals,action)
lambdabar=0.8;
r=0.1;
alpha=500;
delta=1;
p_0=0.05; %optimal to screen
T=5;
p1=0.01015; %prevalence to stop at 
l1=lambda_L(p1,delta);
odefunc=@(t,x,ll) [f(x(1,:),ll*(x(1,:)>p1*1.001)+l1*(x(1,:)<=p1*1.001),delta);...
	ll*(x(1,:)>p1*1.001)+l1*(x(1,:)<=p1*1.001);...
	exp(-r*t).*(C(ll*(x(1,:)>p1*1.001)+l1*(x(1,:)<=p1*1.001),delta,x(1,:))+alpha*x(1,:))];
x0=[p_0;0;0];
[t1,y1]=ode15s(@(t,x) odefunc(t,x,lambdabar),[0,T],x0);
[t2,y2]=ode15s(@(t,x) odefunc(t,x,l1),[0,T],x0);
figure(9)
clf
subplot(3,1,1)
plot(t1,y1(:,2),'k',t2,y2(:,2),'k:')
ylabel('cumulative screening')
legend('optimal','const','Location','NorthOutside','Orientation','Horizontal');
subplot(3,1,2)
plot(t1,y1(:,1),'k',t2,y2(:,1),'k:')
axis([0 T 0 p_0])
ylabel('p','Rotation',0)
subplot(3,1,3)
plot(t1,y1(:,3),'k',t2,y2(:,3),'k:')
xlabel('time [years]')
ylabel('\int_0^t e^{-rt} costs')
set(findall(gcf,'Type','line','linewidth',0.5),'linewidth',1)
set(findall(gcf,'Type','text','fontsize',10),'fontsize',12)
savefig('responsefig');

function savefig(fname)
saveas(gcf,[fname,'.fig'])
print(['-f',num2str(gcf)],'-depsc2',[fname,'.eps'])
system(['epstopdf ',fname,'.eps']);

