text:
#!/usr/bin/env python2

'''
Description:    Download KNMI zipped ascii data and create a csv file with
                station information.
Author:         Ronald van Haren, NLeSC (r.vanharen@esciencecenter.nl)
Created:        -
Last Modified:  -
License:        Apache 2.0
Notes:          -
'''

from lxml.html import parse
import csv
import urllib2
from lxml import html
import numbers
import json
import os
import utils
from numpy import vstack
from numpy import concatenate as npconcatenate
import argparse
from download_synop.ncdf import knmi_ncdf
import pandas
import re
import glob
from numpy import sort
from numpy import concatenate
import collections
import zipfile
from datetime import datetime
from datetime import date
from datetime import timedelta
from numpy import zeros

class get_knmi_reference_data:
    '''
    description
    '''
    def __init__(self, opts):
        #self.csvfile = opts.csvfile
        self.outputdir = opts.outputdir
        self.keep = opts.keep
        self.check_output_dir()
        self.station = opts.stationid
        self.get_station_ids()
        print (self.stationids)
        self.download_station_data()
        #self.get_station_locations()

    def get_station_ids(self):
        '''
        get all stationids from the KNMI website
        '''
        url = 'http://projects.knmi.nl/klimatologie/metadata/index.html'
        page = parse(url)
        url_metadata = page.xpath(".//table/tr/td/a/@href")
        station_name_id = [c.text for c in page.xpath(".//table/tr/td/a")]
        stationids = [s.split()[0] for s in station_name_id]
        bad_chars = '(){}<>'
        self.rgx = re.compile('[%s]' % bad_chars)
        station_names = [re.sub(self.rgx, '', " ".join(s.split()[1:])) for s in station_name_id]
        if len(self.station)==0:
            self.stationids = [stationids[idx] + ' - ' + station_names[idx] for
                              idx in range(0,len(stationids))]
        else:
          idx = stationids.index(self.station)
          self.stationids = [stationids[idx] + ' - ' + station_names[idx]]

    def download_station_data(self):
        ''''
        download zip files containing csv station data
        (complete time series for all KNMI stations)
        '''
        url = 'http://www.knmi.nl/nederland-nu/klimatologie/uurgegevens'
        page = parse(url)
        # find location of stations on web page
        num_stations = len(page.xpath("/html/body/main/div[2]/div"))
        station_elements = [page.xpath("/html/body/main/div[2]/div["+str(idx)+"]/div/div/div[2]/table/thead/tr/th") for idx in range(0,num_stations)]
        station_names = [re.sub(self.rgx, '', x[0].text) if len(x)>0 else 'ndf'
                         for x in station_elements]
        for stationid in self.stationids:
            div_id = str(station_names.index(stationid))
            relpaths = page.xpath("/html/body/main/div[2]/div["+div_id+"]/div/div/div[2]/table/tbody/tr/td/a/@href")
            for path in relpaths:
                try:
                    fullpath = "http:" + path
                    request = urllib2.urlopen(fullpath)
                    filename = os.path.basename(path)
                    outputfile = os.path.join(self.outputdir, filename)
                    if self.keep:
                        if os.path.exists(outputfile):
                            # check if filesize is not null
                            if os.path.getsize(outputfile) > 0:
                                # file exists and is not null, continue next iteration
                                continue
                            else:
                                # file exists but is null, so remove and redownload
                                os.remove(outputfile)
                    elif os.path.exists(outputfile):
                        os.remove(outputfile)
                    #save
                    output = open(outputfile, "w")
                    output.write(request.read())
                    output.close()
                except urllib2.HTTPError:
                      print "Error downloading file " + fullpath
            # get station location
            lat, lon ,elev, name = self.get_station_location(stationid)
            data = self.read_knmi_data(stationid)
            knmi_ncdf(stationid, lat, lon, elev, data, self.outputdir)

    def get_station_location(self, stationid):
        '''
        write station name, id and location to csv file
        '''
        # get station names for stationids
        url = 'http://projects.knmi.nl/klimatologie/metadata/index.html'
        page = parse(url)
        url_metadata = page.xpath(".//table/tr/td/a/@href")
        station_name_id = [c.text for c in page.xpath(".//table/tr/td/a")]
        station_ids = [s.split()[0] for s in station_name_id]
        station_names = [" ".join(s.split()[1:]) for s in station_name_id]
        #for idx, stationid in enumerate(station_id):
        idx = station_ids.index(stationid.split()[0])
        station_name = station_names[idx]
        station_url = os.path.join(os.path.split(url)[0],
                                   url_metadata[idx])
        page = parse(station_url)
        rows = [c.text for c in page.xpath(".//table/tr/td")]
        idx_position = rows.index('Positie:') + 1
        idx_startdate = rows.index('Startdatum:') + 1
        lat, lon = rows[idx_position].encode('UTF-8').replace(
            '\xc2\xb0','').replace(' N.B. ', ',').replace(
                'O.L.','').strip().split(',')
        lat,lon = self.latlon_conversion(lat,lon)
        idx_elevation = rows.index('Terreinhoogte:') + 1
        elevation = float(rows[idx_elevation].encode('UTF-8').split(' ')[0].replace(',','.'))
        return lat, lon, elevation, station_name

    def latlon_conversion(self, lat, lon):
        '''
        conversion of GPS position to lat/lon decimals
            example string for lat and lon input: "52 11'"
        '''
        # latitude conversion
        latd = lat.replace("'","").split()
        lat = float(latd[0]) + float(latd[1])/60
        # longitude conversion
        lond = lon.replace("'","").split()
        lon = float(lond[0]) + float(lond[1])/60
        return lat,lon

    def check_output_dir(self):
        '''
        check if outputdir exists and create if not
        '''
        if not os.path.exists(self.outputdir):
            os.makedirs(self.outputdir)

    def read_knmi_data(self, stationid):
        '''
        Calculate or load KNMI reference data:
            pickled file exists -> load
            pickled file doesn't exist -> calculate
        '''
        # generate filename of KNMI station
        filenames = sort(glob.glob(os.path.join(self.outputdir, 'uurgeg_' + str(stationid.split()[0]) + '*.zip' )))
        # load all csv files in list of dictionaries
        dicts = [load_knmi_data(filename).csvdata for filename in filenames]
        # merge all dictionaries in a super dictionary
        knmi_data = collections.defaultdict(list)
        for idx in range(0,len(dicts)):
          try:
            knmi_data = dict((k, npconcatenate((knmi_data.get(k), dicts[idx].get(k)))) for k in set(knmi_data.keys() + dicts[idx].keys()))
          except ValueError:
            # cannot concatenate empty arrays
            knmi_data = dict((k, dicts[idx].get(k)) for k in dicts[idx].keys())
        # return dictionary with all variables/time steps
        return knmi_data

class load_knmi_data:
    def __init__(self, filename):
        self.filename = filename
        self.load_file()
        self.process_reference_data()

    def load_file(self):
        '''
        function description
        '''
        # load the zip file
        zipf = zipfile.ZipFile(self.filename)
        # name of csv name in zip file
        txtname = os.path.splitext(os.path.basename(self.filename))[0] + '.txt'
        # read the data in the txt file
        header_found = False
        header_row = 0
        with zipfile.ZipFile.open(zipf, txtname) as f:
          for line in f:
            if '# STN' in line:
              break
            if len(line.strip())>0:
              header_row += 1
        self.csvdata = pandas.read_csv(self.filename, engine='c',
                                       header=header_row,
                                       skipinitialspace=True).to_dict('list')
        # TODO: use parse_dates=['YYYYMMDD']

    def process_reference_data(self):
        '''
        process the reference csv data
        '''
        ## Convert time to datetime.datetime object
        # hours should be 0-23 instead of 1-24
        self.csvdata['HH'] = [item if item!=24 else 0 for item in
                              self.csvdata['HH']]
        # combine date and hour into one string
        dstring = [str(self.csvdata['YYYYMMDD'][idx]) +
                   str(self.csvdata['HH'][idx]).zfill(2) for idx in
                   range(0,len(self.csvdata['YYYYMMDD']))]
        # create datetime object
        self.csvdata['datetime'] = [datetime.strptime(
            str(item), ('%Y%m%d%H')) for item in dstring]
        # Correct conversion of the datestring
        # the date of the night -> HH=24 is HH=0 on the next day!
        self.csvdata['datetime'] = [c+ timedelta(days=1) if c.hour==0 else
                                    c for c in self.csvdata['datetime']]
        # rain is (-1 for <0.05 mm), set to 0
        self.csvdata['RH'] = [0 if item == -1 else
                                  item for item in self.csvdata['RH']]
        # process all variables that need to be divided by 10
        for variable in ['T10', 'T', 'RH', 'FF', 'TD']:
            # T10: temperature at 10 cm height, divide by 10 to convert to degC
            # T: temperature at 1.50 m height, divide by 10 to convert to degC
            # RH: rain
            # FF: wind speed
            self.csvdata[variable] = [round(0.1 * item, 1) if item != -999 else
                                item for item in self.csvdata[variable]]
        # multiply pressure to convert to Pascal
        self.csvdata['P'] = [round(float(10 * item), 3) if
                             item != -999 else item for item in
                             self.csvdata['P']]
        # SWD
        self.csvdata['Q'] = [round(float(10000 * item)/3600, 5) if
                             item != -999 else item for item in
                             self.csvdata['Q']]

reason: from numpy


text:
# Uses a learnable embedding for the countries which is fed as input to both encoders (at the respective dense layers).

# Copyright 2020 (c) Cognizant Digital Business, Evolutionary AI. All rights reserved. Issued under the Apache 2.0 License.

import os
import urllib.request

from keras.layers.core import RepeatVector

# Suppress noisy Tensorflow debug logging
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'

# noinspection PyPep8Naming
import keras.backend as K
import numpy as np
import pandas as pd
from keras.callbacks import EarlyStopping
from keras.constraints import Constraint
from keras.layers import Concatenate
from keras.layers import Dense
from keras.layers import Embedding
from keras.layers import Input
from keras.layers import LSTM
from keras.layers import Reshape
from keras.layers import Lambda
from keras.models import Model

import ongoing.predictors.base as base
from ongoing.predictors.base import BasePredictor

# See https://github.com/OxCGRT/covid-policy-tracker
DATA_URL = "https://raw.githubusercontent.com/OxCGRT/covid-policy-tracker/master/data/OxCGRT_latest.csv"

ROOT_DIR = os.path.dirname(os.path.abspath(__file__))
DATA_PATH = os.path.join(ROOT_DIR, os.pardir, os.pardir, 'data')
TEMPERATURE_DATA_FILE_PATH = os.path.join(DATA_PATH, "temperature_data.csv")
TEMPERATURE_COLUMN = 'temp,C'
HOLIDAY_COLUMN = 'Holiday'

NB_LOOKBACK_DAYS = 21
NB_TEST_DAYS = 14
WINDOW_SIZE = 7
US_PREFIX = "United States / "
NUM_TRIALS = 10
NUM_EPOCHS = 1000
LSTM_SIZE = 32
EMBED_SIZE = 4
NPI_DELAY = 0
TEMP_SCALE = 20.  # divide temperature values by 20 so they're roughly in the range 0-2
AVG_EARTH_TEMP = 16./TEMP_SCALE  # average temperature on earth (used to predict for locations where temperature data is missing)
DISCARD_DATA_BFR = np.datetime64("2020-02-01")

class Positive(Constraint):

    def __call__(self, w):
        return K.abs(w)

# Functions to be used for lambda layers in model
def _combine_r_and_d(x):
    r, d = x
    return r * (1. - d)

class tempGeoLSTMPredictor(BasePredictor):
    """
    A class that computes a fitness for Prescriptor candidates.
    """

    def __init__(self, path_to_model_weights=None, path_to_geos=None,
                 use_embedding=True, embed_size=EMBED_SIZE, lstm_size=LSTM_SIZE, nb_lookback_days=NB_LOOKBACK_DAYS,
                 nb_test_days=NB_TEST_DAYS, window_size=WINDOW_SIZE, npi_delay=NPI_DELAY,
                 num_trials=NUM_TRIALS, num_epochs=NUM_EPOCHS, seed=base.SEED):

        super().__init__(seed=seed)
        self.use_embedding=use_embedding
        self.embed_size = embed_size
        self.lstm_size = lstm_size
        self.nb_lookback_days = nb_lookback_days
        self.nb_test_days = nb_test_days
        self.window_size = window_size
        self.npi_delay = npi_delay
        self.num_trials = num_trials
        self.num_epochs = num_epochs

        if path_to_model_weights:
            # Load model weights
            nb_context = 4  # time series of new cases rate, death rate, temperature, and holidays are used
            nb_action = len(base.NPI_COLUMNS)
            if self.use_embedding:
                with open(path_to_geos, 'r') as f:
                    self.geos = [g.rstrip() for g in f.readlines()]
            else:
                self.geos = None
            self.predictor, _ = self._construct_model(nb_context=nb_context,
                                                      nb_action=nb_action,
                                                      embed_size=self.embed_size,
                                                      lstm_size=self.lstm_size,
                                                      nb_lookback_days=self.nb_lookback_days)
            self.predictor.load_weights(path_to_model_weights)

        # read and preprocess temperature data
        self.temp_df = pd.read_csv(TEMPERATURE_DATA_FILE_PATH,
                                   parse_dates=['Date'],
                                   encoding="ISO-8859-1",
                                   dtype={"RegionName": str,
                                          "RegionCode": str},
                                   error_bad_lines=False)
        self.temp_df["GeoID"] = np.where(self.temp_df["RegionName"].isnull(),
                                         self.temp_df["CountryName"],
                                         self.temp_df["CountryName"] + ' / ' + self.temp_df["RegionName"])
        self.temp_df[TEMPERATURE_COLUMN] = self.temp_df[TEMPERATURE_COLUMN]/TEMP_SCALE

        self.country_samples = None  # will be set when fit() or predict() are called

    def predict(self, data=None, start_date=None, end_date=None):
        if self.train_df is None:
            raise Exception("train_df must be defined before calling predict()")

        if data is None:
            data = self.test_df
        if start_date is None:
            start_date = pd.to_datetime(data.Date.min(), format='%Y-%m-%d')
        if end_date is None:
            end_date = pd.to_datetime(data.Date.max(), format='%Y-%m-%d')

        train_df = pd.merge(self.train_df, self.temp_df, on=['CountryName', 'RegionName', 'GeoID', 'Date'], how='left')
        train_df[TEMPERATURE_COLUMN] = train_df[TEMPERATURE_COLUMN].fillna(AVG_EARTH_TEMP)
        train_df[HOLIDAY_COLUMN] = train_df[HOLIDAY_COLUMN].fillna(0)
        self.country_samples = self._create_country_samples(train_df,
                                                            list(self.train_df.GeoID.unique()),
                                                            self.nb_lookback_days,
                                                            self.npi_delay,
                                                            self.nb_test_days)

        nb_days = (end_date - start_date).days + 1

        # Prepare the output
        forecast = {"GeoID": [],
                    "CountryName": [],
                    "RegionName": [],
                    "Date": [],
                    "PredictedDailyTotalCases": [],
                    "PredictedDailyNewCases": [],
                    "PredictedDailyTotalDeaths": [],
                    "PredictedDailyNewDeaths": []}

        # For each requested geo
        geos = data.GeoID.unique()
        for g in geos:
            if self.use_embedding and g not in self.geos:
                # the model was not trained for this geo: return zeroes
                print("WARNING: The model was not trained for {}".format(g))
                pred_total_cases = [0] * nb_days
                pred_new_cases = [0] * nb_days
                pred_total_deaths = [0] * nb_days
                pred_new_deaths = [0] * nb_days
                geo_start_date = start_date
            else:
                cdf = self.train_df[self.train_df.GeoID == g]

                if (len(cdf) == 0) or (g not in self.country_samples.keys()):
                    # we don't have historical data for this geo: return zeroes
                    print("WARNING: No historical data for {}".format(g))
                    pred_total_cases = [0] * nb_days
                    pred_new_cases = [0] * nb_days
                    pred_total_deaths = [0] * nb_days
                    pred_new_deaths = [0] * nb_days
                    geo_start_date = start_date
                else:
                    last_known_date = cdf.Date.max()
                    # Start predicting from start_date, unless there's a gap since last known date
                    geo_start_date = min(last_known_date + np.timedelta64(1, 'D'), start_date)
                    npis_gdf = data[(data.Date >= geo_start_date - pd.Timedelta(days=self.npi_delay)) & (data.Date <= end_date - pd.Timedelta(days=self.npi_delay))]
                    if geo_start_date.year == end_date.year:
                        temp_gdf = self.temp_df[(self.temp_df.Date >= geo_start_date.replace(year=2020)) & (self.temp_df.Date <= end_date.replace(year=2020))]
                    else:
                        temp_gdf1 = self.temp_df[(self.temp_df.Date >= geo_start_date.replace(year=2020)) & (self.temp_df.Date <= pd.to_datetime('2020-12-31', format='%Y-%m-%d'))]
                        temp_gdf2 = self.temp_df[(self.temp_df.Date >= pd.to_datetime('2020-01-01', format='%Y-%m-%d')) & (self.temp_df.Date <= end_date.replace(year=2020))]
                        temp_gdf = pd.concat([temp_gdf1, temp_gdf2])

                    pred_total_cases, pred_new_cases, pred_total_deaths, pred_new_deaths = self._get_new_cases_preds(cdf, g, npis_gdf, temp_gdf)

            # Append forecast data to results to return
            country = data[data.GeoID == g].iloc[0].CountryName
            region = data[data.GeoID == g].iloc[0].RegionName
            for i, (ptot_cases, pnew_cases, ptot_deaths, pnew_deaths) in enumerate(zip(pred_total_cases, pred_new_cases, pred_total_deaths, pred_new_deaths)):
                forecast["GeoID"].append(g)
                forecast["CountryName"].append(country)
                forecast["RegionName"].append(region)
                current_date = geo_start_date + pd.offsets.Day(i)
                forecast["Date"].append(current_date)
                forecast["PredictedDailyTotalCases"].append(ptot_cases)
                forecast["PredictedDailyNewCases"].append(pnew_cases)
                forecast["PredictedDailyTotalDeaths"].append(ptot_deaths)
                forecast["PredictedDailyNewDeaths"].append(pnew_deaths)

        forecast_df = pd.DataFrame.from_dict(forecast)
        # Return only the requested predictions
        return forecast_df[(forecast_df.Date >= start_date) & (forecast_df.Date <= end_date)]

    def _get_new_cases_preds(self, c_df, g, npis_df, temp_df):
        cdf = c_df[c_df.ConfirmedCases.notnull()]
        initial_context_input = self.country_samples[g]['X_test_context'][-1]
        initial_action_input = self.country_samples[g]['X_test_action'][-1]
        if self.use_embedding:
            country_id = np.array([self.geos.index(g)])
        else:
            country_id = np.array([0.])
        # Predictions with passed npis
        cnpis_df = npis_df[npis_df.GeoID == g]
        npis_sequence = np.array(cnpis_df[base.NPI_COLUMNS])
        ctemp_df = temp_df[temp_df.GeoID == g]
        if ctemp_df.empty:
            print("WARNING: No temperature data available for {}".format(g))
            temp_sequence = AVG_EARTH_TEMP*np.ones(npis_sequence.shape[0])
            holiday_sequence = np.zeros(npis_sequence.shape[0])
        else:
            temp_sequence = np.array(ctemp_df[TEMPERATURE_COLUMN])
            holiday_sequence = np.array(ctemp_df[HOLIDAY_COLUMN])
        # Get the predictions with the passed NPIs
        preds = self._roll_out_predictions(self.predictor,
                                           initial_context_input,
                                           initial_action_input,
                                           country_id,
                                           npis_sequence,
                                           temp_sequence,
                                           holiday_sequence)
        preds_cases = preds[:,0]
        preds_deaths = preds[:,1]
        # Gather info to convert to total cases
        prev_confirmed_cases = np.array(cdf.ConfirmedCases)
        prev_new_cases = np.array(cdf.NewCases)
        initial_total_cases = prev_confirmed_cases[-1]
        pop_size = np.array(cdf.Population)[-1]  # Population size doesn't change over time
        prev_confirmed_deaths = np.array(cdf.ConfirmedDeaths)
        prev_new_deaths = np.array(cdf.NewDeaths)
        initial_total_deaths = prev_confirmed_deaths[-1]

        # Compute predictor's forecast
        pred_total_cases, pred_new_cases = base.convert_ratios_to_total_cases(
            preds_cases,
            self.window_size,
            prev_new_cases,
            initial_total_cases,
            pop_size)

        # Compute predictor's deaths forecast
        pred_total_deaths, pred_new_deaths = base.convert_ratios_to_total_deaths(
            preds_deaths,
            self.window_size,
            prev_new_deaths,
            initial_total_deaths)

        return pred_total_cases, pred_new_cases, pred_total_deaths, pred_new_deaths

    @staticmethod
    def _create_country_samples(df: pd.DataFrame, geos: list, nb_lookback_days: int, npi_delay: int, nb_test_days: int) -> dict:
        """
        For each country, creates numpy arrays for Keras
        :param df: a Pandas DataFrame with historical data for countries (the "Oxford" dataset)
        :param geos: a list of geo names
        :return: a dictionary of train and test sets, for each specified country
        """
        context_column = ['PredictionRatio', 'DeathRatio', TEMPERATURE_COLUMN, HOLIDAY_COLUMN]
        action_columns = base.NPI_COLUMNS
        outcome_column = ['PredictionRatio', 'DeathRatio']
        country_samples = {}
        for i, g in enumerate(geos):
            cdf = df[df.GeoID == g]
            cdf = cdf[cdf.ConfirmedCases.notnull()]
            context_data = np.array(cdf[context_column])
            action_data = np.array(cdf[action_columns])
            outcome_data = np.array(cdf[outcome_column])
            context_samples = []
            action_samples = []
            outcome_samples = []
            nb_total_days = outcome_data.shape[0]
            for d in range(nb_lookback_days+npi_delay, nb_total_days):
                context_samples.append(context_data[d - nb_lookback_days: d])
                action_samples.append(action_data[d-npi_delay - nb_lookback_days: d-npi_delay])
                outcome_samples.append(outcome_data[d])
            if len(outcome_samples) > 0:
                X_context = np.stack(context_samples, axis=0)
                X_action = np.stack(action_samples, axis=0)
                X_country = i*np.ones(X_context.shape[0])
                y = np.stack(outcome_samples, axis=0)
                country_samples[g] = {
                    'X_context': X_context,
                    'X_action': X_action,
                    'X_country': X_country,
                    'y': y,
                    'X_test_context': X_context[-nb_test_days:],
                    'X_test_action': X_action[-nb_test_days:],
                    'X_test_country': X_country[-nb_test_days:],
                    'y_test': y[-nb_test_days:],
                }
        return country_samples

    # Function for performing roll outs into the future
    def _roll_out_predictions(self, predictor, initial_context_input, initial_action_input, country_id, future_action_sequence, future_temperature_sequence, future_holiday_sequence):
        nb_roll_out_days = future_action_sequence.shape[0]
        pred_output = np.zeros((nb_roll_out_days, 2))
        context_input = np.expand_dims(np.copy(initial_context_input), axis=0)
        action_input = np.expand_dims(np.copy(initial_action_input), axis=0)
        country_input = np.expand_dims(np.copy(country_id), axis=0)
        for d in range(nb_roll_out_days):
            action_input[:, :-1] = action_input[:, 1:]
            # Use the passed actions
            action_sequence = future_action_sequence[d]
            action_input[:, -1] = action_sequence
            if self.use_embedding:
                inputs = [context_input, action_input, country_input]
            else:
                inputs = [context_input, action_input]

            pred = predictor.predict(inputs)
            pred_output[d] = pred[-1]
            context_input[:, :-1] = context_input[:, 1:]
            context_input[:, -1, 0:2] = pred[-1]
            context_input[:, -1, 2] = future_temperature_sequence[d]
            context_input[:, -1, 3] = future_holiday_sequence[d]
        return pred_output

    def fit(self):
        if self.train_df is None:
            raise Exception("train_df must be defined bfr calling predict()")

        # merge the two dataframes (keep only rows where new cases rate and temperature are available)
        train_df = pd.merge(self.train_df, self.temp_df, on=['CountryName', 'RegionName', 'GeoID', 'Date'], how='inner')
        train_df = train_df[train_df.Date >= DISCARD_DATA_BFR]
        self.country_samples = self._create_country_samples(train_df,
                                                            list(train_df.GeoID.unique()),
                                                            self.nb_lookback_days,
                                                            self.npi_delay,
                                                            self.nb_test_days)
        self.geos = list(self.country_samples.keys())

        # Aggregate data for training
        all_X_context_list = [self.country_samples[c]['X_context']
                              for c in self.country_samples]
        all_X_action_list = [self.country_samples[c]['X_action']
                             for c in self.country_samples]
        all_X_country_list = [self.country_samples[c]['X_country']
                              for c in self.country_samples]
        all_y_list = [self.country_samples[c]['y']
                      for c in self.country_samples]
        X_context = np.concatenate(all_X_context_list)
        X_action = np.concatenate(all_X_action_list)
        X_country = np.concatenate(all_X_country_list)
        y = np.concatenate(all_y_list)

        # Clip outliers
        MIN_VALUE = 0.
        MAX_VALUE = 2.
        X_context = np.clip(X_context, MIN_VALUE, MAX_VALUE)
        y = np.clip(y, MIN_VALUE, MAX_VALUE)

        X_context, X_action, X_country, y = self._permute_data(X_context, X_action, X_country, y)
        self.predictor, training_model = self._construct_model(nb_context=X_context.shape[-1],
                                                      nb_action=X_action.shape[-1],
                                                      embed_size=self.embed_size,
                                                      lstm_size=self.lstm_size,
                                                      nb_lookback_days=self.nb_lookback_days)
        history = self._train_model(training_model, X_context, X_action, X_country, y, epochs=self.num_epochs, verbose=0)
        top_epoch = np.argmin(history.history['val_loss'])
        train_loss = history.history['loss'][top_epoch]
        val_loss = history.history['val_loss'][top_epoch]
        print('Train Loss:', train_loss)
        print('Val Loss:', val_loss)

    # Shuffling data prior to train/val split
    def _permute_data(self, X_context, X_action, X_country, y):
        p = np.random.permutation(y.shape[0])
        X_context = X_context[p]
        X_action = X_action[p]
        X_country = X_country[p]
        y = y[p]
        return X_context, X_action, X_country, y

    # Construct model
    def _construct_model(self, nb_context, nb_action, embed_size=10, lstm_size=32, nb_lookback_days=21):
        if self.use_embedding:
            # Create country embedding
            country_id = Input(shape=(1,),
                               name='country_id')
            emb = Embedding(len(self.geos), embed_size)(country_id)
            emb = Reshape((embed_size,))(emb)

        # Create context encoder
        context_input = Input(shape=(nb_lookback_days, nb_context),
                              name='context_input')
        x = LSTM(lstm_size,
                 return_sequences=False,
                 name='context_lstm')(context_input)
        if self.use_embedding:
            x = Concatenate(axis=1)([x, emb])  # concatenate the output of the LSTM with the country embedding prior to the dense layer
        context_output = Dense(units=1,
                               activation='softplus',
                               name='context_dense')(x)

        # Create action encoder
        # Every aspect is monotonic and nonnegative except final bias
        action_input = Input(shape=(nb_lookback_days, nb_action),
                             name='action_input')
        x = LSTM(units=lstm_size,
                 kernel_constraint=Positive(),
                 recurrent_constraint=Positive(),
                 bias_constraint=Positive(),
                 return_sequences=False,
                 name='action_lstm')(action_input)
        if self.use_embedding:
            x = Concatenate(axis=1)([x, emb])  # concatenate the output of the LSTM with the country embedding prior to the dense layer
        action_output = Dense(units=2,
                              activation='sigmoid',
                              # kernel_constraint=Positive(),
                              name='action_dense')(x)

        # Create prediction model
        model_output = Lambda(_combine_r_and_d, name='prediction')(
            [context_output, action_output])
        if self.use_embedding:
            inputs = [context_input, action_input, country_id]
        else:
            inputs = [context_input, action_input]

        model = Model(inputs=inputs,
                      outputs=[model_output])
        model.compile(loss='mae', optimizer='adam')

        # Create training model, which includes loss to measure
        # variance of action_output predictions
        training_model = Model(inputs=inputs,
                               outputs=[model_output])
        training_model.compile(loss='mae',
                               optimizer='adam')

        return model, training_model

    # Train model
    def _train_model(self, training_model, X_context, X_action, X_country, y, epochs=1, verbose=0):
        early_stopping = EarlyStopping(patience=20,
                                       restore_best_weights=True)
        if self.use_embedding:
            inputs = [X_context, X_action, X_country]
        else:
            inputs = [X_context, X_action]
        history = training_model.fit(inputs, [y],
                                     epochs=epochs,
                                     batch_size=32,
                                     validation_split=0.1,
                                     callbacks=[early_stopping],
                                     verbose=verbose)
        return history

    def save_model(self, path_to_weights, path_to_country_list=None):
        self.predictor.save_weights(path_to_weights)
        if self.use_embedding and path_to_country_list is not None:
            with open(path_to_country_list, 'w') as f:
                f.writelines("{}\n".format(g) for g in self.geos)

if __name__ == '__main__':
    # Run all test cases
    # model = tempGeoLSTMPredictor('./ongoing/predictors/tempgeolstm/models/model.h5', './ongoing/predictors/tempgeolstm/models/countries.txt')
    model = tempGeoLSTMPredictor(use_embedding=False)
    model.evaluate()
    # model.choose_train_test_split(start_date=pd.to_datetime('2021-01-01', format='%Y-%m-%d'),
    #                               end_date=pd.to_datetime('2021-01-01', format='%Y-%m-%d'),
    #                               update_data=True)
    # model.fit()
    # model.save_model('./ongoing/predictors/tempgeolstm/models/model.h5', './ongoing/predictors/tempgeolstm/models/countries.txt')

ext: py

max_stars_repo_name: bradyneal/covid-xprize-comp

max_issues_repo_name: bradyneal/covid-xprize-comp

max_forks_repo_name: bradyneal/covid-xprize-comp

reason: import numpy


text:
import siconos.numerics as sn
import os.path
import scipy.sparse
import scipy.sparse.linalg
import numpy as np
from siconos.tests_setup import working_dir

# this could be changed at install time
data_dir = working_dir + '/data/'
# solvers = [sn.SICONOS_GLOBAL_FRICTION_3D_NSGS, sn.SICONOS_GLOBAL_FRICTION_3D_NSN_AC]
solvers = (sn.SICONOS_GLOBAL_FRICTION_3D_NSGS,)
solvers_reduced1 = (sn.SICONOS_FRICTION_3D_NSGS, sn.SICONOS_FRICTION_3D_NSN_AC)
solvers_reduced2 = (sn.SICONOS_FRICTION_3D_NSN_AC,)  # sn.SICONOS_FRICTION_3D_NSN_FB)
solvers_reduced3 = (sn.SICONOS_FRICTION_3D_NSGS,)

def condensed_from_global(fcp):
    # spsolve expect the indices to be cint aka 32 bits int
    # Hence, we do some magic to cirumvent those issues.
    Mcoo = scipy.sparse.coo_matrix(fcp.M)
    Mcsc = scipy.sparse.csc_matrix(Mcoo)
    Hcoo = scipy.sparse.coo_matrix(fcp.H)
    Hcsc = scipy.sparse.csc_matrix(Hcoo)
    WW = Hcsc.T.dot(scipy.sparse.linalg.spsolve(Mcsc, Hcsc))
    qprime = fcp.H.T.dot(scipy.sparse.linalg.spsolve(Mcsc, fcp.q))

    fcp_reduced = sn.FrictionContactProblem()
    fcp_reduced.dimension = fcp.dimension
    fcp_reduced.numberOfContacts = fcp.numberOfContacts

    # this is a hack to deal with the inability of fc3d solvers to work with
    # sparse matrices
    _, Wsbm = sn.SBM_from_csparse(3, WW)
    fcp_reduced.M = Wsbm
    fcp_reduced.mu = fcp.mu
    fcp_reduced.q = fcp.b + qprime
    return fcp_reduced

def solve_reduced(fcp, solver_reduced):
    SO_reduced = sn.SolverOptions(solver_reduced)
    SO_reduced.iparam[0] =100000
    SO_reduced.dparam[0] = np.sqrt(fcp.numberOfContacts)*1e-9
    size_reaction = fcp.numberOfContacts * 3
    reaction_reduced = np.zeros((size_reaction,))
    velocities_reduced = np.zeros((size_reaction,))

    return sn.fc3d_driver(fcp, reaction_reduced, velocities_reduced, SO_reduced)

def solve_global(fcp, solver):
    SO = sn.SolverOptions(solver)
    SO.dparam[0] = np.sqrt(fcp.numberOfContacts)*1e-9

    n, m = fcp.H.shape
    size_reaction = m
    reaction = np.zeros((size_reaction,))
    velocities = np.zeros((size_reaction,))
    global_velocities = np.zeros((n,))
    return sn.gfc3d_driver(fcp, reaction, velocities, global_velocities, SO)

def test_gfc3d():
    data_files = ('LMGC_GFC3D_CubeH8.hdf5', 'LMGC_GlobalFrictionContactProblem00046.hdf5')
    mark_as_failed = False

    for d in data_files:
        full_path = data_dir + d
        if os.path.isfile(full_path):

            fcp = sn.globalFrictionContact_fclib_read(full_path)
            for s in solvers:
                res = solve_global(fcp, s)
                if res:
                    print('Solver {:} on problem {:} failed with info = {:}'.format(sn.solver_options_id_to_name(s), d, res))
                    mark_as_failed = True
                else:
                    print('Solver {:} on problem {:} is ok'.format(sn.solver_options_id_to_name(s), d))

            fcp_reduced = condensed_from_global(fcp)
            for s in solvers_reduced1:
                res = solve_reduced(fcp_reduced, s)
                if res:
                    print('Solver {:} on problem {:} in reduced form failed with info = {:}'.format(sn.solver_options_id_to_name(s), d, res))
                    mark_as_failed = True
                else:
                    print('Solver {:} on problem {:} is ok'.format(sn.solver_options_id_to_name(s), d))

    assert mark_as_failed is False

def test_fc3d():
    data_files = ['Capsules-i125-1213.hdf5']
    mark_as_failed = False

    for d in data_files:
        full_path = data_dir + d
        if os.path.isfile(full_path):
            sn.numerics_set_verbose(1)
            fcp = sn.frictionContact_fclib_read(full_path)
            for s in solvers_reduced3:
                res = solve_reduced(fcp, s)
                if res:
                    print('Solver {:} on problem {:} failed with info = {:}'.format(sn.solver_options_id_to_name(s), d, res))
                    mark_as_failed = True
                else:
                    print('Solver {:} on problem {:} is ok'.format(sn.solver_options_id_to_name(s), d))

    assert mark_as_failed is False

max_stars_repo_name: ljktest/siconos

max_issues_repo_name: ljktest/siconos

max_forks_repo_name: ljktest/siconos

reason: import numpy,import scipy


text:
import time
import json
import warnings

import numpy as np
import pandas as pd

from pysr import pysr, best, best_callable
from baseline_utils import nostdout, processDataFiles, relativeErr

warnings.filterwarnings("ignore")

def generate_results(file_path, save_path):

    results_data = []
    test_set = processDataFiles([file_path])
    test_set = test_set.strip().split("\n")

    for idx, sample in enumerate(test_set):
        t = json.loads(sample)
        X = np.array(t["X"])
        y = np.array(t["Y"])
        X_test = np.array(t["XT"])
        y_test = np.array(t["YT"])
        raw_eqn = t["EQ"]
        skeleton = t["Skeleton"]

        try:
            start_time = time.time()
            with nostdout():
                model = pysr(X, y, npop=100, niterations=10,
                             binary_operators=[
                                 "plus", "mult", "div", "pow", "sub"
                             ],
                             unary_operators=[
                                 "neg", "exp", "log_abs", "sqrt_abs", "sin"
                             ],
                             temp_equation_file=True,
                             verbosity=0
                             )
            equation_pred = best(model)
            model_pred = best_callable(model)
            pred_y = model_pred(X_test)
            train_time = time.time() - start_time
            err = relativeErr(y_test, pred_y)
        except ValueError:
            equation_pred = "NA"
            pred_y = "NA"
            err = "NA"
            train_time = "NA"

        results_data.append({
            "test_index": idx,
            "true_equation": raw_eqn,
            "true_skeleton": skeleton,
            "predicted_equation": equation_pred,
            "predicted_y": pred_y,
            "rel_err": err,
            "dsr_time": train_time
        })

    results = pd.DataFrame(results_data)
    results.to_csv(save_path, index=False)

if __name__ == "__main__":

    data_path = "./data/0_1_0_12062021_083325.json"
    save_path = "./results/0_1_0_12062021_083325_pysr.csv"

    generate_results(data_path, save_path)

max_stars_repo_name: t-taniai/symbolicgpt

max_issues_repo_name: t-taniai/symbolicgpt

max_forks_repo_name: t-taniai/symbolicgpt