<jupyter_start><jupyter_text>Simulate aircraft trajectory<jupyter_code>dt = .01 # 1/100 second
one_over_dt = int(1./dt)
Ts = 100
T = int(Ts/dt) # 10000
x = np.zeros((T, 5)) # [p, V, \psi]
x[0, :] = np.array([100, 100, 100, 10, .1]) # initial state

V = 100 # aircraft airspeed
u = np.zeros((T, 3)) # [\gamma, \phi, dV]
u[20*one_over_dt:500*one_over_dt, 0] = .02
u[10*one_over_dt:12*one_over_dt, 1] = .5
u[50*one_over_dt:52*one_over_dt, 1] = -.3
u[130*one_over_dt:152*one_over_dt, 1] = -.1

for i in range(T-1):
    p = x[i, :3]
    V = x[i, 3]
    psi = x[i, 4]
    gamma, phi, d_V = u[i, :]
    
    d_p = V * np.array([
        np.sin(psi)*np.cos(gamma),
        np.cos(psi)*np.cos(gamma),
        np.sin(gamma)
    ])
    d_psi = (g*np.tan(phi)/V)
    x[i+1, :3] = p + dt * d_p
    x[i+1, 3] = V + dt * d_V
    x[i+1, 4] = psi + dt * d_psi
    
p = x[:, :3]
p = p[::one_over_dt, :]
noise = np.random.multivariate_normal([0, 0, 0], 50*np.eye(3), size=Ts)
p_wiggle = p + noise
# longitudinal view

plt.plot(p[:,0],p[:,1])
plt.scatter(p_wiggle[:,0],p_wiggle[:,1], s=4, alpha=.5)
plt.show()
# vertical profile

plt.plot(p[:,2])
plt.scatter(np.arange(p_wiggle.shape[0]),p_wiggle[:,2], s=4, alpha=.5)
plt.show()<jupyter_output><empty_output><jupyter_text>Smooth with various values of $\lambda_1$ and $\lambda_2$<jupyter_code>N = Ts+20
a = np.zeros(N)
a[:-20:2] = 1
Phat = np.zeros((N, 3))
Phat[:-20:2, :] = p_wiggle[::2, :]

lambda_1s = (1e-6,1e4,1e10)
lambda_2s = (1e-6,1e4,1e7)

fig, axs = plt.subplots(nrows=len(lambda_1s), ncols=len(lambda_2s), sharex=True, sharey=True)
fig.subplots_adjust(left=.08, bottom=.05, right=.99, top=.97)
for i, lambda_1 in enumerate(lambda_1s):
    for j, lambda_2 in enumerate(lambda_2s): 
        P = smooth(a, Phat, N, lambda_1, lambda_2)
        axs[i, j].scatter(p_wiggle[::2, 0], p_wiggle[::2, 1], s=3, alpha=.6, c='gray')
        axs[i, j].plot(P[:, 0], P[:, 1], linewidth=1, c='black', alpha=.8)
        axs[i, j].xaxis.set_major_formatter(plt.NullFormatter())
        axs[i, j].yaxis.set_major_formatter(plt.NullFormatter())
        axs[i,j].set_yticks([])
        axs[i,j].set_xticks([])
z = 40
axs[0, 0].set_ylabel('$1 \\times 10^{-6}$')
axs[1, 0].set_ylabel('$1 \\times 10^{4}$')
axs[2, 0].set_ylabel('$1 \\times 10^{10}$')
axs[-1, 0].set_xlabel('$1 \\times 10^{-6}$')
axs[-1, 1].set_xlabel('$1 \\times 10^{4}$')
axs[-1, 2].set_xlabel('$1 \\times 10^{7}$')
fig.set_size_inches(6, 6)
plt.show()<jupyter_output><empty_output><jupyter_text>Test out-of-sample-validation<jupyter_code>def train_test_split(data, fraction_val):
    data_shuffled = shuffle(data)
    n_val = int(data.shape[0]*fraction_val)

    return data_shuffled[n_val:], data_shuffled[:n_val]

def convert_dataset(data, L):
    Phat = np.zeros((L, 3))
    a = np.zeros(L)
    for i in range(data.shape[0]):
        t, p = int(data[i, 0]), data[i, 1:]
        if t < L:
            Phat[t, :] += p
            a[t] += 1

    return Phat, a
def val_loss(a, P, Phat):
    return np.linalg.norm(np.diag(a).dot(P-Phat), 2)

def smooth_cross_validation(traj, fraction_val=.4, L=225):
    N = max(L, int(np.max(traj[:, 0])))

    best_loss, best_lambda_1, best_lambda_2 = float("inf"), 1., 1.
    train, test = train_test_split(traj, fraction_val)
    Phat_train, a_train = convert_dataset(train, N)
    Phat_test, a_test = convert_dataset(test, N)
    for log_lambda_1 in [-6, 4, 10]:
        for log_lambda_2 in [-6, 4, 10]:
            loss = 0.0
            P = smooth(a_train, Phat_train, N, 10**log_lambda_1, 10**log_lambda_2)
            loss = val_loss(a_test, P, Phat_test) / 50.
            if loss < best_loss:
                best_loss = loss
                best_lambda_1, best_lambda_2 = 10**log_lambda_1, 10**log_lambda_2
            print ("\\num{%.3f} & " % loss, end="")
        print ("\\\\")
    print ("best:", best_lambda_1, best_lambda_2)
    Phat, a = convert_dataset(traj, N)
    return smooth(a, Phat, N, best_lambda_1, best_lambda_2)
traj = np.c_[np.arange(p_wiggle.shape[0])*2, p_wiggle]
smooth_cross_validation(traj, L=int(np.max(traj[0, :])+1))
pass<jupyter_output>\num{1.675} & \num{1.198} & \num{2.724} & \\
\num{1.138} & \num{1.143} & \num{2.724} & \\
\num{4.999} & \num{4.999} & \num{4.999} & \\
best: 10000 1e-06
{'fraction_non_alphanumeric': 0.2057923239934071, 'fraction_numerical': 0.0661643513068048, 'mean_word_length': 3.411214953271028, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '537094', 'license': 'MIT', 'n_tokens_mistral': 2090, 'n_tokens_neox': 1913, 'n_words': 357, 'path': 'TrajectoryReconstruction.ipynb', 'repo_name': 'sisl/terminal-airspace-models', 'size': '4247'}