text:
# Time series analysis of 2020 worldwide coronavirus outbreak
#
# At first data were taken from Johns Hopkins CSSE repository
# https://github.com/CSSEGISandData/COVID-19
#
# but since this data source for non-understandable reasons on March 23 totally
# changed their data format, without reporting confirmed cases any more
# I changed to datahub
# https://github.com/datasets/covid-19/tree/master/data.

module C

using  Dates, DSP, PyPlot, DifferentialEquations, Polynomials

include("read_data.jl")

function smooth_log(x)
# make log-transform and smooth data
    lx = log10.(x)
    b = (1/2)*ones(2)     # flat (rectangular) window -> moving average
    lx = filtfilt(b,lx)
    10 .^ lx, lx
end

function local_fit(x,y,di)
# local slope of graph y(x) from regression in running window
# window width: 2*di+1
  len = length(x)
  slope = zeros(len)
  for i=1:len
    i1 = max(1,i-di)
    i2 = min(len,i+di)
    #p = polyfit(x[i1:i2],y[i1:i2],1)
    p = fit(x[i1:i2],y[i1:i2],1)
    cp = coeffs(p)
    if length(cp) == 1
        slope[i] = 0.0
    else
      slope[i] = coeffs(p)[2]
    end
  end
  slope, log(10.) .* slope
end

# doubling time when slope is taken from log10-transformed data
doubling_time(slope) = log(2.)/log(10.)/slope

# **********************************************************
function nlys_Asia()
   di = 3

   t_cn,I_cn,R_cn,D_cn = covid_ts("China",         Date(2020,1,22))
   #(D_cn, C_cn, R_cn, dates_cn) =  read_datahub("China",Date(2020,1,22))
   #println(R_cn)
   #t_cn = 1:length(C_cn)
   #I_cn = C_cn #  .- R_cn .- D_cn
   I_cn, lI_cn = smooth_log(I_cn)
   #println(t_cn)
   #println(lI_cn)
   slope_cn,growth_cn = local_fit(t_cn,lI_cn,di)

   di = 2
   t_jp,I_jp,R_jp,D_jp = covid_ts("Japan",         Date(2020,1,24))
   #(D_jp, C_jp, R_jp, dates_jp) =  read_datahub("Japan",Date(2020,1,24))
   #t_jp = 1:length(C_jp)
   #I_jp = C_jp .- R_jp .- D_jp
   I_jp, lI_jp = smooth_log(I_jp)
   slope_jp,growth_jp = local_fit(t_jp,lI_jp,di)

   t_kr,I_kr,R_kr,D_kr = covid_ts("Korea, South",  Date(2020,2,17))
   #(D_kr, C_kr, R_kr, dates_kr) =  read_datahub("Korea, South",Date(2020,2,17))
   #t_kr = 1:length(C_kr)
   #I_kr = C_kr .- R_kr .- D_kr
   I_kr, lI_kr = smooth_log(I_kr)
   slope_kr,growth_kr = local_fit(t_kr,lI_kr,di)

   t_ir,I_ir,R_ir,D_ir = covid_ts("Iran",          Date(2020,2,22))
   #(D_ir, C_ir, R_ir, dates_ir) =  read_datahub("Iran",Date(2020,2,19))
   #t_ir = 1:length(C_ir)
   #I_ir = C_ir .- R_ir .- D_ir
   I_ir, lI_ir = smooth_log(I_ir)
   slope_ir,growth_ir = local_fit(t_ir,lI_ir,di)

  fig=figure("Outbreak Asia",figsize=(12,8))
  clf()
  subplot(3,4,1)         # China
      plot(t_cn,I_cn,"r*")
      #plot(dates_nl,I_nl,"r*")
      ylabel("active cases",fontsize=12)

      title("China,"* today)
  ax1=subplot(3,4,5)
      #semilogy(dates_nl,I_nl,"r*", label=country_it, markersize=8)
      semilogy(t_cn,I_cn,"r*")
      #semilogy(t_nl,10 .^ lI_fit_nl,"k")
      #ax1.xaxis_date()
      #ax1.xaxis.set_major_formatter(majorformatter)
      #ax1.xaxis.set_major_locator(majorlocator)
      #fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      #xlim(Date(2020,2,21),Date(2020,3,26))
      #ax1.xaxis.grid()
      ax1.yaxis.grid()
      ylim(5,5e5)
      ylabel("active cases",labelpad=16,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_cn,growth_cn,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      #ax.yaxis.grid(which="major", color="r", linestyle="-", linewidth=2)
      ax.yaxis.grid()
     # annotate("growth 0.107",xy=[20,2],
     # xycoords="axes fraction"
     # )
      #ax.xaxis.set_major_formatter(majorformatter)
      #fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=20)

  subplot(3,4,2)          # Japan
      plot(t_jp,I_jp,"r*")
      #ylabel("infected Italy")
      title("Japan,"*today)
  ax=subplot(3,4,6)
      semilogy(t_jp,I_jp,"r*")
      #semilogy(t_no,10 .^ lI_fit_no,"k")
      ax.yaxis.grid()
      ylim(5,1e5)
  ax=subplot(3,4,10)
      plot(t_jp,growth_jp,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,3)              # Korea, South
      plot(t_kr,I_kr,"r*")
      #ylabel("infected Spain")
      title("Korea, South,"* today)
  ax=subplot(3,4,7)
      semilogy(t_kr,I_kr,"r*")
     # semilogy(t_ch,10 .^ lI_fit_ch,"k")
      ax.yaxis.grid()
      ylim(5,1e5)
      #ylabel("log10 infected Spain")
  ax= subplot(3,4,11)
      plot(t_kr,growth_kr,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,4)              # Iran
      plot(t_ir,I_ir,"r*")
      #ylabel("infected")
      title("Iran,"* today)
  ax=subplot(3,4,8)
      semilogy(t_ir,I_ir,"r*")
     # semilogy(t_at,10 .^ lI_fit_at,"k")
      ax.yaxis.grid()
      ylim(5,1e5)
      #ylabel("log10 infected")
  ax = subplot(3,4,12)
      plot(t_ir,growth_ir,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")
  PyPlot.tight_layout()
  savefig(folder*"outbreak_Asia.png")

end

# **********************************************************
function nlys_Europe()
   di = 2

   t_de,I_de,R_de,D_de = covid_ts("Germany",       Date(2020,2,23))
   #(D_de, C_de, R_de, dates_de) =  read_datahub("Germany",Date(2020,2,24))
   #t_de = 1:length(C_de)
   #I_de = C_de #.- R_de .- D_de
   I_de, lI_de = smooth_log(I_de)
   slope_de,growth_de = local_fit(t_de,lI_de,di)

   t_it,I_it,R_it,D_it = covid_ts("Italy",         Date(2020,2,21))
   #(D_it, C_it, R_it, dates_it) =  read_datahub("Italy",Date(2020,2,21))
   #t_it = 1:length(C_it)
   #I_it = C_it #.- R_it .- D_it
   I_it, lI_it = smooth_log(I_it)
   slope_it,growth_it = local_fit(t_it,lI_it,di)

   t_fr,I_fr,R_fr,D_fr = covid_ts("France",        Date(2020,2,27))
   #(D_fr, C_fr, R_fr, dates_fr) =  read_datahub("France",Date(2020,2,25))
   #t_fr = 1:length(C_fr)
   #I_fr = C_fr #.- R_fr .- D_fr
   I_fr, lI_fr = smooth_log(I_fr)
   slope_fr,growth_fr = local_fit(t_fr,lI_fr,di)

   t_es,I_es,R_es,D_es = covid_ts("Spain",         Date(2020,2,25))
   #(D_es, C_es, R_es, dates_es) =  read_datahub("Spain",Date(2020,2,25))
   #t_es = 1:length(C_es)
   #I_es = C_es #.- R_es .- D_es
   I_es, lI_es = smooth_log(I_es)
   slope_es,growth_es = local_fit(t_es,lI_es,di)

  fig=figure("Outbreak Europe",figsize=(12,8))
  clf()
  subplot(3,4,1)         # Germany
      plot(t_de,I_de,"r*")
      #plot(dates_nl,I_nl,"r*")
      ylabel("active cases",fontsize=12)

      title("Germany,"* today)
  ax1=subplot(3,4,5)
      #semilogy(dates_nl,I_nl,"r*", label=country_it, markersize=8)
      semilogy(t_de,I_de,"r*")
      #semilogy(t_nl,10 .^ lI_fit_nl,"k")
      #ax1.xaxis_date()
      #ax1.xaxis.set_major_formatter(majorformatter)
      #ax1.xaxis.set_major_locator(majorlocator)
      #fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      #xlim(Date(2020,2,21),Date(2020,3,26))
      #ax1.xaxis.grid()
      ax1.yaxis.grid()
      ylim(5,2e5)
      ylabel("active cases",labelpad=16,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_de,growth_de,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      #ax.yaxis.grid(which="major", color="r", linestyle="-", linewidth=2)
      ax.yaxis.grid()
     # annotate("growth 0.107",xy=[20,2],
     # xycoords="axes fraction"
     # )
      #ax.xaxis.set_major_formatter(majorformatter)
      #fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=20)

  subplot(3,4,2)          # Italy
      plot(t_it,I_it,"r*")
      #ylabel("infected Italy")
      title("Italy,"*today)
  ax=subplot(3,4,6)
      semilogy(t_it,I_it,"r*")
      #semilogy(t_no,10 .^ lI_fit_no,"k")
      ax.yaxis.grid()
      ylim(5,2e5)
  ax=subplot(3,4,10)
      plot(t_it,growth_it,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,3)              # France
      plot(t_fr,I_fr,"r*")
      #ylabel("infected Spain")
      title("France,"* today)
  ax=subplot(3,4,7)
      semilogy(t_fr,I_fr,"r*")
     # semilogy(t_ch,10 .^ lI_fit_ch,"k")
      ax.yaxis.grid()
      ylim(5,2e5)
      #ylabel("log10 infected Spain")
  ax= subplot(3,4,11)
      plot(t_fr,growth_fr,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,4)              # Spain
      plot(t_es,I_es,"r*")
      #ylabel("infected")
      title("Spain,"* today)
  ax=subplot(3,4,8)
      semilogy(t_es,I_es,"r*")
     # semilogy(t_at,10 .^ lI_fit_at,"k")
      ax.yaxis.grid()
      ylim(5,2e5)
      #ylabel("log10 infected")
  ax = subplot(3,4,12)
      plot(t_es,growth_es,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")
  PyPlot.tight_layout()
  savefig(folder*"outbreak_Europe.png")

end

# **********************************************************
# **********************************************************
# **********************************************************
function nlys_World()
   di = 2

   t_us,I_us,R_us,D_us = covid_ts("US",            Date(2020,2,22))
   I_us, lI_us = smooth_log(I_us)
   slope_us,growth_us = local_fit(t_us,lI_us,di)

   t_ca,I_ca,R_ca,D_ca = covid_ts("Canada",        Date(2020,2,27))
   I_ca, lI_ca = smooth_log(I_ca)
   slope_ca,growth_ca = local_fit(t_ca,lI_ca,di)

   t_uk,I_uk,R_uk,D_uk = covid_ts("United Kingdom",Date(2020,2,27))
   I_uk, lI_uk = smooth_log(I_uk)
   slope_uk,growth_uk = local_fit(t_uk,lI_uk,di)

   t_il,I_il,R_il,D_il = covid_ts("Israel",        Date(2020,2,27))
   I_il, lI_il = smooth_log(I_il)
   slope_il,growth_il = local_fit(t_il,lI_il,di)

  fig=figure("Outbreak World",figsize=(12,8))
  clf()
  subplot(3,4,1)         # US
      plot(t_us,I_us,"r*")
      #plot(dates_nl,I_nl,"r*")
      ylabel("active cases",fontsize=12)

      title("US,"* today)
  ax1=subplot(3,4,5)
      semilogy(t_us,I_us,"r*")
      ax1.yaxis.grid()
      ylim(5,5e6)
      ylabel("active cases",labelpad=16,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_us,growth_us,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=20)

  subplot(3,4,2)          # Canada
      plot(t_ca,I_ca,"r*")
      title("Canada,"*today)
  ax=subplot(3,4,6)
      semilogy(t_ca,I_ca,"r*")
      ax.yaxis.grid()
      ylim(5,5e4)
  ax=subplot(3,4,10)
      plot(t_ca,growth_ca,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,3)              # United Kingdom
      plot(t_uk,I_uk,"r*")
      title("United Kingdom,"* today)
  ax=subplot(3,4,7)
      semilogy(t_uk,I_uk,"r*")
      ax.yaxis.grid()
      ylim(5,5e5)
  ax= subplot(3,4,11)
      plot(t_uk,growth_uk,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,4)              # Israel
      plot(t_il,I_il,"r*")
      title("Israel,"* today)
  ax=subplot(3,4,8)
      semilogy(t_il,I_il,"r*")
      ax.yaxis.grid()
      ylim(5,5e4)
  ax = subplot(3,4,12)
      plot(t_il,growth_il,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
  PyPlot.tight_layout()
  savefig(folder*"outbreak_World.png")

end

# **********************************************************
# **********************************************************
# **********************************************************
function nlys_SouthAmerica()
   di = 2

   t_br,I_br,R_br,D_br = covid_ts("Brazil",        Date(2020,2,22))
   I_br, lI_br = smooth_log(I_br)
   slope_br,growth_br = local_fit(t_br,lI_br,di)

   t_cl,I_cl,R_cl,D_cl = covid_ts("Chile",        Date(2020,2,27))
   I_cl, lI_cl = smooth_log(I_cl)
   slope_cl,growth_cl = local_fit(t_cl,lI_cl,di)

   t_ar,I_ar,R_ar,D_ar = covid_ts("Argentina",Date(2020,2,27))
   I_ar, lI_ar = smooth_log(I_ar)
   slope_ar,growth_ar = local_fit(t_ar,lI_ar,di)

   t_ec,I_ec,R_ec,D_ec = covid_ts("Ecuador",        Date(2020,2,27))
   I_ec, lI_ec = smooth_log(I_ec)
   slope_ec,growth_ec = local_fit(t_ec,lI_ec,di)

  fig=figure("Outbreak South America",figsize=(12,8))
  clf()
  subplot(3,4,1)         # Brazil
      plot(t_br,I_br,"r*")
      ylabel("active cases",fontsize=12)
      title("Brazil,"* today)
  ax1=subplot(3,4,5)
      semilogy(t_br,I_br,"r*")
      ax1.yaxis.grid()
      ylim(5,5e5)
      ylabel("active cases",labelpad=16,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_br,growth_br,"r*-")
      ylim(0,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=20)

  subplot(3,4,2)          # Chile
      plot(t_cl,I_cl,"r*")
      title("Chile,"*today)
  ax=subplot(3,4,6)
      semilogy(t_cl,I_cl,"r*")
      ax.yaxis.grid()
      ylim(5,5e4)
  ax=subplot(3,4,10)
      plot(t_cl,growth_cl,"r*-")
      ylim(0,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,3)              # Argentina
      plot(t_ar,I_ar,"r*")
      title("Argentina,"* today)
  ax=subplot(3,4,7)
      semilogy(t_ar,I_ar,"r*")
      ax.yaxis.grid()
      ylim(5,1e5)
  ax= subplot(3,4,11)
      plot(t_ar,growth_ar,"r*-")
      ylim(0,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,4)              # Ecuador
      plot(t_ec,I_ec,"r*")
      title("Ecuador,"* today)
  ax=subplot(3,4,8)
      semilogy(t_ec,I_ec,"r*")
      ax.yaxis.grid()
      ylim(5,5e4)
  ax = subplot(3,4,12)
      plot(t_ec,growth_ec,"r*-")
      ylim(0,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
  PyPlot.tight_layout()
  savefig(folder*"outbreak_SouthAmerica.png")

end

# **********************************************************
# **********************************************************
# **********************************************************
function nlys_Africa()
   di = 2

   t_eg,I_eg,R_eg,D_eg = covid_ts("Egypt",        Date(2020,3,8))
   I_eg, lI_eg = smooth_log(I_eg)
   slope_eg,growth_eg = local_fit(t_eg,lI_eg,di)

   t_sa,I_sa,R_sa,D_sa = covid_ts("South Africa",        Date(2020,3,8))
   I_sa, lI_sa = smooth_log(I_sa)
   slope_sa,growth_sa = local_fit(t_sa,lI_sa,di)

   t_al,I_al,R_al,D_al = covid_ts("Algeria",Date(2020,3,8))
   I_al, lI_al = smooth_log(I_al)
   slope_al,growth_al = local_fit(t_al,lI_al,di)

   t_ca,I_ca,R_ca,D_ca = covid_ts("Cameroon",   Date(2020,3,12))
   I_ca, lI_ca = smooth_log(I_ca)
   slope_ca,growth_ca = local_fit(t_ca,lI_ca,di)

  fig=figure("Africa",figsize=(12,8))
  clf()
  subplot(3,4,1)         # Egypt
      plot(t_eg,I_eg,"r*")
      ylabel("active cases",fontsize=12)
      title("Egypt,"* today)
  ax1=subplot(3,4,5)
      semilogy(t_eg,I_eg,"r*")
      ax1.yaxis.grid()
      ylim(5,5e5)
      ylabel("active cases",labelpad=16,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_eg,growth_eg,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=20)

  subplot(3,4,2)          # South Africa
      plot(t_sa,I_sa,"r*")
      title("South Africa,"*today)
  ax=subplot(3,4,6)
      semilogy(t_sa,I_sa,"r*")
      ax.yaxis.grid()
      ylim(5,5e4)
  ax=subplot(3,4,10)
      plot(t_sa,growth_sa,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,3)              # Algeria
      plot(t_al,I_al,"r*")
      title("Algeria,"* today)
  ax=subplot(3,4,7)
      semilogy(t_al,I_al,"r*")
      ax.yaxis.grid()
      ylim(5,1e5)
  ax= subplot(3,4,11)
      plot(t_al,growth_al,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,4)              # Cameroon
      plot(t_ca,I_ca,"r*")
      title("Cameroon,"* today)
  ax=subplot(3,4,8)
      semilogy(t_ca,I_ca,"r*")
      ax.yaxis.grid()
      ylim(5,1e4)
  ax = subplot(3,4,12)
      plot(t_ca,growth_ca,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
  PyPlot.tight_layout()
  savefig(folder*"outbreak_Africa.png")

end

# **********************************************************
# **********************************************************
# **********************************************************
function nlys_Russia()
   di = 2

   t_ru,I_ru,R_ru,D_ru = covid_ts("Russia",        Date(2020,2,29))
   I_ru, lI_ru = smooth_log(I_ru)
   slope_ru,growth_ru = local_fit(t_ru,lI_ru,di)

   t_cz,I_cz,R_cz,D_cz = covid_ts("Czechia",        Date(2020,3,12))
   I_cz, lI_cz = smooth_log(I_cz)
   slope_cz,growth_cz = local_fit(t_cz,lI_cz,di)

   t_hu,I_hu,R_hu,D_hu = covid_ts("Hungary",        Date(2020,3,5))
   I_hu, lI_hu = smooth_log(I_hu)
   slope_hu,growth_hu = local_fit(t_hu,lI_hu,di)

   t_pl,I_pl,R_pl,D_pl = covid_ts("Poland",Date(2020,2,27))
   I_pl, lI_pl = smooth_log(I_pl)
   slope_pl,growth_pl = local_fit(t_pl,lI_pl,di)

  fig=figure("Outbreak Russia",figsize=(12,8))
  clf()
  subplot(3,4,1)         # Russia
      plot(t_ru,I_ru,"r*")
      ylabel("active cases",fontsize=12)
      title("Russia,"* today)
  ax1=subplot(3,4,5)
      semilogy(t_ru,I_ru,"r*")
      ax1.yaxis.grid()
      ylim(5,5e5)
      ylabel("active cases",labelpad=16,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_ru,growth_ru,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=20)

  subplot(3,4,2)          # Czechia
      plot(t_cz,I_cz,"r*")
      title("Czechia,"*today)
  ax=subplot(3,4,6)
      semilogy(t_cz,I_cz,"r*")
      ax.yaxis.grid()
      ylim(5,5e4)
  ax=subplot(3,4,10)
      plot(t_cz,growth_cz,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,3)              # Hungary
      plot(t_hu,I_hu,"r*")
      title("Hungary,"* today)
  ax=subplot(3,4,7)
      semilogy(t_hu,I_hu,"r*")
      ax.yaxis.grid()
      ylim(5,1e5)
  ax= subplot(3,4,11)
      plot(t_hu,growth_hu,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,4)              # Poland
      plot(t_pl,I_pl,"r*")
      title("Poland,"* today)
  ax=subplot(3,4,8)
      semilogy(t_pl,I_pl,"r*")
      ax.yaxis.grid()
      ylim(5,1e5)
  ax = subplot(3,4,12)
      plot(t_pl,growth_pl,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
  PyPlot.tight_layout()
  savefig(folder*"outbreak_Russia.png")

end

# **********************************************************
# **********************************************************
# **********************************************************
function nlys_Australia()
   di = 2

   t_au,I_au,R_au,D_au = covid_ts("Australia",        Date(2020,2,29))
   I_au, lI_au = smooth_log(I_au)
   slope_au,growth_au = local_fit(t_au,lI_au,di)

   t_nz,I_nz,R_nz,D_nz = covid_ts("New Zealand",        Date(2020,3,12))
   I_nz, lI_nz = smooth_log(I_nz)
   slope_nz,growth_nz = local_fit(t_nz,lI_nz,di)

   t_ph,I_ph,R_ph,D_ph = covid_ts("Philippines",        Date(2020,3,5))
   I_ph, lI_ph = smooth_log(I_ph)
   slope_ph,growth_ph = local_fit(t_ph,lI_ph,di)

   t_ml,I_ml,R_ml,D_ml = covid_ts("Malaysia",Date(2020,2,27))
   I_ml, lI_ml = smooth_log(I_ml)
   slope_ml,growth_ml = local_fit(t_ml,lI_ml,di)

  fig=figure("Outbreak Australia",figsize=(12,8))
  clf()
  subplot(3,4,1)         # Australia
      plot(t_au,I_au,"r*")
      ylabel("active cases",fontsize=12)
      title("Australia,"* today)
  ax1=subplot(3,4,5)
      semilogy(t_au,I_au,"r*")
      ax1.yaxis.grid()
      ylim(5,5e5)
      ylabel("active cases",labelpad=16,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_au,growth_au,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=20)

  subplot(3,4,2)          # New Zealand
      plot(t_nz,I_nz,"r*")
      title("New Zealand,"*today)
  ax=subplot(3,4,6)
      semilogy(t_nz,I_nz,"r*")
      ax.yaxis.grid()
      ylim(5,5e4)
  ax=subplot(3,4,10)
      plot(t_nz,growth_nz,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,3)              # Philippines
      plot(t_ph,I_ph,"r*")
      title("Philippines,"* today)
  ax=subplot(3,4,7)
      semilogy(t_ph,I_ph,"r*")
      ax.yaxis.grid()
      ylim(5,1e5)
  ax= subplot(3,4,11)
      plot(t_ph,growth_ph,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)

  subplot(3,4,4)              # Malaysia
      plot(t_ml,I_ml,"r*")
      title("Malaysia,"* today)
  ax=subplot(3,4,8)
      semilogy(t_ml,I_ml,"r*")
      ax.yaxis.grid()
      ylim(5,1e4)
  ax = subplot(3,4,12)
      plot(t_ml,growth_ml,"r*-")
      ylim(-0.1,0.5)
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
  PyPlot.tight_layout()
  savefig(folder*"outbreak_Australia.png")

end

# **********************************************************
# **********************************************************
# **********************************************************
# **********************************************************
function nlys_smallEurope()
   di = 2

   t_nl,I_nl,R_nl,D_nl = covid_ts("Netherlands",   Date(2020,2,23))
   #(D_nl, C_nl, R_nl, dates_nl) =  read_datahub("Netherlands",Date(2020,2,23))
   #t_nl = 1:length(C_nl)
   #I_nl = C_nl .- R_nl .- D_nl
   I_nl, lI_nl = smooth_log(I_nl)
   slope_nl,growth_nl = local_fit(t_nl,lI_nl,di)

   t_no,I_no,R_no,D_no = covid_ts("Norway",        Date(2020,2,28))
   #(D_no, C_no, R_no, dates_no) =  read_datahub("Norway",Date(2020,2,28))
   #t_no = 1:length(C_no)
   #I_no = C_no .- R_no .- D_no
   I_no, lI_no = smooth_log(I_no)
   slope_no,growth_no = local_fit(t_no,lI_no,di)

   t_ch,I_ch,R_ch,D_ch = covid_ts("Switzerland",   Date(2020,2,27))
   #(D_ch, C_ch, R_ch, dates_ch) =  read_datahub("Switzerland",Date(2020,2,25))
   #t_ch = 1:length(C_ch)
   #I_ch = C_ch .- R_ch .- D_ch
   I_ch, lI_ch = smooth_log(I_ch)
   slope_ch,growth_ch = local_fit(t_ch,lI_ch,di)

   t_at,I_at,R_at,D_at = covid_ts("Austria",       Date(2020,2,29))
   #(D_at, C_at, R_at, dates_at) =  read_datahub("Austria",Date(2020,2,29))
   #t_at = 1:length(C_at)
   #I_at = C_at .- R_at .- D_at
   I_at, lI_at = smooth_log(I_at)
   slope_at,growth_at = local_fit(t_at,lI_at,di)

   t_dk,I_dk,R_dk,D_dk = covid_ts("Denmark",        Date(2020,2,27))
   #(D_dk, C_dk, R_dk, dates_dk) =  read_datahub("Denmark",Date(2020,2,27))
   #t_dk = 1:length(C_dk)
   #I_dk = C_dk .- R_dk .- D_dk
   I_dk, lI_dk = smooth_log(I_dk)
   slope_dk,growth_dk = local_fit(t_dk,lI_dk,di)

   t_se,I_se,R_se,D_se = covid_ts("Sweden",        Date(2020,2,27))
   #(D_dk, C_dk, R_dk, dates_dk) =  read_datahub("Denmark",Date(2020,2,27))
   #t_dk = 1:length(C_dk)
   #I_dk = C_dk .- R_dk .- D_dk
   I_se, lI_se = smooth_log(I_se)
   slope_se,growth_se = local_fit(t_se,lI_se,di)

   t_tr,I_tr,R_tr,D_tr = covid_ts("Turkey",        Date(2020,2,27))
   #(D_dk, C_dk, R_dk, dates_dk) =  read_datahub("Denmark",Date(2020,2,27))
   #t_dk = 1:length(C_dk)
   #I_dk = C_dk .- R_dk .- D_dk
   I_tr, lI_tr = smooth_log(I_tr)
   slope_tr,growth_tr = local_fit(t_tr,lI_tr,di)

   t_fi,I_fi,R_fi,D_fi = covid_ts("Finland",        Date(2020,2,27))
   #(D_dk, C_dk, R_dk, dates_dk) =  read_datahub("Denmark",Date(2020,2,27))
   #t_dk = 1:length(C_dk)
   #I_dk = C_dk .- R_dk .- D_dk
   I_fi, lI_fi = smooth_log(I_fi)
   slope_fi,growth_fi = local_fit(t_fi,lI_fi,di)

# ********************************************
  majorformatter = matplotlib.dates.DateFormatter("%d.%m.%Y")
  #majorformatter = matplotlib.dates.DateFormatter("%d.%m")
  majorlocator = matplotlib.dates.DayLocator(interval=7)
  fig=figure("Outbreak Small Europe",figsize=(12,8))
  clf()
  subplot(3,4,1)         # Netherlands
      plot(t_nl,I_nl,"r*")
      #plot(dates_nl,I_nl,"r*")
      ylabel("active cases",fontsize=12)

      title("Netherlands,"* today)
  ax1=subplot(3,4,5)
      #semilogy(dates_nl,I_nl,"r*", label=country_it, markersize=8)
      semilogy(t_nl,I_nl,"r*")
      #semilogy(t_nl,10 .^ lI_fit_nl,"k")
      #ax1.xaxis_date()
      #ax1.xaxis.set_major_formatter(majorformatter)
      #ax1.xaxis.set_major_locator(majorlocator)
      #fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      #xlim(Date(2020,2,21),Date(2020,3,26))
      #ax1.xaxis.grid()
      ax1.yaxis.grid()
      ylim(5,5e4)
      ylabel("active cases",labelpad=16,fontsize=12)
   ax = subplot(3,4,9)
      plot(t_nl,growth_nl,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      #ax.yaxis.grid(which="major", color="r", linestyle="-", linewidth=2)
      ax.yaxis.grid()
     # annotate("growth 0.107",xy=[20,2],
     # xycoords="axes fraction"
     # )
      #ax.xaxis.set_major_formatter(majorformatter)
      #fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      xlabel("time (days)",fontsize=12)
      ylabel("growth",fontsize=12,labelpad=20)

  subplot(3,4,2)          # Turkey
      plot(t_tr,I_tr,"r*")
      #ylabel("infected Italy")
      title("Turkey,"*today)
  ax=subplot(3,4,6)
      semilogy(t_tr,I_tr,"r*")
      #semilogy(t_dk,10 .^ lI_fit_dk,"k")
      ax.yaxis.grid()
      ylim(5,1e5)
  ax=subplot(3,4,10)
      plot(t_tr,growth_tr,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,3)              # Switzerland
      plot(t_ch,I_ch,"r*")
      #ylabel("infected Spain")
      title("Switzerland,"* today)
  ax=subplot(3,4,7)
      semilogy(t_ch,I_ch,"r*")
     # semilogy(t_ch,10 .^ lI_fit_ch,"k")
      ax.yaxis.grid()
      ylim(5,5e4)
      #ylabel("log10 infected Spain")
  ax= subplot(3,4,11)
      plot(t_ch,growth_ch,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,4)              # Austria
      plot(t_at,I_at,"r*")
      #ylabel("infected")
      title("Austria,"* today)
  ax=subplot(3,4,8)
      semilogy(t_at,I_at,"r*")
     # semilogy(t_at,10 .^ lI_fit_at,"k")
      ax.yaxis.grid()
      ylim(5,5e4)
      #ylabel("log10 infected")
  ax = subplot(3,4,12)
      plot(t_at,growth_at,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")
  PyPlot.tight_layout()
  savefig(folder*"outbreak_europe1.png")

# ********************************************
  fig=figure("Outbreak Small Europe2",figsize=(12,8))
  clf()

  subplot(3,4,1)          # Sweden
      plot(t_se,I_se,"r*")
      #ylabel("infected Italy")
      title("Sweden,"*today)
  ax=subplot(3,4,5)
      semilogy(t_se,I_se,"r*")
      #semilogy(t_dk,10 .^ lI_fit_dk,"k")
      ax.yaxis.grid()
      ylim(5,1e5)
  ax=subplot(3,4,9)
      plot(t_se,growth_se,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,2)          # Denmark
      plot(t_dk,I_dk,"r*")
      #ylabel("infected Italy")
      title("Denmark,"*today)
  ax=subplot(3,4,6)
      semilogy(t_dk,I_dk,"r*")
      #semilogy(t_dk,10 .^ lI_fit_dk,"k")
      ax.yaxis.grid()
      ylim(5,1e4)
  ax=subplot(3,4,10)
      plot(t_dk,growth_dk,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,3)          # Norway
      plot(t_no,I_no,"r*")
      #ylabel("infected Italy")
      title("Norway,"*today)
  ax=subplot(3,4,7)
      semilogy(t_no,I_no,"r*")
      #semilogy(t_no,10 .^ lI_fit_no,"k")
      ax.yaxis.grid()
      ylim(5,1e4)
  ax=subplot(3,4,11)
      plot(t_no,growth_no,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  subplot(3,4,4)          # Finland
      plot(t_fi,I_fi,"r*")
      #ylabel("infected Italy")
      title("Finland,"*today)
  ax=subplot(3,4,8)
      semilogy(t_fi,I_fi,"r*")
      #semilogy(t_no,10 .^ lI_fit_no,"k")
      ax.yaxis.grid()
      ylim(5,1e4)
  ax=subplot(3,4,12)
      plot(t_fi,growth_fi,"r*-")
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2,0.25])
      ax.yaxis.grid()
      xlabel("time (days)",fontsize=12)
      #ylabel("growth")

  PyPlot.tight_layout()
  savefig(folder*"outbreak_europe2.png")
end

# **********************************************************
# **********************************************************
# **********************************************************
function nlys_ItDe()

  t_it,I_it,R_it,D_it = covid_ts("Italy",Date(2020,2,21))
  #(D_it, C_it, R_it, dates_it) =  read_datahub("Italy",Date(2020,2,21))
  #(D_it, I_it, dates_it) =  read_datahub("Italy",Date(2020,2,21))
  #t_it = 1:length(I_it)
  #I_it = C_it .- R_it .- D_it
  dates_it = Date(2020,2,21) .+ Day.(t_it .- 1)
  I_it, lI_it = smooth_log(I_it)
  #lI_it = log10.(I_it)
  di = 2
  slope_it,growth_it = local_fit(t_it,lI_it,di)

  p_tot_it1 = fit(t_it[3:8],lI_it[3:8],1)  # fit slope 1st phase Italy
  coef = coeffs(p_tot_it1)[2]
  println("slope Italy Phase 1 ",coef, " T12 ", doubling_time(coef))
  #lI_fit_it1 = polyval(p_tot_it1, t_it)
  lI_fit_it1 = p_tot_it1.(t_it)

  p_tot_it2 = fit(t_it[9:22],lI_it[9:22],1)  # fit slope 2nd phase Italy
  coef = coeffs(p_tot_it2)[2]
  println("slope Italy Phase 2 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it2 = p_tot_it2.(t_it)

  p_tot_it3 = fit(t_it[23:35],lI_it[23:35],1)  # fit slope 2nd phase Italy
  coef = coeffs(p_tot_it3)[2]
  println("slope Italy Phase 3 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it3 = p_tot_it3.(t_it)

  # ********************
  t_de,I_de,R_de,D_de = covid_ts("Germany",  Date(2020,2,24))
  #(D_de, C_de, R_de, dates_de) =  read_datahub("Germany",Date(2020,2,24))
  #(D_de, I_de, dates_de) =  read_datahub("Germany",Date(2020,2,24))
  #t_de = 1:length(I_de)
  #I_de = C_de .- R_de .- D_de
  dates_de = Date(2020,2,24) .+ Day.(t_de .- 1 )
  I_de, lI_de = smooth_log(I_de)
  #lI_de = log10.(I_de)
  di = 2
  slope_de,growth_de = local_fit(t_de,lI_de,di)
  t_de = t_de .+ 3 # bring to the same time axis

  p_tot_de1 = fit(t_de[4:11],lI_de[4:11],1)  # fit full slope Germany
  coef = coeffs(p_tot_de1)[2]
  println("slope Germany ",coef, " T12 ", doubling_time(coef))
  lI_fit_de1= p_tot_de1.(t_de)

  p_tot_de2 = fit(t_de[12:25],lI_de[12:25],1)  # fit full slope Germany
  coef = coeffs(p_tot_de2)[2]
  println("slope Germany ",coef, " T12 ", doubling_time(coef))
  lI_fit_de2= p_tot_de2.(t_de)

  p_tot_de3 = fit(t_de[26:34],lI_de[26:34],1)  # fit full slope Germany
  coef = coeffs(p_tot_de3)[2]
  println("slope Germany ",coef, " T12 ", doubling_time(coef))
  lI_fit_de3= p_tot_de3.(t_de)

  p_tot_de4 = fit(t_de[35:end],lI_de[35:end],1)  # fit full slope Germany
  coef = coeffs(p_tot_de4)[2]
  println("slope Germany ",coef, " T12 ", doubling_time(coef))
  lI_fit_de4= p_tot_de4.(t_de)

   majorformatter = matplotlib.dates.DateFormatter("%d.%m.%Y")
   #majorformatter = matplotlib.dates.DateFormatter("%d.%m")
   majorlocator = matplotlib.dates.DayLocator(interval=7)
   lang = "en"
   if lang == "en"
       country_de = "Germany"
       country_it = "Italy"
   else
       country_de = "Deutschland"
       country_it = "Italien"
   end
   fig=figure("Germany_Italy", figsize=(12,6))
   clf()
   subplot(1,2,1)
      #semilogy(t_it[1:12],10 .^ lI_fit_it1[1:12],"k")
      #semilogy(t_it[5:25],10 .^ lI_fit_it2[5:25],"k")
      #semilogy(t_it[19:28],10 .^ lI_fit_it3[19:28],"k")
      #semilogy(t_it,I_it,"b*", label="Italy", markersize=8)
      #semilogy(t_de[1:15],10 .^ lI_fit_de1[1:15],"k")
      #semilogy(t_de[8:end],10 .^ lI_fit_de2[8:end],"k")
      #semilogy(t_de,I_de,"r*",label="Germany",markersize=8)
      semilogy(dates_it[1:12],10 .^ lI_fit_it1[1:12],"k")
      semilogy(dates_it[5:25],10 .^ lI_fit_it2[5:25],"k")
      #semilogy(dates_it[19:37],10 .^ lI_fit_it3[19:37],"k")
      semilogy(dates_it,I_it,"b*", label=country_it, markersize=8)
      #semilogy(dates_it,I_it,"b*", label="Italy", markersize=8)
      semilogy(dates_de[1:15],10 .^ lI_fit_de1[1:15],"k")
      semilogy(dates_de[8:25],10 .^ lI_fit_de2[8:25],"k")
      semilogy(dates_de[23:35],10 .^ lI_fit_de3[23:35],"k")
      #semilogy(dates_de[32:end],10 .^ lI_fit_de4[32:end],"k")
      #semilogy(dates_de,I_de,"r*",label="Germany",markersize=8)
      semilogy(dates_de,I_de,"r*",label=country_de,markersize=8)
      ax1 = gca()
      ax1.xaxis_date()
      ax1.xaxis.set_major_formatter(majorformatter)
      ax1.xaxis.set_major_locator(majorlocator)
      fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      xlim(Date(2020,2,21),today_date)
      ylim(8,2e5)
      ax1.xaxis.grid()
      ax1.yaxis.grid()
      if lang == "en"
       # ylabel("active cases",fontsize=14)
        xlabel("time",fontsize=14)
        #xlabel("time (Feb 21 - March 19)", fontsize=14)
      else
        #xlabel("Zeit", fontsize=14)
        ylabel("aktive Fälle",fontsize=14)
      end
      #xlabel("time (Feb 21 - March 19)", fontsize=14)
      legend(loc="upper left",fancybox="true")
   subplot(1,2,2)
      #plot(t_it,growth_it,"b*-",label="Italy",markersize=8)
      #plot(t_de,growth_de,"r*-",label="Germany",markersize=8)
        plot(dates_it,growth_it,"b*-",label=country_it,markersize=8)
        plot(dates_de,growth_de,"r*-",label=country_de,markersize=8)
      #plot([10,10],[0,0.6],"grey",linewidth=1)
      #ylim(0,0.22)
      #yticks([0,0.05,0.1,0.15,0.2])
      ax1 = gca()
      ax1.xaxis_date()
      ax1.xaxis.set_major_formatter(majorformatter)
      ax1.xaxis.set_major_locator(majorlocator)
      fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      ylim(-0.1,0.5)
      xlim(Date(2020,2,21),today_date)
      #yticks([0,0.05,0.1,0.15,0.2])
      ax1.xaxis.grid()
      ax1.yaxis.grid()
      if lang == "en"
        #xlabel("time",fontsize=14)
        ylabel("growth rate",fontsize=14)
        #xlabel("time (Feb 21 - March 19)", fontsize=14)
      else
        #xlabel("Zeit", fontsize=14)
        ylabel("Wachstumsrate",fontsize=14)
      end
      legend(loc="upper right",fancybox="true")
  PyPlot.tight_layout()
  savefig(folder*"ItalyGermany_March.png")
  #savefig("../figures/ItalienDeutschland_March.png")
end

# **********************************************************
function nlys_ItUS()
  #(D_it, C_it, R_it, dates_it) =  read_datahub("Italy",Date(2020,2,22))
  #(D_it, I_it, dates_it) =  read_datahub("Italy",Date(2020,2,22))
  #t_it = 1:length(I_it)
  t_it,I_it,R_it,D_it = covid_ts("Italy",Date(2020,2,21))
  #I_it = C_it .- R_it .- D_it
  dates_it = Date(2020,2,21) .+ Day.(t_it .- 1)
  I_it, lI_it = smooth_log(I_it)
  #lI_it = log10.(I_it)
  di = 2
  slope_it,growth_it = local_fit(t_it,lI_it,di)

  p_tot_it1 = fit(t_it[2:7],lI_it[2:7],1)  # fit slope 1st phase Italy
  coef = coeffs(p_tot_it1)[2]
  println("slope Italy Phase 1 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it1 = p_tot_it1.(t_it)

  p_tot_it2 = fit(t_it[8:21],lI_it[8:21],1)  # fit slope 2nd phase Italy
  coef = coeffs(p_tot_it2)[2]
  println("slope Italy Phase 2 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it2 = p_tot_it2.(t_it)

  p_tot_it3 = fit(t_it[22:end],lI_it[22:end],1)  # fit slope 2nd phase Italy
  coef = coeffs(p_tot_it3)[2]
  println("slope Italy Phase 3 ",coef, " T12 ", doubling_time(coef))
  lI_fit_it3 = p_tot_it3.(t_it)

  # ********************  US
  t_us,I_us,R_us,D_us = covid_ts("US",            Date(2020,2,22))
  #(D_us, C_us, R_us, dates_us) =  read_datahub("US",Date(2020,2,22))
  #(D_us, I_us, dates_us) =  read_datahub("US",Date(2020,2,22))
  #t_us = 1:length(I_us)
  #I_us = C_us # for the US we just take the confirmed cases
  dates_us = Date(2020,2,22) .+ Day.(t_us .- 1 )
  I_us, lI_us = smooth_log(I_us)
  #lI_us = log10.(I_us)
  p_tot_us = fit(t_us[11:26],lI_us[11:26],1)  # fit full slope USA
  lI_fit_us= p_tot_us.(t_us)
  coef = coeffs(p_tot_us)[2]
  slope_us, growth_us = local_fit(t_us,lI_us,di)
  println("slope US ",coef, " T12 ", doubling_time(coef))

   majorformatter = matplotlib.dates.DateFormatter("%d.%m.%Y")
   #majorformatter = matplotlib.dates.DateFormatter("%d.%m")
   majorlocator = matplotlib.dates.DayLocator(interval=7)

   fig = figure("ItalyUS", figsize=(12,6))
   clf()
   subplot(1,2,1)
      #semilogy(t_it[1:12],10 .^ lI_fit_it1[1:12],"k")
      #semilogy(t_it[5:25],10 .^ lI_fit_it2[5:25],"k")
      #semilogy(t_it[19:28],10 .^ lI_fit_it3[19:28],"k")
      #semilogy(t_it,I_it,"b*", label="Italy", markersize=8)
      #semilogy(t_us,I_us,"r*",label="US",markersize=8)
      #semilogy(t_us[8:27],10 .^ lI_fit_us[8:27],"k")
      semilogy(dates_it[1:11],10 .^ lI_fit_it1[1:11],"k")
      semilogy(dates_it[4:24],10 .^ lI_fit_it2[4:24],"k")
      semilogy(dates_it[18:31],10 .^ lI_fit_it3[18:31],"k")
      semilogy(dates_it,I_it,"b*", label="Italy", markersize=8)
      semilogy(dates_us,I_us,"r*",label="US",markersize=8)
      semilogy(dates_us[8:30],10 .^ lI_fit_us[8:30],"k")
      ax1 = gca()
      ax1.xaxis_date()
      ax1.xaxis.set_major_formatter(majorformatter)
      ax1.xaxis.set_major_locator(majorlocator)
      ax1.xaxis.grid()
      ax1.yaxis.grid()
      fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      xlim(Date(2020,2,21),today_date)
      ylim(8,1e6)
      ylabel("active cases",fontsize=14)
      #xlabel("time (Feb 21 - March 19)", fontsize=14)
      legend(loc="upper left",fancybox="true")
   subplot(1,2,2)
      #plot(t_it,growth_it,"b*-",label="Italy",markersize=8)
      #plot(t_us,growth_us,"r*-",label="US",markersize=8)
      plot(dates_it,growth_it,"b*-",label="Italy",markersize=8)
      plot(dates_us,growth_us,"r*-",label="US",markersize=8)
      #plot([10,10],[0,0.6],"grey",linewidth=1)
      #ylim(0,0.22)
      #yticks([0,0.05,0.1,0.15,0.2])
      ax1 = gca()
      ax1.xaxis_date()
      ax1.xaxis.set_major_formatter(majorformatter)
      ax1.xaxis.set_major_locator(majorlocator)
      ax1.xaxis.grid()
      ax1.yaxis.grid()
      fig.autofmt_xdate(bottom=0.2,rotation=30,ha="right")
      xlim(Date(2020,2,21),today_date)
      ylim(-0.1,0.5)
      #yticks([0,0.05,0.1,0.15,0.2])
      #xlabel("time (days)",fontsize=14)
      #xlabel("time (Feb 21 - March 19)", fontsize=14)
      ylabel("growth rate",fontsize=14)
      legend(loc="upper right",fancybox="true")
  PyPlot.tight_layout()
  savefig(folder*"/ItalyUS_March.png")

end

# ********************************************
function corona_scaling()

   t_it,I_it,R_it,D_it = covid_ts("Italy",Date(2020,2,20))
   I_it, lI_it = smooth_log(I_it)
   lt_it = log10.(t_it)
   p_tot_it = fit(lt_it,lI_it,1)  # fit full slope Germany
   println("loglog It ", coeffs(p_tot_it)[2])
   I_fit_it = 10 .^ p_tot_it.(lt_it)

   #lt_de = log10.(t_de)
   #p_tot_de = fit(lt_de,lI_de,1)  # fit full slope Germany
   #println("loglog It ", coeffs(p_tot_de)[2])
   #I_fit_de = 10 .^ polyval(p_tot_de, lt_de)

   figure("scaling")
   clf()
   #subplot(2,1,1)
   title("Power-law scaling of COVID-19 outbreak in Italy")
   loglog(t_it,I_it,"*-", markersize=8)
   loglog(t_it,I_fit_it,"-k", label="regression line, slope=2.75")
   ax = gca()
   ax.yaxis.grid()
   ax.xaxis.grid()
   xlabel("time (days)",fontsize=12)
   ylabel("active cases (Italy)",fontsize=12)
   legend(loc="upper left",fancybox="true")
   savefig(folder*"Italy_scaling.png")

end

#nlys_country(country, date)
today = " August 19"
today_date = Date(2020,8,20)
folder = "../figures/"

nlys_Asia()
nlys_Europe()
nlys_World()
nlys_SouthAmerica()
nlys_Africa()
nlys_Russia()
nlys_Australia()
nlys_smallEurope()
nlys_ItDe()
#nlys_ItUS()

#corona_scaling()

# Russia, Poland
# Indonesia, Phillippines, Malaysia
# India, Pakistan

end   # module


text:
function _make_radii(a::Number, b::Number, c::Number)
    dab = abs(a - b)
    dbc = abs(b - c)
    dac = abs(a - c)
    rhs = [dab; dbc; dac]
    A = [1 1 0; 0 1 1; 1 0 1]

    radii = A \ rhs
end

"""
    kiss(z1,z2,z3)

Given three points (as complex numbers) find three circles
centered at those points that are mutually tangent.

See: `touch_points`
"""
function kiss(a::Number, b::Number, c::Number)
    rads = _make_radii(a, b, c)
    ctrs = [a, b, c]

    [Circle(ctrs[k], rads[k]) for k = 1:3]
end

export kiss, touch_points

"""
    touch_points(z1,z2,z3)

Given three points (as complex numbers) determine the points of tangency
of three mutually tangent circles centered at those points.

See: `kiss`
"""
function touch_points(a::Number, b::Number, c::Number)
    radii = _make_radii(a, b, c)
    ra, rb, rc = radii

    ab = (b - a) * (ra / (ra + rb)) + a
    ac = (c - a) * (ra / (ra + rc)) + a
    bc = (c - b) * (rb / (rb + rc)) + b
    return [ab; ac; bc]
end

"""
    soddy_radii(r1, r2, r3)
Given three radii (real numbers) of mutually tangent circles,
give the radii of a fourth circle that is mutually tangent to the
three that are given as a 2-tuple of real numbers.

The positive result (first entry) is the radius of a small circle
nestled between the other three, and the negative result (second entry)
is negative the radius of a circle that is tangent to and surrounds the
three.
"""
function soddy_radii(r1, r2, r3)
    k1 = 1 / r1
    k2 = 1 / r2
    k3 = 1 / r3

    common = 2 * sqrt(k1 * k2 + k2 * k3 + k3 * k1)

    k4a = k1 + k2 + k3 + common
    k4b = k1 + k2 + k3 - common

    return 1 / k4a, 1 / k4b
end

soddy_radii(C1::Circle, C2::Circle, C3::Circle) = soddy_radii(C1.r, C2.r, C3.r)

function soddy_center(C1::Circle, C2::Circle, C3::Circle)
    r = soddy_radii(C1, C2, C3)[1]

    D1 = Circle(center(C1), radius(C1) + r)
    D2 = Circle(center(C2), radius(C2) + r)

    zz = collect(D1 ∩ D2)

    d1 = abs(center(C3) - zz[1]) - (radius(C3) + r) |> abs
    d2 = abs(center(C3) - zz[2]) - (radius(C3) + r) |> abs

    d1 < d2 ? zz[1] : zz[2]

end

export soddy

"""
    soddy(C1,C2,C3)
Given three mutually tangent circles, return the circle that is
tangent to all three an nesteled between them.
"""
function soddy(C1::Circle, C2::Circle, C3::Circle)
    r = soddy_radii(C1, C2, C3)[1]
    z = soddy_center(C1, C2, C3)
    return Circle(z, r)
end

ext: jl


text:
# # [Stochastic Gradient Descent](@id SGDTutorial)
#
# This tutorial illustrates how to use the [`stochastic_gradient_descent`](@ref)
# solver and different [`DirectionUpdateRule`](@ref)s in order to introduce
# the average or momentum variant, see [Stochastic Gradient Descent](https://en.wikipedia.org/wiki/Stochastic_gradient_descent).
#
# Computationally we look at a very simple but large scale problem,
# the Riemannian Center of Mass or [Fréchet mean](https://en.wikipedia.org/wiki/Fréchet_mean):
# For given points ``p_i ∈\mathcal M``, ``i=1,…,N`` this optimization problem reads
#
# ```math
# \operatorname*{arg\,min}_{x∈\mathcal M} \frac{1}{2}\sum_{i=1}^{N}
#   \operatorname{d}^2_{\mathcal M}(x,p_i),
# ```
# which of course can be (and is) solved by a gradient descent, see the [introductionary tutorial](@ref Optimize).
# If ``N`` is very large it might be quite expensive to evaluate the complete gradient.
# A remedy is, to evaluate only one of the terms at a time and choose a random order for these.
#
# We first initialize the manifold (see [])
export_folder = joinpath( #src
    @__DIR__, #src
    "..", #src
    "..", #src
    "docs", #src
    "src", #src
    "assets", #src
    "images", #src
    "tutorials", #src
) #src
using Manopt, Manifolds, Random, Colors
# and we define some colors from [Paul Tol](https://personal.sron.nl/~pault/)
black = RGBA{Float64}(colorant"#000000")
TolVibrantOrange = RGBA{Float64}(colorant"#EE7733") # Start
TolVibrantBlue = RGBA{Float64}(colorant"#0077BB") # a path
TolVibrantTeal = RGBA{Float64}(colorant"#009988") # points
nothing #hide
#
# And optain a large data set
n = 5000
σ = π / 12
M = Sphere(2)
x = 1 / sqrt(2) * [1.0, 0.0, 1.0]
Random.seed!(42)
data = [exp(M, x, random_tangent(M, x, Val(:Gaussian), σ)) for i in 1:n]
nothing #hide
# which looks like
#
asymptote_export_S2_signals( #src
    export_folder * "/centerAndLargeData.asy"; #src
    points=[[x], data], #src
    colors=Dict(:points => [TolVibrantOrange, TolVibrantTeal]), #src
    dot_sizes=[2.5, 1.25], #src
    camera_position=(1.0, 0.5, 0.5), #src
) #src
render_asymptote(export_folder * "/centerAndLargeData.asy"; render=2) #src
#md # ```julia
#md # asymptote_export_S2_signals("centerAndLargeData.asy";
#md #     points = [ [x], data],
#md #     colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal]),
#md #     dot_sizes = [2.5, 1.0], camera_position = (1.,.5,.5)
#md # )
#md # render_asymptote("centerAndLargeData.asy"; render = 2)
#md # ```
#md #
#md # ![The data of noisy versions of $x$](../assets/images/tutorials/centerAndLargeData.png)
#
# Note that due to the construction of the points as zero mean tangent vectors, the mean should
# be very close to our initial point `x`.
#
# In order to use the stochastic gradient, we now need a function that returns the vector of gradients.
# There are two ways to define it in `Manopt.jl`: as one function, that returns a vector or a vector of funtions.
#
# The first variant is of course easier to define, but the second is more efficient when only evaluating one of the gradients.
# For the mean we have as a gradient
#
# ```math
#  gradF(x) = \sum_{i=1}^N \operatorname{grad}f_i(x) \quad \text{where} \operatorname{grad}f_i(x) = -\log_x p_i
# ```
#
# Which we define as
F(M, x) = 1 / (2 * n) * sum(map(p -> distance(M, x, p)^2, data))
gradF(M, x) = [grad_distance(M, p, x) for p in data]
gradf = [(M, x) -> grad_distance(M, p, x) for p in data];
# The calls are only slightly different, but notice that accessing the 2nd gradient element
# requires evaluating all logs in the first function.
# So while you can use both `gradF` and `gradf` in the following call, the second one is faster:
@time x_opt1 = stochastic_gradient_descent(M, gradF, x);
# versus
@time x_opt2 = stochastic_gradient_descent(M, gradf, x);
# This result is reasonably close. But we can improve it by using a [`DirectionUpdateRule`](@ref),
# namely:
# On the one hand [`MomentumGradient`](@ref), which requires both the manifold and the initial value,
#    in order to keep track of the iterate and parallel transport the last direction to the current iterate.
#    you can also set a `vector_transport_method`, if `ParallelTransport()` is not
#    available on your manifold. Here we simply do
@time x_opt3 = stochastic_gradient_descent(
    M, gradf, x; direction=MomentumGradient(M, x, StochasticGradient(zero_vector(M, x)))
);
# And on the other hand the [`AverageGradient`](@ref) computes an average of the last `n` gradients, i.e.
@time x_opt4 = stochastic_gradient_descent(
    M, gradf, x; direction=AverageGradient(M, x, 10, StochasticGradient(zero_vector(M, x)))
);
# note that the default [`StoppingCriterion`](@ref) is a fixed number of iterations.
#
# Note that since you can apply both also in case of [`gradient_descent`](@ref),
# i.e. to use [`IdentityUpdateRule`](@ref) and evaluate the classical gradient,
# both constructors have to know that internally the default evaluation of the Stochastic gradient
# (choosing one gradient $\operatorname{grad}f_k$ at random) has to be specified.
#
# For this small example you can of course also use a gradient descent with [`ArmijoLinesearch`](@ref),
# but it will be a little slower usually
@time x_opt5 = gradient_descent(
    M, F, (M, x) -> sum(gradF(M, x)), x; stepsize=ArmijoLinesearch()
);
# but it is for sure faster than the variant above that evaluates the full gradient on every iteration,
# since stochastic gradient descent takes more iterations.
#
# Note that all 5 of couse yield the same result
[distance(M, x, y) for y in [x_opt1, x_opt2, x_opt3, x_opt4, x_opt5]]


text:
"""
``RotationAngle3D`` is an abstract type to be used
as the first argument for functions.
It represents the rotation angle
around X, Y and Z axis.
"""
abstract RotationAngle3D


text:
function SpearmansCorrelation()
    return SpearmansCorrelation(())
end

function SpearmansCorrelation(arg0::RankingAlgorithm)
    return SpearmansCorrelation((RankingAlgorithm,), arg0)
end

function SpearmansCorrelation(arg0::RealMatrix)
    return SpearmansCorrelation((RealMatrix,), arg0)
end

function SpearmansCorrelation(arg0::RealMatrix, arg1::RankingAlgorithm)
    return SpearmansCorrelation((RealMatrix, RankingAlgorithm), arg0, arg1)
end

function compute_correlation_matrix(obj::SpearmansCorrelation, arg0::RealMatrix)
    return jcall(obj, "computeCorrelationMatrix", RealMatrix, (RealMatrix,), arg0)
end

function compute_correlation_matrix(obj::SpearmansCorrelation, arg0::Vector{Vector{jdouble}})
    return jcall(obj, "computeCorrelationMatrix", RealMatrix, (Vector{Vector{jdouble}},), arg0)
end

function correlation(obj::SpearmansCorrelation, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "correlation", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function get_correlation_matrix(obj::SpearmansCorrelation)
    return jcall(obj, "getCorrelationMatrix", RealMatrix, ())
end

function get_rank_correlation(obj::SpearmansCorrelation)
    return jcall(obj, "getRankCorrelation", PearsonsCorrelation, ())
end


text:
let
    srand(100)

    function f_s(x::Vector)
        (x[1] - 5.0)^4
    end

    function rosenbrock_s(x::Vector)
        (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    end

    initial_x = [0.0]
    upper = [100.0]
    lower = [-100.0]
    n_particles = 4

    res = Optim.optimize(f_s, initial_x, method=ParticleSwarm(lower, upper, n_particles),
                         iterations=100)
    @assert norm(Optim.minimizer(res) - [5.0]) < 0.1

    initial_x = [0.0, 0.0]
    lower = [-20., -20.]
    upper = [20., 20.]
    n_particles = 5
    res = Optim.optimize(rosenbrock_s, initial_x, method=ParticleSwarm(lower, upper, n_particles),
                             iterations=300)
    @assert norm(Optim.minimizer(res) - [1.0, 1.0]) < 0.1

    # Add UnconstrainedProblems here; currently they take too many iterations to be
    # feasible
end


text:
# This file contains functions for the Russell DEA model
"""
    RussellDEAModel
An data structure representing a Russell DEA model.
"""
struct RussellDEAModel <: AbstractTechnicalDEAModel
    n::Int64
    m::Int64
    s::Int64
    orient::Symbol
    rts::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    thetaX::Matrix
    thetaY::Matrix
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    dearussell(X, Y)
Compute the Russell model using data envelopment analysis for inputs X and outputs Y.

# Optional Arguments
- `orient=:Input`: chooses the Russell input mode. For the Russell output model choose `:Output`. For the Russell graph model choose `:Graph`.
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `slack=true`: computes input and output slacks.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> dearussell(X, Y, rts = :VRS)
Russell DEA Model
DMUs = 11; Inputs = 2; Outputs = 1
Orientation = Input; Returns to Scale = VRS
───────────────────────────────────────────────
    efficiency     effX1     effX2      slackY1
───────────────────────────────────────────────
1     1.0       1.0       1.0       0.0
2     0.79386   0.587719  1.0       9.88357e-15
3     1.0       1.0       1.0       0.0
4     1.0       1.0       1.0       0.0
5     0.603175  0.277778  0.928571  4.0
6     1.0       1.0       1.0       3.78178e-16
7     1.0       1.0       1.0       0.0
8     1.0       1.0       1.0       0.0
9     1.0       1.0       1.0       0.0
10    0.485746  0.563492  0.408     0.0
11    0.882353  1.0       0.764706  0.0
───────────────────────────────────────────────
```
"""
function dearussell(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    orient::Symbol = :Input, rts::Symbol = :CRS, slack::Bool = true,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix,Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::RussellDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        if orient == :Input || orient == :Output
            optimizer = DEAOptimizer(:LP)
        else
            optimizer = DEAOptimizer(:NLP)
        end
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    thetaXi = zeros(n, m)
    thetaYi = zeros(n, s)
    lambdaeff = spzeros(n, nref)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Create the optimization model
        if orient == :Input
            # Input orientation
            deamodel = newdeamodel(optimizer)

            mno0 = sum(x0 .!= 0)

            @variable(deamodel, lambda[1:nref] >= 0)
            @variable(deamodel, theta[1:m] <= 1)
            @objective(deamodel, Min, 1 / mno0 * sum(theta[t] for t in 1:m if x0[t] != 0 ))

            @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) <= theta[j] * x0[j])
            @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) >= y0[j])
        elseif orient == :Output
            # Output orientation
            deamodel = newdeamodel(optimizer)

            sno0 = sum(y0 .!= 0)

            @variable(deamodel, lambda[1:nref] >= 0)
            @variable(deamodel, theta[1:s] >= 1)
            @objective(deamodel, Max, 1 / sno0 * sum(theta[t] for t in 1:s if y0[t] != 0 ))

            @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) <= x0[j])
            @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) >= theta[j] * y0[j])
        elseif orient == :Graph
            # Graph orientation
            deamodel = newdeamodel(optimizer)
            set_silent(deamodel)

            mno0 = sum(x0 .!= 0)
            sno0 = sum(y0 .!= 0)

            @variable(deamodel, lambda[1:nref] >= 0)
            @variable(deamodel, theta[1:m] <= 1)
            @variable(deamodel, phi[1:s] >= 1)

            @NLobjective(deamodel, Min, 1 / (mno0 + sno0) * (sum(theta[t] for t in 1:m if x0[t] != 0 ) + sum(1/phi[t] for t in 1:s if y0[t] != 0) ))

            @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) == theta[j] * x0[j])
            @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) == phi[j] * y0[j])
        else
            throw(ArgumentError("`orient` must be :Input, :Output or :Graph"));
        end

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            throw(ArgumentError("`rts` must be :CRS or :VRS"));
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        effi[i]  = JuMP.objective_value(deamodel)
        lambdaeff[i,:] = JuMP.value.(lambda)

        if orient == :Input
            thetaXi[i,:] = JuMP.value.(theta)
        elseif orient == :Output
            thetaYi[i,:] = JuMP.value.(theta)
        elseif orient == :Graph
            thetaXi[i,:] = JuMP.value.(theta)
            thetaYi[i,:] = JuMP.value.(phi)
        end

        # Check termination status
        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

    end

    # Get first-stage X and Y targets
    if orient == :Input
        Xtarget = X .* thetaXi
        Ytarget = Y
    elseif orient == :Output
        Xtarget = X
        Ytarget = Y .* thetaYi
    elseif orient == :Graph
        Xtarget = X .* thetaXi
        Ytarget = Y .* thetaYi
    end

    # Compute slacks
    if (slack == true) && (orient != :Graph)

        # Use additive model with Russell efficient X and Y to get slacks
        russellSlacks = deaadd(Xtarget, Ytarget, :Ones, rts = rts, Xref = Xref, Yref = Yref, optimizer = optimizer)
        slackX = slacks(russellSlacks, :X)
        slackY = slacks(russellSlacks, :Y)

        # Get second-stage X and Y targets
        Xtarget = Xtarget - slackX
        Ytarget = Ytarget + slackY
    else
        if typeof(Xtarget) <: AbstractVector    Xtarget = Xtarget[:,:]  end
        if typeof(Ytarget) <: AbstractVector    Ytarget = Ytarget[:,:]  end

        slackX = Array{Float64}(undef, 0, 0)
        slackY = Array{Float64}(undef, 0, 0)
    end

    if orient == :Input
        thetaYi = Array{Float64}(undef, 0, 0)
    elseif orient == :Output
        thetaXi = Array{Float64}(undef, 0, 0)
    end

    return RussellDEAModel(n, m, s, orient, rts, names, effi, thetaXi, thetaYi, slackX, slackY, lambdaeff, Xtarget, Ytarget)

end

function Base.show(io::IO, x::RussellDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Russell DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if x.orient == :Input
            thetaX = x.thetaX
            if hasslacks == true
                show(io, CoefTable(hcat(eff, thetaX, slackY), ["efficiency"; ["effX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
            else
                show(io, CoefTable(hcat(eff, thetaX), ["efficiency"; ["effX$i" for i in 1:m ] ], dmunames))
            end
        elseif x.orient == :Output
            thetaY = x.thetaY
            if hasslacks == true
                show(io, CoefTable(hcat(eff, thetaY, slackX), ["efficiency"; ["effY$i" for i in 1:s ]; ["slackX$i" for i in 1:m ]], dmunames))
            else
                show(io, CoefTable(hcat(eff, thetaY), ["efficiency"; ["effY$i" for i in 1:s ] ], dmunames))
            end
        elseif x.orient == :Graph
            thetaX = x.thetaX
            thetaY = x.thetaY
            show(io, CoefTable(hcat(eff, thetaX, thetaY), ["efficiency"; ["effX$i" for i in 1:m ]; ["effY$i" for i in 1:s ] ], dmunames))
        end
    end

end

function efficiency(model::RussellDEAModel, type::Symbol)::Matrix

    if (type == :X && (model.orient == :Input  || model.orient == :Graph)) return model.thetaX end
    if (type == :Y && (model.orient == :Output || model.orient == :Graph)) return model.thetaY end

    throw(ArgumentError("$(typeof(model)) with orienation $(model.orient) has no efficiency $(type)"));

end


text:
using BusinessDays
using Test

# lets find out minimum and maximum month for easter
function find_easter_min_max()
    y = 1582
    easter_1582 = BusinessDays.easter_date(Dates.Year(1582))
    month_min = Dates.month(easter_1582)
    date_min = easter_1582
    month_max = Dates.month(easter_1582)
    date_max = easter_1582

    while y <= 2100
        e_date = BusinessDays.easter_date(Dates.Year(y))
        m_e_date = Dates.month(e_date)
        month_min = min(month_min, m_e_date)
        month_max = max(month_max, m_e_date)

        if month_min == m_e_date
            date_min = e_date
        end

        if month_max == m_e_date
            date_max = e_date
        end

        y += 1
    end

    return (month_min, month_max)
end

month_min, month_max = find_easter_min_max()

@test month_min >= 3
@test month_max <= 4


text:
@doc raw"""
    velocity_curve(p, c::Curve, backend)

Velocity curve for a given curve `p` of type `c` calculated using differentiation
defined by `backend`.
If $c$ is a function such that $c\colon [0, 1] \to M$, then
`velocity_curve(c)(t)` is a vector tangent to `c` at `t`.
"""
function velocity_curve(M::AbstractManifold, p, backend::AbstractRiemannianDiffBackend)
    return t -> Manifolds.differential(M, p, t, backend)
end

function velocity_curve(M::AbstractManifold, p, backend::ProjectedDifferenceBackend)
    return t -> project(M.M, p(t), (embed(p(t + backend.dt)) - embed(p(t))) / backend.dt)
end

@doc raw"""
    curve_length(p, c::FunctionCurveSpace, a = 0.0, b = 1.0, dtype = Val(:continuous))

Returns length of the curve `c` between parameters `a` and `b`. By default
`a = 0.0` and `b = 1.0`. Calculates velocity using method `dtype`.
"""
function curve_length(
    M::AbstractManifold,
    p,
    a=0.0,
    b=1.0,
    backend=Manifolds.rdifferential_backend(),
)
    v = velocity_curve(M, p, backend)
    #TODO: add quadrature tolerances to curve_length?
    #TODO: add velocity calculation options
    i_val, error = quadgk(t -> norm(M, p(t), v(t)), a, b)
    return i_val
end

"""
    q_function(M::AbstractManifold, p, X)

Rescales given tangent vector `X` at point `p` to square root of its original length.
"""
function q_function(M::AbstractManifold, p, X)
    norm_X = norm(M, p, X)
    return if norm_X ≈ 0
        return zero_vector(M, p)
    else
        return (1 / sqrt(norm_X)) * X
    end
end

"""
    srvf(M, f, backend)

Calculate the Square Root Velocity Function of a given curve `f`.
Parameter `backend` describes type of differentiaion used. Supported values:
* `Val(:continuous)` -- symbolic differentiation resulting in a continuous
curve,
* `Val(:forward_diff)` -- forward finite differentiaion of a discretized curve
`c` based on its grid.
"""
function srvf(M::AbstractManifold, f, backend::AbstractRiemannianDiffBackend)
    vel = velocity_curve(M, f, backend)
    return t -> q_function(M, f(t), vel(t))
end

"""
    transport_srvf(M::AbstractManifold, c_p, c_X, p)

Transports all tangent vectors `c_X` of curve `c_p` to a given
point `p` along shortest geodesics.
`M` is the respective manifold of curves.
"""
function transport_srvf(M::AbstractManifold, c_p, c_X, p)
    c_out = allocate_result(M, transport_srvf, c_X, c_p)
    transport_srvf!(M, c_out, c_p, c_X, p)
    return c_out
end

"""
    tsrvf(M::AbstractManifold, c, p, backend::AbstractRiemannianDiffBackend)

Calculate SRVF of curve `c` from manifold `M` using method `backend` and transport to `p`.
"""
function tsrvf(M::AbstractManifold, c, p, backend::AbstractRiemannianDiffBackend)
    c_X = srvf(M, c, backend)
    return transport_srvf(M, c, c_X, p)
end

function reverse_srvf(M::AbstractManifold, c_X, initial_point)
    c_out = allocate_result(M, reverse_srvf, c_X)
    reverse_srvf!(M, c_out, c_X, initial_point)
    return c_out
end

"""
    srsf(M::CurveWarpingSpace, p)

Square Root Slope Function of curve warping `p`.
"""
function srsf(M::CurveWarpingSpace, p)
    grid = M.knots
    coeffs = map(p, grid)
    grads = [
        (coeffs[i + 1] - coeffs[i]) / (grid[i + 1] - grid[i]) for i in 1:(length(grid) - 1)
    ]
    push!(grads, grads[end])
    vals = map(sqrt, grads)
    return project(CurveWarpingSRSFSpace(grid), vals)
end

"""
    reverse_srsf(c)

Reverse Square Root Slope Function of a given SRSF of a curve warping.
"""
function reverse_srsf(M::CurveWarpingSpace, p)
    grid = M.knots
    steps = diff(grid)
    ys = accumulate(+, [zero(eltype(p)); steps .* p[1:(end - 1)] .^ 2])
    ys ./= ys[end]
    return make_warping(M, ys)
end
function reverse_srsf(M::CurveWarpingSpace, p::Identity{WarpingCompositionOperation})
    return reverse_srsf(M, identity_element(CurveWarpingGroup(M), p))
end


text:
function optimize!(structure::HorizontalStructure, optimizer::AbstractStructuredOptimizer, x₀::AbstractVector)
    # Sanity check
    supports_structure(optimizer, structure) || throw(UnsupportedStructure{typeof(optimizer), typeof(structure)}())
    # Load structure
    load_structure!(optimizer, structure, x₀)
    # Run structure-exploiting optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end

function set_master_optimizer!(structure::HorizontalStructure, optimizer)
    return nothing
end

function set_master_optimizer_attribute!(::HorizontalStructure, ::MOI.AbstractOptimizerAttribute, value)
    return nothing
end

function set_subproblem_optimizer!(structure::HorizontalStructure, optimizer)
    set_optimizer!(scenarioproblems(structure), optimizer)
    return nothing
end

function set_subproblem_optimizer_attribute!(structure::HorizontalStructure, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(scenarioproblems(structure), attr, value)
    return nothing
end

function cache_solution!(stochasticprogram::StochasticProgram{2}, structure::HorizontalStructure{2}, optimizer::MOI.AbstractOptimizer)
    cache = solutioncache(stochasticprogram)
    # Cache main solution
    variables = decision_variables_at_stage(stochasticprogram, 1)
    constraints = decision_constraints_at_stage(stochasticprogram, 1)
    cache[:solution] = SolutionCache(optimizer, variables, constraints)
    # Cache first-stage solution
    cache[:node_solution_1] = SolutionCache(optimizer, variables, constraints)
    # Cache scenario-dependent solutions (Skip if more than 100 scenarios for performance)
    if num_scenarios(stochasticprogram) <= 100
        cache_solution!(cache, scenarioproblems(structure), optimizer, 2)
    end
    return nothing
end


text:
struct CalibratedCI
    est_num::Float64
    est_num_linear::EmpiricalBayesEstimator #LinearEstimator
    est_denom::Float64
    est_denom_linear::EmpiricalBayesEstimator #LinearEstimator
    est_target::Float64
    calibrated_target::Float64
    ci_left::Float64
    ci_right::Float64
    std::Float64
    max_bias::Float64
    f::BinnedMarginalDensityNeighborhood
    ma::MinimaxCalibrator
end

function estimate(L::LinearEstimator, Xs)
    mean(L.(Xs))
end

function estimate(L::LinearEstimator, target::LinearInferenceTarget, Xs)
    # TODO: Add check that target is correct target
    mean(L.(Xs))
end
# provides lightweight CI functionality for linear functionals
function donoho_ci(Xs, ma::MinimaxCalibrator; conf=0.9)
    est = mean(ma.Q.(Xs))
    sd = std(ma.Q.(Xs))/sqrt(length(Xs))
    max_bias = ma.max_bias
    zz =  get_plus_minus(max_bias, sd, conf)
    (est, est-zz, est+zz)
end

function CEB_ci(Xs_train, Xs_test,
            ds::MixingNormalConvolutionProblem,
            target::PosteriorTarget,
            M_max_num::EmpiricalBayesEstimator, #
            M_max_denom::EmpiricalBayesEstimator=M_max_num;
            f_nb::BinnedMarginalDensityNeighborhood = fit(BinnedMarginalDensityNeighborhood, Xs_train, ds),
            C=:auto, conf=0.9, kwargs...)

    m_train = length(Xs_train)
    m_test = length(Xs_test)

    marginal_grid = ds.marginal_grid
    marginal_h = ds.marginal_h

    num_res = estimate(M_max_num, target.num, Xs_train)
    denom_res = estimate(M_max_denom, target.denom, Xs_train)

    #TODO: Check if denominator not sure to be >0... abort or throw warning
    est_target = num_res/denom_res

    calib_target = CalibratedNumerator(target.num, est_target)
            # Test: Use the Donoho calibrator on the learned function

    #TODO : Change to KWarg..

    if C==:auto
        C = f_nb.C_std*(1+f_nb.η_infl) + f_nb.C_bias
    end

    ma = MinimaxCalibrator(ds, f_nb.f, m_test, calib_target; C=C, kwargs...)

    QXs =  ma.(Xs_test)
    sd = std(QXs)/sqrt(m_test)
    max_bias = ma.max_bias
    zz =  get_plus_minus(max_bias, sd, conf)

    zz = zz/denom_res
    calib_target = est_target + mean(QXs)/denom_res
    ci_left = calib_target - zz
    ci_right = calib_target + zz

    don_ci = CalibratedCI(num_res,
                M_max_num,
                denom_res,
                M_max_denom,
                est_target,
                calib_target,
                ci_left,
                ci_right,
                sd/denom_res[1],
                max_bias/denom_res[1],
                f_nb,
                ma)

    return don_ci
end

function CEB_ci(Xs_train, Xs_test,
            ds::MixingNormalConvolutionProblem,
            target::PosteriorTarget;
            C=:auto,
            conf=0.9,
            kwargs...)
    # fix this
    m_train = length(Xs_train)
    m_test = length(Xs_test)

    marginal_grid = ds.marginal_grid
    marginal_h = ds.marginal_h

    M_bd = marginal_h/sqrt(2*pi)

    f_const = BinnedMarginalDensity(M_bd, marginal_grid, marginal_h)
    # Train
    # Estimate the numerator
    M_max_num = MinimaxCalibrator(ds, f_const, m_train, target.num;
                 C=Inf, kwargs...);

    # Estimate the denominator
    M_max_denom = MinimaxCalibrator(ds, f_const, m_train, target.denom;
                C=Inf, kwargs...);

    CEB_ci(Xs_train, Xs_test, ds, target,
                M_max_num,
                M_max_denom;
                C=C, conf=conf, kwargs...)

end

# TODO: Delete this, instead clean up
function CEB_ci_cb(Xs_train, Xs_test,
            ds::MixingNormalConvolutionProblem, target::PosteriorTarget, C_bias;
            C=:auto, conf=0.9, kwargs...)

            m_train = length(Xs_train)
            m_test = length(Xs_test)

            marginal_grid = ds.marginal_grid
            marginal_h = ds.marginal_h

            cb_num = ComteButucea( target.num, length(Xs_train), marginal_grid)
            cb_denom = ComteButucea( target.denom, length(Xs_train), marginal_grid)

            num_res = estimate(Xs_train, ComteButucea, target.num, marginal_grid)
            denom_res = estimate(Xs_train, ComteButucea, target.denom, marginal_grid)

            if denom_res <= 0
                denom_res = 1e-10
                num_res =0
            end
                    #TODO: Check if denominator not sure to be >0... abort or throw warning
            est_target = num_res[1]/denom_res[1]

            calib_target = CalibratedNumerator(target.num, est_target)
                    # Test: Use the Donoho calibrator on the learned function

            #TODO : Change to KWarg..
            f_nb = fit(BinnedMarginalDensityNeighborhood, Xs_train, marginal_grid)
            f_nb.C_bias = C_bias

            if C==:auto
                C = f_nb.C_std + f_nb.C_bias
            end

            ma = MinimaxCalibrator(ds, f_nb.f, m_test, calib_target; C=C, kwargs...)

            QXs =  ma.(Xs_test)
            sd = std(QXs)/sqrt(m_test)
            max_bias = ma.max_bias
            zz =  get_plus_minus(max_bias, sd, conf)

            zz = zz/denom_res[1]
            calib_target = est_target + mean(QXs)/denom_res[1]
            ci_left = calib_target - zz
            ci_right = calib_target + zz

            don_ci = CalibratedCI(num_res[1],
                        cb_num,
                        denom_res[1],
                        cb_denom,
                        est_target,
                        calib_target,
                        ci_left,
                        ci_right,
                        sd/denom_res[1],
                        max_bias/denom_res[1],
                        f_nb,
                        ma)

            return don_ci
end

function StatsBase.confint(ci::CalibratedCI, target::PosteriorTarget)
    (ci.ci_left, ci.ci_right)
end

function estimate(ci::CalibratedCI, target::PosteriorTarget)
    ci.calibrated_target
end

function maxbias(ci::CalibratedCI, target)
    ci.max_bias
end


text:
using YOLO, Test

enable_info() = Base.CoreLogging.disable_logging(Base.CoreLogging.LogLevel(-1)) #Enable printing of @info messages
disable_info() = Base.CoreLogging.disable_logging(Base.CoreLogging.Info) #disable printing of @info messages

disable_info()
datadir = joinpath(dirname(dirname(pathof(YOLO))), "data")
datasetsdir = joinpath(datadir, "datasets")
pretraineddir = joinpath(datadir, "pretrained")

if !isdir(joinpath(datasetsdir, "voc2007", "VOCdevkit", "VOC2007"))
    @info "Downloading VOC dataset for testing..."
    YOLO.download_dataset("voc2007")
end

@testset "Example datasets have been downloaded" begin
    @test isdir(joinpath(datasetsdir, "voc2007", "VOCdevkit", "VOC2007"))
end

@testset "Example pretrained weights files have been downloaded" begin
    @test isfile(joinpath(
        pretraineddir,
        "v2_tiny",
        "voc2007",
        "v2_tiny_voc.weights",
    ))
end

@testset "Loading and running YOLOv2_tiny_voc pretrained model" begin

    num_images = 2

    voc = YOLO.datasets.VOC.populate()
    @test length(voc.image_paths) == 5011
    @test length(voc.label_paths) == 5011
    @test length(voc.objects) == 20

    settings = YOLO.pretrained.v2_tiny_voc.load(minibatch_size = num_images)
    @test settings.num_classes == 20

    vocloaded = YOLO.load(voc, settings, indexes = collect(1:num_images))
    @test size(vocloaded.imstack_mat) == (416, 416, 3, num_images)
    @test length(vocloaded.paddings) == num_images
    @test length(vocloaded.labels) == num_images
    #Checks that the VOC download hasn't changed
    @test vec(sum(vocloaded.imstack_mat, dims = (1, 2, 3))) ≈ [
        140752.47,
        122024.16,
        126477.125,
        114651.555,
        143701.72,
        196821.47,
        179382.5,
        83476.43,
        132382.72,
        166430.16,
    ][1:num_images]

    model = YOLO.v2_tiny.load(settings)
    YOLO.loadWeights!(model, settings)

    res = model(vocloaded.imstack_mat) #run once to deal with compillation overhead
    #accdata = minibatch(inp, out, settings.minibatch_size; xtype = xtype)

    t = @elapsed for i in 1:10
        model(vocloaded.imstack_mat)
    end

    inference_time = (t / 10) / num_images
    inference_rate = 1 / inference_time

    predictions = YOLO.postprocess(res, settings, conf_thresh = 0.3, iou_thresh = 0.3)
    @test length(predictions) == num_images # images
    @test length(predictions[1]) == 0
    @test length(predictions[2]) == 1

    @test isapprox(predictions[2][1].bbox.x, 0.23406872f0, rtol=0.05)
    @test isapprox(predictions[2][1].bbox.y, 0.28794688f0, rtol=0.05)
    @test isapprox(predictions[2][1].bbox.w, 0.8048826f0, rtol=0.05)
    @test isapprox(predictions[2][1].bbox.h, 0.46979588f0, rtol=0.05)
    @test predictions[2][1].class == 7
    @test isapprox(predictions[2][1].conf, 0.5042128f0, rtol=0.05)

    t = @elapsed for i in 1:10
        YOLO.postprocess(res, settings, conf_thresh = 0.3, iou_thresh = 0.3)
    end

    postprocess_time = (t / 10) / num_images
    postprocess_rate = 1 / postprocess_time

    ## Makie Tests
    #disabled because Makie can't be tested on headless CI
    # scene = YOLO.renderResult(vocloaded.imstack_mat[:,:,:,1], predictions[1], settings, save_file = "test.png")
    # @test isfile("test.png")
    # rm("test.png", force=true)

    enable_info()
    @info "YOLO_v2_tiny inference time per image: $(round(inference_time, digits=4)) seconds ($(round(inference_rate, digits=2)) fps)"
    @info "YOLO_v2_tiny postprocess time per image: $(round(postprocess_time, digits=4)) seconds ($(round(postprocess_rate, digits=2)) fps)"
    @info "Total time per image: $(round(inference_time + postprocess_time, digits=4)) seconds ($(round(1/(inference_time + postprocess_time), digits=2)) fps)"

end


text:
function getRateParams()

    ## initalize
    n = zeros(69, 1); # NFkB Activation Module
    i = zeros(34, 1); # IKK Activation Module
    a = zeros(32, 1);
    t = zeros(17, 1);  # TNF production module

    ## ----- IkB:NFkB Module [including A20 synthesis/deg] -----
    #           ikba       ikbb        ikbe
    n[1:3]  = [7e-5      1e-5       1e-6 ]'; # constitutive txn
    n[4:6]	= [8         0.02       0.3  ]'; # inducible txn
    n[7:9]  = [4.159060426         3          4.999962344    ]'; # Hill Coefficient
    n[10:12] = [0         37         37   ]'; # inducible txn delay
    n[13] = 0.005830947	; # mRNA Degradation [estimated in Lee et al. 2018 Processes]
    n[14:15] = [      3e-3       4e-3 ]'; # mRNA Degradation
    n[16:18] = [0.25      0.25       0.25 ]'; # translation rate

    n[19:21] = [0.09      9e-3       0.045]'; # Free IkB Import
    n[23:25] = [0.012     0.012      0.012]'; # Free IkB Export
    n[27:29] = [0.276     0.0276     0.138]'; # IkB:NFkB Import
    n[30:32] = [0.828     0.414      0.414]'; # IkB:NFkB Export
    n[22]   = 5.4;   # Free NFkB import()
    n[26]   = 0.0048;# Free NFkB export

    n[33:35] = [0.12      0.18       0.18 ]'; # IkB deg cytoplasm
    n[36:38] = [0.12      0.18       0.18 ]'; # IkB deg nucleus
    n[39:41] = [6e-5      6e-5       6e-5 ]'; # Bound IkB deg cyt
    n[42:44] = [6e-5      6e-5       6e-5 ]'; # Bound IkB deg nuc

    n[45:47] = [30        30         30   ]'; # IkB:NFkB Asn cyt
    n[48:50] = [30        30         30   ]'; # IkB:NFkB Asn nuc
    n[51:53] = [6e-5      6e-5       6e-5 ]'; # IkB:NFkB Dsn cyt
    n[54:56] = [6e-5      6e-5       6e-5 ]'; # IkB:NFkB Dsn nuc

    n[57:59] = [0.36      0.12       0.18 ]'; # Free IkB + IKK deg
    n[60:62] = [0.36      0.12       0.18 ]'; # IkB:NFkB + IKK deg

    # A20 [NEW]
    n[63]  = 2e-6;    # basal txn rate 2e-6
    n[64]  = 0.4;     # inducible txn for A[20],0.4
    n[65]  = 3;       # Hill Coefficient
    n[66]  = 0.035;   # mRNA Degradation
    n[67]  = 0.25;    # translation rate
    n[68]  = 0.0029;  # protein degradation
    n[69]  = 120;     # promoter shutdown [experiments show ~120min]

    ## ----- IKK Activation Module -----
    i[1] =     0.1698;  # k_b_lps 6
    i[2] =    0.17765;  # k_in_lps 7
    i[3] =    0.26114 ;  # k_out_lps 8
    i[4] =     13.367;  # k_d_lpsen, 9

    # ligand binding
    i[5]   = 0.19 ;   # k_b_tlr4lps, 10
    i[6]   = 2.7  ;   # k_ub_tlr4lps, 11
    i[7]   = 0.19 ;   # k_b_tlr4lpsen, 12
    i[8]   = 2.7  ;   # k_ub_tlr4lpsen, 13

    # generation & degradation
    i[9]   =   0.037479238 ;   # k_g_tlr4, 14 estimated in Lee et al. 2018
    i[10]  =    0.89603;   # k_d_tlr4,  15
    i[11]  =    9036.88927	;   # k_d_tlr4en, 16 estimated in Lee et al. 2018
    i[16]  =     14.414;     # k_deg_TLR4LPS,  21
    i[17]  =    0.42132;    # k_d_tlr4lpsen, 22

    # shutting
    i[12]  =     0.1344  ;   # k_in_tlr4, 17
    i[13]  =     3.6099  ;   # k_out_tlf4, 18
    i[14]  =   0.23513  ;     # k_in_tlr4lps, 19
    i[15]  =   0.041496  ;   # k_out_tlr4lps, 20

    # activation
    i[18]  =   3.2907 ;    # k_f_myd88 TL activated MyD88 -> (MyD88)_6, 23
    i[19]  = 3; # Hill coefficient for MyD88 activation, 24
    i[20]  =   0.057936;      #  km_a_myd88, 25
    i[21]  =    0.27504;     # k_dis_myd88, 0.02, 26
    i[22]  =    0.38509  ;     # k_a_trif  1e+8, 27
    i[23]  =   0.011822   ;     # k_i_trif, 0.01, 28

    # TRAF6
    i[24]  =     7.477 ; # k_a_TRAF6_MyS88s 66
    i[25]  =     3.4132; # k_a_TRAF6_TRIFs 67
    i[26]  =    0.21778;       # k_i_TRAF6  68

    # IKKK
    i[27]  =    0.34317;   # IKKK_off --> IKKK [TRAF6s mediated], 71
    i[28]  =    5e-7;   # IKKK_off --> IKKK [constitutive],72
    i[29]  = 0.25;    # IKKK     --> IKKK_off [constitutive], 73

    # IKK
    i[30]  = 4753.870211;    # IKK_off  --> IKK [IKKK mediated],74 [estimated in Lee et al. 2018]
    i[31]  = 5e-5 ;    # IKK_off  --> IKK [constitutive],75
    i[32]  = 0.02  ;    # IKK      --> IKK_off [constitutive], 76
    i[33]  = 0.028466263 ;    # IKK      --> IKK_i [constitutive]  , 77 [estimated in Lee et al. 2018]
    i[34]  = 0.02  ;    # IKKi     --> IKK_off [constitutive] , 78

    ## the TNFR parameters
    # ----- IKK Activation Module ----- NEW f
    a[1]   = 0.0154; # pd_m_tnf 45' half life of exogenous TNF, 49

    # tnfrm metabolism [synthesis & degradation]
    a[2]   = 2e-7;   # tnfrm synthesis [txn, tsl, localization], 36
    a[3]   = 0.0058; # tnfrm --> deg  -- 120' halflife, 37

    # TNF-Independent C1 Activation
    a[4]   = 1e-5;   # 3tnfrm --> TNFR, 38
    a[5]   = 0.1;    # TNFR   --> 3tnfrm, 39
    a[6]   = 0.023;  # TNFR internalization -- 30' halflife, 40

    a[7]   = 100;    # TNFR + TTR --> C1_off, 41
    a[8]   = 0.1;    # C1_off --> TNFR + TTR, 42
    a[9]   = 30;     # C1_off --> C1, 43
    a[10]  = 2;      # C1     --> C1_off, 44
    a[11]  = 1000;   # C1     --> C1_off [A20 Mediated], 45
    a[12]  = 0.1;   # C1     --> TNFR + TTR, 46
    a[13]  = 0.023;   # C1_off internalization, 47
    a[14]  = 0.023;   # C1 internalization, 48

    # TNF-dependent C1 Activation
    a[15]  = 1100;   # 3tnfrm + tnf --> TNFRtnf, 50
    a[16]  = 1100;  # TNFR + tnf --> TNFRtnf, 51
    a[17]  = 0.021;  # TNFRtnf   --> TNFR + tnf, 52
    a[18]  = 0.023;   # TNFRtnf internalization -- 30' halflife, 53

    a[19]  = 100;   # TNFRtnf + TTR --> C1tnf_off, 54
    a[20]  = 0.1;   # C1tnf_off --> TNFRtnf + TTR, 55
    a[21]  = 30;   # C1tnf_off --> C1tnf, 56
    a[22]  = 2;  # C1tnf     --> C1tnf_off, 57
    a[23]  = 1000;  # C1tnf     --> C1tnf_off [A20 Mediated], 58
    a[24]  = 0.1;   # C1tnf     --> TNFRtnf + TTR, 59
    a[25]  = 0.023;   # C1tnf_off internalization, 60
    a[26]  = 0.023;   # C1tnf internalization, 61

    a[27]  = 0.021;  # C1tnf_off --> C1_off + tnf, 62
    a[28]  = 1100;  # C1_off + tnf --> C1tnf_off, 63
    a[29]  = 0.021;  # C1tnf    --> C1 + tnf, 64
    a[30]  = 1100;  # C1 + tnf --> C1tnf, 65

    # IKKK

    a[31]  = 500;    # IKKK_off --> IKKK [C1 mediated]?500, 69
    a[32]  = 500;  # IKKK_off --> IKKK [C1tnf mediated], 70

    ## TNF part
    t[1] = 0;# 1e-5;#tnf constitutive txn
    t[2] = 1e-5;# tnf induced txn [rescaled in the whole model], 150
    t[3] = 2;# tnf transcription,Hill coefficient, 151
    t[4] = .65;# tnf transcription induction EC50, 152
    t[5] = .02;# tnf transcript deg. max rate, 154
    t[6] = .4;# tnf nascent process rate, 153
    t[7] = .05;# tnf protein syns rate, 156
    t[8] = .07;# tnf deg rate, 157
    t[9] = .07;# tnf secretion rate, 158

    # from private communication with Dr. Zhang Cheng in Prof. Hoffmann's group
    t[10] = 0.0001;  # K_a0
    t[11] = 0.1;     # K_a
    t[12] = 0.4;     # K_i

    t[13] = 833.3; # estimated
    t[14] = 90;  # assumed time needed for full effects of BFA
    t[15] = 0.01; # assumed rate of deactivation of TRAF6* due to A20

    ##
    t[16] = 0.998934823;  # Effect from BFA on translation coefficient
    t[17] = 0.2;  # denominator term in translation inhibition

    return n, i, a, t
end

NP, IP, AP, TP = getRateParams()

mutable struct RateParams
    NP
    IP
    AP
    TP
    IKK_TOTAL
    flag_noTnfFeedback
end

IKK_TOTAL = sum(INIT[33:35])
flag_noTnfFeedback = false

v = RateParams(NP, IP, AP, TP, IKK_TOTAL, flag_noTnfFeedback)

max_stars_repo_name: jiweiqi/NFkB

max_issues_repo_name: jiweiqi/NFkB

max_forks_repo_name: jiweiqi/NFkB


text:
const IDType = UInt32
const NOT_A_VAR = typemax(IDType)
const EMPTY_IDS = IDType[]

struct Monomial <: AbstractMonomial
    ids::Vector{IDType}
end
Monomial() = Monomial(EMPTY_IDS)
degree(x::Monomial) = length(x.ids)
Base.copy(x::Monomial) = Monomial(copy(x.ids))

firstid(m::Monomial) = m.ids[1]
degree(m::Monomial, id) = count(isequal(id), m.ids)
rmid(m::Monomial, id) = Monomial(filter(!isequal(id), m.ids))

const VARNAME_DICT = Dict{IDType, String}(0 => "x", 1 => "y", 2 => "z")
const SUPERSCRIPTS = ['⁰', '¹', '²', '³', '⁴', '⁵', '⁶', '⁷', '⁸', '⁹']
const LOWERSCRIPTS = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉']

function int2superscript(x)
    mapreduce(d->SUPERSCRIPTS[d + 1], *, Iterators.reverse(digits(x)))
end
function int2lowerscript(x)
    mapreduce(d->LOWERSCRIPTS[d + 1], *, Iterators.reverse(digits(x)))
end
function print_single_monomial(io, v, o, star=false)
    iszero(o) && return print(io, 1)
    star && (PRETTY_PRINT[] || print(io, '*'))
    print(io, VARNAME_DICT[v])
    if o > 1
        if PRETTY_PRINT[]
            print(io, int2superscript(o))
        else
            print(io, '^', o)
        end
    elseif o == 1 || o == 0
    else
        error("unreachable")
    end
end
function Base.show(io::IO, x::Monomial)
    isempty(x.ids) && (print(io, '1'); return)

    ids = x.ids
    v = first(ids)
    count = 0
    star = false
    for id in ids
        if id == v
            count += 1
        else
            print_single_monomial(io, v, count, star)
            star = true
            v = id
            count = 1
        end
    end
    if count > 0
        print_single_monomial(io, v, count, star)
    end
    return nothing
end

function Base.:*(y::Monomial, x::Monomial)
    ids = y.ids
    i = j = 1
    n0 = length(ids)
    n1 = length(x.ids)
    r = Monomial(similar(ids, n0 + n1))
    T = eltype(ids)
    M = typemax(T)
    for k in 1:(n0 + n1)
        a = i <= n0 ? ids[i] : M
        b = j <= n1 ? x.ids[j] : M
        if a < b
            i += 1
            r.ids[k] = a
        else
            j += 1
            r.ids[k] = b
        end
    end
    return r
end

function Base.:(/)(y::Monomial, x::Monomial)
    n = Monomial(similar(y.ids, 0))
    i = j = 1
    ids = y.ids
    n0 = length(ids)
    n1 = length(x.ids)
    T = eltype(ids)
    M = typemax(T)
    while (i + j - 2) < (n0 + n1)
        a = i <= n0 ? ids[i] : M
        b = j <= n1 ? x.ids[j] : M
        if a < b
            push!(n.ids, a)
            i += 1
        elseif a == b
            i += 1
            j += 1
        else
            return n, true
        end
    end
    return n, false
end

Base.:(==)(x::Monomial, y::Monomial) = (x === y) || (x.ids == y.ids)

# graded lex
function Base.isless(x::Monomial, y::Monomial)
    dx = degree(x)
    dy = degree(y)
    dx < dy && return true
    dx > dy && return false
    for i in 1:dx
        vx = x.ids[i]
        vy = y.ids[i]
        vx != vy && return vx > vy
    end
    return false
end

function Base.gcd(x::Monomial, y::Monomial)
    g = Monomial(similar(x.ids, 0))
    i = j = 1
    n0 = length(x.ids)
    n1 = length(y.ids)
    n = n0 + n1 + 2
    T = eltype(x.ids)
    M = typemax(T)
    while i + j < n
        xk = i <= n0 ? x.ids[i] : M
        yk = j <= n1 ? y.ids[j] : M
        if xk < yk
            i += 1
        elseif xk > yk
            j += 1
        else
            push!(g.ids, xk)
            i += 1
            j += 1
        end
    end
    return g
end