[HELP Request] How to use multivariate softDTW barycenter?

I have a problem similar to the one asked in 87, i.e my data set looks as follows, where each row belongs to a sequence and each row is a point in time:

According to the documentation (https://tslearn.readthedocs.io/en/stable/variablelength.html?highlight=variable%20lenght), I can do something like this:

# Each row in the multi-dimension array represents a time-series
# Each array in a row represents a feature, the number of arrays in a row are times the features were encountered
X = to_time_series_dataset([
    [[1,1,1], [1,1,1], [1,1,1]],   # time-series of length 3
    [[1,1,1], [1,1,1]],   # time-series of length 2
    [[4,4,4], [4,4,4]],   # # time-series of length 2
])
y = np.array([0, 0, 1])

X.shape
>> (3, 3, 3)

X
>> array([[[ 0.,  0.,  0.],
        [ 1.,  1.,  1.],
        [ 2.,  2.,  2.]],

       [[ 2.,  2.,  2.],
        [ 3.,  3.,  3.],
        [nan, nan, nan]],

       [[ 4.,  4.,  4.],
        [ 5.,  5.,  5.],
        [nan, nan, nan]]])

The time-series conversion seems good but throws an error when I try to run the following:

initial_barycenter = ts_zeros(sz=5)
bar = softdtw_barycenter(X, init=initial_barycenter)

...
...
...
~/anaconda3/envs/pytorch-build/lib/python3.8/site-packages/sklearn/metrics/pairwise.py in euclidean_distances(X, Y, Y_norm_squared, squared, X_norm_squared)
    260     paired_distances : distances betweens pairs of elements of X and Y.
    261     """
--> 262     X, Y = check_pairwise_arrays(X, Y)
    263 
    264     # If norms are passed as float32, they are unused. If arrays are passed as

~/anaconda3/envs/pytorch-build/lib/python3.8/site-packages/sklearn/metrics/pairwise.py in check_pairwise_arrays(X, Y, precomputed, dtype, accept_sparse, force_all_finite, copy)
    151                              (X.shape[0], X.shape[1], Y.shape[0]))
    152     elif X.shape[1] != Y.shape[1]:
--> 153         raise ValueError("Incompatible dimension for X and Y matrices: "
    154                          "X.shape[1] == %d while Y.shape[1] == %d" % (
    155                              X.shape[1], Y.shape[1]))

ValueError: Incompatible dimension for X and Y matrices: X.shape[1] == 1 while Y.shape[1] == 3

What am I doing wrong?

Start with the variablelength documentation example and the softdtw_barycenter call shown in the issue. Compare the shape produced by to_time_series_dataset with the shape of the ts_zeros initialization, then verify the corrected example against the reported dimension error. Done means the multivariate usage and required initialization shape are clearly documented.