clear
set more off

use "$temp/final_round_2", clear

* Get march 26 information
********************************************************************************

preserve
use "$temp/final", clear

* Store value of bottom third on March 26
	local third = 1/3*100
	_pctile nbre_hospit_corona_cum_pc if date == mdy(03,26,2020) & prediction_26_log3!=.  , p(`third') 
	local pc33 = r(r1)

* Compute cumulative distribution function on March
local j =  0
forvalues i = 0(0.25)500{
local j = `j' + 1
local pt_`j' =  `i'
count if nbre_hospit_corona_cum_pc <= `i' &  date == mdy(03,26,2020)
local n_march26_`j' = r(N)
}

count if nbre_hospit_corona_cum_pc <= `pc33' &  date == mdy(03,26,2020)
local safe_march26 = r(N)

restore

********************************************************************************

gen threshold33 = `pc33'

gen nbre_hospit_corona_cum_9_pc = (nbre_hospit_corona_cum_9/population)*100000

* Compute cumulative distribution function on July 9
********************************************************************************

local j =  0
forvalues i = 0(0.25)500{
local j = `j' + 1
local pt_`j' =  `i'
*count if nbre_hospit_corona_cum_pc <= `i' &  date == mdy(03,26,2020)
*local n_march26_`j' = r(N)
count if nbre_hospit_corona_cum_9_pc <= `i' &  date == mdy(07,09,2020)
local n_july09_`j' = r(N)
}

count if nbre_hospit_corona_cum_9_pc <= `pc33' &  date == mdy(07,09,2020)
local safe_july09 = r(N)

* List of departements still at risk
********************************************************************************
tab departement_name if nbre_hospit_corona_cum_9_pc > `pc33' &  date == mdy(07,09,2020)

* Plot cumulative distribution functions
********************************************************************************

local text_july09 = 4
count if nbre_hospit_corona_cum_9_pc <= `text_july09' &  date == mdy(07,09,2020)
local  pos_text_july09 = r(N)

* Store data
********************************************************************************
clear
gen pt = .
gen n_march26 = .
gen n_july09 = .
set obs `j'
forvalues i  = 1(1)`j'{
replace pt = `pt_`i'' if _n == `i'
replace n_march26 = `n_march26_`i'' if _n == `i'
replace n_july09 = `n_july09_`i'' if _n == `i'
}

local text_march26 = 35
summ n_march26 if inrange(pt,`text_march26'-1,`text_march26'+1)
*count if nbre_hospit_corona_cum_pc <= `text_march26' &  date == mdy(03,26,2020)
local  pos_text_march26 = r(mean)

replace pt = log(1+pt)

local a = log(1+`pc33')
local x_july09 = log(1+`text_july09')
local x_march26 = log(1+`text_march26')

local a_text : di %2.0fc `pc33'

local y0 = "0"
local y0_v =  log(1+real("`y0'"))

local y1 = "5"
local y1_v =  log(1+real("`y1'"))

local y2 = "10"
local y2_v =  log(1+real("`y2'"))

local y3 = "50"
local y3_v =  log(1+real("`y3'"))

local y4 = "100"
local y4_v =  log(1+real("`y4'"))

local y5 = "500"
local y5_v =  log(1+real("`y5'"))

local y6 = "20"
local y6_v =  log(1+real("`y6'"))

local y7 = "1"
local y7_v =  log(1+real("`y7'"))

local y8 = "250"
local y8_v =  log(1+real("`y8'"))

* March 26 and June 19
********************************************************************************

twoway scatter n_march26 pt, msymbol(none) connect(line) lcolor(gs9) lwidth(medthick) ||  ///
		scatter n_july09 pt, msymbol(none) connect(line) lcolor(gs6) lwidth(medthick) lpattern(longdash) || ///
		scatteri `safe_march26' 0 `safe_march26' `a', connect(line) msymbol(none) lcolor(black) || ///
		scatteri 0 `a' `safe_march26' `a', connect(line) msymbol(none) lcolor(black) || ///
		scatteri `safe_march26' `a' `safe_july09' `a', connect(line) msymbol(none) lcolor(black) lpattern(dash) || ///
		scatteri `safe_july09' `a' `safe_july09' 0, connect(line) msymbol(none) lcolor(black) lpattern(dash) || ///
		, ///
		text(`pos_text_july09' `x_july09' "Hospitalizations" "by July 9, 2020", placement(nw) justification(right) margin(vsmall) size(small)) ///
		text(`pos_text_march26' `x_march26' "Hospitalizations" "by March 26, 2020", placement(se) justification(left) margin(vsmall) size(small)) ///
		ylabel(0(12)72 `safe_july09' 96, angle(0) nogrid) ///
		xlabel(`a' "`a_text'" `y0_v' "`y0'" `y1_v' "`y1'" `y2_v' "`y2'" `y3_v' "`y3'" `y4_v' "`y4'" `y5_v' "`y5'" `y6_v' "`y6'" `y7_v' "`y7'" `y8_v' "`y8'") ///
		graphregion(color(white)) plotregion(margin(zero)) ///
		xtitle(" " "32-day cumulated hospitalizations per 100,000 inhabitants") ///
		ytitle("Number of départements" " ") ///
		legend(off)

		graph export "$graphs/epidemics_march26_and_july09_distributions.eps", replace
		graph export "$graphs/epidemics_march26_and_july09_distributions.png", replace
