<jupyter_start><jupyter_text>PyData 2018 NYC: End to End Data Science Without Leaving the GPURandy Zwitch Senior Developer Advocate, [OmniSci](https://omnisci.com) (formerly MapD) October 18, 2018Twitter: [@randyzwitch](https://twitter.com/randyzwitch)  Agenda- Personal Introduction- What is OmniSci? (Also...we're hiring!)- GOAI/RAPIDS- Demo  0. Source data stored in OmniSci GPU databaseAs the starting point for our "End to End" journey, the data are stored in OmniSci, a GPU-accelerated relation database and analytics platform. The data are hourly measurements of electricity demand in the PJM Interconnection transmission area from 1993-01-01 to present. (Brief OmniSci demo) [![Foo](omnisci_screenshot.png)](http://172.16.17.32:9092//dashboard/49)  1. Connecting to OmniSci via PythonOmniSci is part of the [GPU Open Analytics Initiative](http://gpuopenanalytics.com//), with the mission of allowing a data scientist to "explore data, train machine learning algorithms, and build applications while primarily staying on the GPU". The foundation of this workflow is the [GPU Dataframe](https://github.com/rapidsai/libgdf), which OmniSci supports returning results in when using [pymapd](https://github.com/omnisci/pymapd).GOAI is an idealized state and very much a work-in-progress. As such, it is _possible_ to do an end-to-end workflow without ever leaving the GPU. However, in the case of this tutorial, I have worked around some gaps in functionality by copying the data from one GPU to another one with the desired properties. Eventually, these edges will be smoothed out, providing a seamless workflow.<jupyter_code>#pymapd native connector to OmniSci
#numpy.dtype warning a cython issue, safe to ignore
import pymapd
import pandas as pd
from pygdf.dataframe import DataFrame
from credentials import password, user, database #user needs to define themselves
#connect using Ibis (which in turn, uses pymapd behind the scenes)
conn = pymapd.connect(host = "localhost", port = 9091, password = password, dbname = database, user = user) 

#show sample of data
#pymapd is DBI-compliant, so you can pass a query string to pandas
#passing through pandas causes a CPU copy, which is ok because I'm using pandas for pretty-printing
pd.read_sql("select * from hrl_load_metered limit 5", conn)<jupyter_output><empty_output><jupyter_text>2. Pulling Data SampleIn order to build a simple model forecasting electricity demand, we need to get the data from OmniSci into Python. While we're pulling the data, we can define a few features like day of week, month, hour and year.<jupyter_code>#Limiting data to more recent years based on domain knowledge
#Definition of 'RTO' area grew considerably in early 2000s, definition stabilized around then
df = conn.select_ipc("""select
datetime_beginning_ept,
extract(DOW from datetime_beginning_ept) as day_of_week,
extract(MONTH from datetime_beginning_ept) as month_,
extract(HOUR from datetime_beginning_ept) as hour_,
extract(YEAR from datetime_beginning_ept) as year_,
mw,
cast(1 as double) as intercept
from hrl_load_metered
where zone = 'RTO' and datetime_beginning_ept >= '2005-07-01 00:00:00'
order by 1
""")

#DataFrame.from_pandas() a work-around for https://github.com/rapidsai/pygdf/issues/286
rto_demand_gdf = DataFrame.from_pandas(df)

#pretty printing
rto_demand_gdf.head().to_pandas()<jupyter_output>/home/mapdadmin/.local/lib/python3.5/site-packages/pyarrow/util.py:35: FutureWarning: pyarrow.frombuffer is deprecated as of 0.9.0, please use py_buffer instead
  warnings.warn(msg, FutureWarning)<jupyter_text>3. Simple Curve Fitting Using Linear RegressionFrom the crossfiltering demonstration using Immerse, as well as the heatmap, we can see that electricity demand is well described just by knowing the time of day and day of the week. We can validate this hypothesis using linear regression. Because the data are already on the GPU, and since the math of linear regression can be represented in basic linear algebra terms, we'll do that instead of using a machine learning library. Instead, we'll use [CuPy](https://cupy.chainer.org/) to get the coefficients and validate the model fit.  3a. One-hot encoding for day of week, month, hour, yearAlthough month and day of the week are stored as numbers, we don't want to use them that way. Rather, we create 1/0 variables using the `gdf.one_hot_encoding()` functionality.<jupyter_code># Do one-hot encoding on day of week and month
rto_demand_gdf_1 = rto_demand_gdf.one_hot_encoding('day_of_week', 'dow', range(0,6))
rto_demand_gdf_2 = rto_demand_gdf_1.one_hot_encoding('month_', 'month', range(1,12))
rto_demand_gdf_encoded = rto_demand_gdf_2.one_hot_encoding('hour_', 'hour', range(0,23))

# pretty printing
rto_demand_gdf_encoded.head().to_pandas()
type(rto_demand_gdf_encoded)<jupyter_output><empty_output><jupyter_text>3b. Split into training, validation setsSince this is a time-series problem, we'll estimate our model based on a given set of years, then test it against future years.<jupyter_code>#On-GPU processing
train = rto_demand_gdf_encoded.query('year_ <= 2016')
val = rto_demand_gdf_encoded.query('year_ > 2016')
import cupy
from cupy import asarray, dot
from cupy.linalg import inv, qr

#list of predictors
indep = list(train.columns)[6:]

#Predictor matrix: convert to cupy.core.core.ndarray (numpy array on GPU)
X_train = asarray(train.loc[:, indep].as_matrix())
X_val = asarray(val.loc[:, indep].as_matrix())

#Target matrix: convert to cupy.core.core.ndarray (numpy array on GPU)
y_train = asarray(train.loc[:, ['mw']].as_matrix())
y_val = asarray(val.loc[:, ['mw']].as_matrix())<jupyter_output><empty_output><jupyter_text>3c. Calculate model parameters using CuPy and Normal Equation<jupyter_code>#use timeit to highlight speed of GPU calculation
#linear regression generally not calculated this way, because taking inverse computationally expensive
%timeit inv(X_train.T.dot(X_train)).dot(X_train.T).dot(y_train)
#Model coefficients (timeit doesn't save results)
coef = inv(X_train.T.dot(X_train)).dot(X_train.T).dot(y_train)
#Predicted values on GPU
yhat_train = X_train.dot(coef)
#Calculate RMSE on GPU
rmse = cupy.sqrt(cupy.sum(cupy.square(yhat_train - y_train)) / len(y_train))
rmse
#Calculate MAPE on GPU
mape = cupy.mean(cupy.abs((y_train - yhat_train) / y_train))
#9.2% error not bad for such simplistic model, validation sample likely worse
mape<jupyter_output><empty_output><jupyter_text>4. Visualizing Model Performancepymapd provides a method `render_vega` that can be used to render data on the GPU (this is the basis of how the Immerse interface works). Additionally, we can upload the data back into OmniSci and build an Immerse dashboard.For this example, we'll just use local CPU rendering since the data are modest in size.<jupyter_code>lc = pd.DataFrame(cupy.concatenate((yhat_train, y_train), axis=1).get(), columns=["forecasted", "actual"])
lc.index = train["datetime_beginning_ept"].to_array()
#Evaluate model
data = lc.loc[(lc.index >= '2016-10-01 00:00:00') & (lc.index < '2016-11-01 00:00:00')]
p = data.plot(figsize=(20,8))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09113924050632911, 'fraction_numerical': 0.01870604781997187, 'mean_word_length': 5.844080846968239, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 27, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '194961', 'license': 'Apache-2.0', 'n_tokens_mistral': 2353, 'n_tokens_neox': 2200, 'n_words': 806, 'path': 'pydatanyc.ipynb', 'repo_name': 'isabella232/pydatanyc2018', 'size': '7111'}