{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5fa5ca58-b746-4b6f-8ac3-8fba4ac2380c",
   "metadata": {},
   "source": [
    "# Linear Algebra\n",
    "\n",
    "## Basics\n",
    "\n",
    "Scipp supports basic linear algebra operations in 3 dimensions using vectors of length 3 and 3x3 matrices (to define rotations).\n",
    "Variables containing vectors or matrices are created using\n",
    "[scipp.vectors](../generated/functions/scipp.vectors.rst#scipp-vectors) (for an array of vectors)\n",
    "[scipp.vector](../generated/functions/scipp.vector.rst#scipp-vector) (for a single vector)\n",
    "[scipp.matrices](../generated/functions/scipp.matrices.rst#scipp-matrices) (for and array of matrices)\n",
    "[scipp.matrix](../generated/functions/scipp.matrix.rst#scipp-matrix) (for a single matrix).\n",
    "Let us consider an example, creating a 1-D variable containg vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb04ef53-39ef-4320-b7e9-2a7fdc61026a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipp as sc\n",
    "import numpy as np\n",
    "\n",
    "vecs = sc.vectors(dims=['x'], unit='m', values=np.arange(2 * 3).reshape(2, 3))\n",
    "vecs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37049758-fb73-4d18-aadb-64ea3b2ca40b",
   "metadata": {},
   "source": [
    "While `vecs` has only a single dimension `x`, the `values` property exposes the underlying vector value as a 2-D array of `float64`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f73f4960-ec75-4576-aa9f-f2b077b9809a",
   "metadata": {},
   "outputs": [],
   "source": [
    "vecs.values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9bc14bc-d100-4bfb-bcfd-fd89ebd743d1",
   "metadata": {},
   "source": [
    "Access to individual vector components `x`, `y`, and `z` is provided using the special `fields` property:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30cb2021-e3e8-4953-9086-e02d70d00f8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "vecs.fields.x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a0b726e-9d83-4474-a4d4-39903b4ae8f1",
   "metadata": {},
   "source": [
    "These properties return a view to the underlying data and thus support setting element values as well as in-place modification:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03017333-7ba0-49c3-88c6-4d4a8bc98381",
   "metadata": {},
   "outputs": [],
   "source": [
    "vecs.fields.x = 123.0 * sc.units.m\n",
    "vecs.fields.y += 0.123 * sc.units.m\n",
    "vecs.values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60e1b618-ec63-4da2-8374-a4655fdb0a66",
   "metadata": {},
   "source": [
    "This works in an identical manner for matrices with the properties `xx`, `xy`, `xz`, `yx`, `yy`, `yz`, `zx`, `zy`, and `zz`.\n",
    "The `values` property allows for use of, e.g., numpy functionality.\n",
    "We may, e.g., compute the inverse using `numpy.linalg.inv`, which is not supported in scipp yet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9907b95e-3192-4931-9b25-62f8039535e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)  # matrices are not singular for this seed\n",
    "mats = sc.matrices(dims=['x'], values=np.random.rand(2, 3, 3))\n",
    "inv = sc.matrices(dims=['x'], values=np.linalg.inv(mats.values))\n",
    "(mats * inv).values.round(decimals=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd597add-0ecd-47a2-a38c-b28b60a28c09",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "We create some random data and use positions to create a plot with `projection='3d'` to illustrate some concepts.\n",
    "Image this as an image sensor with 12x12 pixels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bec94cc4-37b1-4f26-abbc-b4fbd574bbdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "nx = 12\n",
    "ny = 12\n",
    "x = sc.linspace(dim='x', start=-0.1, stop=0.1, num=nx, unit='m')\n",
    "y = sc.linspace(dim='y', start=-0.1, stop=0.1, num=ny, unit='m')\n",
    "sensor = sc.DataArray(\n",
    "    data=sc.array(dims=['x','y'], values=np.random.rand(nx, ny)),\n",
    "    coords={'position':sc.geometry.position(x, y, 0.0*sc.units.m)}\n",
    ")\n",
    "sc.plot(sensor, projection=\"3d\", positions=\"position\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d852b52-613c-4168-88f2-1de4340c77b8",
   "metadata": {},
   "source": [
    "We can use the vector-element access properties in combination with slicing to offset some of the pixels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "207e32b8-3e7f-4a79-8821-e9bdefe48e98",
   "metadata": {},
   "outputs": [],
   "source": [
    "sensor.coords['position']['x',5:].fields.x += 0.1 * sc.units.m\n",
    "sensor.coords['position']['y',5:].fields.y += 0.05 * sc.units.m\n",
    "\n",
    "sc.plot(sensor, projection=\"3d\", positions=\"position\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbf54580-29d3-40bb-96bd-a93f3ff4c26f",
   "metadata": {},
   "source": [
    "We use `sc.matrix` to create a rotation matrix (in this case to rotate by 30 deg around the `y` axis) and apply it to the positions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83fe9f32-26cc-43d0-b6fa-ae761f09b221",
   "metadata": {},
   "outputs": [],
   "source": [
    "rotation = sc.matrix(value=[[ 0.8660254, 0.0000000, 0.5000000],\n",
    "                            [ 0.0000000, 1.0000000, 0.0000000],\n",
    "                            [-0.5000000, 0.0000000, 0.8660254]])\n",
    "\n",
    "sensor.coords['position']['x',5:] = rotation * sensor.coords['position']['x',5:]\n",
    "sc.plot(sensor, projection=\"3d\", positions=\"position\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b6b58d4-6b8d-4d12-a885-41eb0cec192d",
   "metadata": {},
   "source": [
    "More conveniently, we can create our rotations directly from a rotation vector.\n",
    "Scipp provides `from_rotvec` and `as_rotvec` functions as part of the `scipp.spatial.transform` module.\n",
    "Let's reverse our original rotation by applying a -30 deg around the `y` axis.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6784591d-fa75-4da9-a5ef-52ca9e3ae222",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipp.spatial.transform import from_rotvec\n",
    "rotation_back = from_rotvec(sc.vector(value=[0, -30.0, 0], unit=sc.units.deg))\n",
    "sensor.coords['position']['x',5:] = rotation_back * sensor.coords['position']['x',5:]\n",
    "sc.plot(sensor, projection=\"3d\", positions=\"position\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7296c82-a193-48c3-9b4d-c4339403be85",
   "metadata": {},
   "source": [
    "We can compound rotations by multipying our rotation variable.\n",
    "Lets make rotate the same again in the same axis to give an overall 60 deg around the `y` axis, applied to every pixel in our sensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "770ac127-041a-4051-b126-1470231b8cc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "rotation = rotation * rotation\n",
    "sensor.coords['position'] = rotation * sensor.coords['position']\n",
    "sc.plot(sensor, projection=\"3d\", positions=\"position\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ae7b111-f08f-44ba-bb3e-983721bea736",
   "metadata": {},
   "source": [
    "Scipp also provides vector and scalar products of vectors.\n",
    "Lets determine the surface normal vector for our sensor in the new rotation given the knowledge that all points are coplanar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9af721a-3224-47dd-b98c-62ff5e5e7fa2",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = sensor.coords['position']['x', 0]['y', 0]\n",
    "b = sensor.coords['position']['x', -1]['y', 0]\n",
    "c = sensor.coords['position']['x', 0]['y', -1]\n",
    "norm = sc.cross(c - a, c - b)\n",
    "norm /= sc.norm(norm)\n",
    "norm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15f0f7a5-91f7-4de3-b4ec-dbc643cc7e9f",
   "metadata": {},
   "source": [
    "Now lets verify that the normal vector we just calculated corresponds to a rotation of our original normal vector (0, 0, 1) by 60 degrees around `y`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee4f2417-f1f2-45b3-9c9b-479a6d58dd10",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipp.spatial.transform import as_rotvec\n",
    "sc.isclose(rotation * sc.vector(value=[0, 0, 1]), norm)"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}