Feature extraction of timestamp based on RBF + polar coordinates of a clock

rbf_timestamp_features.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

""""""

import matplotlib
matplotlib.use("TkAgg")

from datetime import datetime as dt
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
from matplotlib import cm, mlab
from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots()
#fig.set_tight_layout(True)


# ------------------------------
# --- Coordinate calculation ---
# ------------------------------

# Conversion is mostly based on the following post:
# http://www.cs.cornell.edu/courses/cs1132/2011sp/module1/module1part8/example1.html


class TimeFeaturesXY(object):
    def __init__(self):
        # 2D coord (x,y) for each of the 4 time attributes: [day, hour, minute, am/pm]
        self.dimensions = 8

    @staticmethod
    def polar2xy(r, theta):
        """
        (x,y) are the Cartesian coordinates of polar coordinates (r,theta).
        r is the radial coordinate and theta is the angle, also called phase.
        theta is in degrees.
        """
        rads = theta * np.pi / 180.0
        x = r * np.cos(rads)
        y = r * np.sin(rads)
        return x, y

    @staticmethod
    def date_polar2xy(d=0.0, h=0.0, m=0.0, s=0.0):

        # Extract angle of AM/PM condition
        pm = 2*(h / 24.0 + m / 1440.0 + s / 86400.0)
        pm_angle = 90.0 - pm * (360 / 2.0)
        # Convert to (x, y) coordinate
        x_pm, y_pm = TimeFeaturesXY.polar2xy(r=1.0, theta=pm_angle)

        # Extract angle of day given in range [0, 6]
        d = d + h / 24.0 + m / 1440.0 + s / 86400.0
        d_angle = 90.0 - d * (360 / 7.0)
        # Convert to (x, y) coordinate
        x_d, y_d = TimeFeaturesXY.polar2xy(r=1.0, theta=d_angle)

        # Extract angle of hour
        h = h + m / 60.0 + s / 3600.0
        h_angle = 90.0 - h * (360 / 12.0)
        # Convert to (x, y) coordinate
        x_h, y_h = TimeFeaturesXY.polar2xy(r=1.0, theta=h_angle)

        # Extract angle of minute
        m = m + s / 60.0
        m_angle = 90.0 - m * (360 / 60.0)
        # Convert to (x, y) coordinate
        x_m, y_m = TimeFeaturesXY.polar2xy(r=1.0, theta=m_angle)

        return (x_d, y_d), (x_h, y_h), (x_m, y_m), (x_pm, y_pm)

    def transform(self, date):

        # Extract time attributes
        d, h, m, s = date.weekday(), date.hour, date.minute, date.second

        # Extract cartesian coordinates
        (x_d, y_d), (x_h, y_h), (x_m, y_m), (x_pm, y_pm) = self.date_polar2xy(d, h, m, s)

        return np.asarray([x_d, y_d, x_h, y_h, x_m, y_m, x_pm, y_pm])

    @staticmethod
    def n_equidistant_points(n_d=3, n_h=4, n_m=4, n_pm=2):

        # Day
        delta_d = 7 / float(n_d)
        points_d = [TimeFeaturesXY.date_polar2xy(d=i*delta_d, h=0, m=0, s=0)[0]
                    for i in range(n_d)]

        # Hour
        delta_h = 12 / float(n_h)
        points_h = [TimeFeaturesXY.date_polar2xy(d=0, h=i*delta_h, m=0, s=0)[1]
                    for i in range(n_h)]

        # Minute
        delta_m = 60 / float(n_m)
        points_m = [TimeFeaturesXY.date_polar2xy(d=0, h=0, m=i*delta_m, s=0)[2]
                    for i in range(n_m)]

        # AM/PM (based on the hour)
        delta_pm = 12 / float(n_pm)
        points_pm = [TimeFeaturesXY.date_polar2xy(d=0, h=i*delta_pm, m=0, s=0)[3]
                     for i in range(n_m)]

        return points_d, points_h, points_m, points_pm


# --- Conversion of time attributes ---


def string_to_datetime(s, date_format='%Y-%m-%d %H:%M:%S'):
    return dt.strptime(s, date_format)


def extract_time_features(t):
    date = string_to_datetime(t)
    return date.weekday(), date.hour, date.minute, date.second


# --------------------------
# --- Plotting functions ---
# --------------------------

day_mapping = dict(enumerate(["Monday", "Tuesday", "Wednesday", "Thursday",
                              "Friday", "Saturday", "Sunday"]))


def plot_clock(date='2017-04-15 15:30:30',
               radius_pm=1.6, radius_d=1.3, radius_h=1.0, radius_m=0.7,
               color_pm='k', color_d='b', color_h='r', color_m='g',
               show=False):
    plt.gcf().clear()
    #fig.clear()

    #fig = plt.figure(1)
    plt.axis([-2.5, 2.5, -2.5, 2.5])
    #ax = fig.add_subplot(1, 1, 1)

    # Circles for each time attribute
    Cp = plt.Circle((0, 0), radius=radius_pm, linewidth=4, color=color_pm, fill=False)
    Cd = plt.Circle((0, 0), radius=radius_d, linewidth=4, color=color_d, fill=False)
    Ch = plt.Circle((0, 0), radius=radius_h, linewidth=4, color=color_h, fill=False)
    Cm = plt.Circle((0, 0), radius=radius_m, linewidth=4, color=color_m, fill=False)

    ax.add_patch(Cp)
    ax.add_patch(Cd)
    ax.add_patch(Ch)
    ax.add_patch(Cm)

    # Extract cartesian coordinates
    d, h, m, s = extract_time_features(date)
    (x_d, y_d), (x_h, y_h), (x_m, y_m), (x_pm, y_pm) = TimeFeaturesXY.date_polar2xy(d=d, h=h, m=m, s=s)
    dict_features = {"Day": {"X": x_d, "Y": y_d},
                     "Hour": {"X": x_h, "Y": y_h},
                     "Minute": {"X": x_m, "Y": y_m},
                     "AM/PM": {"X": x_pm, "Y": y_pm}}
    print "Features: %s" % str(dict_features)

    # Adjust to a given radio for the arrows
    x_h, y_h = x_h*radius_h, y_h*radius_h
    x_m, y_m = x_m*radius_m, y_m*radius_m
    x_d, y_d = x_d*radius_d, y_d*radius_d
    x_pm, y_pm = x_pm * radius_pm, y_pm * radius_pm

    # Generate the title with the given time
    #plt.title('Clock - %s %ih %im (%s)' %
    #          (day_mapping[d], h, m, "PM" if pm else "AM"))
    plt.title('Clock - %s' % date)

    # Now plot!
    plt.arrow(0, 0, x_h*0.9, y_h*0.9, head_width=0.03, head_length=0.08,
              fc=color_h, ec=color_h, linewidth=3.0)
    plt.arrow(0, 0, x_m*0.9, y_m*0.9, head_width=0.05, head_length=0.08,
              fc=color_m, ec=color_m, linewidth=3.0)
    plt.arrow(0, 0, x_d*0.9, y_d*0.9, head_width=0.05, head_length=0.1,
              fc=color_d, ec=color_d, linewidth=3.0)
    plt.arrow(0, 0, x_pm * 0.9, y_pm * 0.9, head_width=0.05, head_length=0.1,
              fc=color_pm, ec=color_pm, linewidth=3.0)

    # This is just a trick to obtain legends for the arrows
    plt.plot([0, x_h*0.9], [0, y_h*0.9], color_h,
             linewidth=2.0, label="Hour")
    plt.plot([0, x_m * 0.9], [0, y_m * 0.9], color_m,
             linewidth=2.0, label="Minute")
    plt.plot([0, x_d * 0.9], [0, y_d * 0.9], color_d,
             linewidth=2.0, label="Day")
    plt.plot([0, x_pm * 0.9], [0, y_pm * 0.9], color_pm,
             linewidth=2.0, label="AM/PM")

    plt.legend(loc='upper right', shadow=True)

    if show:
        plt.show()

    return fig



# TODO: more sigmax for 0 and 30 min and more sigmay for 15 and 45 min

def plot_gaussian(mux=0.0, muy=0.0, sigmax=0.1, sigmay=0.1,
                  show=False, cmap=None):
    if cmap is None:
        cmap = cm.autumn

    delta = 0.025
    x_i = np.arange(-2.0, 2.0, delta)
    y_i = np.arange(-2.0, 2.0, delta)
    X, Y = np.meshgrid(x_i, y_i)
    Z = mlab.bivariate_normal(X, Y,
                              sigmax=sigmax, sigmay=sigmay,
                              mux=mux, muy=muy)

    # Create a simple contour plot with labels using default colors.
    CS = plt.contour(X, Y, Z, cmap=cmap)

    if show:
        plt.show()


def generate_range_dates(start='2017-11-02 00:00:00', end='2017-11-09 00:00:00',
                         freq='1S', date_format='%Y-%m-%d %H:%M:%S'):
    start_date, end_date = dt.strptime(start, date_format), dt.strptime(end, date_format)
    yield map(lambda d: str(d), pd.date_range(start_date, end_date, freq=freq).tolist())


class GaussianXY(object):
    def __init__(self, muy, mux, sigmax=0.1, sigmay=0.1, radius=1.0, cmap=None, label='Gaussian'):
        if cmap is None:
            cmap = cm.autumn

        self.mux = mux
        self.muy = muy
        self.sigmax = sigmax
        self.sigmay = sigmay
        self.radius = radius
        self.cmap = cmap
        self.label = label

    def get_distance(self, x, y):
        return np.linalg.norm((self.mux - x, self.muy - y))

        #return np.exp(-self.sigmax*np.linalg.norm((self.mux - x, self.muy -y)))


def create_gaussians(n_d=3, n_h=4, n_m=4, n_pm=2, radius_pm=1.6, radius_d=1.3, radius_h=1.0, radius_m=0.7):
    points_d, points_h, points_m, points_pm = TimeFeaturesXY.n_equidistant_points(n_d=n_d, n_h=n_h, n_m=n_m, n_pm=n_pm)

    gaussians = []

    for i, (mux_d, muy_d) in enumerate(points_d):
        gaussians.append(GaussianXY(mux=mux_d, muy=muy_d, radius=radius_d, cmap=cm.winter, label='GDay%i' % i))

    for i, (mux_h, muy_h) in enumerate(points_h):
        gaussians.append(GaussianXY(mux=mux_h, muy=muy_h, radius=radius_h, cmap=cm.autumn, label='GHour%i' % i))

    for i, (mux_m, muy_m) in enumerate(points_m):
        gaussians.append(GaussianXY(mux=mux_m, muy=muy_m, radius=radius_m, cmap=cm.summer, label='GMin%i' % i))

    for i, (mux_pm, muy_pm) in enumerate(points_pm):
        gaussians.append(GaussianXY(mux=mux_pm, muy=muy_pm, radius=radius_pm, cmap=cm.gist_gray, label='GAm/Pm%i' % i))

    return gaussians


def plot_clock_gaussians(date='2017-04-15 15:30:30', gaussians=None, show=False):
    plot_clock(date, show=False)

    if gaussians is None:
        gaussians = create_gaussians()

    for g in gaussians:
        plot_gaussian(mux=g.mux*g.radius, muy=g.muy*g.radius, show=False, cmap=g.cmap)

    if show:
        plt.show()


def plot_week_timeline(delta_sec=120., show=False):
    import time

    t = time.time()
    n_x = int(24*60*60 / delta_sec)

    x = [dt.fromtimestamp(t + i*delta_sec) for i in range(n_x)]

    x_day, y_day, x_hour, y_hour, x_min, y_min, x_ampm, y_ampm = [], [], [], [], [], [], [], []

    for date in x:
        d, h, m, s = extract_time_features(date.strftime('%Y-%m-%d %H:%M:%S'))
        (x_d, y_d), (x_h, y_h), (x_m, y_m), (x_pm, y_pm) = TimeFeaturesXY.date_polar2xy(d=d, h=h, m=m, s=s)
        x_day.append(x_d)
        y_day.append(y_d)
        x_hour.append(x_h)
        y_hour.append(y_h)
        x_min.append(x_m)
        y_min.append(y_m)
        x_ampm.append(x_pm)
        y_ampm.append(y_pm)

    plt.plot(x, x_day, label='DayX')
    plt.plot(x, y_day, label='DayY')
    plt.plot(x, x_hour, label='HourX')
    plt.plot(x, y_hour, label='HourY')
    plt.plot(x, x_min, label='MinX')
    plt.plot(x, y_min, label='MinY')
    plt.plot(x, x_ampm, label='AM/PM_X')
    plt.plot(x, y_ampm, label='AM/PM_Y')

    plt.legend(loc='upper right', shadow=True)

    if show:
        plt.show()

def print_distance_gaussians(delta_sec=120., show=False):
    import time

    t = time.time()
    n_x = int(24*60*60 / delta_sec)

    x = [dt.fromtimestamp(t + i*delta_sec) for i in range(n_x)]

    gaussians = create_gaussians()

    x_day, y_day, x_hour, y_hour, x_min, y_min, x_ampm, y_ampm = [], [], [], [], [], [], [], []

    for date in x:
        d, h, m, s = extract_time_features(date.strftime('%Y-%m-%d %H:%M:%S'))
        (x_d, y_d), (x_h, y_h), (x_m, y_m), (x_pm, y_pm) = TimeFeaturesXY.date_polar2xy(d=d, h=h, m=m, s=s)

        str_g = ''
        for g in gaussians[3:7]:
            str_g += '%s: %.4f ' % (g.label, g.get_distance(x_h, y_h))

        print str_g

        x_day.append(x_d)
        y_day.append(y_d)
        x_hour.append(x_h)
        y_hour.append(y_h)
        x_min.append(x_m)
        y_min.append(y_m)
        x_ampm.append(x_pm)
        y_ampm.append(y_pm)

    plt.plot(x, x_day, label='DayX')
    plt.plot(x, y_day, label='DayY')
    plt.plot(x, x_hour, label='HourX')
    plt.plot(x, y_hour, label='HourY')
    plt.plot(x, x_min, label='MinX')
    plt.plot(x, y_min, label='MinY')
    plt.plot(x, x_ampm, label='AM/PM_X')
    plt.plot(x, y_ampm, label='AM/PM_Y')

    plt.legend(loc='upper right', shadow=True)

    if show:
        plt.show()


def plot_week_clock(delta_min=15, show=False):
    # TODO: make gif (https://eli.thegreenplace.net/2016/drawing-animated-gifs-with-matplotlib/)

    freq = str(delta_min)+'T'

    date_list = generate_range_dates(freq=freq)

    gaussians = create_gaussians()

    plot_fun = lambda d: plot_clock_gaussians(d, gaussians=gaussians)

    anim = FuncAnimation(fig, plot_fun, frames=date_list.next(), interval=5)

    if show:
        plt.show()

# TODO: list:
# - alpha en las gaussianas, que varie dependiendo de la ponderacion en base a la distancia

if __name__ == '__main__':
    # plot_week_timeline(show=True)
    #plot_week_clock(show=True)
    print_distance_gaussians()