clear

T=20;

%Data
PhD=[2120 2023];
PostdDoc= [3013 1891];
Faculty = [1557 505];
HHMI = [113 38];
NAS = [210 47];
Award=[53 10];
Elite=HHMI+NAS+Award;


%Postdocs
%Elite male recruits female
q_EM=0.31
%Normal male recruits female
q_NM=0.37
%Elite/normal female recruits female
q_EF=0.46
q_NF=0.46

names={'Elite male','Non-elite male','Elite female','Non-elite female'}

figure(867)
pd=bar([1 2 3 4],[q_EM q_NM q_EF q_NF],'w')
set(gca,'XTickLabel',names)
set(gca,'TickDir','out')
set(gca,'Color',[0.95,0.95,0.95])
set(gca,'FontSize',16)
set(pd,'LineWidth',4)
pd=ylabel('Proportion of female postdocs')
set(pd,'FontSize',16)
pd=xlabel('Lab leader')
set(pd,'FontSize',16)

axis([0 5 0 0.5])



%Start with equality!
N_EM(1)=235;
N_NM(1)=800;
N_EF(1)=235;
N_NF(1)=800;

%Start with reality
N_EM(1)=Elite(1);
N_NM(1)=Faculty(1)-Elite(1);
N_EF(1)=Elite(2);
N_NF(1)=Faculty(2)-Elite(2);


%How many more PhDs in elite labs than non-elite.
EliteLabSize=3;
retire=1/6;
totalPostdocplaces=5000;

%Use this if you want to run for lots of values.
erange=[1:0.1:12]
ei=0;
for e=4
ei=ei+1;
    
eliteFactorElite=e;
eliteFactorNormal=e;

for t=1:T
    
    %Recruit postdocs from 50/50 Male, female elite non-elite.
    %EM
    Postdocplaces(1)=N_EM(t)*EliteLabSize;
    %NM
    Postdocplaces(2)=N_NM(t);
    %EF
    Postdocplaces(3)=N_EF(t)*EliteLabSize;
    %NF
    Postdocplaces(4)=N_NF(t);
    
    %Make total number of places constant
    Postdocplaces=round(Postdocplaces*totalPostdocplaces/sum(Postdocplaces));
    
    %Postdocs defined by gender and whether they have been in an elite .
    %Elite male labs
    females_recruited=sum(rand(Postdocplaces(1),1)<q_EM);
    FE(t)=females_recruited;
    ME(t)=Postdocplaces(1)-females_recruited;
    
    %Normal male labs
    females_recruited=sum(rand(Postdocplaces(2),1)<q_NM);
    FN(t)=females_recruited;
    MN(t)=Postdocplaces(2)-females_recruited;
    
    %Elite female labs
    females_recruited=sum(rand(Postdocplaces(3),1)<q_EF);
    FE(t)=FE(t)+females_recruited;
    ME(t)=ME(t)+Postdocplaces(3)-females_recruited;
    
    %Normal female labs
    females_recruited=sum(rand(Postdocplaces(4),1)<q_NF);
    FN(t)=FN(t)+females_recruited;
    MN(t)=MN(t)+Postdocplaces(4)-females_recruited;
    
    %Recruit new faculty, purely on basis of eliteness.
    probF=FE(t)*eliteFactorElite+FN(t);
    probM=ME(t)*eliteFactorElite+MN(t);
    probF=probF/(probF+probM);
 
    store_probFE(t)=probF;
       
    %Assume that becoming elite is total random
    %First retire
    totalEliteplaces=round(N_EF(t)*retire)+round(N_EM(t)*retire);
    %Then recruit
    females_recruited=sum(rand(totalEliteplaces,1)<probF);
    N_EF(t+1)=N_EF(t)-round(N_EF(t)*retire)+females_recruited;
    N_EM(t+1)=N_EM(t)-round(N_EM(t)*retire)+totalEliteplaces-females_recruited;
  
    probF=FE(t)*eliteFactorNormal+FN(t);
    probM=ME(t)*eliteFactorNormal+MN(t);
    probF=probF/(probF+probM);
  
    totalNormalplaces=round(N_NF(t)*retire)+round(N_NM(t)*retire);
    females_recruited=sum(rand(totalNormalplaces,1)<probF);
    N_NF(t+1)=N_NF(t)-round(N_NF(t)*retire)+females_recruited;
    N_NM(t+1)=N_NM(t)-round(N_NM(t)*retire)+totalNormalplaces-females_recruited;

    store_probFN(t)=probF;
   
    
    
end
    ratioE(ei)=N_EM(T)/N_EF(T); 
    ratioN(ei)=N_NM(T)/N_NF(T); 
end

figure(1)
years=2014+[0:T-1]*5
hold off
pd=plot(years,N_NF(1:T),'r')
set(pd,'LineWidth',4)
set(gca,'TickDir','out')
set(gca,'Color',[0.95,0.95,0.95])
set(gca,'FontSize',16)
hold on
pd=plot(years,N_NM(1:T),'b')
set(pd,'LineWidth',4)
pd=plot(years,N_EM(1:T),'b:')
set(pd,'LineWidth',4)
pd=plot(years,N_EF(1:T),'r:')
set(pd,'LineWidth',4)
legend('Normal females','Normal males','Elite males','Elite females')
  
pd=ylabel('Number of faculty')
set(pd,'FontSize',14)
pd=xlabel('Year')
set(pd,'FontSize',14)
axis([2014 2014+T*5 0 1200])


figure(2)
hold off
pd=plot(years,N_NM(1:T)./N_NF(1:T),'k')
set(pd,'LineWidth',4)
set(gca,'TickDir','out')
set(gca,'FontSize',16)
set(gca,'Color',[0.95,0.95,0.95])
hold on
pd=plot(years,N_EM(1:T)./N_EF(1:T),'r')
set(pd,'LineWidth',4)
legend('Non-Elite','Elite')
  
pd=ylabel('Male:Female ratio')
set(pd,'FontSize',14)
pd=xlabel('Year')
set(pd,'FontSize',14)
set(gca,'FontSize',16)
axis([2014 2014+T*5 0 4])

figure(3)
hold off
pd=plot(years,FN(1:T),'r')
set(pd,'LineWidth',4)
hold on
pd=plot(years,MN(1:T),'b')
set(pd,'LineWidth',4)
pd=plot(years,ME(1:T),'b:')
set(pd,'LineWidth',4)
pd=plot(years,FE(1:T),'r:')
set(pd,'LineWidth',4)
legend('Non-elite females','Non-elite males','Non-elite males','Non-elite females')
 



figure(4)
pd=plot(erange,ratioE,'k')
set(pd,'LineWidth',4)
hold on
pd=plot(erange,ratioN,'r')
set(pd,'LineWidth',4)
pd=xlabel('Elite factor')
set(pd,'FontSize',14)
pd=ylabel('Faculty gander ratio')
set(pd,'FontSize',14)

