{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"import matplotlib_inline\n",
"matplotlib_inline.backend_inline.set_matplotlib_formats('png')\n",
"import seaborn as sns\n",
"sns.set_context(\"paper\")\n",
"sns.set_style(\"ticks\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Homework 1\n",
"\n",
"## References\n",
"\n",
"+ Lectures 1 through 3 (inclusive).\n",
"\n",
"## Instructions\n",
"\n",
"+ Type your name and email in the \"Student details\" section below.\n",
"+ Develop the code and generate the figures you need to solve the problems using this notebook.\n",
"+ For the answers that require a mathematical proof or derivation you should type them using latex. If you have never written latex before and you find it exceedingly difficult, we will likely accept handwritten solutions.\n",
"+ The total homework points are 100. Please note that the problems are not weighed equally.\n",
"\n",
"## Student details\n",
"\n",
"+ **First Name:**\n",
"+ **Last Name:**\n",
"+ **Email:**\n",
"+ **Used generative AI to complete this assignment (Yes/No):**\n",
"+ **Which generative AI tool did you use (if applicable)?:**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 1 - Recursion vs Iteration\n",
"\n",
"This problem adjusted from the [Structure and Interpretation of Computer Programs](https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book.html) book.\n",
"In particular from [this section](https://mitp-content-server.mit.edu/books/content/sectbyfn/books_pres_0/6515/sicp.zip/full-text/book/book-Z-H-11.html#%_sec_1.2.1).\n",
"\n",
"Imagine you are working with a programming language that does not have loops.\n",
"This is how you have to think when writing code in `Jax`.\n",
"Let's say we want to write a function that calculates the factorial of a number:\n",
"\n",
"$$\n",
"n! = n \\times (n-1) \\times (n-2) \\times \\dots \\times 1\n",
"$$\n",
"\n",
"The standard recursive definition of the factorial function is:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def factorial(n):\n",
" if n == 0:\n",
" return 1\n",
" else:\n",
" return n * factorial(n-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is how it can be used:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"factorial(5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's unroll what actually happens behind the scenes:\n",
"\n",
"```python\n",
"factorial(5)\n",
"5 * factorial(4)\n",
"5 * (4 * factorial(3))\n",
"5 * (4 * (3 * factorial(2)))\n",
"5 * (4 * (3 * (2 * factorial(1))))\n",
"5 * (4 * (3 * (2 * 1)))\n",
"5 * (4 * (3 * 2))\n",
"5 * (4 * 6)\n",
"5 * 24\n",
"120\n",
"```\n",
"\n",
"You quickly notice, that the amount of intermediate results that are stored in memory grows exponentially with the input.\n",
"This won't work for large inputs, because you will run out of memory.\n",
"But, there is another way to achieve the same result without exploding memory usage.\n",
"We could start by multiplying 1 by 2, then the result with 3, then the result with 4, and so on.\n",
"So, we keep track of a running product that we update.\n",
"We don't need a loop to do this kind of iteration.\n",
"We can do it with recursion:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def fact_iter(product, counter, max_iter):\n",
" if counter > max_iter:\n",
" return product\n",
" else:\n",
" return fact_iter(counter * product, counter + 1, max_iter)\n",
" \n",
"def good_factorial(n):\n",
" return fact_iter(1, 1, n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Check that this works as before:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"good_factorial(5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is how this unrolls:\n",
"\n",
"```python\n",
"factorial(5)\n",
"fact_iter(1, 1, 5)\n",
"fact_iter(1, 2, 5)\n",
"fact_iter(2, 3, 5)\n",
"fact_iter(6, 4, 5)\n",
"fact_iter(24, 5, 5)\n",
"fact_iter(120, 6, 5)\n",
"120\n",
"```\n",
"\n",
"We say that the second approach is *iterative* and the first approach is *recursive*.\n",
"We want to be writing iterative code, because it is more efficient.\n",
"\n",
"Write iterative code that, given $n$, computes the fibonacci number:\n",
"\n",
"$$\n",
"f_n = f_{n-1} + f_{n-2}\n",
"$$\n",
"\n",
"where $f_0 = 0$ and $f_1 = 1$.\n",
"You should not use a loop!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Your code here - Demonstrate that it works"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Here show how your code works for $n=5$ like I did above with the factorial example.*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 2 - The `foldl` function\n",
"\n",
"The `foldl` function is a higher order function that is used to implement iteration.\n",
"It is defined as follows:\n",
"\n",
"$$\n",
"\\text{foldl}(f, z, [x_1, x_2, \\dots, x_n]) = f(f(\\dots f(f(z, x_1), x_2), \\dots), x_n)\n",
"$$\n",
"\n",
"where $f$ is a function that takes two arguments and $z$ is the initial value.\n",
"In words, `foldl` takes a function $f$, an initial value $z$, and a list $[x_1, x_2, \\dots, x_n]$.\n",
"It then applies $f$ to $z$ and the first element of the list, then applies $f$ to the result of the previous application and the second element of the list, and so on.\n",
"\n",
"Implement `foldl` in `Python`. Pay attention to create an iterative implementation."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Your code here - Demonstrate that it works"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use your `foldl` function to implement the `sum` function and the `product` function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Your code here - Demonstrate that it works"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 3 - No Loops in Jax\n",
"\n",
"Use `Jax`'s [`jax.lax.scan`](https://jax.readthedocs.io/en/latest/_autosummary/jax.lax.scan.html) to implement and `jit` a function that returns the Fibonacci sequence up to a given number.\n",
"Don't bother using integer types, just use `float32` for everything."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Your code here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem 4 - Feigenbaum Map\n",
"\n",
"Consider the function:\n",
"\n",
"$$\n",
"f(x; r) = r x (1 - x)\n",
"$$\n",
"\n",
"where $r$ is a parameter.\n",
"One can define dynamics on the real line by iterating this function:\n",
"\n",
"$$\n",
"x_{n+1} = f(x_n; r)\n",
"$$\n",
"\n",
"where $x_n$ is the state at time $n$.\n",
"\n",
"This map exhibits a [period doubling cascade](https://en.wikipedia.org/wiki/Feigenbaum_constants) as $r$ increases.\n",
"\n",
"Write a function in `jax`, call it `logistic_map`, that takes a lot of $r$'s and $x_0$'s as inputs and returns the first $n$ states of the system.\n",
"You should independently vectorize for the $r$'s and the $x_0$'s.\n",
"And you should `jit`.\n",
"Use `jax.lax.scan` to implement the iteration."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Your code here - Demonstrate that it works"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Test your code here:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x0s = jnp.linspace(0, 1, 100)\n",
"rs = jnp.linspace(0, 4, 1_000)\n",
"n = 10_000\n",
"data = logistic_map(rs, x0s, n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Your shape should be `(1000, 100, 10000)`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Discard all but the last iteration:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data = data[:, :, -1:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Make the famous period doubling plot. The plot will take a while and it will take a lot of memory. I suggest you restart your kernel before moving to the next problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(rs,\n",
" data.reshape(data.shape[0], data.shape[1] * data.shape[2]).T,\n",
" '.k',\n",
" ms=0.1,\n",
" alpha=0.5\n",
");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem 5 - Analysis of Nonlinear Dynamical System\n",
"\n",
"Consider the dynamical system:\n",
"\n",
"$$\n",
"\\dot{x_1} = \\mu x_1 + x_2 - x_1^2,\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"\\dot{x_2} = -x_1 + \\mu x_2 + 2 x_1^2.\n",
"$$\n",
"\n",
"Use the random initial conditions:\n",
"\n",
"$$\n",
"x(0) \\sim N\\left(\n",
" \\begin{pmatrix}\n",
" 0 \\\\\n",
" 0\n",
" \\end{pmatrix},\n",
" \\begin{pmatrix}\n",
" \\sigma^2 & 0 \\\\\n",
" 0 & \\sigma^2\n",
" \\end{pmatrix}\n",
"\\right)\n",
"$$\n",
"\n",
"First, write code that solves the differential equation given the initial and the parameter $\\mu$.\n",
"Make sure your code is vectorized with respect to the initial conditions and that it can be `jit`ed."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Your code and evidence that it works here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"+ Use first order sensitivity analysis to compute the mean and covariance matrix of the solution for the time interval $t \\in [0, 10]$.\n",
"+ Implement a simple Monte Carlo procedure to compare the results of the sensitivity analysis.\n",
"+ Do it for three different values of $\\mu$, $\\mu=0, 0.01$, and $0.066$.\n",
"+ Use $\\sigma=0.01$.\n",
"+ Plot the mean in the $x_1 - x_2$ plane for each value of $\\mu$. Compare local sensitivity analysis to Monte Carlo.\n",
"+ Plot the standard deviation of $x_1$ and $x_2$ as a function of time for each value of $\\mu$. Compare local sensitivity analysis to Monte Carlo."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Your code here - You will need several blocks and discussion"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}