#!/bin/bash

# dz_sarscov2_phi_i - Algerian SARS-CoV-2 pandemic data with phi_i from ArXiv:2007.11779v1

# (C) Wikimedia Commons User:Boud, 2020 CC-BY-SA-4.0
# https://creativecommons.org/licenses/by-sa/4.0/
# https://creativecommons.org/2015/10/08/cc-by-sa-4-0-now-one-way-compatible-with-gplv3/

# Refs: see
# * https://arxiv.org/abs/arXiv:2007.11779
# * https://zenodo.org/record/3951152

## Download daily data

## files to download
TMPdays=/tmp/tmp.alldata.tmponly
TMP28=/tmp/tmp.phi.tmponly

## plotting files
TMP0=/tmp/tmp.dz.0
TMPplus=/tmp/tmp.dz.plus
TMPminus=/tmp/tmp.dz.minus
TMPlowest28_1=/tmp/tmp.dz.lowest1
TMPlowest28_2=/tmp/tmp.dz.lowest2

LANG=C
DEBUG=NO
#DEBUG=YES

if [ "x$1" = "xw" ]; then
    wget -O ${TMPdays} https://zenodo.org/record/3951152/files/WP_C19CCTF_SARSCoV2.dat
    wget -O ${TMP28} https://zenodo.org/record/3951152/files/phi_N_28days.dat
fi

## Extract DZ data
day_one=$(grep DZ, ${TMPdays} | awk -F , '{print $1}'|head -n1)
N_days=$(grep DZ, ${TMPdays} | wc |awk '{print $1}')

dz_daily=$(grep DZ, ${TMPdays} |awk -F , '{cnew=$3-cold; cold=$3; print cnew}')
if [ "x${DEBUG}" = "xYES" ]; then
    printf "${dz_daily}" | tr '\n' ' ' ; printf "\n"
fi

## Extract DZ least noisy period
dz_lowest_28_day1=$(grep ^DZ ${TMP28} | awk '{print $NF}')
dz_lowest_28_date1=$(grep DZ, ${TMPdays} | grep -n "${dz_lowest_28_day1}" | \
                            sed -e 's/:/,/' | awk -F , '{print $1}')
dz_lowest_28_date2=$(echo ${dz_lowest_28_date1} |awk '{print $1+27}')
if [ "x${DEBUG}" = "xYES" ]; then
    printf "dz_lowest_28_date1=${dz_lowest_28_date1}\n"
    printf "dz_lowest_28_date2=${dz_lowest_28_date2}\n"
fi

## Create median model (see Methods section of ArXiv:2007.11779)
dz_median_0=$(echo ${dz_daily[@]} | tr ' ' '\n' | \
                    LANG=C awk '{i+=1; c[5]=c[4]; c[4]=c[3]; c[3]=c[2]; c[2]=c[1]; c[1]=$1; \
                         d[4]=c[5]; d[3]=c[4]; d[2]=c[2]; d[1]=c[1]; \
                         asort(d,e); median = (e[2]+e[3])*0.5; \
                         if(i>=5)print median}')

first=$(echo ${dz_median_0}|awk '{print $1}')
last=$(echo ${dz_median_0}|awk '{print $NF}')
dz_median="${first} ${first} ${dz_median_0} ${last} ${last}"

if [ "x${DEBUG}" = "xYES" ]; then
    printf "${dz_median}\n"
fi

## print days and Poisson bounds
echo ${dz_median} | tr ' ' '\n' | awk '{print NR,$1-sqrt($1)}' > ${TMPminus}

echo ${dz_median} | tr ' ' '\n' | awk '{print NR,$1+sqrt($1)}' > ${TMPplus}

echo ${dz_daily} | tr ' ' '\n' | awk '{print NR,$1}' > ${TMP0}

FORMAT=svg
#FORMAT=ps
X1=0
X2=${N_days}
Y1=0
Y2=700
X_label="date starting ${day_one}"
Y_label="number of infections"
title="least noisy 28 days: starting ${dz_lowest_28_day1}  "

X_label="date commencant le ${day_one}"
Y_label="nombre d'infections"
title="28 j. moindre bruit: debut=${dz_lowest_28_day1}"

printf "${dz_lowest_28_date1} ${Y1}\n${dz_lowest_28_date1} ${Y2}\n" > ${TMPlowest28_1}
printf "${dz_lowest_28_date2} ${Y1}\n${dz_lowest_28_date2} ${Y2}\n" > ${TMPlowest28_2}

graph -T${FORMAT} -w 0.75 -W0.001 -X "${X_label}" -Y "${Y_label}" \
      -x ${X1} ${X2} -y ${Y1} ${Y2} ${TMPminus} \
      --reposition 0.0 0.0 1.0 --blankout 0.0 \
      -x ${X1} ${X2} -y ${Y1} ${Y2} ${TMPplus} \
      --reposition 0.0 0.0 1.0 --blankout 0.0 \
      -C --pen-colors "1=blue:2=green:3=red" -m0 -S 16 \
      -x ${X1} ${X2} -y ${Y1} ${Y2} ${TMP0} \
      --reposition 0.0 0.0 1.0 --blankout 0.0 \
      -m3 -W0.001 \
      --frame-color red -L "${title}" --grid-style 1 \
      -x ${X1} ${X2} -y ${Y1} ${Y2} ${TMPlowest28_2} \
      --reposition 0.0 0.0 1.0 --blankout 0.0 \
      -m3 -W0.001 \
      --frame-color black -L "" \
      -x ${X1} ${X2} -y ${Y1} ${Y2} ${TMPlowest28_1} \
      > SARS-CoV-2-DZ.${FORMAT}