Event demo

Events demo¶

This guide shows how to use the low level events generator package to find event sets and event interval sets of basic functions. It also show the usage of autodif in finding event root sensitivity to model parameters.

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from godot.core import num, tempo, events

We define the following function:

$$ f(t) = t * \cos(\omega * t) $$

In [2]:

freq = num.TwoPi
ref_epoch = tempo.Epoch("0.0 TDB")
def f(e):
    t = (e - ref_epoch) / tempo.SecondsInDay
    return t * np.cos(t * freq)

To plot this function we generate some data

In [3]:

ran = tempo.EpochRange(tempo.Epoch("0.0 TDB"), tempo.Epoch("4.0 TDB"))
grid = ran.createGrid(1e-2 * tempo.SecondsInDay)
t = [e.mjd() for e in grid]
x = [f(e) for e in grid]

We use matplotlib to plot the function.

In [4]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
plt.grid()
plt.legend()
plt.show()

We can find all the roots of this function using the event finding module

In [5]:

tol = 1e-9
eps = 1e-3
event_grid = ran.createGrid(0.3 * tempo.SecondsInDay)
ev = events.generateEventSet(f, eps, event_grid, tol)
event_t = [e.mjd() for e in event_grid]
event_x = [f(e) for e in event_grid]
roots = [entry.value().mjd() for entry in ev]

Once we have the roots, we can add them to the plot to visualize them.

In [6]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
plt.plot(roots, [0] * len(roots), 'o', label="roots")
plt.grid()
plt.legend()
plt.show()

We can do the same to generate event interval sets where the function is positive.

In [7]:

f_above_zero = events.generateEventIntervalSet(f, eps, event_grid, tol)
f_ranges = [tempo.EpochRange(entry.start().value(), entry.end().value()) for entry in f_above_zero]

And visualize the intervals on the function.

In [8]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
for ran in f_ranges:
    grid_ev = ran.createGrid(1e-2 * tempo.SecondsInDay)
    t_ev = [e.mjd() for e in grid_ev]
    x_ev = [f(e) for e in grid_ev]
    plt.plot([ran.start().mjd(), ran.end().mjd()], [0, 0], '-ok', linewidth=3)
    plt.plot(t_ev, x_ev, '-k', linewidth=3)
plt.grid()
plt.legend()
plt.show()

To demonstrate how to combine event intervals, we define the corresponding sine function counterpart:

$$ f(t) = t * \sin(\omega * t) $$

In [9]:

def g(e):
    t = (e - ref_epoch) / tempo.SecondsInDay
    return t * np.sin(t * freq)
y = [g(e) for e in grid]

And visualize them together.

In [10]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
plt.plot(t, y, label="g")
plt.grid()
plt.legend()
plt.show()

We can do the same event computation on the second function.

In [11]:

g_above_zero = events.generateEventIntervalSet(g, eps, event_grid, tol)
g_ranges = [tempo.EpochRange(entry.start().value(), entry.end().value()) for entry in g_above_zero]

And visualize the event intervals where each function is positive.

In [12]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
plt.plot(t, y, label="g")
for ran in f_ranges:
    grid_ev = ran.createGrid(1e-2 * tempo.SecondsInDay)
    t_ev = [e.mjd() for e in grid_ev]
    x_ev = [f(e) for e in grid_ev]
    plt.plot([ran.start().mjd(), ran.end().mjd()], [0, 0], '-ok', linewidth=3)
    plt.plot(t_ev, x_ev, '-k', linewidth=3)
for ran in g_ranges:
    grid_ev = ran.createGrid(1e-2 * tempo.SecondsInDay)
    t_ev = [e.mjd() for e in grid_ev]
    y_ev = [g(e) for e in grid_ev]
    plt.plot([ran.start().mjd(), ran.end().mjd()], [0, 0], '-or', linewidth=3)
    plt.plot(t_ev, y_ev, '-r', linewidth=3)
plt.grid()
plt.legend()
plt.show()

We can easily combine event intervals to form unions and intersections between them.

In [13]:

f_or_g_above_zero = events.any([f_above_zero, g_above_zero])
f_or_g_ranges = [tempo.EpochRange(entry.start().value(), entry.end().value()) for entry in f_or_g_above_zero]

We can visualize the union of the two event intervals sets, where either function is positive.

In [14]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
plt.plot(t, y, label="g")
for ran in f_or_g_ranges:
    grid_ev = ran.createGrid(1e-2 * tempo.SecondsInDay)
    t_ev = [e.mjd() for e in grid_ev]
    x_or_y_ev = [max(f(e), g(e)) for e in grid_ev]
    plt.plot([ran.start().mjd(), ran.end().mjd()], [0, 0], '-ok', linewidth=3)
    plt.plot(t_ev, x_or_y_ev, '-k', linewidth=3)
plt.grid()
plt.legend()
plt.show()

Similarly, we can also compute the intersection of the two.

In [15]:

f_and_g_above_zero = events.all([f_above_zero, g_above_zero])
f_and_g_ranges = [tempo.EpochRange(entry.start().value(), entry.end().value()) for entry in f_and_g_above_zero]

And visualize when BOTH functions are positive.

In [16]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
plt.plot(t, y, label="g")
for ran in f_and_g_ranges:
    grid_ev = ran.createGrid(1e-2 * tempo.SecondsInDay)
    t_ev = [e.mjd() for e in grid_ev]
    x_and_y_ev = [min(f(e), g(e)) for e in grid_ev]
    plt.plot([ran.start().mjd(), ran.end().mjd()], [0, 0], '-ok', linewidth=3)
    plt.plot(t_ev, x_and_y_ev, '-k', linewidth=3)
plt.grid()
plt.legend()
plt.show()

Event finding using partial derivatives¶

In the second part of the guide we show how to work with functions with dependent model parameters and how to compute the sensitivites of the roots.

In [17]:

from godot.core import autodif as ad

We re-define our first function, but now using autodif variables. We can also define this function in a more generic way, as later we will se.

In [18]:

freq = ad.Scalar(num.TwoPi, "freq")
def f2(e):
    time_from_ref_epoch = (e - ref_epoch) / tempo.SecondsInDay
    return time_from_ref_epoch * ad.cos(time_from_ref_epoch * freq)

We repeat the same excersise as before to plot it.

In [19]:

ran = tempo.EpochRange(tempo.Epoch("0.0 TDB"), tempo.Epoch("4.0 TDB"))
grid = ran.createGrid(1e-2 * tempo.SecondsInDay)
t = [e.mjd() for e in grid]
x = [f2(e).value() for e in grid]
x_old = [f2(e).value() for e in grid]

In [20]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
plt.grid()
plt.legend()
plt.show()

Now we do root finding for event intervals, but this time using autodif. We will compute the length of the third interval where the function is positive. We can see that apart from finding the width of the interval, we also obtain the sensitivity to the frequency model parameter.

In [21]:

eis = events.generateXEventIntervalSet(f2, eps, event_grid, tol)
idx = 2
ran = tempo.XEpochRange(eis[idx].start().value(), eis[idx].end().value())
w = ran.width() / tempo.SecondsInDay
print(w)

         Value |          freq
  5.000000e-01 | -7.957747e-02

In [22]:

print(w.gradient()["freq"])

[-0.07957747]

We visualize the third interval of interest.

In [23]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
grid_ev = ran.createGrid(1e-2 * tempo.SecondsInDay)
t_ev = [e.mjd() for e in grid_ev]
x_ev = [f2(e).value() for e in grid_ev]
plt.plot([ran.start().mjd(), ran.end().mjd()], [0, 0], '-ok', linewidth=4)
plt.plot(t_ev, x_ev, '-k', linewidth=4)
plt.grid()
plt.legend()
plt.show()

What we aim to do as a final excersive is to use a simple gradient descent algorithm to find the frequency model parameter that changes the third interval to take 0.6 units.

In [24]:

import godot.core.autodif.bridge as br

We create functions to compute the values and gradients of the event interval duration.

In [25]:

idx = 2

def f(e, freq):
    time_from_ref_epoch = (e - ref_epoch) / tempo.SecondsInDay
    return time_from_ref_epoch * br.cos(time_from_ref_epoch * freq)

def func(freq, target):
    eis = events.generateEventIntervalSet(lambda e: f(e, freq), eps, event_grid, tol)
    ran = tempo.EpochRange(eis[idx].start().value(), eis[idx].end().value())
    return ran.width() - target

def grad(freq):
    freq_x = ad.Scalar(freq, "freq")
    eis = events.generateXEventIntervalSet(lambda e: f(e, freq_x), eps, event_grid, tol)
    ran = tempo.XEpochRange(eis[idx].start().value(), eis[idx].end().value())
    return ran.width().gradient()["freq"]

We write a simple gradient descent algorithm for one dimension.

In [26]:

def gradient_descent(x0, func, grad, tol):
    x = x0
    it = 0
    max_it = 10
    while np.abs(func(x)) > tol and it < max_it:
        dx = -func(x) / grad(x)
        x = x + dx
        print(func(x) / tempo.SecondsInDay)
    return x

We run the gradient descent algorithm and target 0.6 days duration. As output, we have the optimal frequency to target the duration.

In [27]:

freq0 = num.TwoPi
target = 0.6 * tempo.SecondsInDay
tol = 1e-6
freq_opt = gradient_descent(freq0, lambda freq: func(freq, target), grad, tol)
print(freq_opt)

0.02499999999999924
0.0009615384615392065
1.5360039317244998e-06
3.09800843415023e-12
[5.23598776]

freq0 = num.TwoPi target = 0.6 * tempo.SecondsInDay tol = 1e-9 freq_opt = gradient_descent(freq0, lambda freq: func(freq, target), grad, tol) print(freq_opt)

In [28]:

freq = freq_opt
x = [f(e, freq) for e in grid]
eis = events.generateXEventIntervalSet(f2, eps, event_grid, tol)
ran = tempo.XEpochRange(eis[idx].start().value(), eis[idx].end().value())
w = ran.width() / tempo.SecondsInDay
print(w)

         Value
  6.000000e-01

freq = freq_opt x = [f(e, freq) for e in grid] eis = events.generateXEventIntervalSet(f2, eps, event_grid, tol) ran = tempo.XEpochRange(eis[idx].start().value(), eis[idx].end().value()) w = ran.width() / tempo.SecondsInDay print(w)

In [29]:

plt.figure(figsize=(10, 5))
plt.xlabel("Time (day)")
plt.ylabel("Function value (1)")
plt.plot(t, x, label="f")
plt.plot(t, x_old, '--b', label="f_old", linewidth=0.2)
grid_ev = ran.createGrid(1e-2 * tempo.SecondsInDay)
t_ev = [e.mjd() for e in grid_ev]
x_ev = [f2(e).value() for e in grid_ev]
plt.plot([ran.start().mjd(), ran.end().mjd()], [0, 0], '-ok', linewidth=4)
plt.plot(t_ev, x_ev, '-k', linewidth=4)
plt.grid()
plt.legend()
plt.show()

In [ ]: