Maths with Sympy

docs/_data/navigation.yaml

@@ -84,6 +84,8 @@ pages:
     link: /python/jupyter-notebooks
   - name: Encrypting Your Input
     link: /python/encrypting
+  - name: Maths with SymPy
+    link: /python/sympy
 - name: "2015"
   link: /2015/
   mainnav: True

docs/python/sympy.md

@@ -1,64 +1,240 @@
 ---
-title: SymPy
+title: Maths with SymPy
 main_img:
-  name: Jupyter Notebook
-  link: /assets/images/aoc2017-jn.png
+  name: SymPy
+  link: /assets/images/sympylogo.png
 tags: 
-  - name: Jupyter
-    link: https://jupyter.org/
-  - name: Jupyter Project Documentation
-    link: https://docs.jupyter.org/en/latest/index.html
-  - name: Jupyter Lab Documentation
-    link: https://jupyterlab.readthedocs.io/en/stable/
-  - name: Installing Jupyter
-    link: https://jupyterlab.readthedocs.io/en/stable/getting_started/installation.html
-  - name: Notebook Tips and Tricks
-    link: https://www.dataquest.io/blog/jupyter-notebook-tips-tricks-shortcuts/
-  - name: Anaconda
-    link: https://learning.anaconda.cloud/get-started-with-anaconda?source=install
-  - name: SciPy Docker Compose
-    link: https://gist.github.com/derailed-dash/b1d9eb511e336ba837da234518c09842
-  - name: Google Colab
-    link: https://colab.research.google.com/
-  - name: Anaconda Cloud
-    link: https://nb.anaconda.cloud/
+  - name: SymPy
+    link: https://www.sympy.org/en/index.html
+  - name: SymPy Tutorial
+    link: https://www.tutorialspoint.com/sympy/index.htm
 ---
+<script
+  src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"
+  type="text/javascript">
+</script> 
 ## Page Contents
 
 Work in progress...
 
-- [What are Jupyter Notebooks?](#what-are-jupyter-notebooks)
-- [A Few Benefits of Notebooks](#a-few-benefits-of-notebooks)
-- [Ideal Scenarios for Using Notebooks](#ideal-scenarios-for-using-notebooks)
+- [What is SymPy?](#what-is-sympy)
+- [Basic Examples](#basic-examples)
+  - [Simplifying and Expanding](#simplifying-and-expanding)
+  - [Assigning Variables and Evaluating](#assigning-variables-and-evaluating)
+  - [Differentiation and Integration](#differentiation-and-integration)
+  - [Solving Quadratics](#solving-quadratics)
+  - [Including Only Real Solutions](#)
+  - [Evaluating to Obtain Numeric Values](#evaluating-to-obtain-numeric-values)
+- [AoC Examples](#aoc-examples)
 
-## What are Jupyter Notebooks?
+## What is SymPy?
 
-Integrated Development Environments (IDEs) have been the go-to choice for developers seeking a dedicated space to write, test, and debug their code. However, as data science and machine learning began gaining traction, there was increasing need for a more interactive and data-centric environment. Enter **Jupyter Notebooks**: a web-based application that has redefined the way we interact with code and data, the way we document our code, and the way we share it.
+SymPy is so cool!  I came across when I was trying to solve some equations as part of 2023 AoC puzzles. It is computer alegebra systems (CAS): a library that allows us to perfrom algebra and mathematical computations. 
 
-Unlike traditional IDEs, **Jupyter Notebooks allow for code, data, and multimedia to coexist in a shared space**. With its cell-based structure, users can write and execute code in small chunks, making it easier to test ideas and see results in real-time. This attribute alone sets it apart from the linear approach of traditional coding environments, making Jupyter Notebooks a beloved tool among data scientists and analysts.
+You can use it for things like:
 
-## A Few Benefits of Notebooks
+- Simplifying, rearranging, factoring and expanding algebraic expressions and equations.
+- Solving equations of various types, including linear, polynomial, and differential.
+- Performing calculus, including differentiation and integration.
+- Plotting graphs.
+- Generating LaTeX output.
 
-- **Interactive computing:** Real-time feedback on code execution, facilitating a more iterative and exploratory approach to problem-solving. Code is organised into _cells_. Cells can be executed and re-executed.  The output of those cells (such as the values assigned to variables) are available to all subsequent cells in the Notebook.  This is useful, since we can run a cell, and then experiment with changes to the cells that come _after_ the cell.
-- **Documentation and sharing:** You have the ability to intertwine code, text, and multimedia. You can place markdown or HTML content alongside your code. This greatly improves your ability to document your solutions and share insights.  And you can structure your Notebook using chapters.
-- **Images and visualisation**: You can include images and visualisations, equations, animations and video, alongside your code.
-- **Reproducibility:** Encourages reproducible research by packaging code, data, and narrative into a single document.
-- **Extension ecosystem:** A vast array of extensions and libraries are available, enhancing functionality and integration with other tools.
-- **Run shell commands:** You can run shell commands from within the Notebook! This is useful for running file system commands, or - for example - to pre-install packages with `pip`.
-- **Easy to share:** Since all the code, documentation and output is in a single file, this makes sharing (and collaborating) incredibly easy.
-- **Portable:** Jupyter Notebooks are inherently portable, and can be run on a wide variety of platforms.
+## Basic Examples
 
-Here's an example of a Notebook, with chapter structure, and with an introduction written in markdown:
+### Simplifying and Expanding
 
-![Notebook Example](/assets/images/notebook-md.png)
+We need to define variables before we can use them in SymPy. We do this by creating `symbols`.
 
-And here's an example, where we're dynamically generating an image with code, and showing it after the cell:
+One interesting observation here is that when you evaluate a SymPy expression a Jupyter notebook, it is rendered in _pretty-printed_ (mathematical) format. But when we print using `print()`, it is printed in conventional Python style.
 
-![Markdown, Cells, and Graphical Output](/assets/images/notebook-generating-image.png)
+If you want to explicitly pretty-print from a Python notebook, you can use `IPython.display`.
 
-## Ideal Scenarios for Using Notebooks
+Alternatively, you can optionally run `init_printing()` to enable the best printer for your environment. Or run `init_session()` to automatically import everything in SymPy, create some common Symbols, setup plotting, and run `init_printing()`.
 
-- **Data Analysis and Visualization:** Whether it's cleaning data, performing statistical analysis, or creating visualizations, Jupyter Notebooks provide an interactive environment that fosters insight discovery.
-- **Machine Learning and Model Development:** Building and testing models become streamlined with the ability to execute code incrementally and visualize results on the fly.
-- **Educational Purposes:** An excellent tool for teaching coding and data science concepts, given its interactive nature and the ability to annotate code with rich text.
-- **Collaborative Research:** Facilitates collaborative efforts by allowing multiple users to interact with the notebook and share their findings seamlessly.
+```python
+import sympy
+from sympy import latex
+from IPython.display import display, Markdown
+
+sympy.init_session() # set up default behavour and printing - not essential
+
+a, b, x, y = sympy.symbols("a b x y") # define symbols (variables)
+a = (x + 1)**2 # set a to be an expresssion
+b = x**2 + 2*x + 1 # set b to be an equivalent expression
+display(Markdown(f"$a = {latex(a)}$")) # pretty-print a
+print(f"{latex('a = ')}{latex(a)}") # If we just want the raw latex
+
+display(Markdown(f"$b = {latex(b)}$")) # pretty-print b
+
+# Expand a...
+display(Markdown("Expanding $a$..."))
+expanded = sympy.expand(a)
+display(Markdown(f"$a = {latex(expanded)}$"))
+
+# Show a-b without simplifying...
+a_minus_b = a-b
+display(Markdown(f"$a-b = {latex(a_minus_b)}$"))
+
+# Simplify...
+display(Markdown("Simplifying..."))
+display(Markdown(f"$a-b$ $= {latex(sympy.simplify(a_minus_b))}$"))
+```
+
+![Simplifying and expanding](/assets/images/sympy_simplify_and_expand.png)
+
+We can also test if two expressions are the same, e.g.
+
+```python
+a.equals(b)
+```
+
+### Assigning Variables and Evaluating
+
+You can't simply assign values to the Python variable. If you did so, they would cease to be a symbol. Instead, you need to use `subs()` to associate a value with a SymPy symbol.
+
+```python
+x, y = sympy.symbols("x y")
+expr = x+y
+result = expr.subs({x: 2, y: 5})
+print(f"{expr=}, {result=}")
+```
+
+![SymPy variables](/assets/images/sympy_variables.png){:style="width:480px"}
+
+### Differentiation and Integration
+
+```python
+import sympy
+
+x = sympy.symbols("x")
+expr = x**2
+
+deriv = sympy.diff(expr)
+print(deriv)
+
+integral = sympy.integrate(deriv)
+print(integral)
+```
+
+![SymPy basics](/assets/images/sympy_basics.png)
+
+### Solving Quadratics
+
+Here, I use `solve()` to determine the values of $x$. I've specified `dict=True` such that $x$ can be retrieved by passing its symbol as a dictionary key to the returned `solutions` dictionary.
+
+```python
+expr = x**2 + 3*x - 10
+solutions = sympy.solve(expr, x, dict=True)
+solutions
+
+for solution in solutions:
+    print(f"x={solution[x]}")
+```
+
+![SymPy quadratic](/assets/images/sympy_quadratic.png){:style="width:480px"}
+
+Note that the expressions still need to be written in valid Python. So we can write `3*x`, but we can't write `3x`.
+
+### Including Only Real Solutions
+
+What if you want to ignore _complex_ solutions and only include _real_ solutions? In this case, you can use `solveset()`, and specify the domain as `domain=sympy.S.Reals`. (Whereas `domain=S.Complexes` is the default.)
+
+```python
+expressions = []
+expressions.append(x**2 - 9) # there are only real solutions
+expressions.append(x**2 + 9) # there are only complex solutions
+
+for expr in expressions:
+    display(expr)
+    solutions = sympy.solveset(expr, x)
+    print(f"There are {len(solutions)} solutions for x...")
+    display(solutions)
+
+display(expressions[-1])
+solutions = sympy.solveset(expressions[-1], x, domain=sympy.S.Reals)
+print(f"There are {len(solutions)} real solutions for x.")
+```
+
+### Evaluating to Obtain Numeric Values
+
+We can evaluate to obtain the float value of a symbol, and display to an arbitrary level of precision: 
+
+```python
+expr = sympy.sqrt(8)
+display(expr)
+display(expr.evalf(4)) # to 4 digits of precisions
+```
+
+![SymPy Numeric Evaluation](/assets/images/sympy_numeric_eval.png){:style="width:480px"}
+
+## AoC Examples
+
+### Boat Racing Quadratic
+
+In this example taken from [2023 Day 6](https://colab.research.google.com/github/derailed-dash/Advent-of-Code/blob/master/src/AoC_2023/Dazbo's_Advent_of_Code_2023.ipynb#scrollTo=2jT8WTv-Cfgo){:target="_blank"}, I'm solving a quadratic:
+
+$$
+h^2 - th + d = 0
+$$
+
+```python
+def solve_part2_sympy_quadratic(data):
+    """ h^2 - th + d = 0 """
+    race_duration = int("".join(x for x in data[0].split(":")[1].split()))
+    distance = int("".join(x for x in data[1].split(":")[1].split()))
+    logger.debug(f"{race_duration=}, {distance=}")  
+
+    # solve using quadratic with SymPy       
+    h = sympy.symbols("h", real=True)
+    equation = sympy.Eq(h**2 - race_duration*h + distance, 0)
+    solutions = sympy.solve(equation, h, dict=True)
+    answers = [solution[h].evalf() for solution in solutions] # there should be two
+```
+
+### Trajectory Intersection
+
+In this example taken from [2023 Day 24](https://colab.research.google.com/github/derailed-dash/Advent-of-Code/blob/master/src/AoC_2023/Dazbo's_Advent_of_Code_2023.ipynb#scrollTo=9agx4CcYCfhR){:target="_blank"}, I have a series of equations that are necessary to find several unknown variables. Specifically:
+
+$$
+\begin{align}
+t &= \frac{x_{r} - x_{h}}{v_{x_{h}} - v_{x_{r}}} = \frac{y_{r} - y_{h}}{v_{y_{h}} - v_{y_{r}}} = \frac{z_{r} - z_{h}}{v_{z_{h}} - v_{z_{r}}} \\
+\notag \\
+(x_{r} - x_{h})(v_{y_{h}} - v_{y_{r}}) &= (y_{r} - y_{h})(v_{x_{h}} - v_{x_{r}}) \\
+(y_{r} - y_{h})(v_{z_{h}} - v_{z_{r}}) &= (z_{r} - z_{h})(v_{y_{h}} - v_{y_{r}}) \\
+(z_{r} - z_{h})(v_{x_{h}} - v_{x_{r}}) &= (x_{r} - x_{h})(v_{z_{h}} - v_{z_{r}}) \\
+\end{align}
+$$
+
+And here's the code:
+
+```python
+def solve_part2(data: list[str]):
+    """ 
+    Determine the sum of the rock's (x,y,z) coordinate at t=0, for a rock that will hit every hailstone
+    in our input data. The rock has constant velocity and is not affected by collisions.
+    """
+    stones = parse_stones(data)
+    logger.debug(f"We have {len(stones)} stones.")
+
+    # define SymPy rock symbols - these are our unknowns representing:
+    # initial rock location (xr, yr, zr)
+    # rock velocity (vxr, vyr, vzr)
+    xr, yr, zr, vxr, vyr, vzr = sympy.symbols("xr yr zr vxr vyr vzr")
+
+    equations = [] # we assemble a set of equations that must be true
+    for stone in stones[:10]: # we don't need ALL the stones to find a solution. We need just enough.
+        x, y, z = stone.posn
+        vx, vy, vz = stone.velocity
+        equations.append(sympy.Eq((xr-x)*(vy-vyr), (yr-y)*(vx-vxr)))
+        equations.append(sympy.Eq((yr-y)*(vz-vzr), (zr-z)*(vy-vyr)))
+
+    try:
+        solutions = sympy.solve(equations)[0] # SymPy does the hard work
+    except sympy.core.sympify.SympifyError as e:
+        logger.error(f"Could not find a solution: {e}")
+        return None
+
+    logger.debug(solutions)
+    x, y, z = solutions[xr], solutions[yr], solutions[zr]
+    logger.info(f"{x=},{y=},{z=}")
+```