linsolve_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:48:15.534405Z",
     "start_time": "2020-06-13T00:48:15.253820Z"
    }
   },
   "outputs": [],
   "source": [
    "import ast, linsolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# `linsolve` Quickstart Guide\n",
    "\n",
    "`linsolve` is a module providing high-level tools for (linearizing and) performing a least-squares optimization on a system of equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:29:02.687347Z",
     "start_time": "2020-06-13T00:29:02.245024Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import linsolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "## Solving Linear Least-Squares Systems of Equations with `LinearSolver`\n",
    "\n",
    "In `linsolve`, data is expressed as python dictionary object where equations are passed in as a dictionary where which maps strings containing all unknowns (the RHS of the system) to numbers (the LHS of the system). The strings are parsed according to python syntax and all variables are automatically identified (as long as they are valid python variable names) and solved for.\n",
    "\n",
    "For background, see <https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)>\n",
    "\n",
    "As an example, here is how you would the following system\n",
    "\n",
    "$3 x + 4 y = 11 \\\\ -x -3y = -7$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:29:02.782442Z",
     "start_time": "2020-06-13T00:29:02.688957Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'x': 1.0, 'y': 2.0}\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import linsolve\n",
    "\n",
    "#Set up dictionary of equations\n",
    "data = {'3*x+4*y': 11.0, '-1*x-3*y': -7.0}\n",
    "\n",
    "#Build solver\n",
    "ls = linsolve.LinearSolver(data)\n",
    "\n",
    "#Execute solver\n",
    "sol = ls.solve()\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "### Array Systems\n",
    "\n",
    "Often, it its useful to solve the same system of equations for many different \"measurements\" on the LHS of the system of equations. In this case, the dictionary of equations can map strings to $N$-dimensional numpy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.891288Z",
     "start_time": "2020-06-13T00:06:10.548Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "data = {'3*x+4*y': np.array([11.0,22.0]), '-1*x-3*y': np.array([-7.0,-14.0])}\n",
    "ls = linsolve.LinearSolver(data)\n",
    "sol = ls.solve()\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "### Constants\n",
    "\n",
    "While all terms that evaluate to valid python variables in the dictionary keys are assumed to be variables, it is possible to designate some of them as constants. They are passed `LinearSolver` via its `kwargs`. Constants can be expressed either as single numbers or as arrays that can broadcast according to `numpy` rules to the shape of the data, allowing the equations to change between solves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.891963Z",
     "start_time": "2020-06-13T00:06:10.551Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "data = {'3*x+a*y': np.array([11.0,22.0]), 'b*x-3*y': np.array([-7.0,-14.0])}\n",
    "ls = linsolve.LinearSolver(data, a=4, b=np.array([-1,-2]))\n",
    "sol = ls.solve()\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "Any variable names that are not assigned values are assumed to be parameters in need of solutions.  However, `linsolve` being a linear solver, if more than one variable per term is unknown (i.e. if the expression is non-linear), `linsolve` will raise an `AssertionError`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.892560Z",
     "start_time": "2020-06-13T00:06:10.552Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    ls = linsolve.LinearSolver({'a*x':6., 'b*x':8.}) # raises an AssertionError because 'a' and 'b' are not defined.\n",
    "except(AssertionError):\n",
    "    print('AssertionError encountered.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "### Complex Conjugation \n",
    "\n",
    "Complex conjugation of variables is accomplished via trailing underscore in the variable name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.893232Z",
     "start_time": "2020-06-13T00:06:10.554Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "data = {'x': 1.0+1.0j, 'y_': 1.0+1.0j}\n",
    "ls = linsolve.LinearSolver(data)\n",
    "sol = ls.solve()\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "### Weights\n",
    "\n",
    "When a system of equations is overdetermined (as is usually the case in data analysis), `linsolve` lets you apply different weights to each equation. The default behavior is to weight all equations equally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.893948Z",
     "start_time": "2020-06-13T00:06:10.557Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "data = {'1*x': 2.0, 'x': 1.0}\n",
    "ls = linsolve.LinearSolver(data)\n",
    "sol = ls.solve()\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "However, you can provide weights through the `wgts` keyword argument. Weights are multiplied by the data, so in standard inverse covariance weighting they can be thought of as $1/\\sigma^2$ where $\\sigma$ is the noise on that equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.894732Z",
     "start_time": "2020-06-13T00:06:10.559Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "data = {'1*x': 2.0, 'x': 1.0}\n",
    "wgts = {'1*x': 1.0, 'x': .5}\n",
    "ls = linsolve.LinearSolver(data, wgts=wgts)\n",
    "sol = ls.solve()\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "This solution, $x = 5/3$ is twice $(1/W^2)$ as close to $x=2$ as $x=1$. This is because it minimizes $\\chi^2$, defined as\n",
    "\n",
    "$\\chi^2 = \\sum_i \\frac{(d_i - m_i)^2}{\\sigma_i^2} = \\sum_i W_i^2 (d_i - m_i)^2$\n",
    "\n",
    "where $i$ indexes over equations, $d_i$ is the data (RHS), $m_i$ is the model (LHS evaluated at `sol`), and $W_i$ are the weights. We can calculate $\\chi^2$ with this handy function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.895376Z",
     "start_time": "2020-06-13T00:06:10.561Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "print(ls.chisq(sol))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "Weighting allows for the implementation of inverse variance weighting, but assumes that the noise in each equation/measurement is independent (i.e. it does not allow for inverse *covariance* weighting)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "## Linearizing Systems of Equations \n",
    "\n",
    "`linsolve` is also able to solve very basic nonlinear systems of equations using either `LogProductSolver` or `LinProductSolver`.\n",
    "\n",
    "For background, see <https://en.wikipedia.org/wiki/Non-linear_least_squares>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "### `LogProductSolver` \n",
    "\n",
    "`LogProductSolver` can tackle a special class of non-linear equations where all expressions are products. For example, we can solve\n",
    "\n",
    "$ab = 2 \\\\ bc = 1 \\\\ ac = 2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.896074Z",
     "start_time": "2020-06-13T00:06:10.562Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "data = {'a*b': 2.0, 'b*c': 1.0, 'a*c': 2.0}\n",
    "lps = linsolve.LogProductSolver(data)\n",
    "sol = lps.solve()\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "This is acheived by reducing the problem to a linear system of equations by taking the logarithm of both sides. In the case of an overdetermined system of equation, this can produce biased results that do not always minimize $\\chi^2$. However, it is still useful as a starting guess for `LinProductSolver`.\n",
    "\n",
    "Just as with `LinearSolver`, `LogProductSolver` can take array systems, constants passed as keyword arguments, complex conjugated variables, and equation weighting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "### `LinProductSolver`\n",
    "\n",
    "`LinProductSolver` is a more general linearizer of equations. It can tackle equations that are sums (or differences) of products of variables (no quotients, powers, exponents, etc.). Parentheses are not allowed; the user is expected to expand manually. For example, `LinProductSolver` can tackle equations of the form:\n",
    "\n",
    "$10ab + 3bc - 2ac = 32 \\\\ 2aa + 3bb - 2cc = -4 \\\\ -aa + 2bc + ac = 14$\n",
    "\n",
    "\n",
    "It uses the [Gauss-Newton algorithm](https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm) to find a solution that minimzes $\\chi^2$ (thought it may not be unique). This is an second-order iterative algorithm whose success depends on a decent initial guess `sol0`, formatted the same way as all `linsolve` solutions. The equation solved at each step is a Taylor expansion about `sol0`, which for the first example equation above gives\n",
    "\n",
    "$10(a_0+\\Delta a)(b_0+\\Delta b) + 3(b_0 + \\Delta b)(c_0 + \\Delta c) - 2(a_0 + \\Delta a)(c_0 + \\Delta c) = 32 \\\\ 10a\\Delta b + 10b\\Delta a + 3b\\Delta c + 3c \\Delta b - 2a\\Delta c - 2c\\Delta a + \\mathcal{O}(\\Delta^2) = 32 - 10a_0b_0 + 3b_0c_0 - 2a_0c_0$ \n",
    "\n",
    "`LinProductSolver` then solves for all the $\\Delta$ terms and updates `sol0` to `sol` accordingly. Just as with `LinearSolver`, `LinProductSolver` can take array systems, constants passed as keyword arguments, complex conjugated variables, and equation weighting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.896715Z",
     "start_time": "2020-06-13T00:06:10.567Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "data = {'10*a*b + 3*b*c - 2*a*c': 32, \n",
    "        '2*a*a +3*b*b -2*c*c': -4,\n",
    "        '-a*a + 2*b*c + a*c': 14}\n",
    "sol0 = {'a': .9, 'b': 2.1, 'c': 3.2}\n",
    "lps = linsolve.LinProductSolver(data, sol0)\n",
    "meta, sol = lps.solve_iteratively()\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "`LinProductSolver.solve_iteratively()` produces both a solution dictionary in the standard `linsolve` format, but also produces `meta`, a metadata dictionary that contains:\n",
    "\n",
    "* `conv_crit`: A convergence criterion (default `1e-10`) below which to stop iterating. Converegence is measured $\\ell^2$-norm of the change in the solution of all the variables since the previous iteration divided by the $\\ell^2$-norm of the solution itself.\n",
    "* `chisq`: $\\chi^2$ as calculated above for the final iteration\n",
    "* `iter`: the number of iterations ran. The iterative solver runs until until either convergence or `maxiter` is hit (default 50)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-13T00:06:10.897422Z",
     "start_time": "2020-06-13T00:06:10.568Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "print(meta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "#### One Caveat About Variable Names in `LinProductSolver`\n",
    "\n",
    "One should avoid having variables that differ from each other by only a starting `d`, e.g. `dram` and `ram`. Internally, `LinProductSolver` performs a Taylor expansion, treating the variables in `sol0` as constants and creating a corresponding set of variables representing differential adjustments that have a `d` prepended."
   ]
  }
 ],
 "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.7.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}