How to use the symengine.Symbol function in symengine

# -*- coding: utf-8 -*-

"""
Integration test of jitcode_restricted_lyap and jitcode_transversal_lyap by comparing them to each other.
"""

from itertools import combinations
from jitcode import jitcode_restricted_lyap, y, jitcode_transversal_lyap
import numpy as np
from scipy.stats import sem
from symengine import Symbol

a = -0.025794
b =  0.01
c =  0.02
k = Symbol("k")

scenarios = [
	{
	"f":
		[
			y(0)*(a-y(0))*(y(0)-1.0) - y(3) + k*(y(1)-y(0)),
			y(1)*(a-y(1))*(y(1)-1.0) - y(4) + k*(y(2)-y(1)),
			y(2)*(a-y(2))*(y(2)-1.0) - y(5) + k*(y(0)-y(2)),
			b*y(0) - c*y(3),
			b*y(1) - c*y(4),
			b*y(2) - c*y(5),
		],
	"vectors":
		[
			[1.,1.,1.,0.,0.,0.],
			[0.,0.,0.,1.,1.,1.]

"""
Tests whether things works independent of where symbols are imported from.
"""

import jitcode
import jitcode.sympy_symbols
import sympy
import symengine

symengine_manually = [
		symengine.Symbol("t",real=True),
		symengine.Function("y",real=True),
		symengine.cos,
	]

sympy_manually = [
		sympy.Symbol("t",real=True),
		sympy.Function("y",real=True),
		sympy.cos,
	]

jitcode_provisions = [
		jitcode.t,
		jitcode.y,
		symengine.cos,
	]

T = 5

sech = lambda x: 2/(exp(x)+exp(-x))
eps = 1e-5
abs = lambda x: sqrt(x**2+eps**2)

ε = [ 0.03966, 0.03184, 0.02847 ]
ν = [ 1, 2.033, 3.066 ]
μ = [ 0.16115668456085775, 0.14093420256851111, 0.11465065353644151 ]
ybar_0 = 0
τ = 1.7735
ζ = [ 0.017940997406325931, 0.015689701773967984, 0.012763648066925721 ]

anchors_past = Symbol("anchors_past")
difference = Symbol("difference")
factor_μ = Symbol("factor_μ")
factor_ζ = Symbol("factor_ζ")

ydot = [ y(i) for i in range(3,6) ]
y_tot      = sum( current_y(i)                for i in range(  3) )
ydot_tot   = sum( current_y(i)                for i in range(3,6) )
y_past     = sum( past_y (t-τ,i,anchors_past) for i in range(  3) )
ydot_past  = sum( past_y (t-τ,i,anchors_past) for i in range(3,6) )
yddot_past = sum( past_dy(t-τ,i,anchors_past) for i in range(3,6) )

helpers = {
		( anchors_past, anchors(t-τ) ),
		( difference, ybar_0-y_past ),
		( factor_μ, sech(difference)**2 * (yddot_past + 2*ydot_past**2*tanh(difference)) ),
		( factor_ζ, 2*abs(y_tot)*ydot_tot ),
	}

a = 0.165
	b = 0.2
	c = 10.0
	k = 0.01
	
	# get adjacency matrix of a small-world network
	A = small_world_network(
		number_of_nodes = N,
		nearest_neighbours = 20,
		rewiring_probability = 0.1
		)
	
	# generate differential equations
	# -------------------------------
	
	sum_z = symengine.Symbol("sum_z")
	helpers = [( sum_z, sum( y(3*j+2) for j in range(N) ) )]
	
	def f():
		for i in range(N):
			coupling_sum = sum( y(3*j)-y(3*i) for j in range(N) if A[i,j] )
			coupling_term = k * symengine.sin(t) * coupling_sum
			yield -ω[i] * y(3*i+1) - y(3*i+2) + coupling_term
			yield  ω[i] * y(3*i) + a*y(3*i+1)
			coupling_term_2 = k * (sum_z-N*y(3*i+2))
			yield b + y(3*i+2) * (y(3*i) - c) + coupling_term_2
	
	# integrate
	# ---------
	
	initial_state = np.random.random(3*N)

sech = lambda x: 2/(exp(x)+exp(-x))
eps = 1e-5
abs = lambda x: sqrt(x**2+eps**2)

ε = [ 0.03966, 0.03184, 0.02847 ]
ν = [ 1, 2.033, 3.066 ]
μ = [ 0.16115668456085775, 0.14093420256851111, 0.11465065353644151 ]
ybar_0 = 0
τ = 1.7735
ζ = [ 0.017940997406325931, 0.015689701773967984, 0.012763648066925721 ]

anchors_past = Symbol("anchors_past")
difference = Symbol("difference")
factor_μ = Symbol("factor_μ")
factor_ζ = Symbol("factor_ζ")

ydot = [ y(i) for i in range(3,6) ]
y_tot      = sum( current_y(i)                for i in range(  3) )
ydot_tot   = sum( current_y(i)                for i in range(3,6) )
y_past     = sum( past_y (t-τ,i,anchors_past) for i in range(  3) )
ydot_past  = sum( past_y (t-τ,i,anchors_past) for i in range(3,6) )
yddot_past = sum( past_dy(t-τ,i,anchors_past) for i in range(3,6) )

helpers = {
		( anchors_past, anchors(t-τ) ),
		( difference, ybar_0-y_past ),
		( factor_μ, sech(difference)**2 * (yddot_past + 2*ydot_past**2*tanh(difference)) ),
		( factor_ζ, 2*abs(y_tot)*ydot_tot ),
	}

f = { y(i):ydot[i] for i in range(3) }

symbols = (additional_helper(i) for i in count())
				)
			helpers_wc.extend(symengine.sympify(_cse[0]))
			f_sym_wc = symengine.sympify(_cse[1][0])
		
		arguments = [
			("self", "dde_integrator * const"),
			("t", "double const"),
			("y", "double", self.n),
			]
		helper_i = 0
		anchor_i = 0
		converted_helpers = []
		self.past_calls = 0
		self.substitutions = {
				control_par: symengine.Symbol("self->parameter_"+control_par.name)
				for control_par in self.control_pars
			}
		
		def finalise(expression):
			expression = expression.subs(self.substitutions)
			self.past_calls += count_calls(expression,anchors)
			return expression
		
		if helpers_wc:
			get_helper = symengine.Function("get_f_helper")
			set_helper = symengine.Function("set_f_helper")
			
			get_anchor = symengine.Function("get_f_anchor_helper")
			set_anchor = symengine.Function("set_f_anchor_helper")
			
			for helper in helpers_wc:

def _prepare_lambdas(self):
		if self.callback_functions:
			raise NotImplementedError("Callbacks do not work with lambdification. You must use the C backend.")
		
		if not hasattr(self,"_lambda_subs") or not hasattr(self,"_lambda_args"):
			if self.helpers:
				warn("Lambdification handles helpers by plugging them in. This may be very inefficient")
			
			self._lambda_subs = list(reversed(self.helpers))
			self._lambda_args = [t]
			for i in range(self.n):
				symbol = symengine.Symbol("dummy_argument_%i"%i)
				self._lambda_args.append(symbol)
				self._lambda_subs.append((y(i),symbol))
			self._lambda_args.extend(self.control_pars)

from symengine import var, sympify, function_symbol, Symbol
import sympy
s = open("expr.txt").read()
print "Converting to SymPy..."
e = sympy.sympify(s)
print "Converting to SymEngine..."
ce = sympify(e)
print "    Done."
print "SymPy subs:"
t1 = clock()
f = e.subs(sympy.Function("q5")(sympy.Symbol("t")), sympy.Symbol("sq5"))
t2 = clock()
print "Total time:", t2-t1, "s"
print "SymEngine subs:"
t1 = clock()
cf = ce.subs(function_symbol("q5", Symbol("t")), Symbol("sq5"))
t2 = clock()
print "Total time:", t2-t1, "s"

print "SymPy diff:"
t1 = clock()
g = f.diff(sympy.Symbol("sq5"))
t2 = clock()
print "Total time:", t2-t1, "s"
print "SymEngine diff:"
t1 = clock()
cg = cf.diff(Symbol("sq5"))
t2 = clock()
print "Total time:", t2-t1, "s"

import numpy as np

sech = lambda x: 2/(exp(x)+exp(-x))
eps = 1e-5
abs = lambda x: sqrt(x**2+eps**2)

ε = [ 0.03966, 0.03184, 0.02847 ]
ν = [ 1, 2.033, 3.066 ]
μ = [ 0.16115668456085775, 0.14093420256851111, 0.11465065353644151 ]
ybar_0 = 0
τ = 1.7735
ζ = [ 0.017940997406325931, 0.015689701773967984, 0.012763648066925721 ]

anchors_past = Symbol("anchors_past")
difference = Symbol("difference")
factor_μ = Symbol("factor_μ")
factor_ζ = Symbol("factor_ζ")

ydot = [ y(i) for i in range(3,6) ]
y_tot      = sum( current_y(i)                for i in range(  3) )
ydot_tot   = sum( current_y(i)                for i in range(3,6) )
y_past     = sum( past_y (t-τ,i,anchors_past) for i in range(  3) )
ydot_past  = sum( past_y (t-τ,i,anchors_past) for i in range(3,6) )
yddot_past = sum( past_dy(t-τ,i,anchors_past) for i in range(3,6) )

helpers = {
		( anchors_past, anchors(t-τ) ),
		( difference, ybar_0-y_past ),
		( factor_μ, sech(difference)**2 * (yddot_past + 2*ydot_past**2*tanh(difference)) ),
		( factor_ζ, 2*abs(y_tot)*ydot_tot ),
	}

from jitcdde.past import Past, Anchor
import jitcdde._python_core as python_core
from jitcxde_common import jitcxde, checker
from jitcxde_common.helpers import sort_helpers, sympify_helpers, find_dependent_helpers
from jitcxde_common.symbolic import collect_arguments, count_calls, replace_function
from jitcxde_common.transversal import GroupHandler
import chspy

_default_min_step = 1e-10

#sigmoid = lambda x: 1/(1+np.exp(-x))
#sigmoid = lambda x: 1 if x>0 else 0
sigmoid = lambda x: (np.tanh(x)+1)/2

#: the symbol for time for defining the differential equation. You may just as well define the an analogous symbol directly with SymEngine or SymPy, but using this function is the best way to get the most of future versions of JiTCDDE, in particular avoiding incompatibilities. You can import a SymPy variant from the submodule `sympy_symbols` instead (see `SymPy vs. SymEngine`_ for details).
t = symengine.Symbol("t", real=True)

def y(index,time=t):
	"""
	the function representing the DDE’s past and present states used for defining the differential equation. The first integer argument denotes the component. The second, optional argument is a symbolic expression denoting the time. This automatically expands to using `current_y`, `past_y`, and `anchors`; so do not be surprised when you look at the output and it is different than what you entered or expected. You can import a SymPy variant from the submodule `sympy_symbols` instead (see `SymPy vs. SymEngine`_ for details).
	"""
	if time == t:
		return current_y(index)
	else:
		return past_y(time, index, anchors(time))

def dy(index,time):
	"""
	This is the function representing the DDE’s past derivative used for defining the differential equation. The first integer argument denotes the component. The second, argument denotes the time. If you use this, you get a neutral DDE which may make addressing initial discontinuities more difficult. Do not use this to get the current derivative. Instead compute it using your dynamical equations. This will not work with tangential Lyapunov exponents.

	**This feature is experimental.**