#!/usr/bin/env python
# coding: utf-8

# # Symbolic Computation
# Symbolic computation deals with symbols, representing them exactly, instead of numerical approximations (floating point). 
# 
# We will start with the following [borrowed](https://docs.sympy.org/latest/tutorial/intro.html) tutorial to introduce the concepts of SymPy. Devito uses SymPy heavily and builds upon it in its DSL. 

# In[2]:


import math

math.sqrt(3)


# In[3]:


math.sqrt(8)


# $\sqrt(8) = 2\sqrt(2)$, but it's hard to see that here

# In[4]:


import sympy
sympy.sqrt(3)


# SymPy can even simplify symbolic computations

# In[5]:


sympy.sqrt(8)


# In[6]:


from sympy import symbols
x, y = symbols('x y')
expr = x + 2*y
expr


# Note that simply adding two symbols creates an expression. Now let's play around with it. 

# In[7]:


expr + 1


# In[8]:


expr - x


# Note that `expr - x` was not `x + 2y -x`

# In[9]:


x*expr


# In[10]:


from sympy import expand, factor
expanded_expr = expand(x*expr)
expanded_expr


# In[11]:


factor(expanded_expr)


# In[12]:


from sympy import diff, sin, exp

diff(sin(x)*exp(x), x)


# In[13]:


from sympy import limit

limit(sin(x)/x, x, 0)


# ### Exercise
# 
# Solve $x^2 - 2 = 0$ using sympy.solve

# In[13]:


# Type solution here
from sympy import solve


# ## Pretty printing

# In[14]:


from sympy import init_printing, Integral, sqrt

init_printing(use_latex='mathjax')


# In[15]:


Integral(sqrt(1/x), x)


# In[16]:


from sympy import latex

latex(Integral(sqrt(1/x), x))


# More symbols.
# Exercise: fix the following piece of code

# In[18]:


# NBVAL_SKIP
# The following piece of code is supposed to fail as it is
# The exercise is to fix the code

expr2 = x + 2*y +3*z


# ### Exercise 
# 
# Solve $x + 2*y + 3*z$ for $x$

# In[19]:


# Solution here
from sympy import solve


# Difference between symbol name and python variable name

# In[20]:


x, y = symbols("y z")


# In[21]:


x


# In[22]:


y


# In[23]:


# NBVAL_SKIP
# The following code will error until the code in cell 16 above is
# fixed
z


# Symbol names can be more than one character long

# In[24]:


crazy = symbols('unrelated')

crazy + 1


# In[25]:


x = symbols("x")
expr = x + 1
x = 2


# What happens when I print expr now? Does it print 3?

# In[26]:


print(expr)


# How do we get 3?

# In[27]:


x = symbols("x")
expr = x + 1
expr.subs(x, 2)


# ## Equalities

# In[28]:


x + 1 == 4


# In[29]:


from sympy import Eq

Eq(x + 1, 4)


# Suppose we want to ask whether $(x + 1)^2 = x^2 + 2x + 1$

# In[30]:


(x + 1)**2 == x**2 + 2*x + 1


# In[31]:


from sympy import simplify

a = (x + 1)**2
b = x**2 + 2*x + 1

simplify(a-b)


# ### Exercise 
# Write a function that takes two expressions as input, and returns a tuple of two booleans. The first if they are equal symbolically, and the second if they are equal mathematically.

# ## More operations

# In[32]:


z = symbols("z")
expr = x**3 + 4*x*y - z
expr.subs([(x, 2), (y, 4), (z, 0)])


# In[33]:


from sympy import sympify

str_expr = "x**2 + 3*x - 1/2"
expr = sympify(str_expr)
expr


# In[34]:


expr.subs(x, 2)


# In[35]:


expr = sqrt(8)


# In[36]:


expr


# In[37]:


expr.evalf()


# In[38]:


from sympy import pi

pi.evalf(100)


# In[39]:


from sympy import cos

expr = cos(2*x)
expr.evalf(subs={x: 2.4})


# ### Exercise
# 

# In[41]:


from IPython.core.display import Image 
Image(filename='figures/comic.png') 


# Write a function that takes a symbolic expression (like pi), and determines the first place where 789 appears.
# Tip: Use the string representation of the number. Python starts counting at 0, but the decimal point offsets this

# ## Solving an ODE

# In[42]:


from sympy import Function

f, g = symbols('f g', cls=Function)
f(x)


# In[43]:


f(x).diff()


# In[44]:


diffeq = Eq(f(x).diff(x, x) - 2*f(x).diff(x) + f(x), sin(x))
diffeq


# In[45]:


from sympy import dsolve

dsolve(diffeq, f(x))


# ## Finite Differences

# In[48]:


f = Function('f')
dfdx = f(x).diff(x)
dfdx.as_finite_difference()


# In[49]:


from sympy import Symbol

d2fdx2 = f(x).diff(x, 2)
h = Symbol('h')
d2fdx2.as_finite_difference(h)


# Now that we have seen some relevant features of vanilla SymPy, let's move on to Devito, which could be seen as SymPy finite differences on steroids!