3G4Q3cii.ipynb

2014Q3.ipynb

{
  "cells": [
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "import numpy as np",
      "execution_count": 1,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## Getting to the answer ASAP"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "First lets define the left subdivision matrix and the bezier matrix, such that for a 1d curve q(t) = TMLG"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "L = (1/8)*np.array([[8, 0, 0, 0],[4,4,0,0],[2,4,2,0],[1,3,3,1]])\nM = np.array([[-1,3,-3,1],[3,-6,3,0],[-3,3,0,0],[1,0,0,0]])\nprint('L: ')\nprint(L)\nprint('\\n')\nprint('M: ')\nprint(M)",
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": "L: \n[[1.    0.    0.    0.   ]\n [0.5   0.5   0.    0.   ]\n [0.25  0.5   0.25  0.   ]\n [0.125 0.375 0.375 0.125]]\n\n\nM: \n[[-1  3 -3  1]\n [ 3 -6  3  0]\n [-3  3  0  0]\n [ 1  0  0  0]]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Now then, the centre point is at point s = 1 and t = 1\nTherefore at this point $T = S$ and additionally\n$SML = TML$"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "S = np.array([1,1,1,1])\nT = +S",
      "execution_count": 3,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "When we combine it into a doubly subdivided curve \nthe z values of left subdivision looks like\n$SML(Q_ZL^TM^TT^T)$\n\nLets evaluate it at the point (1.0,1.0) which we known is the corner of the patch, and thus the centre point"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "Q_z = np.array([[0,1,1,0],[1,2,2,1],[1,2,2,1],[0,1,1,0]])\nprint('Z at the centre is: ', S@M@L@(Q_z@L.T@M.T@T.T))",
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Z at the centre is:  1.5\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Great, we have the correct answer. \n\nIn fact an even easier way would have been to ignore the subdivision matrix and simply evaluate at $S=0.5,t=0.5$"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "S_temp = np.array([0.5**3,0.5**2,0.5,1])\nT_temp = +S_temp\nprint('Z at the centre is: ', S_temp@M@(Q_z@M.T@T_temp.T))",
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Z at the centre is:  1.5\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## But why does the crib answer work?\n\nConsider the form of $SM$"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "print('SM: ', S@M)",
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": "SM:  [0 0 0 1]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "If we were to multiply this by L we would get just the last row of L"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "print('SML: ', S@M@L)\nprint('L[-1,:]: ', L[-1,:])",
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": "SML:  [0.125 0.375 0.375 0.125]\nL[-1,:]:  [0.125 0.375 0.375 0.125]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "So what we now need to consider is the form of $SMLQ_z(TML)^T$\n\nLets relabel the last row in L as x to make it more clear it is just a vector"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "x = L[-1,:]\nprint('x: ', x)",
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": "x:  [0.125 0.375 0.375 0.125]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "So now what we really need to think about is the form of $xQz^T$, specifically we need to make this into a form of some weight matrix $W$ times the z values $Q_z$.\n\nTo do this, consider this double sided operation on a new matrix:    \n\n$x\\begin{bmatrix} 1 & 0 & 0 & 0\\\\ 0& 0 & 0 & 0 \\\\ 0& 0 & 0 & 0\\\\0& 0 & 0 & 0 \\end{bmatrix}x^T$"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "mat = np.zeros((4,4))\nmat[0,0] = 1\nprint('x@mat@x.T: ', x@mat@x.T)",
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "text": "x@mat@x.T:  0.015625\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "We will label a matrix of this form that is zero except for a 1 at position (i,j) as $S_{(i,j)}$\n\nIn fact we can write any 2d matrix in terms of sums of scalar coefficients $\\alpha_{(i,j)}$ and these matrices $S_{(i,j)}$:\n\n$Q_z = \\sum_j^n\\sum_i^m\\alpha_{(i,j)}S_{(i,j)}$\n\nThe $\\alpha_{(i,j)}$ are just $ Qz_{(i,j)} $\n\nThe upshot of this is that if we evaluate the response of each of these $S_{(i,j)}$ matrices to our double sided operation with x we will get the amount that each coefficient is of $Q_z$ contributes to the final result, basically a weighting..."
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "weight_matrix = np.zeros((4,4))\nfor i in range(4):\n    for j in range(4):\n        Smatrix = np.zeros((4,4))\n        Smatrix[i,j] = 1\n        w = x@Smatrix@x.T\n        weight_matrix[i,j] = w\n\nprint('Calculated weight matrix')\nprint(weight_matrix)",
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Calculated weight matrix\n[[0.015625 0.046875 0.046875 0.015625]\n [0.046875 0.140625 0.140625 0.046875]\n [0.046875 0.140625 0.140625 0.046875]\n [0.015625 0.046875 0.046875 0.015625]]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "This weight matrix multiplies element-wise with the z values and you sum them all up to give you an answer"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "print('Crib weight matrix')\ncrib_matrix = (1/64)*np.array([[1,3,3,1],[3,9,9,3],[3,9,9,3],[1,3,3,1]])\nprint(crib_matrix)\nprint('\\n')\nprint('Is our matrix the same at the one from the crib: ')\nprint(crib_matrix==weight_matrix)",
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Crib weight matrix\n[[0.015625 0.046875 0.046875 0.015625]\n [0.046875 0.140625 0.140625 0.046875]\n [0.046875 0.140625 0.140625 0.046875]\n [0.015625 0.046875 0.046875 0.015625]]\n\n\nIs our matrix the same at the one from the crib: \n[[ True  True  True  True]\n [ True  True  True  True]\n [ True  True  True  True]\n [ True  True  True  True]]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "print('Z value from weight matrix method ', np.sum(weight_matrix*Q_z))",
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Z value from weight matrix method  1.5\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Another fun fact, this weight matrix is the same as the outer product of the vector with itself"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "np.outer(x,x)",
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 13,
          "data": {
            "text/plain": "array([[0.015625, 0.046875, 0.046875, 0.015625],\n       [0.046875, 0.140625, 0.140625, 0.046875],\n       [0.046875, 0.140625, 0.140625, 0.046875],\n       [0.015625, 0.046875, 0.046875, 0.015625]])"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "_draft": {
      "nbviewer_url": "https://gist.github.com/c980a798ce36e96a7a4dac254b2da1da"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.2",
      "nbconvert_exporter": "python",
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "version": 3,
        "name": "ipython"
      }
    },
    "gist": {
      "id": "c980a798ce36e96a7a4dac254b2da1da",
      "data": {
        "description": "3G4Q3cii.ipynb",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}