{
 "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 several transformations.\n",
    "Variables containing vectors or transformations are created using\n",
    "\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",
    "- functions in [scipp.spatial](../generated/modules/scipp.spatial.rst) (for transformations), e.g.,\n",
    "    - [scipp.spatial.rotations](../generated/modules/scipp.spatial.rotations.rst) (for and array of rotations)\n",
    "    - [scipp.spatial.linear_transform](../generated/modules/scipp.spatial.linear_transform.rst) (for a single linear transform, i.e. rotation + scaling)\n",
    "    \n",
    "Let us consider an example, creating a 1-D variable containing vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb04ef53-39ef-4320-b7e9-2a7fdc61026a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipp as sc\n",
    "import scipp.spatial\n",
    "import numpy as np\n",
    "import plopp as pp\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.spatial.linear_transforms(dims=['x'], values=np.random.rand(2, 3, 3))\n",
    "inv = sc.spatial.linear_transforms(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 three-dimensional scatter plot 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.spatial.as_vectors(x, y, 0.0 * sc.units.m)},\n",
    ")\n",
    "pp.scatter3d(sensor, pos=\"position\", pixel_size=0.01)"
   ]
  },
  {
   "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",
    "pp.scatter3d(sensor, pos=\"position\", pixel_size=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbf54580-29d3-40bb-96bd-a93f3ff4c26f",
   "metadata": {},
   "source": [
    "We use `sc.spatial.linear_transform` 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.spatial.linear_transform(\n",
    "    value=[\n",
    "        [0.8660254, 0.0000000, 0.5000000],\n",
    "        [0.0000000, 1.0000000, 0.0000000],\n",
    "        [-0.5000000, 0.0000000, 0.8660254],\n",
    "    ]\n",
    ")\n",
    "\n",
    "sensor.coords['position']['x', 5:] = rotation * sensor.coords['position']['x', 5:]\n",
    "pp.scatter3d(sensor, pos=\"position\", pixel_size=0.01)"
   ]
  },
  {
   "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 the `rotations_from_rotvecs` function as part of the `scipp.spatial` 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 import rotations_from_rotvecs\n",
    "\n",
    "rotation_back = rotations_from_rotvecs(\n",
    "    rotation_vectors=sc.vector(value=[0, -30.0, 0], unit=sc.units.deg)\n",
    ")\n",
    "sensor.coords['position']['x', 5:] = rotation_back * sensor.coords['position']['x', 5:]\n",
    "pp.scatter3d(sensor, pos=\"position\", pixel_size=0.01)"
   ]
  },
  {
   "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",
    "pp.scatter3d(sensor, pos=\"position\", pixel_size=0.01)"
   ]
  },
  {
   "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": [
    "sc.isclose(rotation * sc.vector(value=[0, 0, 1]), norm)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}