Source code for pysb.export.matlab

"""
A class for converting a PySB model to a set of ordinary differential
equations for integration in MATLAB.

Note that for use in MATLAB, the name of the ``.m`` file must match the name of
the exported MATLAB class (e.g., ``robertson.m`` for the example below).

For information on how to use the model exporters, see the documentation
for :py:mod:`pysb.export`.

Output for the Robertson example model
======================================

Information on the form and usage of the generated MATLAB class is contained in
the documentation for the MATLAB model, as shown in the following example for
``pysb.examples.robertson``::

    classdef robertson
        % A simple three-species chemical kinetics system known as "Robertson's
        % example", as presented in:
        % 
        % H. H. Robertson, The solution of a set of reaction rate equations, in Numerical
        % Analysis: An Introduction, J. Walsh, ed., Academic Press, 1966, pp. 178-182.
        % 
        % A class implementing the ordinary differential equations
        % for the robertson model.
        %
        % Save as robertson.m.
        %
        % Generated by pysb.export.matlab.MatlabExporter.
        %
        % Properties
        % ----------
        % observables : struct
        %     A struct containing the names of the observables from the
        %     PySB model as field names. Each field in the struct
        %     maps the observable name to a matrix with two rows:
        %     the first row specifies the indices of the species
        %     associated with the observable, and the second row
        %     specifies the coefficients associated with the species.
        %     For any given timecourse of model species resulting from
        %     integration, the timecourse for an observable can be
        %     retrieved using the get_observable method, described
        %     below.
        %
        % parameters : struct
        %     A struct containing the names of the parameters from the
        %     PySB model as field names. The nominal values are set by
        %     the constructor and their values can be overriden
        %     explicitly once an instance has been created.
        %
        % Methods
        % -------
        % robertson.odes(tspan, y0)
        %     The right-hand side function for the ODEs of the model,
        %     for use with MATLAB ODE solvers (see Examples).
        %
        % robertson.get_initial_values()
        %     Returns a vector of initial values for all species,
        %     specified in the order that they occur in the original
        %     PySB model (i.e., in the order found in model.species).
        %     Non-zero initial conditions are specified using the
        %     named parameters included as properties of the instance.
        %     Hence initial conditions other than the defaults can be
        %     used by assigning a value to the named parameter and then
        %     calling this method. The vector returned by the method
        %     is used for integration by passing it to the MATLAB
        %     solver as the y0 argument.
        %
        % robertson.get_observables(y)
        %     Given a matrix of timecourses for all model species
        %     (i.e., resulting from an integration of the model),
        %     get the trajectories corresponding to the observables.
        %     Timecourses are returned as a struct which can be
        %     indexed by observable name.
        %
        % Examples
        % --------
        % Example integration using default initial and parameter
        % values:
        %
        % >> m = robertson();
        % >> tspan = [0 100];
        % >> [t y] = ode15s(@m.odes, tspan, m.get_initial_values());
        %
        % Retrieving the observables:
        %
        % >> y_obs = m.get_observables(y)
        %
        properties
            observables
            parameters
        end

        methods
            function self = robertson()
                % Assign default parameter values
                self.parameters = struct( ...
                    'k1', 0.040000000000000001, ...
                    'k2', 30000000, ...
                    'k3', 10000, ...
                    'A_0', 1, ...
                    'B_0', 0, ...
                    'C_0', 0);

                % Define species indices (first row) and coefficients
                % (second row) of named observables
                self.observables = struct( ...
                    'A_total', [1; 1], ...
                    'B_total', [2; 1], ...
                    'C_total', [3; 1]);
            end

            function initial_values = get_initial_values(self)
                % Return the vector of initial conditions for all
                % species based on the values of the parameters
                % as currently defined in the instance.

                initial_values = zeros(1,3);
                initial_values(1) = self.parameters.A_0; % A()
                initial_values(2) = self.parameters.B_0; % B()
                initial_values(3) = self.parameters.C_0; % C()
            end

            function y = odes(self, tspan, y0)
                % Right hand side function for the ODEs

                % Shorthand for the struct of model parameters
                p = self.parameters;

                % A();
                y(1,1) = -p.k1*y0(1) + p.k3*y0(2)*y0(3);
                % B();
                y(2,1) = p.k1*y0(1) - p.k2*power(y0(2), 2) - p.k3*y0(2)*y0(3);
                % C();
                y(3,1) = p.k2*power(y0(2), 2);
            end

            function y_obs = get_observables(self, y)
                % Retrieve the trajectories for the model observables
                % from a matrix of the trajectories of all model
                % species.

                % Initialize the struct of observable timecourses
                % that we will return
                y_obs = struct();

                % Iterate over the observables;
                observable_names = fieldnames(self.observables);
                for i = 1:numel(observable_names)
                    obs_matrix = self.observables.(observable_names{i});
                    if isempty(obs_matrix)
                        y_obs.(observable_names{i}) = zeros(size(y, 1), 1);
                        continue
                    end
                    species = obs_matrix(1, :);
                    coefficients = obs_matrix(2, :);
                    y_obs.(observable_names{i}) = ...
                                    y(:, species) * coefficients';
                end
            end
        end
    end
"""

import pysb
import pysb.bng
import sympy
import re
from io import StringIO
from pysb.export import Exporter, pad, ExpressionsNotSupported, \
    CompartmentsNotSupported


[docs]class MatlabExporter(Exporter): """A class for returning the ODEs for a given PySB model for use in MATLAB. Inherits from :py:class:`pysb.export.Exporter`, which implements basic functionality for all exporters. """
[docs] def export(self): """Generate a MATLAB class definition containing the ODEs for the PySB model associated with the exporter. Returns ------- string String containing the MATLAB code for an implementation of the model's ODEs. """ if self.model.expressions: raise ExpressionsNotSupported() if self.model.compartments: raise CompartmentsNotSupported() output = StringIO() pysb.bng.generate_equations(self.model) docstring = '' if self.docstring: docstring += self.docstring.replace('\n', '\n % ') # Substitute underscores for any dots in the model name model_name = self.model.name.replace('.', '_') # -- Parameters and Initial conditions ------- # Declare the list of parameters as a struct params_str = 'self.parameters = struct( ...\n'+' '*16 params_str_list = [] for i, p in enumerate(self.model.parameters): # Add parameter to struct along with nominal value cur_p_str = "'%s', %.17g" % (_fix_underscores(p.name), p.value) # Decide whether to continue or terminate the struct declaration: if i == len(self.model.parameters) - 1: cur_p_str += ');' # terminate else: cur_p_str += ', ...' # continue params_str_list.append(cur_p_str) # Format and indent the params struct declaration params_str += ('\n'+' '*16).join(params_str_list) # Fill in an array of the initial conditions based on the named # parameter values initial_values_str = ('initial_values = zeros(1,%d);\n'+' '*12) % \ len(self.model.species) initial_values_str += ('\n'+' '*12).join( ['initial_values(%d) = self.parameters.%s; %% %s' % (i+1, _fix_underscores(ic.value.name), ic.pattern) for i, ic in enumerate(self.model.initials)]) # -- Build observables declaration -- observables_str = 'self.observables = struct( ...\n'+' '*16 observables_str_list = [] for i, obs in enumerate(self.model.observables): # Associate species and coefficient lists with observable names, # changing from zero- to one-based indexing cur_obs_str = "'%s', [%s; %s]" % \ (_fix_underscores(obs.name), ' '.join([str(sp+1) for sp in obs.species]), ' '.join([str(c) for c in obs.coefficients])) # Decide whether to continue or terminate the struct declaration: if i == len(self.model.observables) - 1: cur_obs_str += ');' # terminate else: cur_obs_str += ', ...' # continue observables_str_list.append(cur_obs_str) # Format and indent the observables struct declaration observables_str += ('\n'+' '*16).join(observables_str_list) # -- Build ODEs ------- # Build a stringified list of species species_list = ['%% %s;' % s for i, s in enumerate(self.model.species)] # Build the ODEs as strings from the model.odes array odes_list = ['y(%d,1) = %s;' % (i+1, sympy.ccode(self.model.odes[i])) for i in range(len(self.model.odes))] # Zip the ODEs and species string lists and then flatten them # (results in the interleaving of the two lists) odes_species_list = [item for sublist in zip(species_list, odes_list) for item in sublist] # Flatten to a string and add correct indentation odes_str = ('\n'+' '*12).join(odes_species_list) # Change species names from, e.g., '__s(0)' to 'y0(1)' (note change # from zero-based indexing to 1-based indexing) odes_str = re.sub(r'__s(\d+)', \ lambda m: 'y0(%s)' % (int(m.group(1))+1), odes_str) # Change C code 'pow' function to MATLAB 'power' function odes_str = re.sub(r'pow\(', 'power(', odes_str) # Prepend 'p.' to named parameters and fix any underscores for i, p in enumerate(self.model.parameters): odes_str = re.sub(r'\b(%s)\b' % p.name, 'p.%s' % _fix_underscores(p.name), odes_str) # -- Build final output -- output.write(pad(r""" classdef %(model_name)s %% %(docstring)s %% A class implementing the ordinary differential equations %% for the %(model_name)s model. %% %% Save as %(model_name)s.m. %% %% Generated by pysb.export.matlab.MatlabExporter. %% %% Properties %% ---------- %% observables : struct %% A struct containing the names of the observables from the %% PySB model as field names. Each field in the struct %% maps the observable name to a matrix with two rows: %% the first row specifies the indices of the species %% associated with the observable, and the second row %% specifies the coefficients associated with the species. %% For any given timecourse of model species resulting from %% integration, the timecourse for an observable can be %% retrieved using the get_observable method, described %% below. %% %% parameters : struct %% A struct containing the names of the parameters from the %% PySB model as field names. The nominal values are set by %% the constructor and their values can be overriden %% explicitly once an instance has been created. %% %% Methods %% ------- %% %(model_name)s.odes(tspan, y0) %% The right-hand side function for the ODEs of the model, %% for use with MATLAB ODE solvers (see Examples). %% %% %(model_name)s.get_initial_values() %% Returns a vector of initial values for all species, %% specified in the order that they occur in the original %% PySB model (i.e., in the order found in model.species). %% Non-zero initial conditions are specified using the %% named parameters included as properties of the instance. %% Hence initial conditions other than the defaults can be %% used by assigning a value to the named parameter and then %% calling this method. The vector returned by the method %% is used for integration by passing it to the MATLAB %% solver as the y0 argument. %% %% %(model_name)s.get_observables(y) %% Given a matrix of timecourses for all model species %% (i.e., resulting from an integration of the model), %% get the trajectories corresponding to the observables. %% Timecourses are returned as a struct which can be %% indexed by observable name. %% %% Examples %% -------- %% Example integration using default initial and parameter %% values: %% %% >> m = %(model_name)s(); %% >> tspan = [0 100]; %% >> [t y] = ode15s(@m.odes, tspan, m.get_initial_values()); %% %% Retrieving the observables: %% %% >> y_obs = m.get_observables(y) %% properties observables parameters end methods function self = %(model_name)s() %% Assign default parameter values %(params_str)s %% Define species indices (first row) and coefficients %% (second row) of named observables %(observables_str)s end function initial_values = get_initial_values(self) %% Return the vector of initial conditions for all %% species based on the values of the parameters %% as currently defined in the instance. %(initial_values_str)s end function y = odes(self, tspan, y0) %% Right hand side function for the ODEs %% Shorthand for the struct of model parameters p = self.parameters; %(odes_str)s end function y_obs = get_observables(self, y) %% Retrieve the trajectories for the model observables %% from a matrix of the trajectories of all model %% species. %% Initialize the struct of observable timecourses %% that we will return y_obs = struct(); %% Iterate over the observables; observable_names = fieldnames(self.observables); for i = 1:numel(observable_names) obs_matrix = self.observables.(observable_names{i}); if isempty(obs_matrix) y_obs.(observable_names{i}) = zeros(size(y, 1), 1); continue end species = obs_matrix(1, :); coefficients = obs_matrix(2, :); y_obs.(observable_names{i}) = ... y(:, species) * coefficients'; end end end end """, 0) % {'docstring': docstring, 'model_name': model_name, 'params_str':params_str, 'initial_values_str': initial_values_str, 'observables_str': observables_str, 'params_str': params_str, 'odes_str': odes_str}) return output.getvalue()
def _fix_underscores(name): if name.startswith('_'): return 'X' + name else: return name