<jupyter_start><jupyter_text>Section 2 $$\begin{aligned}C(X_r, X_s) =& \frac{\sqrt{\frac{1}{n}\sum_{i=1}^n\left[x_r(i) - x_s(i)\right]^2}}{\sqrt{\frac{1}{n}\sum_{i=1}^nx_r(i)^2} + \sqrt{\frac{1}{n}\sum_{i=1}^nx_s(i)^2}} \\=& \frac{\sqrt{\phantom{\frac{1}{n}}\sum_{i=1}^n\left[x_r(i) - x_s(i)\right]^2}}{\sqrt{\phantom{\frac{1}{n}}\sum_{i=1}^nx_r(i)^2} + \sqrt{\phantom{\frac{1}{n}}\sum_{i=1}^nx_s(i)^2}}\end{aligned}$$  Benchmark problem Two parallel plots have the same trend. But, they end up showing maximally inconsistent by TIC.<jupyter_code>consistent_fig()

n = 100
x_r = np.ones(n)
x_s_hflip = -x_r
tic = swm.tic(x_r, x_s_hflip)
plt.plot(x_r, label='$X_r$', color='red', lw=1)
plt.plot(x_s_hflip, label='$X_s$', color='blue', lw=1)
plt.title(f"TIC={tic}")
sns.despine()
plt.legend();<jupyter_output><empty_output><jupyter_text>This applies generally -- TIC can't adapt to translations. And, worse than that, the variance adjustments in the denominator make it so translations of the same magnitude but different signs have different scores.<jupyter_code>consistent_fig()


np.random.seed(100)
n = 11
c = 1
x_r = np.sqrt(np.arange(n))
plt.plot(x_r, label='$x_r$', color='red', lw=1)

for marker, translation in zip("^xov", [4, 2, -2, -4]):
    x_s = x_r + translation
    tic = swm.tic(x_r, x_s)
    plt.plot(x_s, label='$x_{s_1}$ ' + f'(TIC={tic:0.3f})', color='blue', marker=marker, lw=1)

sns.despine()
plt.legend();<jupyter_output><empty_output><jupyter_text>Trend Problem Even if the trends are maximally inconsistent, TIC can't pick it up, showing a (relatively) better score,<jupyter_code>consistent_fig()


np.random.seed(100)
n = 100
c = 1
x_r = np.arange(3)+1
x_s = -(x_r - 4)

tic = swm.tic(x_r, x_s)

plt.plot(x_r, label='$x_{r}$', marker='o', color='red', lw=1)
plt.plot(x_s, label='$x_{s}$', marker='o', color='blue', lw=1)
plt.ylim(0, 4)
sns.despine()
plt.title(f"TIC = {tic:0.3f}")
plt.legend();
consistent_fig()


np.random.seed(100)
n = 100
c = 1
x_r_prime = np.diff(x_r)
x_s_prime = np.diff(x_s)

tic = swm.tic(x_r_prime, x_s_prime)

plt.plot(x_r_prime, label='$x_{r}$', marker='o', color='red', lw=1)
plt.plot(x_s_prime, label='$x_{s}$', marker='o', color='blue', lw=1)
sns.despine()
plt.title(f"TIC = {tic:0.3f}")
plt.legend();<jupyter_output><empty_output><jupyter_text>Simulation Model Validation Method The metric they derived (2) $$C(X_r, X_s) =\left\{	\begin{array}{ll}		\frac{\sqrt{\sum_{i=1}^n \left[x_r(i) - x_s(i)\right]^2}}{\sqrt{\sum_{i=1}^n x_r(i)^2}}  & \mbox{if } \sqrt{\sum_{i=1}^n x_r(i)^2} \ne 0 \\		\sqrt{\sum_{i=1}^n x_s(i)^2} & \mbox{if } \sqrt{\sum_{i=1}^n x_r(i)^2}  = 0	\end{array}\right.$$ confuses me a little. If the sum of squares for $x_r$ is zero as in the lower equation's condition, then it's just a line at 0. I'm not sure why you would want to compare the simulation output to what is essentially no action. (Maybe some robotics domain knowledge that I'm missing?) Also, the $\xi$ is essentially acting as a temperature parameter in,$$D(X_r, X_s) = e^{-\xi C(X_r, X_s)}$$and,$$T(X_r, X_s) = e^{-\xi\lambda \delta(X'_r, X'_s)}$$to squash the new metric to $(0, 1)$. They say it is,> set by domain expert according to the purpose of simulation applicationso it's not as useful as TIC in that now there is a hyper-parameter to select. It's also a bit annoying in that their squashing function flipped the interpretability: Now close to 1 is a good match; before it mean maximally different.  Minor error There are multiple errors in their equation (4) and (5). They have, $$\begin{aligned}X'_r &= \langle x'_r(2), x'_r(2), \dots, x'_r(n-1)\rangle = \left\langle \frac{x_r(2) - x_r(1)}{\Delta t} + \frac{x_r(3) - x_r(2)}{\Delta t} + \cdots + \frac{x_r(n) - x_r(n-1)}{\Delta t} \right\rangle \\X'_s &= \langle x'_s(2), x'_s(2) \dots, x'_s(n-1)\rangle = \left\langle \frac{x_s(2) - x_s(1)}{\Delta t} + \frac{x_s(3) - x_s(2)}{\Delta t} + \cdots + \frac{x_s(n) - x_s(n-1)}{\Delta t} \right\rangle \\\end{aligned}$$it should be,$$\begin{aligned}X'_r &= \langle x'_r(2), x'_r(3), \dots, x'_r(n)\rangle = \left\langle \frac{x_r(2) - x_r(1)}{\Delta t} + \frac{x_r(3) - x_r(2)}{\Delta t} + \cdots + \frac{x_r(n) - x_r(n-1)}{\Delta t} \right\rangle \\X'_s &= \langle x'_s(2), x'_s(3), \dots, x'_s(n)\rangle = \left\langle \frac{x_s(2) - x_s(1)}{\Delta t} + \frac{x_s(3) - x_s(2)}{\Delta t} + \cdots + \frac{x_s(n) - x_s(n-1)}{\Delta t} \right\rangle \\\end{aligned}$$ So basically, take the rate of change for each observation in the reference and simulation, use that for the "trends" series, too.   Integration of position conformity and trend conformity This section assumes a set of simulation outputs, $$O_s = \{X_{s1}, X_{s2}, \dots, X_{sp} \}$$which they use to construct a matrix,$$A = \begin{bmatrix}a_{11} & a_{12} \\a_{21} & a_{22} \\\vdots & \vdots \\a_{p1} & a_{p2} \\\end{bmatrix}$$with, $$\begin{aligned}a_{i1} &= D(X_r, X_{si}) \\a_{i2} &= T(X_r, X_{si})\end{aligned}$$ So the first column is the improved and squashed TIC on position and the second is on trend. They then translate A into S by,$$s_{ij} = \frac{a_{ij}}{\bar{a}_{\mathbf{\bullet} j}} $$where $$a_{\mathbf{\bullet} j} = \sum_{k=1}^p \frac{a_{kj}}{p}$$So, both columns are now in column-mean units. $$S = \begin{bmatrix}s_{11} & s_{12} \\s_{21} & s_{22} \\\vdots & \vdots \\s_{p1} & s_{p2} \\\end{bmatrix}$$ **Here's where the presentation gets a little weird for me.** They define,$$y = k_1s_1 + k_2s_2$$which you can use to select the best simulation. It's a linear combination of the squashed/improved TIC metric on the position and trend, respectively. They further define, $$k_1^2 + k_2^2 = 1$$so that $k$, the weights, represent a unit-vector. That would be a requirement of PCA. (I'm not sure why it would be mean centered (via division) versus mean-shifted so it was zero).But, here's the thing -- this whole section just derives PCA. See Replication section below.  Application They give the expected A matrix via Table 2 and 3.<jupyter_code>expected_D_col = np.array([0.729, 0.733, 0.663, 0.813, 0.840, 0.724, 0.69, 0.555, 0.8])
expected_T_col = np.array([1, 0.05, 1, 0.517, 0.730, 0.652, 0.368, .185, .192])
expected_A = np.vstack([expected_D_col, expected_T_col]).T
expected_A.T<jupyter_output><empty_output><jupyter_text>But they don't give the data used for Fig 3. I think this is *close* to correct.<jupyter_code>t = np.arange(21)

x_sqrt = np.sqrt(t)
X_r = x_sqrt
two_spikes = x_sqrt.copy() + 0.7
two_spikes[6] = 5.7
two_spikes[11] = 5.7

simulations = []
simulations.append(('a', x_sqrt+1))
simulations.append(('b', two_spikes))
simulations.append(('c', x_sqrt+1.2))
simulations.append(('d', t * 5/20))
simulations.append(('e', t**(1/3)+1))
simulations.append(('f', t**(1/4)+2))
simulations.append(('g', 2 + np.ones_like(x_sqrt)))
simulations.append(('e', 2 + -t**(1/8)+1))
simulations.append(('f', x_sqrt + np.sin(2*np.pi * t/20 * 2.25)))

fig = plt.figure(figsize=(15, 15))

sns.set_style('whitegrid')

shared_ax = None
i = 1
for l, X_s in simulations:
    if i == 1:
        ax = fig.add_subplot(3, 3, i)
        shared_ax = ax
        plt.xlim((0, 20))
        plt.ylim((0, 6))
    else:
        ax = fig.add_subplot(3, 3, i, sharex=ax, sharey=ax)
    
    
    plt.plot(t, X_r, marker='o', color='red', lw=1)
    plt.plot(t, X_s, marker='d', color='blue', lw=1)
    
    plt.title(l)
    sns.despine()
    i += 1<jupyter_output><empty_output><jupyter_text>The also don't give their $\xi$. Using the deduced data set, it seems like $\xi=1$,<jupyter_code>consistent_fig()

ξ_space = np.linspace(0.01, 2, 2001)
scores = []
for ξ in ξ_space:
    A = swm.make_matrix_A(x_sqrt, [x_s for _, x_s in simulations], ξ)
    ϵ = np.abs(expected_A - A).sum()
    scores.append(ϵ)
    
plt.plot(ξ_space, scores)
plt.xlabel(r'$\xi$')
plt.ylabel('Absolute error compared to paper matrix A');<jupyter_output><empty_output><jupyter_text>Replication The "paydirt" of all of this is supposed to be equation (18) which they have as,$$y = 0.0223s_1 + 0.9998s_2$$ What is this? It the projection of matrix $S$ onto the first principle component of S. I'm not sure why they felt compelled to derive the whole thing. The only difference is how they center it. The most obvious way to show this? Using the data from above,<jupyter_code>ξ = 1
A = swm.make_matrix_A(x_sqrt, [x_s for _, x_s in simulations], ξ)
A
S = swm.make_matrix_S(x_sqrt, [x_s for _, x_s in simulations], ξ)
S<jupyter_output><empty_output><jupyter_text>Here is the PCA according to sklearn,<jupyter_code>from sklearn.decomposition import PCA

# Create a PCA estimator.
# (We only care about the 1st principle component)
pca = PCA(1)

# Fit it on S
pca.fit(S)

# Recover the principle components.
# These are the k_1 and k_2 in their paper.
k_ = swm.corrected_components(pca.components_[0])
k_<jupyter_output><empty_output><jupyter_text>That's it.  Or, "by hand,"<jupyter_code># Get covariance matrix as per (12)
V = np.cov(S.T)

# Get eigen values and vectors
# That's what they get to in Equation (16) 
λ, k = np.linalg.eig(V)

# Get the biggest eigenvalue
i = np.argmax(λ)

# Select the associated eigenvector
k = k[:, 1]
k<jupyter_output><empty_output><jupyter_text>This is different only in sign. But, it's an eigenvector, so it's a unit vector,<jupyter_code>(k**2).sum()<jupyter_output><empty_output><jupyter_text>and as such, it's only the vector that matters. If both of the weights are negative, you can negatate them, it doesn't change the mathematical object in this context.<jupyter_code>consistent_fig()

plt.scatter([0], [0], marker=10)
plt.quiver(*k, scale=3, lw=1, color='red')
plt.quiver(*k_, scale=3, lw=1, color='blue');
consistent_fig()
ax = plt.gca()

ax.scatter(S[:, 0], S[:, 1]);

for i, label in enumerate(c for c, _ in simulations):
    ax.annotate(label, S[i, :] + np.array([0.01, 0.01]))
y = S @ k_
y
consistent_fig()
ax = plt.gca()

y = S @ k_

ax.scatter(y, y);

for i, label in enumerate(c for c, _ in simulations):
    ax.annotate(label, np.array([y[i], y[i]]) + np.array([0.0, -0.1]))
df = pd.DataFrame({'y': y}, index=list('abcdefghi'))
df['result'] = df['y'].rank(ascending=False)
df.T<jupyter_output><empty_output><jupyter_text>Or, using the replicated package,<jupyter_code>y = swm.compute_y(x_sqrt, [x_s for _, x_s in simulations], ξ)
y
y.argmax() # Best model by row index in S
rankdata(-y)  # Best is highest, so negate<jupyter_output><empty_output><jupyter_text>Follow Up Questions I mentioned I wasn't sure how the $\xi$ parameter would affect things. As a temperature parameter in isolation, *it still preserves relative position*.<jupyter_code>fig = plt.figure(figsize=(15, 15))

sns.set_style('whitegrid')


shared_ax = None
i = 1
for l, X_s in simulations:
    plt.xlim((0, 2))
    plt.ylim((0, 1))
    
    x = np.linspace(0.1, 2.0)
    y = []
    for ξ in x:
        y.append(swm.squashed_tic_improved(x_sqrt, X_s, ξ))
    
    plt.plot(x, y, marker='o', lw=1, label=l)
    
    sns.despine()
    i += 1
plt.legend();<jupyter_output><empty_output><jupyter_text>**and, when done with the PCA, relative rank is preserved as well**. So ξ doesn't matter all that much as long as you select it consistently.<jupyter_code>sims = [sim for _, sim in simulations]
    
rows = []
for ξ in x:
    rows.append(rankdata(-swm.compute_y(x_sqrt, sims, ξ)))
fig = plt.figure(figsize=(15, 15))
sns.set_style('whitegrid')
plt.plot(rows);<jupyter_output><empty_output><jupyter_text>Phone Call Follow Up I refactored a bit so now there are four new functions, - `equation_8`- `position_metric`- `trend_metric`- `mixed_metric` `equation_8` implements  $$\delta(X_r, X_s) =\left\{	\begin{array}{ll}		\frac{\sqrt{\sum_{i=1}^n \left(\left(x_r(i+1)-X_r(i)\right) - \left(x_s(i+1)-X_s(i)\right)\right)^2}}{\sqrt{\sum_{i=1}^n \left(x_s(i+1)-X_s(i)\right)^2}}  & \mbox{if } \sqrt{\sum_{i=1}^n x'_r(i)^2} \ne 0 \\		\sqrt{\sum_{i=1}^n \left(x_s(i+1)-X_s(i)\right)^2} & \mbox{if } \sqrt{\sum_{i=1}^n x'_r(i)^2}  = 0	\end{array}\right.$$and makes use of `np.diff(x)`,<jupyter_code>np.diff([1, 2, 4, 3])<jupyter_output><empty_output><jupyter_text>Computed on the following graph (e)<jupyter_code>consistent_fig()

t = np.arange(21)
x_sqrt = np.sqrt(t)
X_s = t**(1/3)+1
plt.plot(t, X_r, marker='o', color='red', lw=1)
plt.plot(t, X_s, marker='d', color='blue', lw=1);
swm.equation_8(X_r, X_s)
d = swm.position_metric(X_r, X_s, ξ=1)
d
t = swm.trend_metric(X_r, X_s, ξ=1)
t<jupyter_output><empty_output><jupyter_text>I picked (e) to illustrate that, without the PCA method, this no longer is in [0, 1), it's now in [0, 2).<jupyter_code>d + t<jupyter_output><empty_output><jupyter_text>Mixed metric mixes them as, $$w_{trend}D(X_r, X_s) + (1-w_{trend})T(X_r, X_s)$$which keeps it in [0, 1)<jupyter_code>swm.mixed_metric(0, X_r, X_s, ξ=1) == t
swm.mixed_metric(1, X_r, X_s, ξ=1) == d
swm.mixed_metric(0.5, X_r, X_s, ξ=1) == (d + t) / 2<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18387498056289847, 'fraction_numerical': 0.0299331363707044, 'mean_word_length': 4.69915817456801, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 74, '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': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '224851', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5538, 'n_tokens_neox': 5222, 'n_words': 1356, 'path': 'thiel_app/clean_replication.ipynb', 'repo_name': 'zlite/PX4_flight_review', 'size': '12862'}