{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "3b80a41b",
            "cell_type": "markdown",
            "source": "<div class='alert alert-warning'>\n\nSciPy's interactive examples with Jupyterlite are experimental and may not always work as expected. Execution of cells containing imports may result in large downloads (up to 60MB of content for the first import from SciPy). Load times when importing from SciPy may take roughly 10-20 seconds. If you notice any problems, feel free to open an [issue](https://github.com/scipy/scipy/issues/new/choose).\n\n</div>",
            "metadata": {}
        },
        {
            "id": "aa25cbd8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy.spatial.transform import Rotation as R",
            "outputs": []
        },
        {
            "id": "c9af00fb",
            "cell_type": "markdown",
            "source": "Here we run the baseline Kabsch algorithm to best align two sets of\nvectors, where there is noise on the last two vector measurements of\nthe ``b`` set:\n",
            "metadata": {}
        },
        {
            "id": "fb16f7ec",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "a = [[0, 1, 0], [0, 1, 1], [0, 1, 1]]\nb = [[1, 0, 0], [1, 1.1, 0], [1, 0.9, 0]]\nrot, rssd, sens = R.align_vectors(a, b, return_sensitivity=True)\nrot.as_matrix()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0., 0., 1.],\n       [1., 0., 0.],\n       [0., 1., 0.]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "d7561a4c",
            "cell_type": "markdown",
            "source": "When we apply the rotation to ``b``, we get vectors close to ``a``:\n",
            "metadata": {}
        },
        {
            "id": "7cb9d050",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rot.apply(b)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0. , 1. , 0. ],\n       [0. , 1. , 1.1],\n       [0. , 1. , 0.9]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a860231a",
            "cell_type": "markdown",
            "source": "The error for the first vector is 0, and for the last two the error is\nmagnitude 0.1. The `rssd` is the square root of the sum of the\nweighted squared errors, and the default weights are all 1, so in this\ncase the `rssd` is calculated as\n``sqrt(1 * 0**2 + 1 * 0.1**2 + 1 * (-0.1)**2) = 0.141421356237308``\n",
            "metadata": {}
        },
        {
            "id": "8a96f122",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "a - rot.apply(b)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 0., 0.,  0. ],\n       [ 0., 0., -0.1],\n       [ 0., 0.,  0.1]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "e4982b75",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.sqrt(np.sum(np.ones(3) @ (a - rot.apply(b))**2))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.141421356237308"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b4c31550",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rssd",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.141421356237308"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "4466744a",
            "cell_type": "markdown",
            "source": "The sensitivity matrix for this example is as follows:\n",
            "metadata": {}
        },
        {
            "id": "9d080b13",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "sens",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0.2, 0. , 0.],\n       [0. , 1.5, 1.],\n       [0. , 1. , 1.]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "705e8685",
            "cell_type": "markdown",
            "source": "Special case 1: Find a minimum rotation between single vectors:\n",
            "metadata": {}
        },
        {
            "id": "8a4d1297",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "a = [1, 0, 0]\nb = [0, 1, 0]\nrot, _ = R.align_vectors(a, b)\nrot.as_matrix()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0., 1., 0.],\n       [-1., 0., 0.],\n       [0., 0., 1.]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "aae068de",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rot.apply(b)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([1., 0., 0.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "831a9fb4",
            "cell_type": "markdown",
            "source": "Special case 2: One infinite weight. Here we find a rotation between\nprimary and secondary vectors that can align exactly:\n",
            "metadata": {}
        },
        {
            "id": "7e7b8377",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "a = [[0, 1, 0], [0, 1, 1]]\nb = [[1, 0, 0], [1, 1, 0]]\nrot, _ = R.align_vectors(a, b, weights=[np.inf, 1])\nrot.as_matrix()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0., 0., 1.],\n       [1., 0., 0.],\n       [0., 1., 0.]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "52201727",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rot.apply(b)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0., 1., 0.],\n       [0., 1., 1.]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b96e887d",
            "cell_type": "markdown",
            "source": "Here the secondary vectors must be best-fit:\n",
            "metadata": {}
        },
        {
            "id": "bdb84dc7",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "a = [[0, 1, 0], [0, 1, 1]]\nb = [[1, 0, 0], [1, 2, 0]]\nrot, _ = R.align_vectors(a, b, weights=[np.inf, 1])\nrot.as_matrix()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0., 0., 1.],\n       [1., 0., 0.],\n       [0., 1., 0.]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "0adf04f5",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rot.apply(b)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0., 1., 0.],\n       [0., 1., 2.]])"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}