BaseContinuousSystemModel

classdef BaseContinuousSystemModel < BaseSystemModel
%BASECONTINUOUSSYSTEMMODEL  Abstract base class for a continuous time system
%   BaseContinuousSystemModel is a base class that encapsulates a
%   noisy system with continuous time evolution.  The class is
%   abstract, and cannot be directly instantiated.
%
%   It models the following stochastic differential equation:
%     dX = f(X,t) dt + A dW
%   with f a function handle, and A is a matrix. This equation is
%   often written as
%     dx/dt = f(x,t) + w    with w = white noise with covariance A * A'
%
%   BaseContinuousSystemModel methods:
%      simulate     - simulate a noisy measurement.
%      observerData - given a Gaussian prediction for the state and an actual 
%                     measurement, returns covariance data based on
%                     observation model.  

    properties
        f;  % f(x,t)
        A;
        J;  % Jacobian of f with respect to x, if known
        nSteps = 1;  % number of steps in integration between times
    end

    methods
        function sysModel = BaseContinuousSystemModel(f, A, varargin)
            if nargin == 0
                return;
            end
            
            if isa(f, 'BaseContinuousSystemModel')
                assert(nargin == 1, 'Too many arguments');
                
                sysModel.f = f.f;
                sysModel.A = f.A;
                sysModel.J = f.J;
                sysModel.nSteps = f.nSteps;
            else
                sysModel.f = f;
                sysModel.A = A;

                if ~isempty(varargin)
                    % now parse option/value pairs
                    for k = 1:2:length(varargin)
                        sysModel.(varargin{k}) = varargin{k+1};
                    end                    
                end
            end
        end

        function newState = simulate(this, oldState, T, curTime)
        % Simulate time evolution of the system for T units of time
        % returns a state vector
        % We use the Euler-Maruyama method, which is simple, but
        % (unfortunately) quite inaccurate if the time interval is not
        % sufficiently small.
            switch nargin
              case {0, 1}
                error('Too few arguments');
              case 2
                T = 1;
                u = 0;
                curTime = 0;
              case 3
                u = 0;
                curTime = 0;
              case 4
                curTime = 0;
            end
            
            nSteps = this.nSteps;
            dt = T/nSteps;
            
            A = this.A;
            nNoise = size(A, 2); % dimension of noise vector

            newState = oldState;
            
            for tau = 1:nSteps;
                newState = newState + this.f(newState, curTime) * dt ...
                    + A * sqrt(dt) * randn(nNoise,1);
                curTime = curTime + dt;
            end;
        end
    end

end