<jupyter_start><jupyter_text>[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/catboost/tutorials/blob/master/example_usages/COVID_19.ipynb)This is an example on how to predict COVID-19 spread. The goal is to predict Confirmed cases and Fatalities for future dates based on data for previous dates.  Table of contents* [General description](general_description)* [Install and import necessary packages ](install_and_import)* [Get main data from Kaggle](get_main_data)* [Feature engineering](features)  * [Time delay embedding features](time_delay_embedding_features)  * [Day feature](day_feature)  * [WeekDay feature](week_day_feature)  * [Days since Xth Confirmed case and Xth Fatality features](days_since_features)  * ['Distance to Origin' feature](distance_to_origin_feature)  * [Functions for merging external data](functions_for_merging)  * [Country Area](country_area_feature)  * [Country population features](country_population_features)  * [Country Population density feature](country_population_density_feature)  * [Country Smoking rate feature](country_smoking_rate_feature)  * [Country hospital beds per 1000 people](country_hospital_beds_feature)  * [Country Health Expenditure](country_health_expenditure_feature)* [Prepare data for training](prepare_for_training)* [Models training](models_training)* [Feature importance in the models](feature_importance)* [Create predictions for eval and test data](predictions)* [Plots with predictions](plots_with_predictions)  General description The main data to train on is taken from https://www.kaggle.com/c/covid19-global-forecasting-week-1/data.Regions are specified using two data fields: ``Province/State`` and ``Country/Region``.This data is also enriched with region coordinates in ``Lat`` and ``Long`` fields.Train data contains cumulative numbers of Confirmed cases and Fatalities for the different regions for the number of consecutive dates (contained in ``Date`` field) from 2020-01-22 to 2020-03-24. Test data is for dates from 2020-03-12 to 2020-04-23. Test data does not have Confirmed cases and Fatalities.We will also add additional data for features from other sources:* Country Land area from [World bank data](https://data.worldbank.org/indicator/AG.LND.TOTL.K2)* Country Smoking rate from [World bank data](https://data.worldbank.org/indicator/SH.PRV.SMOK)* Country Population Age and Sex distribution from [UN World Population Prospects](https://population.un.org/wpp/Download/Standard/CSV/)* Country Hospital beds per 1000 people from [World Bank data](https://data.worldbank.org/indicator/sh.med.beds.zs)* Country Health expenditure per Capita, PPP (current international $) from [World Bank data](https://data.worldbank.org/indicator/SH.XPD.CHEX.PP.CD) We will construct models that predict new cases and new fatalities (in the logarithmic scale) for the next day based on time delay embedding features for the last 30 days plus additional features from sources described above. For forecasting further into the future (more than a single day) we will use an incremental (by days) approach: we will use predicted values of new cases and new fatalities for days with already calculated predictions as input data for time delay embedding features for the next prediction day and then use the same procedure for the new next prediction day and so on.  Install and import necessary packages<jupyter_code>!pip install kaggle -U
!pip install pandas -U
!pip install catboost -U

import os
import pathlib
import re

import numpy as np
import pandas as pd
import sklearn.preprocessing

import geopy.distance

import catboost as cb<jupyter_output>Requirement already up-to-date: kaggle in /usr/local/lib/python3.6/dist-packages (1.5.6)
Requirement already satisfied, skipping upgrade: certifi in /usr/local/lib/python3.6/dist-packages (from kaggle) (2020.4.5.1)
Requirement already satisfied, skipping upgrade: six>=1.10 in /usr/local/lib/python3.6/dist-packages (from kaggle) (1.12.0)
Requirement already satisfied, skipping upgrade: urllib3<1.25,>=1.21.1 in /usr/local/lib/python3.6/dist-packages (from kaggle) (1.24.3)
Requirement already satisfied, skipping upgrade: python-dateutil in /usr/local/lib/python3.6/dist-packages (from kaggle) (2.8.1)
Requirement already satisfied, skipping upgrade: requests in /usr/local/lib/python3.6/dist-packages (from kaggle) (2.21.0)
Requirement already satisfied, skipping upgrade: tqdm in /usr/local/lib/python3.6/dist-packages (from kaggle) (4.38.0)
Requirement already satisfied, skipping upgrade: python-slugify in /usr/local/lib/python3.6/dist-packages (from kaggle) (4.0.0)
Requirement already satisf[...]<jupyter_text>Get main data from Kaggle As we download from Kaggle we need to provide Kaggle authetification cridentials. We will use the method suggested [here](https://gist.github.com/jayspeidell/d10b84b8d3da52df723beacc5b15cb27gistcomment-2814834). Replace KAGGLE_USERNAME and KAGGLE_KEY with your own.<jupyter_code>os.environ['KAGGLE_USERNAME'] = "XXX"
os.environ['KAGGLE_KEY'] = "YYY"<jupyter_output><empty_output><jupyter_text>Let's download the data now.<jupyter_code>kaggle_data_name = 'covid19-global-forecasting-week-1'
data_folder = os.path.join(pathlib.Path.home(), 'datasets', kaggle_data_name)

!mkdir -p ${data_folder}
!kaggle competitions download -c {kaggle_data_name} --path {data_folder}<jupyter_output>Warning: Looks like you're using an outdated API Version, please consider updating (server 1.5.6 / client 1.5.4)
Downloading submission.csv to /root/datasets/covid19-global-forecasting-week-1

  0% 0.00/108k [00:00<?, ?B/s]
100% 108k/108k [00:00<00:00, 23.0MB/s]
Downloading train.csv to /root/datasets/covid19-global-forecasting-week-1
  0% 0.00/947k [00:00<?, ?B/s]
100% 947k/947k [00:00<00:00, 62.2MB/s]
Downloading test.csv to /root/datasets/covid19-global-forecasting-week-1
  0% 0.00/540k [00:00<?, ?B/s]
100% 540k/540k [00:00<00:00, 165MB/s]<jupyter_text>Let's check that data is ok.<jupyter_code>!ls -l {data_folder}<jupyter_output>total 1600
-rw-r--r-- 1 root root 111051 Apr 19 14:10 submission.csv
-rw-r--r-- 1 root root 552851 Apr 19 14:10 test.csv
-rw-r--r-- 1 root root 969831 Apr 19 14:10 train.csv<jupyter_text>Let's load data to pandas.DataFrame.<jupyter_code>original_train_df = pd.read_csv(os.path.join(data_folder,'train.csv'), parse_dates=['Date'])
original_train_df.head()
original_test_df = pd.read_csv(os.path.join(data_folder,'test.csv'), parse_dates=['Date'])
original_test_df.head()<jupyter_output><empty_output><jupyter_text>We will concatenate original train and test dataframes for easier time delay embedding features calculations.Original train and test data intersect by dates. We will use data from original train in the concatenated dataset for such dates.<jupyter_code>last_original_train_date = original_train_df['Date'].max()
print ('last_original_train_date = ', last_original_train_date)

print ('original_test_df.shape =', original_test_df.shape)

original_test_wo_train_df = original_test_df.drop(
    index=original_test_df[original_test_df['Date'] <= last_original_train_date].index
)

print ('original_test_wo_train_df.shape =', original_test_wo_train_df.shape)

# recreate index bacause we will need unique index values later 
main_df = pd.concat([original_train_df, original_test_wo_train_df], ignore_index=True)<jupyter_output>last_original_train_date =  2020-03-24 00:00:00
original_test_df.shape = (12212, 6)
original_test_wo_train_df.shape = (8520, 6)<jupyter_text>There are special `Provinces` that are in fact data from Cruise ships with `Country` information. Make this ships  ``Countries`` and passengers origin country a ``Province/State`` (for proper hierarchy).<jupyter_code>from_cruise_ships = main_df['Province/State'].isin(['From Diamond Princess', 'Grand Princess'])
main_df.loc[from_cruise_ships, ['Province/State','Country/Region']] = main_df.loc[from_cruise_ships, ['Country/Region','Province/State']].values
main_df[main_df['Country/Region'].isin(['From Diamond Princess', 'Grand Princess'])]<jupyter_output><empty_output><jupyter_text>Feature engineering  Time delay embedding features ConfirmedCases and Fatalities are transformed to the logarithmic scale ( log(1+x) to be precise ). Add log(1+x) - transformed number of new cases and new fatalities. Logarithmic scale is natural for this task as  the growth of this values is typically exponential.Then add their historic values as (LogNewCases_prev_day_X, LogNewFataities_prev_day_X, X = 1 to 30) features.Some locations contain broken data (Confirmed Cases or Fatalities contain decreasing values), remove them from data.<jupyter_code>main_df.sort_values(by='Date', inplace=True)

location_columns = ['Country/Region','Province/State']

# row selection functions in pandas do not like nans
for column in location_columns:
    main_df[column].fillna('', inplace=True)


days_history_size = 30


def is_cumulative(increment_series):
    for v in increment_series:
        if (not np.isnan(v)) and (v < 0):
            return False
    return True

print ('data size before removing bad data = ', len(main_df))

for field in ['LogNewConfirmedCases', 'LogNewFatalities']:
    main_df[field] = np.nan

    for prev_day in range(1, days_history_size + 1):
        main_df[field + '_prev_day_%s' % prev_day] = np.nan


for location_name, location_df in main_df.groupby(['Country/Region', 'Province/State']):
    for field in ['ConfirmedCases', 'Fatalities']:
        new_values = location_df[field].values
        new_values[1:] -= new_values[:-1]
        if not is_cumulative(new_values):
            print ('%s for %s, %s is not valid cumulative series, drop it' % ((field,) + location_name))
            main_df.drop(index=location_df.index, inplace=True)
            break
        log_new_values = np.log1p(new_values)
        main_df.loc[location_df.index, 'LogNew' + field] = log_new_values

        for prev_day in range(1, days_history_size + 1):
            main_df.loc[location_df.index[prev_day:], 'LogNew%s_prev_day_%s' % (field, prev_day)] = (
                log_new_values[:-prev_day]
            )

print ('data size after removing bad data = ', len(main_df))
main_df.head()<jupyter_output>data size before removing bad data =  26412
ConfirmedCases for Australia, Northern Territory is not valid cumulative series, drop it
ConfirmedCases for Australia, Queensland is not valid cumulative series, drop it
ConfirmedCases for Azerbaijan,  is not valid cumulative series, drop it
ConfirmedCases for Bahrain,  is not valid cumulative series, drop it
Fatalities for Canada, Quebec is not valid cumulative series, drop it
ConfirmedCases for China, Guizhou is not valid cumulative series, drop it
ConfirmedCases for France, France is not valid cumulative series, drop it
ConfirmedCases for France, Saint Barthelemy is not valid cumulative series, drop it
ConfirmedCases for Grand Princess, US is not valid cumulative series, drop it
ConfirmedCases for Guyana,  is not valid cumulative series, drop it
Fatalities for Iceland,  is not valid cumulative series, drop it
Fatalities for India,  is not valid cumulative series, drop it
ConfirmedCases for Japan,  is not valid cumulative series, drop it
Fa[...]<jupyter_text>Day feature Now, add simple linear Day feature<jupyter_code>first_date = min(main_df['Date'])
    
main_df['Day'] = (main_df['Date'] - first_date).dt.days.astype('int32')<jupyter_output><empty_output><jupyter_text>WeekDay feature Let's add WeekDay feature<jupyter_code>main_df['WeekDay'] = main_df['Date'].transform(lambda d: d.weekday())<jupyter_output><empty_output><jupyter_text>Days since Xth Confirmed case and Xth Fatality features Add Days since Xth Confirmed case and Xth Fatality (where X = 1, 10, 100) features<jupyter_code>thresholds = [1, 10, 100]

for threshold in thresholds:
    main_df['Days_since_ConfirmedCases=%s' % threshold] = np.nan
    main_df['Days_since_Fatalities=%s' % threshold] = np.nan

for location_name, location_df in main_df.groupby(['Country/Region', 'Province/State']):
    for field in ['ConfirmedCases', 'Fatalities']:
        for threshold in thresholds:
            first_day = location_df['Day'].loc[location_df[field] >= threshold].min()
            if not np.isnan(first_day):
                main_df.loc[location_df.index, 'Days_since_%s=%s' % (field, threshold)] = (
                    location_df['Day'].transform(lambda day: -1 if (day < first_day) else (day - first_day))
                )<jupyter_output><empty_output><jupyter_text>'Distance to Origin' feature. Chinese province of Hubei is the original epicentre of the pandemic, so add a feature that represents distance from the region of the sample to the epidemic's origin.<jupyter_code>def get_hubei_coords(df):
    for index, row in df.iterrows():
        if row['Province/State'] == 'Hubei':
            return (row['Lat'], row['Long'])

    raise Exception('Hubei not found in data')


origin_coords = get_hubei_coords(main_df)

main_df['Distance_to_origin'] = main_df.apply(
    lambda row: geopy.distance.distance((row['Lat'], row['Long']), origin_coords).km,
    axis='columns'
)<jupyter_output><empty_output><jupyter_text>Functions for merging external data Let's add some data from other sources.But first, we will need a couple of helper functions with common functionality.World bank data, UN World Population Project data and the main dataset use different names for some of the countries so we will need to remap country names in external datasets for correct merging.<jupyter_code>def merge_with_column_drop(left_df, right_df, right_df_column='Country'):
    df = pd.merge(
        left=left_df,
        right=right_df,
        how='left',
        left_on='Country/Region',
        right_on=right_df_column
    )
    df.drop(columns=right_df_column, inplace=True)
    
    return df

def remap_country_name_from_world_bank_to_main_df_name(country):
    return {
        'Bahamas, The': 'The Bahamas',
        'Brunei Darussalam': 'Brunei',
        'Congo, Rep.': 'Congo (Brazzaville)',
        'Congo, Dem. Rep.': 'Congo (Kinshasa)',
        'Czech Republic': 'Czechia',
        'Egypt, Arab Rep.': 'Egypt',
        'Iran, Islamic Rep.': 'Iran',
        'Korea, Rep.': 'Korea, South',
        'Kyrgyz Republic': 'Kyrgyzstan',
        'Russian Federation': 'Russia',
        'Slovak Republic': 'Slovakia',
        'St. Lucia': 'Saint Lucia',
        'St. Vincent and the Grenadines': 'Saint Vincent and the Grenadines',
        'United States': 'US',
        'Venezuela, RB': 'Venezuela',
    }.get(country, country)

def remap_country_name_from_un_wpp_to_main_df_name(country):
    return {
        'Bahamas': 'The Bahamas',
        'Bolivia (Plurinational State of)': 'Bolivia',
        'Brunei Darussalam': 'Brunei',
        'China, Taiwan Province of China': 'Taiwan*',
        'Congo' : 'Congo (Brazzaville)',
        'Côte d\'Ivoire': 'Cote d\'Ivoire',
        'Democratic Republic of the Congo': 'Congo (Kinshasa)',
        'Gambia': 'The Gambia',
        'Iran (Islamic Republic of)': 'Iran',
        'Republic of Korea': 'Korea, South',
        'Republic of Moldova': 'Moldova',
        'Réunion': 'Reunion',
        'Russian Federation': 'Russia',
        'United Republic of Tanzania': 'Tanzania',
        'United States of America': 'US',
        'Venezuela (Bolivarian Republic of)': 'Venezuela',
        'Viet Nam': 'Vietnam'
    }.get(country, country)

# for read_csv
world_bank_converters={'Country Name': remap_country_name_from_world_bank_to_main_df_name}
un_wpp_converters={'Location': remap_country_name_from_un_wpp_to_main_df_name}<jupyter_output><empty_output><jupyter_text>Country Area Let's add Country Area feature. We will use data from World Bank for that.<jupyter_code>!mkdir -p area
os.chdir('area')
!wget -nc http://api.worldbank.org/v2/en/indicator/AG.LND.TOTL.K2?downloadformat=csv -O area.zip
!unzip -n area.zip  
!ls -l .
!head API_AG.LND.TOTL.K2_DS2_en_csv_v2_937279.csv

area_df = pd.read_csv('./API_AG.LND.TOTL.K2_DS2_en_csv_v2_937279.csv', skiprows=4, converters=world_bank_converters)
os.chdir('..')
area_df.head()<jupyter_output><empty_output><jupyter_text>We will use the last available data.<jupyter_code>year_columns = [str(year) for year in range(1960, 2020)]

area_df['CountryArea'] = area_df[year_columns].apply(
    lambda row: row[row.last_valid_index()] if row.last_valid_index() else np.nan,
    axis='columns'    
)<jupyter_output><empty_output><jupyter_text>Drop unneeded columns<jupyter_code>area_df = area_df[['Country Name', 'CountryArea']]
area_df.head()<jupyter_output><empty_output><jupyter_text>Now, merge this data with main data by country.<jupyter_code>main_df = merge_with_column_drop(
    main_df,
    area_df,
    right_df_column='Country Name'
)
main_df.head()<jupyter_output><empty_output><jupyter_text>Country population features Now, let's add country population data from UN.<jupyter_code>!mkdir -p population
os.chdir('population')
!wget -nc https://population.un.org/wpp/Download/Files/1_Indicators%20\(Standard\)/CSV_FILES/WPP2019_PopulationByAgeSex_Medium.csv

!head ./WPP2019_PopulationByAgeSex_Medium.csv
population_df = pd.read_csv(
    'WPP2019_PopulationByAgeSex_Medium.csv',
     usecols=['Location', 'Time', 'AgeGrp', 'PopMale', 'PopFemale', 'PopTotal'],
    parse_dates=['Time'],
    converters=un_wpp_converters
)
os.chdir('..')
population_df.head()<jupyter_output><empty_output><jupyter_text>We will use only recent (no earlier than 2014) data<jupyter_code>population_df = population_df.loc[
    (population_df['Time'] >= pd.Timestamp(2014,1,1))
    & (population_df['Time'] <= pd.Timestamp(2019,1,1))
]
population_df.head()<jupyter_output><empty_output><jupyter_text>Aggregate data by 20-year groups and also compute cumulative population.<jupyter_code>aggegated_population_df = pd.DataFrame()

for (location, time), group_df in population_df.groupby(['Location', 'Time']):
    # by ['Pop_0-20', 'Pop_20-40', 'Pop_40-60', 'Pop_60-80', 'Pop_80+']
    pop_by_age_groups = [0] * 5 

    pop_male = 0
    pop_female = 0

    for _, row in group_df.iterrows():
        age_grp_start = int( re.split(r'[\-\+]', row['AgeGrp'])[0] )
        pop_by_age_groups[min(age_grp_start // 20, 4)] += row['PopMale'] + row['PopFemale']
        pop_male += row['PopMale']
        pop_female += row['PopFemale']

    ""
    aggegated_population_df = aggegated_population_df.append(
        {
            'Location': location,
            'Time': time,
            'CountryPop_0-20': pop_by_age_groups[0],
            'CountryPop_20-40': pop_by_age_groups[1],
            'CountryPop_40-60': pop_by_age_groups[2],
            'CountryPop_60-80': pop_by_age_groups[3],
            'CountryPop_80+': pop_by_age_groups[4],
            'CountryPopMale': pop_male,
            'CountryPopFemale': pop_female,
            'CountryPopTotal': pop_male + pop_female
        },
        ignore_index=True        
    )

aggegated_population_df.head()<jupyter_output><empty_output><jupyter_text>Select the most recent data.<jupyter_code>aggegated_population_df = aggegated_population_df.sort_values('Time').drop_duplicates(['Location'], keep='last')
aggegated_population_df.head()<jupyter_output><empty_output><jupyter_text>Now, drop Time and join with main data<jupyter_code>aggegated_population_df.drop(columns='Time', inplace=True)

main_df = merge_with_column_drop(
    main_df,
    aggegated_population_df,
    right_df_column='Location'
)
main_df.head()<jupyter_output><empty_output><jupyter_text>Country Population density feature With both Area and Population data we can now compute Population density feature.<jupyter_code>main_df['CountryPopDensity'] = main_df['CountryPopTotal'] / main_df['CountryArea']<jupyter_output><empty_output><jupyter_text>Country Smoking rate feature Let's also add Smoking rate by Country data because smoking can influence the severity of respiratory diseases.<jupyter_code>!mkdir -p smoking
os.chdir('smoking')
!wget -nc http://api.worldbank.org/v2/en/indicator/SH.PRV.SMOK?downloadformat=csv -O smoking.zip
!unzip -n smoking.zip  
!ls -l .
!head API_SH.PRV.SMOK_DS2_en_csv_v2_937317.csv
smoking_df = pd.read_csv('./API_SH.PRV.SMOK_DS2_en_csv_v2_937317.csv', skiprows=4)
os.chdir('..')
smoking_df.head()<jupyter_output><empty_output><jupyter_text>We will use the last available data from recent years (2010 and later).<jupyter_code>recent_year_columns = [str(year) for year in range(2010, 2020)]

smoking_df['CountrySmokingRate'] = smoking_df[recent_year_columns].apply(
    lambda row: row[row.last_valid_index()] if row.last_valid_index() else np.nan,
    axis='columns'    
)<jupyter_output><empty_output><jupyter_text>Drop unneeded columns<jupyter_code>smoking_df = smoking_df[['Country Name', 'CountrySmokingRate']]
smoking_df.head()<jupyter_output><empty_output><jupyter_text>Now, merge this data with main data by country.<jupyter_code>main_df = merge_with_column_drop(
    main_df,
    smoking_df,
    right_df_column='Country Name'
)
main_df.head()<jupyter_output><empty_output><jupyter_text>Country hospital beds per 1000 people Let's add some information about countries' health system.We will add the number of hospital beds per 1000 people first.<jupyter_code>!mkdir -p hospital_beds
os.chdir('hospital_beds')
!wget -nc http://api.worldbank.org/v2/en/indicator/SH.MED.BEDS.ZS?downloadformat=csv -O hospital_beds.zip
!unzip -n hospital_beds.zip  
!ls -l .
!head API_SH.MED.BEDS.ZS_DS2_en_csv_v2_935968.csv
hospital_beds_df = pd.read_csv('./API_SH.MED.BEDS.ZS_DS2_en_csv_v2_935968.csv', skiprows=4)
os.chdir('..')
hospital_beds_df.head()<jupyter_output><empty_output><jupyter_text>We will use the last available data from recent years (2010 and later).<jupyter_code>recent_year_columns = [str(year) for year in range(2010, 2020)]

hospital_beds_df['CountryHospitalBedsRate'] = hospital_beds_df[recent_year_columns].apply(
    lambda row: row[row.last_valid_index()] if row.last_valid_index() else np.nan,
    axis='columns'    
)<jupyter_output><empty_output><jupyter_text>Drop unneeded columns<jupyter_code>hospital_beds_df = hospital_beds_df[['Country Name', 'CountryHospitalBedsRate']]
hospital_beds_df.head()<jupyter_output><empty_output><jupyter_text>Now, merge this data with main data by country.<jupyter_code>main_df = merge_with_column_drop(
    main_df,
    hospital_beds_df,
    right_df_column='Country Name'
)
main_df.head()<jupyter_output><empty_output><jupyter_text>Country Health Expenditure Now, let's add Health expenditure per capita, PPP (in current international $) feature.<jupyter_code>!mkdir -p health_expenditure
os.chdir('health_expenditure')
!wget -nc http://api.worldbank.org/v2/en/indicator/SH.XPD.CHEX.PP.CD?downloadformat=csv -O health_expenditure.zip
!unzip -n health_expenditure.zip  
!ls -l .
!head API_SH.XPD.CHEX.PP.CD_DS2_en_csv_v2_939222.csv
health_expenditure_df = pd.read_csv('./API_SH.XPD.CHEX.PP.CD_DS2_en_csv_v2_939222.csv', skiprows=4)
os.chdir('..')
health_expenditure_df.head()<jupyter_output><empty_output><jupyter_text>We will use the last available data from recent years (2010 and later).<jupyter_code>recent_year_columns = [str(year) for year in range(2010, 2020)]

health_expenditure_df['CountryHealthExpenditurePerCapitaPPP'] = health_expenditure_df[recent_year_columns].apply(
    lambda row: row[row.last_valid_index()] if row.last_valid_index() else np.nan,
    axis='columns'    
)<jupyter_output><empty_output><jupyter_text>Drop unneeded columns<jupyter_code>health_expenditure_df = health_expenditure_df[['Country Name', 'CountryHealthExpenditurePerCapitaPPP']]
health_expenditure_df.head()<jupyter_output><empty_output><jupyter_text>Now, merge this data with main data by country.<jupyter_code>main_df = merge_with_column_drop(
    main_df,
    health_expenditure_df,
    right_df_column='Country Name'
)
main_df.head()<jupyter_output><empty_output><jupyter_text>Prepare data for training Split data to train, validation and test by date.Construct features (by removing label and non-features data from the original dataframe) and labels for train, validation and test parts.<jupyter_code>last_train_date = pd.Timestamp(2020, 3, 11)
last_eval_date = pd.Timestamp(2020, 3, 24)
last_test_date = pd.Timestamp(2020, 4, 23)

train_df = main_df[main_df['Date'] <= last_train_date].copy()
eval_df = main_df[(main_df['Date'] > last_train_date) & (main_df['Date'] <= last_eval_date)].copy()
test_df = main_df[main_df['Date'] > last_eval_date].copy()

# return features_df, labels
def preprocess_df(df):
    labels = df[['LogNewConfirmedCases', 'LogNewFatalities']].copy()
    features_df = df.drop(
        columns=['Id', 'ForecastId', 'ConfirmedCases', 'LogNewConfirmedCases', 'Fatalities', 'LogNewFatalities', 'Date']
    ).copy()
    
    return features_df, labels

train_features_df, train_labels = preprocess_df(train_df)
eval_features_df, eval_labels = preprocess_df(eval_df)
test_features_df, _ = preprocess_df(test_df)
train_features_df.head()<jupyter_output><empty_output><jupyter_text>Categorical features<jupyter_code>cat_features = ['Province/State', 'Country/Region']<jupyter_output><empty_output><jupyter_text>Models training Now, construct models that predict LogNewConfirmedCases and LogNewFatalities for the next day.We will use RMSE as an objective and evaluation metric for training, that will correspond to RMSLE (Root Mean Square Logarithmic Error) metric for NewConfirmedCases and NewFatalities.$$RMSLE(labels, predictions) = \sqrt{\frac{1}{n}\sum_{i=1}^{n}{\Big(log(labels_i + 1) - log(predictions_i + 1)\Big)^2}}$$Such metric is natural for this task (model should weight errors between 1000 and 1100 and 100000 and 110000 approximately the same).For prediction values in the original scale model prediction output should be transformed back to the original scale with the reverse transformation (can be preformed using numpy.expm1 function).We'll train models for 1000 iterations.<jupyter_code>iterations = 1000<jupyter_output><empty_output><jupyter_text>CatBoost allows Categorical features values can be used as is in DataFrame.<jupyter_code>print ('catboost version', cb.__version__)

catboost_models = {}

for prediction_name in ['LogNewConfirmedCases', 'LogNewFatalities']:
    model = cb.CatBoostRegressor(
        has_time=True,
        iterations=iterations
    )

    model.fit(
        train_features_df,
        train_labels[prediction_name],
        eval_set=(eval_features_df, eval_labels[prediction_name]),
        cat_features=cat_features,
        verbose=100  # print metrics each 100 iterations
    )
    
    catboost_models[prediction_name] = model
    print ('CatBoost: prediction of %s: RMSLE on test = %s' % (prediction_name, model.evals_result_['validation']['RMSE'][-1]))<jupyter_output>catboost version 0.22
Learning rate set to 0.074929
0:	learn: 0.8444130	test: 2.2561020	best: 2.2561020 (0)	total: 79.4ms	remaining: 1m 19s
100:	learn: 0.2849727	test: 0.8850558	best: 0.8850558 (100)	total: 2.2s	remaining: 19.6s
200:	learn: 0.2518107	test: 0.8737440	best: 0.8737323 (198)	total: 4.25s	remaining: 16.9s
300:	learn: 0.2288005	test: 0.8734222	best: 0.8722002 (273)	total: 6.27s	remaining: 14.6s
400:	learn: 0.2078918	test: 0.8766954	best: 0.8720236 (318)	total: 8.37s	remaining: 12.5s
500:	learn: 0.1899614	test: 0.8815886	best: 0.8720236 (318)	total: 10.5s	remaining: 10.4s
600:	learn: 0.1765715	test: 0.8839863	best: 0.8720236 (318)	total: 12.6s	remaining: 8.36s
700:	learn: 0.1650007	test: 0.8869404	best: 0.8720236 (318)	total: 14.7s	remaining: 6.26s
800:	learn: 0.1547044	test: 0.8882405	best: 0.8720236 (318)	total: 16.8s	remaining: 4.17s
900:	learn: 0.1458812	test: 0.8884831	best: 0.8720236 (318)	total: 18.9s	remaining: 2.08s
999:	learn: 0.1379156	test: 0.8917756	best: 0.87202[...]<jupyter_text>Feature importance in the models Let's look at feature importances.For CatBoost use [PredictionValuesChange](https://catboost.ai/docs/concepts/fstr.htmlfstr__regular-feature-importance).We'll print only 25 most important features for clarity.<jupyter_code>for prediction_name in ['LogNewConfirmedCases', 'LogNewFatalities']:
    print ('\nCatBoost: prediction of %s. Feature importance. Type=PredictionValuesChange' % prediction_name)
    print (
        catboost_models[prediction_name].get_feature_importance(
            type=cb.EFstrType.PredictionValuesChange,
            prettified=True
        ).head(25).to_string()
    )<jupyter_output>CatBoost: prediction of LogNewConfirmedCases. Feature importance. Type=PredictionValuesChange
                          Feature Id  Importances
0    LogNewConfirmedCases_prev_day_1    22.109213
1    LogNewConfirmedCases_prev_day_2    17.254841
2       Days_since_ConfirmedCases=10    11.929176
3        Days_since_ConfirmedCases=1    11.340144
4    LogNewConfirmedCases_prev_day_3     6.377358
5                                Day     2.969185
6           Days_since_Fatalities=10     2.720522
7      Days_since_ConfirmedCases=100     1.865486
8    LogNewConfirmedCases_prev_day_4     1.732866
9                 Distance_to_origin     1.259510
10         Days_since_Fatalities=100     1.179583
11  LogNewConfirmedCases_prev_day_10     1.139283
12  LogNewConfirmedCases_prev_day_24     1.134655
13   LogNewConfirmedCases_prev_day_5     1.125969
14                  CountryPop_60-80     0.963053
15       LogNewFatalities_prev_day_3     0.914145
16                  CountryPop_40-60     0.747576
17   [...]<jupyter_text>Create predictions for eval and test data. For forecasting further into the future (more than a single day) we will use an incremental (by days) approach:1. Predict Log New Confirmed Cases and Fatalities for the next unknown day.2. Confirmed Cases and Fatalities for this day are incremented from the previous day based on predicted New Confirmed Cases and Fatalities.3. Use predicted Log New Confirmed Cases and Fatalities for this day to initialize time delay embedding features corresponding to this day for future days.This procedure is then repeated for the new next prediction day and so on. Common function for eval and test:<jupyter_code>def predict_for_dataset(df, features_df, prev_day_df, first_date, last_date, update_features_data):
    df['PredictedLogNewConfirmedCases'] = np.nan
    df['PredictedLogNewFatalities'] = np.nan
    df['PredictedConfirmedCases'] = np.nan
    df['PredictedFatalities'] = np.nan

    for day in pd.date_range(first_date, last_date):
        day_df = df[df['Date'] == day]
        day_features_pool = cb.Pool(features_df.loc[day_df.index], cat_features=cat_features)

        # predict LogNew* data
        for prediction_type in ['LogNewConfirmedCases', 'LogNewFatalities']:
            # prediction is imprecise and can produce negative values, clip them
            df.loc[day_df.index, 'Predicted' + prediction_type] = np.maximum(
                catboost_models[prediction_type].predict(day_features_pool),
                0.0
            )

        day_predictions_df = df.loc[day_df.index][
            location_columns + ['PredictedLogNewConfirmedCases', 'PredictedLogNewFatalities']
        ]

        # update Predicted ConfirmedCases and Fatalities
        for field in ['ConfirmedCases', 'Fatalities']:
            prev_day_field = field if day == first_eval_date else ('Predicted' + field)
            merged_df = day_predictions_df.merge(
                right=prev_day_df[location_columns + [prev_day_field]],
                how='inner',
                on=location_columns
            )

            df.loc[day_df.index, 'Predicted' + field] = merged_df.apply(
                lambda row: row[prev_day_field] + np.rint(np.expm1(row['PredictedLogNew' + field])),
                axis='columns'
            ).values

        if update_features_data:
            # fill time delay embedding features based on this day for next days
            for next_day in pd.date_range(day + pd.Timedelta(days=1), last_date):
                next_day_features_df = features_df[df['Date'] == next_day]

                merged_df = next_day_features_df[location_columns].merge(
                    right=day_predictions_df,
                    how='inner',
                    on=location_columns
                )

                prev_day_idx = (next_day - day).days
                for prediction_type in ['LogNewConfirmedCases', 'LogNewFatalities']:
                    features_df.loc[next_day_features_df.index, prediction_type + '_prev_day_%s' % prev_day_idx] = (
                        merged_df['Predicted' + prediction_type].values
                    )

        # select by day_df.index again to get Predicted* columns
        prev_day_df = df.loc[day_df.index]<jupyter_output><empty_output><jupyter_text>Create predictions for eval data<jupyter_code>prev_day_df = train_df.loc[train_df['Date'] == last_train_date] 
first_eval_date = last_train_date + pd.Timedelta(days=1)

predict_for_dataset(eval_df, eval_features_df, prev_day_df, first_eval_date, last_eval_date, update_features_data=False)

eval_df.head()<jupyter_output><empty_output><jupyter_text>Create predictions for test data.<jupyter_code>prev_day_df = eval_df.loc[eval_df['Date'] == last_eval_date] 
first_test_date = last_eval_date + pd.Timedelta(days=1)

predict_for_dataset(test_df, test_features_df, prev_day_df, first_test_date, last_test_date, update_features_data=True)

test_df.head()<jupyter_output><empty_output><jupyter_text>Plots with predictions Let's plot some example graphs with predictions.Plots will use data from reconcatenated main_df.<jupyter_code>main_df = pd.concat([train_df, eval_df, test_df])
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib


def plot_graph(country_region, province_state, field, log_scale=True):
    location_df = main_df.loc[
        (main_df['Country/Region'] == country_region) & (main_df['Province/State'] == province_state)
    ]

    if province_state:
        title = '%s for %s, %s' % (field, country_region, province_state)
    else:
        title = '%s for %s' % (field, country_region)
    
    plt.figure(figsize = (16, 10))
    plt.suptitle(title, fontsize=14)
    if log_scale:
        plt.yscale('log')

    for sub_field in [field, 'Predicted' + field]:
        plt.plot(location_df['Date'], location_df[sub_field], label=sub_field)


    # add vertical lines splitting train, test and eval parts of the graph
    ax = plt.gca()

    # the x coords of this transformation are data, and the
    # y coord are axes
    transform_for_text = matplotlib.transforms.blended_transform_factory(ax.transData, ax.transAxes)

    plt.axvline(x=last_train_date, color='#000000')
    plt.text(first_eval_date, 0.95, 'eval', transform = transform_for_text)
    plt.axvline(x=last_eval_date, color='#000000')
    plt.text(first_test_date, 0.95, 'test', transform = transform_for_text)


    plt.legend()
    plt.show()

plot_graph('US', 'Kansas', 'ConfirmedCases')
plot_graph('US', 'Kansas', 'Fatalities')
plot_graph('China', 'Hubei', 'ConfirmedCases')
plot_graph('China', 'Hubei', 'Fatalities')
plot_graph('Germany', '', 'ConfirmedCases', log_scale=False)
plot_graph('Germany', '', 'Fatalities', log_scale=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10795826106735663, 'fraction_numerical': 0.03468766881806033, 'mean_word_length': 4.747998038895244, 'pattern_counts': {'":': 0, '<': 183, '<?xml version=': 0, '>': 176, 'https://': 11, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '44', 'dir': 'jupyter-structured-clean-dedup', 'id': '452282', 'license': 'Apache-2.0', 'n_tokens_mistral': 12622, 'n_tokens_neox': 11545, 'n_words': 2826, 'path': 'catboost/tutorials/example_usages/COVID_19.ipynb', 'repo_name': 'jochenater/catboost', 'size': '35171'}