Source code for function

# ---
# Filename:
# Description:
# Author: Subhasis Ray
# Maintainer:
# Created: Tue Sep  9 17:59:50 2014 (+0530)
# Version:
# Last-Updated: Sun Dec 20 00:02:50 2015 (-0500)
#           By: subha
#     Update #: 4
# URL:
# Keywords:
# Compatibility:

# Commentary:

# Change log:
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 3, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program; see the file COPYING.  If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
# Floor, Boston, MA 02110-1301, USA.

# Code:

import numpy as np
import sys
import matplotlib.pyplot as plt

import moose

simtime = 1.0

[docs]def example(): """Function objects can be used to evaluate expressions with arbitrary number of variables and constants. We can assign expression of the form:: f(c0, c1, ..., cM, x0, x1, ..., xN, y0,..., yP ) where `c_i`'s are constants and `x_i`'s and `y_i`'s are variables. The constants must be defined before setting the expression and variables are connected via messages. The constants can have any name, but the variable names must be of the form x{i} or y{i} where i is increasing integer starting from 0. The `x_i`'s are field elements and you have to set their number first (function.x.num = N). Then you can connect any source field sending out double to the 'input' destination field of the `x[i]`. The `y_i`'s are useful when the required variable is a value field and is not available as a source field. In that case you connect the `requestOut` source field of the function element to the `get{Field}` destination field on the target element. The `y_i`'s are automatically added on connecting. Thus, if you call:: moose.connect(function, 'requestOut', a, 'getSomeField') moose.connect(function, 'requestOut', b, 'getSomeField') then ``a.someField`` will be assigned to ``y0`` and ``b.someField`` will be assigned to ``y1``. In this example we evaluate the expression: ``z = c0 * exp(c1 * x0) * cos(y0)`` with x0 ranging from -1 to +1 and y0 ranging from -pi to +pi. These values are stored in two stimulus tables called xtab and ytab respectively, so that at each timestep the next values of x0 and y0 are assigned to the function. Along with the value of the expression itself we also compute its derivative with respect to y0 and its derivative with respect to time (rate). The former uses a five-point stencil for the numerical differentiation and has a glitch at y=0. The latter uses backward difference divided by dt. Unlike Func class, the number of variables and constants are unlimited in Function and you can set all the variables via messages. """ demo = moose.Neutral('/model') function = moose.Function('/model/function') function.c['c0'] = 1.0 function.c['c1'] = 2.0 function.x.num = 1 function.expr = 'c0 * exp(c1*x0) * cos(y0) + sin(t)' # mode 0 - evaluate function value, derivative and rate # mode 1 - just evaluate function value, # mode 2 - evaluate derivative, # mode 3 - evaluate rate function.mode = 0 function.independent = 'y0' nsteps = 1000 xarr = np.linspace(0.0, 1.0, nsteps) # Stimulus tables allow you to store sequences of numbers which # are delivered via the 'output' message at each time step. This # is a placeholder and in real scenario you will be using any # sourceFinfo that sends out a double value. input_x = moose.StimulusTable('/xtab') input_x.vector = xarr input_x.startTime = 0.0 input_x.stepPosition = xarr[0] input_x.stopTime = simtime moose.connect(input_x, 'output', function.x[0], 'input') yarr = np.linspace(-np.pi, np.pi, nsteps) input_y = moose.StimulusTable('/ytab') input_y.vector = yarr input_y.startTime = 0.0 input_y.stepPosition = yarr[0] input_y.stopTime = simtime moose.connect(function, 'requestOut', input_y, 'getOutputValue') # data recording result = moose.Table('/ztab') moose.connect(result, 'requestOut', function, 'getValue') derivative = moose.Table('/zprime') moose.connect(derivative, 'requestOut', function, 'getDerivative') rate = moose.Table('/dz_by_dt') moose.connect(rate, 'requestOut', function, 'getRate') x_rec = moose.Table('/xrec') moose.connect(x_rec, 'requestOut', input_x, 'getOutputValue') y_rec = moose.Table('/yrec') moose.connect(y_rec, 'requestOut', input_y, 'getOutputValue') dt = simtime/nsteps for ii in range(32): moose.setClock(ii, dt) moose.reinit() moose.start(simtime) # Uncomment the following lines and the import matplotlib.pyplot as plt on top # of this file to display the plot. plt.subplot(3,1,1) plt.plot(x_rec.vector, result.vector, 'r-', label='z = {}'.format(function.expr)) z = function.c['c0'] * np.exp(function.c['c1'] * xarr) * np.cos(yarr) + np.sin(np.arange(len(xarr)) * dt) plt.plot(xarr, z, 'b--', label='numpy computed') plt.xlabel('x') plt.ylabel('z') plt.legend() plt.subplot(3,1,2) plt.plot(y_rec.vector, derivative.vector, 'r-', label='dz/dy0') # derivatives computed by putting x values in the analytical formula dzdy = function.c['c0'] * np.exp(function.c['c1'] * xarr) * (- np.sin(yarr)) plt.plot(yarr, dzdy, 'b--', label='numpy computed') plt.xlabel('y') plt.ylabel('dz/dy') plt.legend() plt.subplot(3,1,3) # *** BEWARE *** The first two entries are spurious. Entry 0 is # *** from reinit sending out the defaults. Entry 2 is because # *** there is no lastValue for computing real forward difference. plt.plot(np.arange(2, len(rate.vector), 1) * dt, rate.vector[2:], 'r-', label='dz/dt') dzdt = np.diff(z)/dt plt.plot(np.arange(0, len(dzdt), 1.0) * dt, dzdt, 'b--', label='numpy computed') plt.xlabel('t') plt.ylabel('dz/dt') plt.legend() plt.tight_layout()
if __name__ == '__main__': example() # # ends here