def test_nsolve():
    # onedimensional
    from sympy import Symbol, sin, pi
    x = Symbol('x')
    assert nsolve(sin(x), 2) - pi.evalf() < 1e-16
    assert nsolve(Eq(2*x, 2), x, -10) == nsolve(2*x - 2, -10)
    # multidimensional
    x1 = Symbol('x1')
    x2 = Symbol('x2')
    f1 = 3 * x1**2 - 2 * x2**2 - 1
    f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
    f = Matrix((f1, f2)).T
    F = lambdify((x1, x2), f.T, modules='mpmath')
    for x0 in [(-1, 1), (1, -2), (4, 4), (-4, -4)]:
        x = nsolve(f, (x1, x2), x0, tol=1.e-8)
        assert mnorm(F(*x),1) <= 1.e-10
    # The Chinese mathematician Zhu Shijie was the very first to solve this
    # nonlinear system 700 years ago (z was added to make it 3-dimensional)
    x = Symbol('x')
    y = Symbol('y')
    z = Symbol('z')
    f1 = -x + 2*y
    f2 = (x**2 + x*(y**2 - 2) - 4*y)  /  (x + 4)
    f3 = sqrt(x**2 + y**2)*z
    f = Matrix((f1, f2, f3)).T
    F = lambdify((x,  y,  z), f.T, modules='mpmath')
    def getroot(x0):
        root = nsolve((f1,  f2,  f3), (x,  y,  z), x0)
        assert mnorm(F(*root),1) <= 1.e-8
        return root
    assert map(round,  getroot((1,  1,  1))) == [2.0,  1.0,  0.0]

def velocity_field(psi): #takes a symbolic function and returns two lambda functions
#to evaluate the derivatives in both x and y.
   global w
   if velocity_components:
      u = lambdify((x,y), eval(x_velocity), modules='numpy')
      v = lambdify((x,y), eval(y_velocity), modules='numpy')
   else:
      if is_complex_potential:
         print "Complex potential, w(z) given"
         #define u, v symbolically as the imaginary part of the derivatives
         u = lambdify((x, y), sympy.im(psi.diff(y)), modules='numpy')
         v = lambdify((x, y), -sympy.im(psi.diff(x)), modules='numpy')
      else:
         #define u,v as the derivatives 
         print "Stream function, psi given"
         u = sympy.lambdify((x, y), psi.diff(y), 'numpy')
         v = sympy.lambdify((x, y), -psi.diff(x), 'numpy')
   if (branch_cuts): # If it's indicated that there are branch cuts in the mapping,
                      # then we need to return vectorized numpy functions to evaluate
                      # everything numerically, instead of symbolically 
                      # This of course results in a SIGNIFICANT time increase
                      #   (I don't know how to handle more than the primitive root
                      #   (symbolically in Sympy
      return np.vectorize(u),np.vectorize(v)
   else:
       # If there are no branch cuts, then return the symbolic lambda functions (MUCH faster)
      return u,v

def test_msolve():
    x1 = Symbol('x1')
    x2 = Symbol('x2')
    f1 = 3 * x1**2 - 2 * x2**2 - 1
    f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
    f = Matrix((f1, f2)).T
    F = lambdify((x1, x2), f.T)
    # numeric.newton is tested in this example too
    for x0 in [(-1., 1.), (1., -2.), (4., 4.), (-4., -4.)]:
        x = msolve((x1, x2), f, x0, tol=1.e-8)
        assert maxnorm(F(*x)) <= 1.e-11
    # The Chinese mathematician Zhu Shijie was the very first to solve this
    # nonlinear system 700 years ago (z was added to make it 3-dimensional)
    x = Symbol('x')
    y = Symbol('y')
    z = Symbol('z')
    f1 = -x + 2*y
    f2 = (x**2 + x*(y**2 - 2) - 4*y)  /  (x + 4)
    f3 = sqrt(x**2 + y**2)*z
    f = Matrix((f1, f2, f3)).T
    F = lambdify((x,  y,  z), f.T)
    def getroot(x0):
        root = msolve((x,  y,  z),  (f1,  f2,  f3),  x0)
        assert maxnorm(F(*root)) <= 1.e-8
        return root
    assert map(round,  getroot((1.,  1.,  1.))) == [2.0,  1.0,  0.0]

    def compile_function(self, module="numpy"):
        import sympy
        from sympy.utilities.lambdify import lambdify

        expr = sympy.sympify(self._str_expression)

        rvars = sympy.symbols([s.name for s in expr.free_symbols], real=True)
        real_expr = expr.subs({orig: real_ for (orig, real_) in zip(expr.free_symbols, rvars)})
        # just replace with the assumption that all our variables are real
        expr = real_expr

        eval_expr = expr.evalf()
        # Extract parameters
        parameters = [symbol for symbol in expr.free_symbols if symbol.name != "x"]
        parameters.sort(key=lambda x: x.name)  # to have a reliable order
        # Extract x
        x, = [symbol for symbol in expr.free_symbols if symbol.name == "x"]
        # Create compiled function
        self._f = lambdify([x] + parameters, eval_expr, modules=module, dummify=False)
        parnames = [symbol.name for symbol in parameters]
        self._parameter_strings = parnames
        for parameter in parameters:
            grad_expr = sympy.diff(eval_expr, parameter)
            setattr(
                self,
                "_f_grad_%s" % parameter.name,
                lambdify([x] + parameters, grad_expr.evalf(), modules=module, dummify=False),
            )

            setattr(
                self,
                "grad_%s" % parameter.name,
                _fill_function_args(getattr(self, "_f_grad_%s" % parameter.name)).__get__(self, Expression),
            )

 def twin_function_expr(self, value):
     if not value:
         self.twin_function = None
         self.twin_inverse_function = None
         self._twin_function_expr = ""
         self._twin_inverse_sympy = None
         return
     expr = sympy.sympify(value)
     if len(expr.free_symbols) > 1:
         raise ValueError("The expression must contain only one variable.")
     elif len(expr.free_symbols) == 0:
         raise ValueError("The expression must contain one variable, "
                          "it contains none.")
     x = tuple(expr.free_symbols)[0]
     self.twin_function = lambdify(x, expr.evalf())
     self._twin_function_expr = value
     if not self.twin_inverse_function:
         y = sympy.Symbol(x.name + "2")
         try:
             inv = sympy.solveset(sympy.Eq(y, expr), x)
             self._twin_inverse_sympy = lambdify(y, inv)
             self._twin_inverse_function = None
         except BaseException:
             # Not all may have a suitable solution.
             self._twin_inverse_function = None
             self._twin_inverse_sympy = None
             _logger.warning(
                 "The function {} is not invertible. Setting the value of "
                 "{} will raise an AttributeError unless you set manually "
                 "``twin_inverse_function_expr``. Otherwise, set the "
                 "value of its twin parameter instead.".format(value, self))

def find_CI(quartile, variable, critical, search = "left", grid_left = None, grid_right = None, precision = 0.001):
    '''
    We take {regions} points in our domain and test every point for its likelihood value,
    then we look at points where values of ll ratio criteria cross critical value of chi^2 distribution
    It means, we should choose {regions} in such way that at least one point lays in the confidence interval
    at the first step.
    '''
    regions = 20
    solution = solutions[quartile]
    ll_hat = -0.5 * commons.chi2(solution)
    critical_val = chi2.ppf(critical, df=1)
    other1 = variable_to_others[variable][0]
    other2 = variable_to_others[variable][1]
    variable_bounds = (grid_left if not grid_left is None else commons.bounds[quartile][variable_index[variable]][0],
                       grid_right if not grid_right is None else commons.bounds[quartile][variable_index[variable]][1])
    grid = np.linspace(start = variable_bounds[0], stop = variable_bounds[1], num = regions)
    other_bounds = (commons.bounds[quartile][variable_index[other1]], commons.bounds[quartile][variable_index[other2]])
    ll_prev_criteria = None
    #print("Critical value is {}".format(critical_val))
    for i, v in enumerate(grid):
        #print("Check {} = {}".format(variable, v))
        ys = sp.Matrix([variable_ref[other1], variable_ref[other2]])
        pll_func_sym = commons.chi2_sym.subs(variable_ref[variable], v)
        pll_func = lambda args: lambdify(ys, pll_func_sym, 'numpy')(*args)[0, 0]
        pll_func_jacob = lambda args: lambdify(ys, pll_func_sym.jacobian(ys), 'numpy')(*args)
        result = opt.minimize(pll_func, [solution[variable_index[other1]], solution[variable_index[other2]]],
                          method='TNC',
                          jac=pll_func_jacob,
                          bounds=other_bounds)
        pll_val = -.5 * result.fun
        ll_criteria = 2 * (ll_hat - pll_val)
        #print("\tProfile LL ratio value = {}".format(ll_criteria))
        if not ll_prev_criteria is None:
            l = grid[i - 1]
            r = grid[i]
            if ll_prev_criteria > critical_val and ll_criteria < critical_val and search == "left":
                #print("Critical point lays between {} and {}\n".format(l, r))
                if (r - l) < precision:
                    return l, r
                else:
                    return find_CI(quartile, variable, critical,
                                   search = search,
                                   grid_left = l,
                                   grid_right = r,
                                   precision = precision)
            elif ll_prev_criteria < critical and ll_criteria > critical and search == "right":
                #print("Critical point lays between {} and {}\n".format(l, r))
                if (r - l) < precision:
                    return l, r
                else:
                    return find_CI(quartile, variable, critical,
                                   search = search,
                                   grid_left = l,
                                   grid_right = r,
                                   precision = precision)
        ll_prev_criteria = ll_criteria

def msolve(args, f, x0, tol=None, maxsteps=None, verbose=False, norm=None,
           modules=['mpmath', 'sympy']):
    """
    Solves a nonlinear equation system numerically.

    f is a vector function of symbolic expressions representing the system.
    args are the variables.
    x0 is a starting vector close to a solution.

    Be careful with x0, not using floats might give unexpected results.

    Use modules to specify which modules should be used to evaluate the
    function and the Jacobian matrix. Make sure to use a module that supports
    matrices. For more information on the syntax, please see the docstring
    of lambdify.

    Currently only fully determined systems are supported.

    >>> from sympy import Symbol, Matrix
    >>> x1 = Symbol('x1')
    >>> x2 = Symbol('x2')
    >>> f1 = 3 * x1**2 - 2 * x2**2 - 1
    >>> f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
    >>> msolve((x1, x2), (f1, f2), (-1., 1.))
    [-1.19287309935246]
    [ 1.27844411169911]
    """
    if isinstance(f,  (list,  tuple)):
        f = Matrix(f).T
    if len(args) != f.cols:
        raise NotImplementedError('need exactly as many variables as equations')
    if verbose:
        print 'f(x):'
        print f
    # derive Jacobian
    J = f.jacobian(args)
    if verbose:
        print 'J(x):'
        print J
    # create functions
    f = lambdify(args, f.T, modules)
    J = lambdify(args, J, modules)
    # solve system using Newton's method
    kwargs = {}
    if tol:
        kwargs['tol'] = tol
    if maxsteps:
        kwargs['maxsteps'] = maxsteps
    kwargs['verbose'] = verbose
    if norm:
        kwargs['norm'] = norm
    x = newton(f, x0, J, **kwargs)
    return x

    def __init__(self, state_vector, next_state_func, observation_vector, observation_func,
                 x0, initial_cov, process_noise_cov, measurement_noise_cov):
        """Create the filter

        :param state_vector: sympy matrix `x`, the full state of the filter
        :param next_state_func: sympy matrix, `x_next` as a function of `x`
            - The dynamics jacobian will be generated internally
        """
        dynamics_jacobian = self.form_dynamics_jacobian(state_vector, next_state_func)
        observation_jacobian = self.form_observation_jacobian(state_vector, observation_vector, observation_func)
        self.x_dim = dynamics_jacobian.shape[0]
        self.observation_dim = observation_vector.shape[0]

        #
        # Input validation
        # Before the arguments even get off the ship, we validate them

        # observation jacobian asserts
        assert observation_jacobian.shape[1] == self.x_dim, "Observation jacobian has the wrong dimension to map to state"
        assert observation_jacobian.shape[0] == self.observation_dim, "Observation jacobian is not the right shape"

        # measurement covariance asserts
        assert measurement_noise_cov.shape[0] == measurement_noise_cov.shape[1], \
            "Measurement noise is not the same shape as symbolic measurement"
        assert measurement_noise_cov.shape[0] == observation_vector.shape[0], "Measurement noise is not square"

        # process covariance asserts
        assert process_noise_cov.shape[0] == process_noise_cov.shape[1], "Process noise is not square"
        assert process_noise_cov.shape[0] == state_vector.shape[0], "Process noise is not the same shape as symbolic state"

        # initial state asserts
        assert x0.shape[0] == state_vector.shape[0], "Supplied initial state is not the same shape as symbolic state"

        self.make_A = lambdify(
            state_vector,
            dynamics_jacobian,
            'numpy'
        )

        self.make_H = lambdify(
            state_vector,
            observation_jacobian,
            'numpy'
        )

        #
        # State/Covariance initialization
        #
        self.x = x0
        self.P = initial_cov
        self.Q = process_noise_cov
        self.R = measurement_noise_cov

def test_1d():
    B = gen_BrownianMotion()
    Bs = implemented_function("B", B)
    t = sympy.Symbol('t')
    expr = 3*sympy.exp(Bs(t)) + 4
    expected = 3*np.exp(B.y)+4
    ee_vec = lambdify(t, expr, "numpy")
    assert_almost_equal(ee_vec(B.x), expected)
    # with any arbitrary symbol
    b = sympy.Symbol('b')
    expr = 3*sympy.exp(Bs(b)) + 4
    ee_vec = lambdify(b, expr, "numpy")
    assert_almost_equal(ee_vec(B.x), expected)

def make_callable(model, variables, mode=None):
    """
    Take a SymPy object and create a callable function.

    Parameters
    ----------
    model, SymPy object
        Abstract representation of function
    variables, list
        Input variables, ordered in the way the return function will expect
    mode, ['numpy', 'numexpr', 'sympy'], optional
        Method to use when 'compiling' the function. SymPy mode is
        slow and should only be used for debugging. If Numexpr is installed,
        it can offer speed-ups when calling the energy function many
        times on multi-core CPUs.

    Returns
    -------
    Function that takes arguments in the same order as 'variables'
    and returns the energy according to 'model'.

    Examples
    --------
    None yet.
    """
    energy = None
    if mode is None:
        # no mode specified; use numexpr if available, otherwise numpy
        # Note: numexpr support appears to break in multi-component situations
        # See numexpr#167 on GitHub for details
        # For now, default to numpy until a solution/workaround is available
        #if _NUMEXPR:
        #    mode = 'numexpr'
        #else:
        #    mode = 'numpy'
        mode = 'numpy'

    if mode == 'sympy':
        energy = lambda *vs: model.subs(zip(variables, vs)).evalf()
    elif mode == 'numpy':
        logical_np = [{'And': np.logical_and, 'Or': np.logical_or}, 'numpy']
        energy = lambdify(tuple(variables), model, dummify=True,
                          modules=logical_np, printer=NumPyPrinter)
    elif mode == 'numexpr':
        energy = lambdify(tuple(variables), model, dummify=True,
                          modules='numexpr', printer=SpecialNumExprPrinter)
    else:
        energy = lambdify(tuple(variables), model, dummify=True,
                          modules=mode)

    return energy

def make_callable(model, variables, mode=None):
    """
    Take a SymPy object and create a callable function.

    Parameters
    ----------
    model, SymPy object
        Abstract representation of function
    variables, list
        Input variables, ordered in the way the return function will expect
    mode, ['numpy', 'numba', 'sympy'], optional
        Method to use when 'compiling' the function. SymPy mode is
        slow and should only be used for debugging. If Numba is installed,
        it can offer speed-ups when calling the energy function many
        times on multi-core CPUs.

    Returns
    -------
    Function that takes arguments in the same order as 'variables'
    and returns the energy according to 'model'.

    Examples
    --------
    None yet.
    """
    energy = None
    if mode is None:
        # no mode specified; use numba if available, otherwise numpy
        if _NUMBA:
            mode = 'numba'
        else:
            mode = 'numpy'

    if mode == 'sympy':
        energy = lambda *vs: model.subs(zip(variables, vs)).evalf()
    elif mode == 'numpy':
        energy = lambdify(tuple(variables), model, dummify=True,
                          modules=[{'where': np.where}, 'numpy'], printer=NumPyPrinter)
    elif mode == 'numba':
        variables = tuple(variables)
        varsig = 'float64({})'.format(','.join(['float64'] * len(variables)))
        energy = lambdify(variables, model, dummify=True,
                          modules=[{'where': nbwhere}, 'numpy'], printer=NumPyPrinter)
        # target=parallel seems to incur too much overhead on small arrays
        energy = _NUMBA.vectorize([varsig], nopython=True, target='cpu')(energy)
    else:
        energy = lambdify(tuple(variables), model, dummify=True,
                          modules=mode)

    return energy

 def num_eval(self, symbol=None, values=None, uncertainties=None):
     if symbol == None:
         return self.value(), self.uncertainty()
     
     valfunc = lambdify(symbol, self.value(), 'numpy')
     if uncertainties == None:        
         uncertaintyfunc = lambdify(symbol, self.uncertainty(), 'numpy')
         return valfunc(values), uncertaintyfunc(values)
     else:
         u_sym = sympy.symbols('u_sym')
         usquared = self.uncertainty_squared() + \
                 (self.formula.diff(symbol).subs(self.values) * u_sym)**2
         u = sympy.sqrt(usquared)
         uncertaintyfunc = lambdify((symbol, u_sym), u, 'numpy')
         return valfunc(values), uncertaintyfunc(values, uncertainties)

def test_step_function():
    # test step function
    # step function is a function of t
    s = step_function([0,4,5],[2,4,6])
    tval = np.array([-0.1,0,3.9,4,4.1,5.1])
    lam = lambdify(t, s)
    assert_array_equal(lam(tval), [0, 2, 2, 4, 4, 6])
    s = step_function([0,4,5],[4,2,1])
    lam = lambdify(t, s)
    assert_array_equal(lam(tval), [0, 4, 4, 2, 2, 1])
    # Name default
    assert_false(re.match(r'step\d+\(t\)$', str(s)) is None)
    # Name reloaded
    s = step_function([0,4,5],[4,2,1], name='goodie_goodie_yum_yum')
    assert_equal(str(s), 'goodie_goodie_yum_yum(t)')

def test_implemented_function():
    # Here we check if the default returned functions are anonymous - in
    # the sense that we can have more than one function with the same name
    f = implemented_function('f', lambda x: 2*x)
    g = implemented_function('f', lambda x: np.sqrt(x))
    l1 = lambdify(x, f(x))
    l2 = lambdify(x, g(x))
    assert_equal(str(f(x)), str(g(x)))
    assert_equal(l1(3), 6)
    assert_equal(l2(3), np.sqrt(3))
    # check that we can pass in a sympy function as input
    func = sympy.Function('myfunc')
    assert_false(hasattr(func, '_imp_'))
    f = implemented_function(func, lambda x: 2*x)
    assert_true(hasattr(func, '_imp_'))

def test_blocks():
    on_off = [[1,2],[3,4]]
    tval = np.array([0.4,1.4,2.4,3.4])
    b = blocks(on_off)
    lam = lambdify(t, b)
    assert_array_equal(lam(tval), [0, 1, 0, 1])
    b = blocks(on_off, amplitudes=[3,5])
    lam = lambdify(t, b)
    assert_array_equal(lam(tval), [0, 3, 0, 5])
    # Check what happens with names
    # Default is from step function
    assert_false(re.match(r'step\d+\(t\)$', str(b)) is None)
    # Can pass in another
    b = blocks(on_off, name='funky_chicken')
    assert_equal(str(b), 'funky_chicken(t)')

def getRelError(variables,func,showFunc = False):
	"""getError generates a function to calculate the error of a function. I.E. a function that
	gives you the error in its result given the error in its input. Output function is numpy ready. 
	Output function will take twice as many args as variables, one for the var and one for the error in that var.

	arguments
	variables      : a list of sympy symbols in func
	errorVariables : list of sympy ymbols representing the error in each value of variables. must be the same length as variables
	func           : a function containing your variables that you want the error of"""
	ErrorFunc = 0 # need to set function to start value
	erVars    = {}
	for i in range(len(variables)): #run through all variables in the function
		v                  = variables[i]
		dv                 =  Symbol('d'+str(v), positive = True)
		erVars['d'+str(v)] = dv
		D                  = (diff(func,v)*dv)**2
		ErrorFunc          += D

	ErrorFunc = sqrt(ErrorFunc)/func
	if showFunc:
		pprint(ErrorFunc)
		
	variables.extend(erVars.values()) #create a list of all sympy symbols involved
	func      = lambdify(tuple(variables), ErrorFunc ,"numpy") #convert ErrorFunc to a numpy rteady python lambda
	return(func)

def get_error(expr):
    """ calculates error of an expression possibly containing quantities

    Returns: tuple of error and error formula (as Expr object)
    """
    integrand = 0
    error_formula = 0
    for varToDiff in expr.free_symbols:
        if varToDiff.error is not None:
            differential = diff(expr,varToDiff)
            error_formula += ( Symbol(varToDiff.name+"_err",positive=True) * differential )**2
            diffFunction = lambdify(differential.free_symbols,differential, modules="numpy")

            diffValues = []
            for var in differential.free_symbols:
                diffValues.append(var.value)

            integrand += ( varToDiff.error*diffFunction(*diffValues) )**2
    if isinstance(integrand,np.ndarray):
        if (integrand==0).all():
            return (None,None)
    elif integrand == 0:
        return (None,None)

    return (np.sqrt(integrand),sym_sqrt (error_formula))

def get_ode_fcn_floating_dp(g_, a1_, L1_, m1_, I1_, a2_, L2_, m2_, I2_):
    """ Returns a function object that can be called with fnc(t,x, u),
    where x = [q1, q2, q3, q4, qd1, qd2, qd3, qd4], and tau = [tau1, tau2, tau3, tau4], is the torques at each 
    joint respectively. The function implements the ode for the floating inverted
    double pendulum.

    For faster return of already constructed models, the set of parameters
    are checked against a pickled dict.
    """
    

    params = (g_, a1_, L1_, m1_, I1_, a2_, L2_, m2_, I2_)

    
    tau1, tau2, tau3, tau4 = sy.symbols('tau1, tau2, tau3, tau4')
    q1, q2, q3, q4, qd1, qd2, qd3, qd4 = symbols('q1, q2, q3, q4, qd1, qd2, qd3, qd4', real=True)

    s1 = sy.sin(q1); s1_ = Symbol('s1')
    c1 = sy.cos(q1); c1_ = Symbol('c1')
    s2 = sy.sin(q2); s2_ = Symbol('s2')
    c2 = sy.cos(q2); c2_ = Symbol('c2')
    s12 = sy.sin(q1+q2); s12_ = Symbol('s12') 
    c12 = sy.cos(q1+q2); c12_ = Symbol('c12') 



    odes = get_symbolic_ode_floating_dp(params)

    # Substitute functions for faster evaluation
    odes = odes.subs(s1, s1_).subs(c1,c1_).subs(s2,s2_).subs(c2,c2_).subs(s12,s12_).subs(c12,c12_)

    lmb = lambdify( (q1, q2, q3, q4, q1dot, q2dot, q3dot, q4dot, s1_, c1_, s2_, c2_, s12_, c12_, tau1, tau2, tau3, tau4), odes)
    return partial(lambda_ode, lambdafunc=lmb)

def fit(x_data, y_data, fit_function, parameters, weighted=None):

	args = [x_data]
	args.extend(parameters)
	np_func = lambdify(tuple(args), fit_function, "numpy")

	start_params = []
	for p in parameters:
		if p.value == None:
			start_params.append(np.float_(1))
		else:
			if isinstance(p.value,np.ndarray):
				raise ValueError("fit parameter '%s' is a data set." % p.name)
			else:
				start_params.append(p.value)

	if weighted is False:
		errors = None
	else:
		errors = y_data.error

	if weighted is True and y_data.error is None:
		raise RuntimeError("can't perform weighted fit because error of '%s' is not set." % y_data.name)
	params_opt, params_covar = curve_fit (np_func,x_data.value,y_data.value,sigma=errors,p0=start_params)
	params_err = np.sqrt(np.diag(params_covar))

	return (params_opt,params_err)

    def __init__(self, firing_rate, record=False, **kwargs):
        
        if isinstance(firing_rate, str):
            self.firing_rate_string = str(firing_rate)
            self.closure = lambdify(sym_t,symp.parse_expr(self.firing_rate_string, local_dict={'Heaviside':lambda x: Piecewise((0,x<0), (1,x>0),(.5,True))}))
        elif hasattr(firing_rate, "__call__"):
            self.closure = firing_rate
        else:
            self.firing_rate_string = str(firing_rate)
            self.closure = lambdify(sym_t,symp.parse_expr(self.firing_rate_string))

        self.record = record
        self.type = "external"
        
        # Additional metadata:
        self.metadata = kwargs