Source code for tvb.datatypes.equations

# -*- coding: utf-8 -*-
#
#
# TheVirtualBrain-Scientific Package. This package holds all simulators, and
# analysers necessary to run brain-simulations. You can use it stand alone or
# in conjunction with TheVirtualBrain-Framework Package. See content of the
# documentation-folder for more details. See also http://www.thevirtualbrain.org
#
# (c) 2012-2023, Baycrest Centre for Geriatric Care ("Baycrest") and others
#
# 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 of the License, 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 MERCHANTABILITY or FITNESS FOR A
# PARTICULAR PURPOSE.  See the GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along with this
# program.  If not, see <http://www.gnu.org/licenses/>.
#
#
#   CITATION:
# When using The Virtual Brain for scientific publications, please cite it as explained here:
# https://www.thevirtualbrain.org/tvb/zwei/neuroscience-publications
#
#

"""

The Equation datatypes.

.. moduleauthor:: Stuart A. Knock <Stuart@tvb.invalid>

"""
import numpy
from scipy.special import gamma as sp_gamma
from tvb.basic.neotraits.api import HasTraits, Attr, Final, TupleEnum, SubformEnum
from tvb.simulator.backend.ref import RefBase


# In how many points should the equation be evaluated for the plot. Increasing this will
# give smoother results at the cost of some performance
DEFAULT_PLOT_GRANULARITY = 1024


[docs] class EquationsEnum(SubformEnum): """ Superclass of all enums that have equations as values """
# class Equation(basic.MapAsJson, core.Type): # todo: handle the MapAsJson functionality
[docs] class Equation(HasTraits): """Base class for Equation data types.""" equation = Attr( field_type=str, label="Equation as a string", doc=""" the equation as it should be interpreted by numexpr""") # todo: transform these parameters into plain declarative attrs parameters = Attr( field_type=dict, label="Parameters in a dictionary.", default=lambda: {}, doc="""Should be a list of the parameters and their meaning, Traits should be able to take defaults and sensible ranges from any traited information that was provided.""")
[docs] def summary_info(self): """ Gather scientifically interesting summary information from an instance of this datatype. """ return { "Equation type": self.__class__.__name__, "equation": self.equation, "parameters": self.parameters }
[docs] def evaluate(self, var): """ Generate a discrete representation of the equation for the space represented by ``var``. The argument ``var`` can represent a distance, or effective distance, for each node in a simulation. Or a time, or in principle any arbitrary `` space ``. ``var`` can be a single number, a numpy.ndarray or a ?scipy.sparse_matrix? TODO: think this last one is true, need to check as we need it for LocalConnectivity... """ ns = {'var': var} ns.update(self.parameters) return RefBase.evaluate(self.equation, ns)
[docs] def get_series_data(self, min_range=0, max_range=100, step=None): """ NOTE: The symbol from the equation which varies should be named: var Returns the series data needed for plotting this equation. """ if step is None: step = float(max_range - min_range) / DEFAULT_PLOT_GRANULARITY var = numpy.arange(min_range, max_range+step, step) var = var[numpy.newaxis, :] y = self.evaluate(var) result = list(zip(var.flat, y.flat)) return result, False
[docs] class TemporalApplicableEquation(Equation): """ Abstract class introduced just for filtering what equations to be displayed in UI, for setting the temporal component in Stimulus on region and surface. """
[docs] class FiniteSupportEquation(TemporalApplicableEquation): """ Equations that decay to zero as the variable moves away from zero. It is necessary to restrict spatial equation evaluated on a surface to this class, are . The main purpose of this class is to facilitate filtering in the UI, for patters on surface (stimuli surface and localConnectivity). """
[docs] class SpatialApplicableEquation(Equation): """ Abstract class introduced just for filtering what equations to be displayed in UI, for setting model parameters on the Surface level. """
[docs] class DiscreteEquation(FiniteSupportEquation): """ A special case for 'discrete' spaces, such as the regions, where each point in the space is effectively just assigned a value. """ equation = Attr( field_type=str, label="Discrete Equation", default="var", # locked=True, doc="""The equation defines a function of :math:`x`""")
[docs] class Linear(TemporalApplicableEquation): """ A linear equation. """ equation = Final( label="Linear Equation", default="a * var + b", # locked=True, doc=""":math:`result = a * x + b`""") parameters = Attr( field_type=dict, label="Linear Parameters", default=lambda: {"a": 1.0, "b": 0.0})
[docs] class RescaleInterval(Equation): """ Direct rescaling of an interval from [oldMin..oldMax] into [newMin..newMax] """ equation = Final( label="Rescale Equation", default="(var - oldMin) * (newMax - newMin) / (oldMax - oldMin) + newMin", doc=""":math:`result = (x - oldMin) * (newMax - newMin) / (oldMax - oldMin) + newMin`""") parameters = Attr( field_type=dict, label="Rescale Interval Parameters", default=lambda: {"newMin": 0.0, "newMax": 1.0, "oldMin": 0.0, "oldMax": 1.0})
[docs] class Absolute(Equation): """ Absolute value """ equation = Final( label="Absolute Equation", default="abs(var)", doc=""":math:`result = abs(x)`""")
[docs] class Identity(Equation): """ Identity value """ equation = Final( label="Identity Equation", default="var", doc=""":math:`result = x`""")
[docs] class Logarithm(Equation): """ Logarithm Equation """ equation = Final( label="Logarithm Equation", default="log(var)", doc=""":math:`result = log(x)`""")
[docs] class Gaussian(SpatialApplicableEquation, FiniteSupportEquation): """ A Gaussian equation. offset: parameter to extend the behaviour of this function when spatializing model parameters. """ equation = Final( label="Gaussian Equation", default="(amp * exp(-((var-midpoint)**2 / (2.0 * sigma**2))))+offset", doc=""":math:`(amp \\exp\\left(-\\left(\\left(x-midpoint\\right)^2 / \\left(2.0 \\sigma^2\\right)\\right)\\right)) + offset`""") parameters = Attr( field_type=dict, label="Gaussian Parameters", default=lambda: {"amp": 1.0, "sigma": 1.0, "midpoint": 0.0, "offset": 0.0})
[docs] class DoubleGaussian(FiniteSupportEquation): """ A Mexican-hat function approximated by the difference of Gaussians functions. """ equation = Final( label="Double Gaussian Equation", default="""(amp_1 * exp(-((var-midpoint_1)**2 / (2.0 * sigma_1**2)))) - """ """(amp_2 * exp(-((var-midpoint_2)**2 / (2.0 * sigma_2**2))))""", doc=""":math:`amp_1 \\exp\\left(-\\left((x-midpoint_1)^2 / \\left(2.0 \\sigma_1^2\\right)\\right)\\right) - amp_2 \\exp\\left(-\\left((x-midpoint_2)^2 / \\left(2.0 \\sigma_2^2\\right)\\right)\\right)`""") parameters = Attr( field_type=dict, label="Double Gaussian Parameters", default=lambda: {"amp_1": 0.5, "sigma_1": 20.0, "midpoint_1": 0.0, "amp_2": 1.0, "sigma_2": 10.0, "midpoint_2": 0.0})
[docs] class Sigmoid(SpatialApplicableEquation, FiniteSupportEquation): """ A Sigmoid equation. offset: parameter to extend the behaviour of this function when spatializing model parameters. """ equation = Final( label="Sigmoid Equation", default="(amp / (1.0 + exp(-1.8137993642342178 * (radius-var)/sigma))) + offset", doc=""":math:`(amp / (1.0 + \\exp(-\\pi/\\sqrt(3.0) (radius-x)/\\sigma))) + offset`""") parameters = Attr( field_type=dict, label="Sigmoid Parameters", default=lambda: {"amp": 1.0, "radius": 5.0, "sigma": 1.0, "offset": 0.0})
[docs] class GeneralizedSigmoid(TemporalApplicableEquation): """ A General Sigmoid equation. """ equation = Final( label="Generalized Sigmoid Equation", default="low + (high - low) / (1.0 + exp(-1.8137993642342178 * (var-midpoint)/sigma))", doc=""":math:`low + (high - low) / (1.0 + \\exp(-\\pi/\\sqrt(3.0) (x-midpoint)/\\sigma))`""") parameters = Attr( field_type=dict, label="Sigmoid Parameters", default=lambda: {"low": 0.0, "high": 1.0, "midpoint": 1.0, "sigma": 0.3})
[docs] class Sinusoid(TemporalApplicableEquation): """ A Sinusoid equation. """ equation = Final( label="Sinusoid Equation", default="amp * sin(6.283185307179586 * frequency * var)", doc=""":math:`amp \\sin(2.0 \\pi frequency x)` """) parameters = Attr( field_type=dict, label="Sinusoid Parameters", default=lambda: {"amp": 1.0, "frequency": 0.01})
[docs] class Cosine(TemporalApplicableEquation): """ A Cosine equation. """ equation = Final( label="Cosine Equation", default="amp * cos(6.283185307179586 * frequency * var)", doc=""":math:`amp \\cos(2.0 \\pi frequency x)` """) parameters = Attr( field_type=dict, label="Cosine Parameters", default=lambda: {"amp": 1.0, "frequency": 0.01})
[docs] class Alpha(TemporalApplicableEquation): """ An Alpha function belonging to the Exponential function family. """ equation = Final( label="Alpha Equation", default="""where((var-onset) > 0, (alpha * beta) / (beta - alpha) * (exp(-alpha * (var-onset)) """ """ - exp(-beta * (var-onset))), 0.0 * var)""", doc=""":math:`(\\alpha * \\beta) / (\\beta - \\alpha) * (\\exp(-\\alpha * (x-onset)) - \\exp(-\\beta * (x-onset)))` for :math:`(x-onset) > 0`""") parameters = Attr( field_type=dict, label="Alpha Parameters", default=lambda: {"onset": 0.5, "alpha": 13.0, "beta": 42.0})
[docs] class PulseTrain(TemporalApplicableEquation): """ A pulse train , offset with respect to the time axis. **Parameters**: * :math:`\\tau` : pulse width or pulse duration * :math:`T` : pulse repetition period * onset : time of the first pulse * amp : amplitude of the pulse """ equation = Final( label="Pulse Train", default="where((var>onset)&(((var-onset) % T) < tau), amp, 0)", doc=""":math:`\\left\\{{\\begin{array}{rl}amp,&{\\text{if }} ((var-onset) \\mod T) < \\tau \\space and \\space var > onset\\\\0, &{\\text{otherwise }}\\end{array}}\\right.`""") # onset is in milliseconds # T and tau are in milliseconds as well parameters = Attr( field_type=dict, default=lambda: {"T": 42.0, "tau": 13.0, "amp": 1.0, "onset": 30.0}, label="Pulse Train Parameters")
[docs] class HRFKernelEquation(Equation): """Base class for hemodynamic response functions."""
[docs] class Gamma(HRFKernelEquation): """ A Gamma function for the bold monitor. It belongs to the family of Exponential functions. **Parameters**: * :math:`\\tau` : Exponential time constant of the gamma function [seconds]. * :math:`n` : The phase delay of the gamma function. * :math: `factorial` : (n-1)!. numexpr does not support factorial yet. * :math: `a` : Amplitude factor after normalization. **Reference**: .. [B_1996] Geoffrey M. Boynton, Stephen A. Engel, Gary H. Glover and David J. Heeger (1996). Linear Systems Analysis of Functional Magnetic Resonance Imaging in Human V1. J Neurosci 16: 4207-4221 .. note:: might be filtered from the equations used in Stimulus and Local Connectivity. """ # TODO: Introduce a time delay in the equation (shifts the hrf onset) # """:math:`h(t) = \frac{(\frac{t-\delta}{\tau})^{(n-1)} e^{-(\frac{t-\delta}{\tau})}}{\tau(n-1)!}""" # delta = 2.05 seconds -- Additional delay in seconds from the onset of the # time-series to the beginning of the gamma hrf. # delay cannot be negative or greater than the hrf duration. equation = Final( label="Gamma Equation", default="((var / tau) ** (n - 1) * exp(-(var / tau)) )/ (tau * factorial)", doc=""":math:`h(var) = \\frac{(\\frac{var}{\\tau})^{(n-1)}\\exp{-(\\frac{var}{\\tau})}}{\\tau(n-1)!}`.""") parameters = Attr( field_type=dict, label="Gamma Parameters", default=lambda: {"tau": 1.08, "n": 3.0, "factorial": 2.0, "a": 0.1})
[docs] def evaluate(self, var): """ Generate a discrete representation of the equation for the space represented by ``var``. .. note: numexpr doesn't support factorial yet """ # compute the factorial n = int(self.parameters["n"]) product = 1 for i in range(n - 1): product *= i + 1 self.parameters["factorial"] = product _pattern = super().evaluate(var) _pattern /= max(_pattern) _pattern *= self.parameters["a"] return _pattern
[docs] class DoubleExponential(HRFKernelEquation): """ A difference of two exponential functions to define a kernel for the bold monitor. **Parameters** : * :math:`\\tau_1`: Time constant of the second exponential function [s] * :math:`\\tau_2`: Time constant of the first exponential function [s]. * :math:`f_1` : Frequency of the first sine function [Hz]. * :math:`f_2` : Frequency of the second sine function [Hz]. * :math:`amp_1`: Amplitude of the first exponential function. * :math:`amp_2`: Amplitude of the second exponential function. * :math:`a` : Amplitude factor after normalization. **Reference**: .. [P_2000] Alex Polonsky, Randolph Blake, Jochen Braun and David J. Heeger (2000). Neuronal activity in human primary visual cortex correlates with perception during binocular rivalry. Nature Neuroscience 3: 1153-1159 """ equation = Final( label="Double Exponential Equation", default="((amp_1 * exp(-var/tau_1) * sin(2.*pi*f_1*var)) - (amp_2 * exp(-var/ tau_2) * sin(2.*pi*f_2*var)))", doc=""":math:`h(var) = amp_1\\exp(\\frac{-var}{\tau_1}) \\sin(2\\cdot\\pi f_1 \\cdot var) - amp_2\\cdot \\exp(-\\frac{var} {\\tau_2})*\\sin(2\\pi f_2 var)`.""") parameters = Attr( field_type=dict, label="Double Exponential Parameters", default=lambda: {"tau_1": 7.22, "f_1": 0.03, "amp_1": 0.1, "tau_2": 7.4, "f_2": 0.12, "amp_2": 0.1, "a": 0.1, "pi": numpy.pi})
[docs] def evaluate(self, var): """ Generate a discrete representation of the equation for the space represented by ``var``. """ _pattern = super().evaluate(var) _pattern /= max(_pattern) _pattern *= self.parameters["a"] return _pattern
[docs] class FirstOrderVolterra(HRFKernelEquation): """ Integral form of the first Volterra kernel of the three used in the Ballon Windekessel model for computing the Bold signal. This function describes a damped Oscillator. **Parameters** : * :math:`\\tau_s`: Dimensionless? exponential decay parameter. * :math:`\\tau_f`: Dimensionless? oscillatory parameter. * :math:`k_1` : First Volterra kernel coefficient. * :math:`V_0` : Resting blood volume fraction. **References** : .. [F_2000] Friston, K., Mechelli, A., Turner, R., and Price, C., *Nonlinear Responses in fMRI: The Balloon Model, Volterra Kernels, and Other Hemodynamics*, NeuroImage, 12, 466 - 477, 2000. """ equation = Final( label="First Order Volterra Kernel", default="""1/3. * exp(-0.5*(var / tau_s)) * (sin(sqrt(1./tau_f - 1./(4.*tau_s**2)) * var)) / """ """(sqrt(1./tau_f - 1./(4.*tau_s**2)))""", doc=""":math:`G(t - t^{\\prime}) = e^{\\frac{1}{2} \\left(\\frac{t - t^{\\prime}}{\\tau_s} \\right)} \\frac{\\sin\\left((t - t^{\\prime}) \\sqrt{\\frac{1}{\\tau_f} - \\frac{1}{4 \\tau_s^2}}\\right)} {\\sqrt{\\frac{1}{\\tau_f} - \\frac{1}{4 \\tau_s^2}}} \\; \\; \\; \\; \\; \\; for \\; \\; \\; t \\geq t^{\\prime} = 0 \\; \\; \\; \\; \\; \\; for \\; \\; \\; t < t^{\\prime}`.""") parameters = Attr( field_type=dict, label="Mixture of Gammas Parameters", default=lambda: {"tau_s": 0.8, "tau_f": 0.4, "k_1": 5.6, "V_0": 0.02})
[docs] class MixtureOfGammas(HRFKernelEquation): """ A mixture of two gamma distributions to create a kernel similar to the one used in SPM. >> import scipy.stats as sp_stats >> import numpy >> t = numpy.linspace(1,20,100) >> a1, a2 = 6., 10. >> lambda = 1. >> c = 0.5 >> hrf = sp_stats.gamma.pdf(t, a1, lambda) - c * sp_stats.gamma.pdf(t, a2, lambda) gamma.pdf(x, a, theta) = (lambda*x)**(a-1) * exp(-lambda*x) / gamma(a) a : shape parameter theta: 1 / lambda : scale parameter **References**: .. [G_1999] Glover, G. *Deconvolution of Impulse Response in Event-Related BOLD fMRI*. NeuroImage 9, 416-429, 1999. **Parameters**: * :math:`a_{1}` : shape parameter first gamma pdf. * :math:`a_{2}` : shape parameter second gamma pdf. * :math:`\\lambda` : scale parameter first gamma pdf. Default values are based on [G_1999]_: * :math:`a_{1} - 1 = n_{1} = 5.0` * :math:`a_{2} - 1 = n_{2} = 12.0` * :math:`c \\equiv a_{2} = 0.4` Alternative values :math:`a_{2}=10` and :math:`c=0.5` NOTE: gamma_a_1 and gamma_a_2 are placeholders, the true values are computed before evaluating the expression, because numexpr does not support certain functions. NOTE: [G_1999]_ used a different analytical function that can be approximated by this difference of gamma pdfs """ equation = Final( label="Mixture of Gammas", default="(l * var)**(a_1-1) * exp(-l*var) / gamma_a_1 - c * (l*var)**(a_2-1) * exp(-l*var) / gamma_a_2", doc=""":math:`\\frac{\\lambda \\,t^{a_{1} - 1} \\,\\, \\exp^{-\\lambda \\,t}}{\\Gamma(a_{1})} - 0.5 \\frac{\\lambda \\,t^{a_{2} - 1} \\,\\, \\exp^{-\\lambda \\,t}}{\\Gamma(a_{2})}`.""") parameters = Attr( field_type=dict, label="Double Exponential Parameters", default=lambda: {"a_1": 6.0, "a_2": 13.0, "l": 1.0, "c": 0.4, "gamma_a_1": 1.0, "gamma_a_2": 1.0})
[docs] def evaluate(self, var): """ Generate a discrete representation of the equation for the space represented by ``var``. .. note: numexpr doesn't support gamma function """ # get gamma functions self.parameters["gamma_a_1"] = sp_gamma(self.parameters["a_1"]) self.parameters["gamma_a_2"] = sp_gamma(self.parameters["a_2"]) return super().evaluate(var)