Reconstruction of planar shapes from functional principal components.

[aa89113]

Hello,
I am trying to use SRVF framework and fs.fdacurve with planar curves.
I was able to compute the functional principal components, using the following code:

obj = fs.fdacurve(beta,N=M)
obj.karcher_mean(rotation=False)
obj.srvf_align(rotation=False)
obj.karcher_cov()
obj.shape_pca(no=5)

However, I am not sure how can i reconstruct the original planar functions using functional principal scores stored in obj.coef.
The goal would be to reconstruct the original shapes and after compute a reconstruction error metric
Can you please help?

[jdtuck]

You need project back on the basis stored in obj.U, that is the singular vectors. The coefficients are the scores.

[aa89113]

@jdtuck many thanks for your reply.
Performing

((obj.coef.T) @ (obj.U.T))  # n_curve_examples x n_points_per_curve

will project the scores obj.coef on the eigenfunctions. However, this seems to be in the q (?) domain.
I am interested on reconstructing the original dataset from principal components. This would be comparable to obj.beta, but reconstructed from functional scores.
My end goal would be to compare the reconstructed data w/ original dataset.

[jdtuck]

Correct then transform them back to curves using q_to_curve

[aa89113]

@jdtuck many thanks again for your help.

I did the following, but the output does not resemble the original curve, at all ...

reconst_curves_q = ((obj.coef.T) @ (obj.U.T))  # n_curve_examples x n_points_per_curve
reconst_curve_q_first =reconst_curves_q[0].reshape(2, -1)  # 2 x n__points_per_curve
reconst_curve_first = cf.q_to_curve(reconst_curve_q_first)

My example is very similar to:
https://fdasrsf-python.readthedocs.io/en/latest/curve_example.html

Can you please help me?

[jdtuck]

Two things the PCA is done on a tangent space to the Hilbert Sphere, so you need to project back down to the sphere, then convert back to curves (q_to_curve). Have you read the corresponding papers? See the function sample_shapes for an example on how to do this.

fdasrsf_python/fdasrsf/curve_stats.py

Line 360 in 0c5b248

def sample_shapes(self, no=3, numSamp=10):

[aa89113]

I tried to adapt your sample_shapes function to reconstruct original curves from functional scores/coefficients.
See below how it looks ...

The main change was the following line:
v = v + obj.coef[m, i] * np.sqrt(s[m]) * np.vstack((U[0:T, m], U[T : 2 * T, m]))

Instead of multiplying which eigenvector by a random number, I am using the functional scores.
However, the produced output still does not look great.
Can you please tell me what I am doing wrong?

def reconstruct_from_coef(obj, no=3, ):
    n, T = obj.q_mean.shape
    modes = ["O", "C"]
    mode = [i for i, x in enumerate(modes) if x == obj.mode]
    if len(mode) == 0:
        mode = 0
    else:
        mode = mode[0]

    U, s, V = np.linalg.svd(obj.C)

    if mode == 0:
        N = 2
    else:
        N = 10

    epsilon = 1.0 / (N - 1)

    numSamp = obj.coef.shape[1]
    print(f'{numSamp=}')
    samples = np.empty(numSamp, dtype=object)
    for i in range(0, numSamp):
        v = np.zeros((2, T))
        for m in range(0, no):
            # v = v + np.random.randn() * np.sqrt(s[m]) * np.vstack((U[0:T, m], U[T : 2 * T, m]))
            v = v + obj.coef[m, i] * np.sqrt(s[m]) * np.vstack((U[0:T, m], U[T : 2 * T, m]))


        q1 = obj.q_mean
        for j in range(0, N - 1):
            normv = np.sqrt(cf.innerprod_q2(v, v))

            if normv < 1e-4:
                q2 = obj.q_mean
            else:
                q2 = np.cos(epsilon * normv) * q1 + np.sin(epsilon * normv) * v / normv
                if mode == 1:
                    q2 = cf.project_curve(q2)

            # Parallel translate tangent vector
            basis2 = cf.find_basis_normal(q2)
            v = cf.parallel_translate(v, q1, q2, basis2, mode)

            q1 = q2

        samples[i] = cf.q_to_curve(q2)

    obj.samples = samples
    return

[jdtuck]

First,

it should be this

v = v + obj.coef[m, i] * np.vstack((U[0:T, m], U[T : 2 * T, m]))

Also you do not need to parallel translate, really only should have to do this

q1 = obj.q_mean
normv = np.sqrt(cf.innerprod_q2(v, v))
q2 = np.cos(epsilon * normv) * q1 + np.sin(epsilon * normv) * v / normv
if mode == 1:
    q2 = cf.project_curve(q2)

samples[i] = cf.q_to_curve(q2)

[jdtuck]

Also please provide an example, other than it doesn't look right

[aa89113]

This has helped a lot.
I think I only have issues on the scale and translation.
The reconstructed output seems to be much smaller and in a different position vs original data.
For some reason, the scaling seems to be fixed by passing the scale into q_to_curve and then multiply samples[i] by obj.len_q[i]:

samples[i] = cf.q_to_curve(q2, scale=obj.len_q[i]) * obj.len_q[i]

Can you please confirm whether this is correct?

I guess the translation difference is expected ...

[jdtuck]

Correct. The curves are standardized to unit length and centered at the origin before analysis. The centers are stored in cent and the scales are stored in len.

[jdtuck]

You should just pass obj.len into the q_to_curve