Tutorial 8 - Using A Custom Model Fit
=====================================

SAILS provides implementations of several algorithms for fitting autoregressive
models but it is straightforward to create a custom class which implements a
new model fit or uses one from another package.

This tutorial will outline how to create a custom model fit class using the
Vector Autoregression class from ``statsmodels.tsa``. We start by importing SAILS
and creating a simulated time series to model.
 
.. code-block:: python

   import sails
   
   # Create a simulated signal
   sample_rate = 100
   siggen = sails.Baccala2001_fig2()
   X = siggen.generate_signal(sample_rate=sample_rate,num_samples=1000)

We then fit a model using Ordinary Least Squared as implemented in SAILS.

.. code-block:: python

   # Fit an autoregressive model with order 3
   sails_model = sails.OLSLinearModel.fit_model(X,np.arange(4))

Now we will create a new model fit class based on
``sails.modelfit.AbstractLinearModel``. This is a base class which contains a
number of methods and properties to store and compute information on a model.
The ``AbstractLinearModel`` is not usable on its own as the fit_model method is
not implemented. When classes such as ``OLSLinearModel`` are defined in SAILS,
they inherit from ``AbstractLinearModel`` to define the helper functions before a
specific ``fit_model`` method is defined. We can do the same to define a custom
model fit class using an external package. We will first create a new class
which inherits from ``AbstractLinearModel`` and then define a ``fit_model`` method
which computes the model fit and stores the outputs in a standard form.

Here is our custom model fit class, each section is described in the comments in the code.

.. code-block:: python

   # Define a new class inheriting from the SAILS base model
   class TSALinearModel( sails.AbstractLinearModel ):
   
       # Define a fit_model method using the python @classmethod decorator
       @classmethod
       def fit_model( cls, data, delay_vect):
   
           # Some sanity checking of the input matrix We make sure that the input
           # data is in 2d format or 3d format with a single epoch.
           # statsmodels.tsa doesn't currently support fitting multitrial data
           if data.ndim == 3:
               if data.shape[2] > 1:
                   raise ValueError('This class is only implemented for single-trial data')
               # Take first trial if we have 3d data
               data = data[:,:,0]
   
           # Create object - classmethods act as object constructors. cls points
           # to TSALinearModel and ret is a specific instance of TSALinearModel
           # though it is currently empty.
           ret = cls()
   
           # The rest of this method will populate ret with information based on
           # our model fit
   
           # Import the model fit function
           from statsmodels.tsa.api import VAR
   
           # Initialise and fit model - we use a simple VAR with default options.
           # Note that we return the model and results to ret.tsa_model.  This
           # means that the statsmodels.tsa.api.VAR object will be stored and
           # returned with ret. Later we can use this to access the statsmodels
           # model and results though the SAILS object.
           ret.tsa_model = VAR(data.T) # SAILS assumes channels first, TSA samples first
           model_order = len(delay_vect) - 1 # delay_vect includes a leading zero
           ret.tsa_results = ret.tsa_model.fit(model_order)
   
           # The method must assign the following values for SAILS metrics to work properly
           ret.maxorder = model_order
           ret.delay_vect = np.arange(model_order)
           ret.parameters = np.concatenate((-np.eye(data.shape[0])[:,:,None],
                                             ret.tsa_results.coefs.transpose((1,2,0))), axis=2)
           ret.data_cov = sails.find_cov(data.T,data.T)
           ret.resid_cov = sails.find_cov(ret.tsa_results.resid.T,ret.tsa_results.resid.T)
   
           # Return fitted model within an instance of a TSALinearModel
           return ret

It is crucial that the ``fit_model`` class returns an instance of our the
overall class. This instance must contain the following information. Other
functions in SAILS assume that these are stored in a fitted model class with
specific formats and names.

 - ``maxorder``: the model order of the fit
 - ``delay_vect``: the vector of delays used in the model fit
 - ``parameters``: the fitted autoregressive parameters of shape `[num_channels x num_channels x model_order]` with a leading identity
 - ``data_cov``: the covariance matrix of the fitted  data
 - ``resid_cov``: the covariance matrix of the residuls of the fit

Other data can be added in as well (we store ``tsa_model`` and ``tsa_results``
in the example here) but these five must be defined within the returned class.

We can now fit a model using our new class

.. code-block:: python

   tsa_model = TSALinearModel.fit_model(X,np.arange(4))

Finally, we compute connectivity metrics from each model fit and plot a comparison

.. code-block:: python

   freq_vect = np.linspace(0,sample_rate/2)

   sails_metrics = sails.FourierMvarMetrics.initialise(sails_model,sample_rate,freq_vect)
   tsa_metrics = sails.FourierMvarMetrics.initialise(tsa_model,sample_rate,freq_vect)

   PDC = np.concatenate( (sails_metrics.partial_directed_coherence,tsa_metrics.partial_directed_coherence),axis=3)

   sails.plotting.plot_vector(PDC,freq_vect,line_labels=['SAILS','TSA'],diag=True,x_label='Frequency (Hz'))

We see that the partial directed coherence from the two models is nearly identical.

.. image:: tutorial8_1.png