{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Inverse problems"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"css_file = 'https://raw.githubusercontent.com/ngcm/training-public/master/ipython_notebook_styles/ngcmstyle.css'\n",
"HTML(url=css_file)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So far we've looked at a variety of tests applied to *working, correct* code. All these tests have shown that the code is behaving as expected, that are mental model of the problem is being reproduced by the concrete computational model implemented.\n",
"\n",
"However, what if one of the tests failed? How could we go from a failing test to understanding what's wrong with the code? If we want to be formal about it, this falls into the category of [inverse problems](http://en.wikipedia.org/wiki/Inverse_problem), and in general [inverse problems are hard](http://en.wikipedia.org/wiki/Inverse_problem#Mathematical_considerations)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Examples"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take the n-body code again:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy\n",
"from matplotlib import pyplot\n",
"%matplotlib inline\n",
"from matplotlib import rcParams\n",
"rcParams['font.family'] = 'serif'\n",
"rcParams['font.size'] = 16"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# The Computer Language Benchmarks Game\n",
"# http://benchmarksgame.alioth.debian.org/\n",
"#\n",
"# originally by Kevin Carson\n",
"# modified by Tupteq, Fredrik Johansson, and Daniel Nanz\n",
"# modified by Maciej Fijalkowski\n",
"# 2to3\n",
"\n",
"import sys \n",
"\n",
"def combinations(l):\n",
" result = []\n",
" for x in range(len(l) - 1):\n",
" ls = l[x+1:]\n",
" for y in ls:\n",
" result.append((l[x],y))\n",
" return result\n",
"\n",
"PI = 3.14159265358979323\n",
"SOLAR_MASS = 4 * PI * PI\n",
"DAYS_PER_YEAR = 365.24\n",
"\n",
"BODIES = {\n",
" 'sun': ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], SOLAR_MASS),\n",
"\n",
" 'jupiter': ([4.84143144246472090e+00,\n",
" -1.16032004402742839e+00,\n",
" -1.03622044471123109e-01],\n",
" [1.66007664274403694e-03 * DAYS_PER_YEAR,\n",
" 7.69901118419740425e-03 * DAYS_PER_YEAR,\n",
" -6.90460016972063023e-05 * DAYS_PER_YEAR],\n",
" 9.54791938424326609e-04 * SOLAR_MASS),\n",
"\n",
" 'saturn': ([8.34336671824457987e+00,\n",
" 4.12479856412430479e+00,\n",
" -4.03523417114321381e-01],\n",
" [-2.76742510726862411e-03 * DAYS_PER_YEAR,\n",
" 4.99852801234917238e-03 * DAYS_PER_YEAR,\n",
" 2.30417297573763929e-05 * DAYS_PER_YEAR],\n",
" 2.85885980666130812e-04 * SOLAR_MASS),\n",
"\n",
" 'uranus': ([1.28943695621391310e+01,\n",
" -1.51111514016986312e+01,\n",
" -2.23307578892655734e-01],\n",
" [2.96460137564761618e-03 * DAYS_PER_YEAR,\n",
" 2.37847173959480950e-03 * DAYS_PER_YEAR,\n",
" -2.96589568540237556e-05 * DAYS_PER_YEAR],\n",
" 4.36624404335156298e-05 * SOLAR_MASS),\n",
"\n",
" 'neptune': ([1.53796971148509165e+01,\n",
" -2.59193146099879641e+01,\n",
" 1.79258772950371181e-01],\n",
" [2.68067772490389322e-03 * DAYS_PER_YEAR,\n",
" 1.62824170038242295e-03 * DAYS_PER_YEAR,\n",
" -9.51592254519715870e-05 * DAYS_PER_YEAR],\n",
" 5.15138902046611451e-05 * SOLAR_MASS) }\n",
"\n",
"\n",
"SYSTEM = list(BODIES.values())\n",
"PAIRS = combinations(SYSTEM)\n",
"\n",
"\n",
"def advance(dt, n, bodies=SYSTEM, pairs=PAIRS):\n",
"\n",
" for i in range(n):\n",
" for (([x1, y1, z1], v1, m1),\n",
" ([x2, y2, z2], v2, m2)) in pairs:\n",
" dx = x1 - x2\n",
" dy = y1 - y2\n",
" dz = z1 - z2\n",
" mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-1.5))\n",
" b1m = m1 * mag\n",
" b2m = m2 * mag\n",
" v1[0] -= dx * b2m\n",
" v1[1] -= dy * b2m\n",
" v1[2] -= dz * b2m\n",
" v2[0] += dx * b1m\n",
" v2[1] += dy * b1m\n",
" v2[2] += dz * b1m\n",
" for (r, [vx, vy, vz], m) in bodies:\n",
" r[0] += dt * vx\n",
" r[1] += dt * vy\n",
" r[2] += dt * vz"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, let's define a number of `advance` functions with small errors in."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def advance_dx_error(dt, n, bodies=SYSTEM, pairs=PAIRS):\n",
"\n",
" for i in range(n):\n",
" for (([x1, y1, z1], v1, m1),\n",
" ([x2, y2, z2], v2, m2)) in pairs:\n",
" dx = x1 - x2\n",
" dy = y1 - y2\n",
" dz = z1 - z2\n",
" mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-1.5))\n",
" b1m = m1 * mag\n",
" b2m = m2 * mag\n",
" v1[0] -= dx * b2m\n",
" v1[1] -= dx * b2m\n",
" v1[2] -= dx * b2m\n",
" v2[0] += dx * b1m\n",
" v2[1] += dy * b1m\n",
" v2[2] += dz * b1m\n",
" for (r, [vx, vy, vz], m) in bodies:\n",
" r[0] += dt * vx\n",
" r[1] += dt * vy\n",
" r[2] += dt * vz"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def advance_sign_error(dt, n, bodies=SYSTEM, pairs=PAIRS):\n",
"\n",
" for i in range(n):\n",
" for (([x1, y1, z1], v1, m1),\n",
" ([x2, y2, z2], v2, m2)) in pairs:\n",
" dx = x1 - x2\n",
" dy = y1 - y2\n",
" dz = z1 - z2\n",
" mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-1.5))\n",
" b1m = m1 * mag\n",
" b2m = m2 * mag\n",
" v1[0] += dx * b2m\n",
" v1[1] += dy * b2m\n",
" v1[2] += dz * b2m\n",
" v2[0] += dx * b1m\n",
" v2[1] += dy * b1m\n",
" v2[2] += dz * b1m\n",
" for (r, [vx, vy, vz], m) in bodies:\n",
" r[0] += dt * vx\n",
" r[1] += dt * vy\n",
" r[2] += dt * vz"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def advance_power_error(dt, n, bodies=SYSTEM, pairs=PAIRS):\n",
"\n",
" for i in range(n):\n",
" for (([x1, y1, z1], v1, m1),\n",
" ([x2, y2, z2], v2, m2)) in pairs:\n",
" dx = x1 - x2\n",
" dy = y1 - y2\n",
" dz = z1 - z2\n",
" mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-0.5))\n",
" b1m = m1 * mag\n",
" b2m = m2 * mag\n",
" v1[0] -= dx * b2m\n",
" v1[1] -= dy * b2m\n",
" v1[2] -= dz * b2m\n",
" v2[0] += dx * b1m\n",
" v2[1] += dy * b1m\n",
" v2[2] += dz * b1m\n",
" for (r, [vx, vy, vz], m) in bodies:\n",
" r[0] += dt * vx\n",
" r[1] += dt * vy\n",
" r[2] += dt * vz"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Which of our previous tests passes, and which fails? Could we go from a failing test to understanding which parts of the codes are wrong?"
]
}
],
"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.4.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}