Example scripts of shift registration methods and interpolation of FDataGrid

Readme.md

Examples of fda.registration

Shift Registration

• shift_registration1.py: Periodic registration
• shift_registration2.py: Non-periodic registration
• landmark_shift.py: Align using landmarks

Interpolation

• curve_interpolation.py: Interpolation of curves R -> R^n
• surface_interpolation.py: Surfaces R^2 -> R^n
• manifold_interpolation.py: Manifolds R^m -> R^n

Other

• landmark_registration.py
• shift_basis.py

curve_interpolation.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

from fda.grid import FDataGrid
import mpl_toolkits.mplot3d

# Data parameters
nsamples = 3  # Number of samples
nscat = 15
nfine = 120  # Number of points per sample

phase_sd = .1  # Standard deviation of phase variation
amp_sd = .1  # Standard deviation of amplitude variation
error_sd = .0  # Standard deviation of gaussian noise
period = 1  # Period of the sine used to generate the data
xlim= (0,1)

# Interpolation parameters
order = 3
s = 0
monotone = False

# Dimension of data of last example
n = 5
nscat_n = n*5 # Number of samples last example


def noise_sin(t, nsamples, period=1, phase_sd=0, amp_sd=0, error_sd=0):

    phase_variation = np.outer(np.random.normal(0, phase_sd,  nsamples),
                               np.ones(len(t)))

    error = np.random.normal(0, error_sd, (nsamples, len(t)))

    amp = np.diag(np.random.normal(1, amp_sd, nsamples))

    return amp @ np.sin(2*np.pi*t/period + phase_variation) + error


if __name__ == '__main__':

    # Seaborn style
    plt.style.use('seaborn')

    # Fixing random state for reproducibility
    np.random.seed(98765)

    # Matrix with times where each sample will be evaluated
    t = np.linspace(xlim[0], xlim[1], nscat)
    T = np.linspace(xlim[0], xlim[1], nfine)

    # First dimension
    data_1 = noise_sin(t, nsamples, period, phase_sd, amp_sd, error_sd)

    # Data R^1->R^1
    fdgrid = FDataGrid(data_1, t, xlim, interpolation_order=order,
                       monotone=monotone, smoothness_parameter=s)

    plt.figure()
    fdgrid.scatter()
    plt.xlabel(r"$t$")
    plt.ylabel(r"$\sin(t)$")
    plt.plot(T,fdgrid(T).T)

    # Data R^1->R^2
    data_2 = noise_sin(t, nsamples, .5*period, phase_sd, amp_sd, error_sd)

    data2_1 = np.dstack((data_1, data_2))
    fdgrid_2_1 = FDataGrid(data2_1, t, xlim, interpolation_order=order,
                           monotone=monotone,  smoothness_parameter=s)

    values = fdgrid_2_1(T)

    plt.figure()
    plt.ylabel(r"$\sin(\frac{1}{2}t)$")
    plt.xlabel(r"$\sin(t)$")
    for i in range(nsamples):
        plt.scatter(data_1[i], data_2[i])
        plt.plot(values[i,:,0], values[i,:,1])

    # Data R^1->R^3
    data_3 = noise_sin(t, nsamples, (1./3)*period, phase_sd, amp_sd, error_sd)

    data3_1 = np.dstack((data_1, data_2, data_3))


    fdgrid_3_1 = FDataGrid(data3_1, t, xlim, interpolation_order=order,
                           monotone=monotone,  smoothness_parameter=s)

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.set_xlabel(r"$\sin(t)$")
    ax.set_ylabel(r"$\sin(\frac{1}{2}t)$")
    ax.set_zlabel(r"$\sin(\frac{1}{3}t)$")

    values = fdgrid_3_1(T)


    for i in range(nsamples):
        ax.plot(values[i,:,0], values[i,:,1], values[i,:,2])
        ax.scatter(data_1[i], data_2[i], data_3[i])

    # Functions R^1->R^n
    f, axarr = plt.subplots(n, n, sharex=True, sharey=True)

    t = np.linspace(xlim[0], xlim[1], nscat_n)
    data = np.dstack([noise_sin(t, nsamples, (1./j)*period, phase_sd, amp_sd, error_sd) for j in range(1,n+1)])

    fdgrid = FDataGrid(data, t, xlim, interpolation_order=order,
                           monotone=monotone,  smoothness_parameter=s)

    values = fdgrid(T)


    for i in range(n):
        for j in range(n):
            #for k in range(nsamples):
            #    axarr[i][j].scatter(data[k,:,i], data[k,:,j])
            if i==0:
                axarr[i][j].title.set_text(r"$\sin(\frac{2 \pi t}{" +str(j+1) +"})$")

            if j==0:
                axarr[i][j].set_ylabel(r"$\sin(\frac{2 \pi t}{" +str(i+1) +"})$")
            axarr[i][j].plot(values[:,:,i].T, values[:,:,j].T, linewidth=1.1)
            axarr[i][j].set_yticklabels([])
            axarr[i][j].set_xticklabels([])

    plt.show()

landmark_registration.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

from fda.grid import FDataGrid
import mpl_toolkits.mplot3d

from fdasrsf.utility_functions import invertGamma

# Data parameters
nsamples = 10  # Number of samples
nscat = 15
nfine = 150  # Number of points per sample

order = 3

sd = 1.3
loc = 4


if __name__ == '__main__':

    plt.style.use('seaborn')

    a = np.empty((nsamples,2))
    a[:,0] = np.random.normal(loc=-loc,scale=sd, size=nsamples)
    a[:,1] = np.random.normal(loc=loc,scale=sd, size=nsamples)


    t = np.linspace(-10,10,nfine)
    data = np.empty((nsamples,nfine))


    for i in range(nsamples):
        data[i] = np.exp(-(a[i,0]-t)**2) + np.exp(-  (a[i,1]-t)**2  )



    x = FDataGrid(data, t, interpolation_order=order)

    x.plot()

    for i in range(nsamples):
        plt.scatter(a[i], [1,1])

    plt.title("Samples $x_i(t)$")

    t1 = np.full((nsamples,1),-10)
    tf = np.full((nsamples,1),10)

    p = np.hstack((t1,a,tf))

    h = FDataGrid(p, (-10,-loc,loc,10), monotone=True, interpolation_order=order)


    plt.figure()

    #h.scatter()

    s = t
    h_v = h(s)
    plt.plot(s, h_v.T)
    h.scatter()
    plt.title("Warping functions $h_i(t)$")


    for i in range(nsamples):
        h_v[i] = 20*invertGamma((h_v[i] + 10)/20) - 10

    plt.figure()
    hI = FDataGrid(h_v, t, interpolation_order=order)
    hI.plot()
    for i in range(nsamples):
        plt.scatter(p[i],hI(p[i])[i])
    plt.title("Inverse warping functions $h_i^{-1}(t)$")

    thI = hI(t)
    x_t = x(t)
    h_t = h(t)


    plt.figure()
    for i in range(nsamples):
        plt.plot(t, x(h(t)).T)
        plt.scatter((-10,-loc,loc,10),(0,1,1,0))

    plt.title("Registered samples $x_i^*(t)=x_i(h_i(t))$")

    plt.figure()
    for i in range(nsamples):
        plt.plot(thI[i], x_t[i])
        plt.scatter((-10,-loc,loc,10),(0,1,1,0))

    plt.title("Plot of $(h_i^{-1}(t), x_i(t))$")







    plt.show()

landmark_shift.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
from scipy import stats

import matplotlib.pyplot as plt

from fda.basis import BSpline
from fda.grid import FDataGrid
from fda.registration import landmark_shift

# Data parameters
nsamples = 12  # Number of samples
nfine = 100  # Number of points per sample
sd = .1 # sd of gaussian curves
shift_sd = .2  # Standard deviation of phase variation
amp_sd = .1  # Standard deviation of amplitude variation
error_sd = .05  # Standard deviation of gaussian noise
xlim = (-1, 1) # Domain range

# Location of the landmark
# Possible values are None, a number or a callable
location = 0.
ext = 'periodic'

# Basis parameters
nknots = 20

# Plot options
width = .8  # Width of the samples curves
samples_color = 'teal'
mean_color = 'black'
c_color = 'maroon'

def noise_gaussian(t, nsamples, sd, shift_sd, amp_sd, error_sd):
    """Noisy Gaussian curve function

    Args:
        t (ndarray): Array of times
        nsamples (float): Number of samples
        period (float): Period of the gaussian function sin(2*pi*t/period)
        shift_sd: Standard deviation of the shift variation normaly distributed
        amp_sd: Standard deviation of the amplitude variation
        error_sd: Standard deviation of the error of the samples

    Returns:
        ndarray with the samples evaluated. Each row is a sample and each
        column is a discrete time and the maximun of the curves
    """

    shift = np.random.normal(0, shift_sd,  nsamples)

    shift_variation = np.outer(shift, np.ones(len(t)))

    error = np.random.normal(0, error_sd, (nsamples, len(t)))

    amp = np.diag(np.random.normal(1, amp_sd, nsamples))

    tsamples = t - shift_variation

    return amp @ stats.norm.pdf(tsamples, 0, sd) + error, shift



if __name__ == '__main__':

    # Matplotlib stylesheet
    plt.style.use('seaborn')

    # Fixing random state for reproducibility
    np.random.seed(987654)

    # Matrix with times where each sample will be evaluated
    t = np.linspace(xlim[0], xlim[1], nfine)

    # Noisy gaussian data, with amplitude variation and gaussian error
    data, shift = noise_gaussian(t, nsamples, sd, shift_sd, amp_sd, error_sd)
    fdgrid = FDataGrid(data, t, xlim, 'Raw data')

    # Real gaussian function
    gaussian = stats.norm.pdf(t, 0, sd)

    # Plot the samples
    plt.figure()
    plt.xlim(xlim)
    l1 = fdgrid.plot(label='samples', c=samples_color, linewidth=width)
    l2 = plt.plot(t, gaussian, label='gaussian', c=c_color, linestyle='dashed')
    l3 = fdgrid.mean().plot(label='mean', c=mean_color)
    plt.legend(handles=[l1[0], l3[0], l2[0]], loc=1)


    # Convert to BSplines
    fd = fdgrid.to_basis(BSpline(xlim, nbasis=nknots+2))
    unregmean = fd.mean()  # Mean of unregistered curves

    # Values at landmarks
    y = fd.evaluate_shifted([0], shift)

    # Plots the smoothed curves
    ax = plt.figure()
    plt.title('Unregistered curves')
    plt.xlim(xlim)
    l1 = fd.plot(label='samples', c=samples_color, linewidth=width)
    l2 = plt.plot(t, gaussian, label='gaussian', c=c_color, linestyle='dashed')
    l3 = unregmean.plot(label='mean', c=mean_color)
    l4 = plt.scatter(shift, y, label='landmarks', c=c_color)
    plt.legend(handles=[l1[0], l3[0], l2[0], l4], loc=1)

    ylim = plt.ylim()

    # Shift registered curves
    regbasis = landmark_shift(fd, shift, location, ext=ext)
    regmean = regbasis.mean()  # Registered mean

    # Values at 0, the new location of landmarks
    yreg = regbasis.evaluate([0])

    # Plots the registered curves
    plt.figure()
    plt.title('Registered curves')
    plt.xlim(xlim)
    plt.ylim(ylim)
    l1 = regbasis.plot(label='samples', c=samples_color, linewidth=width)
    l2 = plt.plot(t, gaussian, label='gaussian', c=c_color, linestyle='dashed')
    l3 = regmean.plot(label='mean', c=mean_color)
    l4 = plt.scatter(nsamples*[0], yreg, label='landmarks', c=c_color)
    plt.legend(handles=[l1[0], l3[0], l2[0], l4], loc=1)


    plt.show()

manifold_interpolation.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

from fda.grid import FDataGrid
import mpl_toolkits.mplot3d

# Data parameters
nsamples = 100  # Number of samples
nscat = 15
nfine = 50  # Number of points per sample

phase_sd = .1  # Standard deviation of phase variation
amp_sd = .1  # Standard deviation of amplitude variation
error_sd = .0  # Standard deviation of gaussian noise
period = 1  # Period of the sine used to generate the data
xlim= (0,1)

# Interpolation parameters
order = 3
s = 0
monotone = False

# Dimension of data of last example
n = 5
nscat_n = n*5 # Number of samples last example


if __name__ == '__main__':

    # Seaborn style
    plt.style.use('seaborn')

    # Fixing random state for reproducibility
    np.random.seed(98765)

    # Matrix with times where each sample will be evaluated
    t = np.linspace(xlim[0], xlim[1], nscat)

    T = np.linspace(xlim[0], xlim[1], nfine)

    X, Y, Z = np.meshgrid(t,t,t)

    R = np.sqrt(X + Y + Z)

    data = np.empty((nsamples, nscat, nscat, nscat,2))

    for i in range(nsamples):
        data[i,...,0] = np.sin(R)
        data[i,...,1] = np.cos(R)


    fdgrid = FDataGrid(data, (t,t,t), interpolation_order=0)

    # Check if interpolated values in nodes are equal
    differences = fdgrid((t,t,t),grid=True) - data

    all_zeros = not np.any(differences)
    if all_zeros:
        print("Coincidence")

    interpolated = fdgrid((T,T,T),grid=True)

    for t in (0,int(nscat/2)):
        for s in (0,int(nscat/2)):
            plt.figure()
            for i in range(nsamples):

                plt.scatter(data[i,t,s,:,0],data[i,t,s,:,1])
                plt.plot(interpolated[i,t,s,:,0],interpolated[i,t,s,:,1])

    plt.show()

shift_basis.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

from fda.basis import *
from fda.grid import FDataGrid

# Data parameters
nsamples = 15  # Number of samples
nfine = 100  # Number of points per sample
phase_sd = .6  # Standard deviation of phase variation
amp_sd = .05  # Standard deviation of amplitude variation
error_sd = .2  # Standard deviation of gaussian noise
period = 1  # Period of the sine used to generate the data
nbasis = 11  # Number of fourier basis elements

# Plot options
width = .8  # Width of the samples curves
samples_color = 'teal'
mean_color = 'black'
curve_color = 'maroon'
ylim = (-1.75, 1.75)
xlim = (0, 1)

iterations = 5  # Number of iterations in the step-by-step figure

if __name__ == '__main__':

    # Seaborn style
    plt.style.use('seaborn')

    # Fixing random state for reproducibility
    np.random.seed(98765)

    # Matrix with times where each sample will be evaluated
    t = np.linspace(xlim[0], xlim[1], nfine)

    # Noisy sine data, with amplitude variation and gaussian error
    data = np.outer(np.ones(nsamples), np.sin(2*np.pi / period * t))
    basis = Fourier(domain_range=xlim, nbasis=nbasis, period=period)
    fd_fourier = FDataGrid(data, t, xlim, 'Raw Data').to_basis(basis)

    # Fourier Basis
    plt.figure()
    fd_fourier.plot()
    plt.title("Fourier basis")

    plt.figure()
    fd_fourier.shift(3).plot() # Probamos una translacion por escalar
    plt.title("Fourier basis scalar shift")

    plt.figure()
    fd_fourier.shift(np.linspace(-0.25,0.25,nsamples)).plot()
    plt.title("Fourier basis default shift")
    plt.xlim(xlim)

    plt.figure()
    fd_fourier.shift(np.linspace(-0.25,0.25,nsamples), ext="slice").plot()
    plt.title("Fourier basis slice")
    plt.xlim(xlim)

    plt.figure()
    fd_fourier.shift(np.linspace(-0.25,0.25,nsamples), ext="const").plot()
    plt.title("Fourier basis const extrapolation")
    plt.xlim(xlim)

    # BSplines Basis
    xlim = (-3,3)
    t = np.linspace(xlim[0], xlim[1], nfine)
    data = np.outer(np.ones(nsamples), np.exp(- t**2))
    basis = BSpline(domain_range=xlim, nbasis = nbasis)
    fd_splines = FDataGrid(data, t, xlim, 'Raw Data').to_basis(basis)


    plt.figure()
    fd_splines.plot()
    plt.title("Bsplines")

    plt.figure()
    fd_splines.shift(5).plot()

    plt.figure()
    plt.title("BSplines Non-periodic extrapolation")
    fd_splines.shift(np.linspace(-0.25,0.25,nsamples)).plot()
    plt.xlim(xlim)

    plt.figure()
    plt.title("BSPlines periodic extrapolation")
    fd_splines.shift(np.linspace(-0.25,0.25,nsamples), ext="periodic").plot()
    plt.xlim(xlim)

    plt.figure()
    plt.title("BSPlines const extrapolation")
    fd_splines.shift(np.linspace(-0.25,0.25,nsamples), ext="const").plot()
    plt.xlim(xlim)

    plt.show()

shift_registration1.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

from fda.basis import FDataBasis, Fourier
from fda.grid import FDataGrid
from fda.registration import shift_registration

# Data parameters
nsamples = 15  # Number of samples
nfine = 100  # Number of points per sample
phase_sd = .6  # Standard deviation of phase variation
amp_sd = .05  # Standard deviation of amplitude variation
error_sd = .2  # Standard deviation of gaussian noise
period = 1  # Period of the sine used to generate the data
nbasis = 11  # Number of fourier basis elements

# Plot options
width = .8  # Width of the samples curves
samples_color = 'teal'
mean_color = 'black'
curve_color = 'maroon'
ylim = (-1.75, 1.75)
xlim = (0, 1)
iterations = 5  # Number of iterations in the step-by-step figure


def noise_sin(t, nsamples, period=1, phase_sd=0, amp_sd=0, error_sd=0):
    """Sine noisy function

    Generates several samples of the function

    ..math::
        x_i(t) = a_i \sin(\frac{2 \pi t}{period} + \delta_i) + \epsilon(t)

    evaluated at the samples points. :math: `a_i, \delta_{i}` and :math:
    `\epsilon(t_{ij})` are normaly distributed.

    Args:
        t (ndarray): Array of times
        nsamples (float): Number of samples
        period (float): Period of the sine function
        phase_sd (float): Standard deviation of the phase variation
        amp_sd (float): Standard deviation of the amplitude variation
        error_sd (float): Standard deviation of the error of the samples

    Returns:
        ndarray: Matrix with the samples evaluated. Each row is a sample and
        each column is a discrete time.
    """

    phase_variation = np.outer(np.random.normal(0, phase_sd,  nsamples),
                               np.ones(len(t)))

    error = np.random.normal(0, error_sd, (nsamples, len(t)))

    amp = np.diag(np.random.normal(1, amp_sd, nsamples))

    return amp @ np.sin(2*np.pi*t/period + phase_variation) + error


if __name__ == '__main__':

    # Seaborn style
    plt.style.use('seaborn')

    # Fixing random state for reproducibility
    np.random.seed(98765)

    # Matrix with times where each sample will be evaluated
    t = np.linspace(xlim[0], xlim[1], nfine)

    # Noisy sine data, with amplitude variation and gaussian error
    data = noise_sin(t, nsamples, period, phase_sd, amp_sd, error_sd)
    fdgrid = FDataGrid(data, t, xlim, 'Raw Data')
    # Real sine function
    sine = np.sin(2*np.pi*t/period)

    # Plot the samples
    plt.figure()
    plt.ylim(ylim)
    plt.xlim(xlim)
    l1 = fdgrid.plot(label='samples', c=samples_color, linewidth=width)
    l2 = plt.plot(t, sine, label='sine', c=curve_color, linestyle='dashed')
    l3 = fdgrid.mean().plot(label='mean', c=mean_color)
    plt.legend(handles=[l1[0], l3[0], l2[0]], loc=1)

    # Curves smoothed with the matrix penalty method
    # Curves smoothed with the matrix penalty method
    fd = fdgrid.to_basis(Fourier(xlim, nbasis, period))
    unregmean = fd.mean()  # Mean of unregistered curves

    # Plots the smoothed curves
    plt.figure()
    plt.title('Unregistered curves')
    plt.ylim(ylim)
    plt.xlim(xlim)
    l1 = fd.plot(label='samples', c=samples_color, linewidth=width)
    l2 = plt.plot(t, sine, label='sine', c=curve_color, linestyle='dashed')
    l3 = unregmean.plot(label='mean', c=mean_color)
    plt.legend(handles=[l1[0], l3[0], l2[0]], loc=1)

    # Shift registered curves
    regbasis = shift_registration(fd)
    regmean = regbasis.mean()  # Registered mean

    # Plots the registered curves
    plt.figure()
    plt.title('Registered curves')
    plt.ylim(ylim)
    plt.xlim(xlim)
    l1 = regbasis.plot(label='samples', c=samples_color,
                       linewidth=width)
    l2 = plt.plot(t, sine, label='sine', c=curve_color, linestyle='dashed')
    l3 = regmean.plot(label='mean', c=mean_color)
    plt.legend(handles=[l1[0], l3[0], l2[0]], loc=1)

    # Plots the process step by step
    f, axarr = plt.subplots(iterations+1, 1, sharex=True, sharey=True)
    axarr[0].title.set_text('Step by step registration')
    plt.xlim(xlim)

    fd.plot(ax=axarr[0], c=samples_color, linewidth=width)
    axarr[0].set_ylabel('Unregistered')

    for i in range(1, iterations+1):
        # tol=0 to realize all the iterations
        regfd = shift_registration(fd, maxiter=i, tol=0.)
        regfd.plot(ax=axarr[i], c=samples_color, linewidth=width)
        axarr[i].set_ylabel('%d iteration%s' % (i, '' if i == 1 else 's'))

    plt.show()

shift_registration2.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
from scipy import stats

import matplotlib.pyplot as plt

from fda.basis import FDataBasis, BSpline
from fda.grid import FDataGrid
from fda.registration import shift_registration

# Data parameters
nsamples = 9  # Number of samples
nfine = 100  # Number of points per sample
sd = .1
shift_sd = .05  # Standard deviation of phase variation
amp_sd = 0 #.1  # Standard deviation of amplitude variation
error_sd = .05  # Standard deviation of gaussian noise
xlim = (-1, 1) # Domain range

# Basis parameters
nbasis = 7  # Number of fourier basis elements
nknots = 20

# Registration parameters
maxiter = 20
tol = 1e-5

# Plot options
width = .8  # Width of the samples curves
samples_color = 'teal'
mean_color = 'black'
curve_color = 'maroon'
ylim = None
iterations = 5  # Number of iterations in the step-by-step figure

def noise_gaussian(t, nsamples, sd, shift_sd, amp_sd, error_sd):
    """Noisy Gaussian curve function

    Args:
        t (ndarray): Array of times
        nsamples (float): Number of samples
        period (float): Period of the gaussian function sin(2*pi*t/period)
        shift_sd: Standard deviation of the shift variation normaly distributed
        amp_sd: Standard deviation of the amplitude variation
        error_sd: Standard deviation of the error of the samples

    Returns:
        ndarray with the samples evaluated. Each row is a sample and each
        column is a discrete time.
    """

    shift_variation = np.outer(np.random.normal(0, shift_sd,  nsamples),
                               np.ones(len(t)))

    error = np.random.normal(0, error_sd, (nsamples, len(t)))

    amp = np.diag(np.random.normal(1, amp_sd, nsamples))

    tsamples = t - shift_variation

    return amp @ stats.norm.pdf(tsamples, 0, sd) + error



if __name__ == '__main__':

    # Matplotlib stylesheet
    plt.style.use('seaborn')

    # Fixing random state for reproducibility
    np.random.seed(1)

    # Matrix with times where each sample will be evaluated
    t = np.linspace(xlim[0], xlim[1], nfine)

    # Noisy gaussian data, with amplitude variation and gaussian error
    data = noise_gaussian(t, nsamples, sd, shift_sd, amp_sd, error_sd)

    # Real gaussian function
    gaussian = stats.norm.pdf(t, 0, sd)

    # Plot the samples
    plt.figure()
    plt.title('Raw data')
    #plt.ylim(ylim)
    plt.xlim(xlim)
    l1 = plt.plot(t, data.T, label='samples', c=samples_color, linewidth=width)
    l2 = plt.plot(t, gaussian, label='gaussian', c=curve_color, linestyle='dashed')
    l3 = plt.plot(t, np.mean(data.T, axis=1), label='mean', c=mean_color)
    plt.legend(handles=[l1[0], l3[0], l2[0]], loc=1)


    knots = np.cos (np.pi * np.arange(nknots + 1) / nknots)
    knots = np.linspace(-1,1,nknots)
    # Curves smoothed with the matrix penalty method
    fd = FDataBasis.from_data(data, t, BSpline(xlim, knots=knots))
    unregmean = fd.mean()  # Mean of unregistered curves

    # Plots the smoothed curves
    plt.figure()
    plt.title('Unregistered curves')
    #plt.ylim(ylim)
    plt.xlim(xlim)
    l1 = fd.plot(label='samples', c=samples_color, linewidth=width)
    l2 = plt.plot(t, gaussian, label='gaussian', c=curve_color, linestyle='dashed')
    l3 = unregmean.plot(label='mean', c=mean_color)
    plt.legend(handles=[l1[0], l3[0], l2[0]], loc=1)

    # Shift registered curves
    regbasis = shift_registration(fd, maxiter=maxiter,tol=tol)
    regmean = regbasis.mean()  # Registered mean


    # Plots the registered curves
    plt.figure()
    plt.title('Registered curves')
    #plt.ylim(ylim)
    plt.xlim(xlim)
    l1 = regbasis.plot(label='samples', c=samples_color,
                       linewidth=width)
    l2 = plt.plot(t, gaussian, label='gaussian', c=curve_color, linestyle='dashed')
    l3 = regmean.plot(label='mean', c=mean_color)
    plt.legend(handles=[l1[0], l3[0], l2[0]], loc=1)

    # Plots the process step by step
    f, axarr = plt.subplots(iterations+1, 1, sharex=True, sharey=True)
    axarr[0].title.set_text('Step by step registration')
    plt.xlim(xlim)

    fd.plot(ax=axarr[0], c=samples_color, linewidth=width)
    axarr[0].set_ylabel('Unregistered')

    for i in range(1, iterations+1):
        # tol=0 to realize all the iterations
        regfd = shift_registration(fd, maxiter=i, tol=0.)
        regfd.plot(ax=axarr[i], c=samples_color, linewidth=width)
        axarr[i].set_ylabel('%d iteration%s' % (i, '' if i == 1 else 's'))

    plt.show()

surface_interpolation.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

from fda.basis import FDataBasis, Fourier
from fda.grid import FDataGrid
from fda.registration import shift_registration, mse_decomposition

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Data parameters
nsamples = 2  # Number of samples
nscat = 10 # Number of points generated per sample
nfine = 75  # Number of points interpolated per sample

xlim = (0,1)

# Interpolation parameters
order = 3
s = 0
derivative = 0

# Dimension of image in the last example
n_image = 4


if __name__ == '__main__':

    # Seaborn style
    plt.style.use('seaborn')

    # Fixing random state for reproducibility
    np.random.seed(98765)

    # Example of function R^2 -> R^1
    data = np.empty((nsamples, nscat, nscat, 1))

    # Points of grid
    t = s = np.linspace(xlim[0], xlim[1], nscat)
    X, Y = np.meshgrid(t, s, indexing='ij')
    # Point of interpolated grid
    T = S = np.linspace(xlim[0], xlim[1], nfine)
    Xn, Yn = np.meshgrid(T, S, indexing='ij')


    # Generates data
    for i in range(nsamples):

        a, b = np.random.normal(loc=1, scale=.2, size=2)
        data[i,:,:,0] = a*np.sin(2*np.pi*Y) + b*np.sin(np.pi*X)

    # Scatter of data
    for _ in range(2):
        fig = plt.figure()
        ax = fig.add_subplot(111, projection='3d')
        for i in range(nsamples):
            ax.scatter(X,Y,data[i,:,:,0])

    # Creation of an FDGrid
    fdgrid = FDataGrid(data,  (t, s), interpolation_order=order, smoothness_parameter=s)


    z = fdgrid((T,S), derivative, grid=True) # Values interpolated

    for i in range(nsamples):
        ax.plot_surface(Xn, Yn, z[i])

    # Example R^2 -> R^n
    data = np.empty((nsamples, nscat, nscat, n_image))

    for j in range(n_image):
        for i in range(nsamples):
            a, b = np.random.normal(loc=1, scale=.2, size=2)
            data[i,:,:,j] = a*np.sin(2*(j+1)*np.pi*Y) + b*np.sin((2*j+2)*np.pi*X)


    fdgrid_2 = FDataGrid(data,  (t, s), interpolation_order=order, smoothness_parameter=s)

    # set up a figure twice as wide as it is tall
    fig = plt.figure(figsize=plt.figaspect(1./n_image))
    plt.title(r"$\mathbb{R}^2 \rightarrow \mathbb{R}^"+str(n_image)+"$")

    Z = fdgrid_2((T,S), derivative, grid=True)

    for i in range(n_image):
        ax = fig.add_subplot(1, n_image, i+1, projection='3d')

        for j in range(nsamples):
            ax.plot_surface(Xn, Yn, Z[j,:,:,i])
            ax.scatter(X, Y, data[j,:,:,i])

    # Alternative plot
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(Z[0,:,:,0], Z[0,:,:,1], Z[0,:,:,2], cmap=cm.coolwarm)
    ax.scatter(data[0,:,:,0], data[0,:,:,1], data[0,:,:,2])

    plt.show()