{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "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": "1",
   "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": "2",
   "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": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "vecs.values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "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": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "vecs.fields.x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "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": "7",
   "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": "8",
   "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": "9",
   "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": "10",
   "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": "11",
   "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, cbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "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": "13",
   "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, cbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "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": "15",
   "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, cbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "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": "17",
   "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, cbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "We can compound rotations by multiplying 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": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "rotation = rotation * rotation\n",
    "sensor.coords['position'] = rotation * sensor.coords['position']\n",
    "pp.scatter3d(sensor, pos=\"position\", pixel_size=0.01, cbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "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": "21",
   "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": "22",
   "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": "23",
   "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
}