EUREKA INVTOP TRIANG SPUR MATRAN.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:99e975d0a4a0b853030687186e6f40e672a8cc3c74b015655614d6fcb52ca01b"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import sys\n",
      "# sys.path.append(\"/Users/becqg/Documents/code/Robinson\")\n",
      "sys.path.append(\"C:\\Users\\Guillaume\\Documents\\projets\\code\\Robinson\")\n",
      "import ER\n",
      "import imp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> EUREKA</h2>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "LR = 2\n",
      "R = numpy.array([2.5625, -1.25])\n",
      "G = numpy.array([-1.25, 0])\n",
      "F = numpy.empty((LR, ))\n",
      "A = numpy.empty((LR, ))\n",
      "(F, A) = ER.EUREKA(LR, R, G, F, A)\n",
      "print(F)\n",
      "print(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[-0.6401249  -0.31225605]\n",
        "[ 1.          0.48780488]\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "V = 2.5625\n",
      "D = -1.25\n",
      "A1 = 1.\n",
      "F1 = -1.25 / 2.5625\n",
      "print(F1)\n",
      "Q = F1 * (-1.25)\n",
      "print(Q)\n",
      "A2 = -D / V\n",
      "print(A2)\n",
      "V += A2 * D\n",
      "F2 = (0 - Q) / V\n",
      "print(F2)\n",
      "F1 = F1 + F2 * A2\n",
      "print(F1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-0.487804878049\n",
        "0.609756097561\n",
        "0.487804878049\n",
        "-0.312256049961\n",
        "-0.64012490242\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "LR = 2\n",
      "R = numpy.array([10., 4])\n",
      "G = numpy.array([2., 0])\n",
      "F = numpy.empty((LR, ))\n",
      "A = numpy.empty((LR, ))\n",
      "(F, A) = ER.EUREKA(LR, R, G, F, A)\n",
      "print(F)\n",
      "print(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.23809524 -0.0952381 ]\n",
        "[ 1.  -0.4]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> INVTOP </h2>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "LR = 5\n",
      "R = numpy.array([10.0, 4.0, -1.0, -4.0, -4.0])\n",
      "F = numpy.empty(LR)\n",
      "A = numpy.empty(LR)\n",
      "G = numpy.array([1., 0., 0., 0., 0.])\n",
      "(F, A) = ER.EUREKA(LR, R, G, F, A)\n",
      "print(F)\n",
      "print(A)\n",
      "\n",
      "V = 10.0\n",
      "D = 4.0\n",
      "A1 = 1.\n",
      "F1 = 1. / V\n",
      "print(F1)\n",
      "Q = F1 * (4.0)\n",
      "print(Q)\n",
      "# IF(LR.EQ.1) RETURN\n",
      "# L = 2\n",
      "A2 = -D / V\n",
      "print(A2)\n",
      "# IF(L.EQ.2) GO TO 2\n",
      "# 2\n",
      "V += A2 * D\n",
      "F2 = (0.0 - Q) / V\n",
      "# L3 = L - 1 = 1\n",
      "# J = 1\n",
      "# K = L - J + 1 = 2\n",
      "F1 += F2 * A2\n",
      "print(\"F1: {0}\".format(F1))\n",
      "print(\"F2: {0}\".format(F2))\n",
      "D = 0.0\n",
      "Q = 0.0\n",
      "# I = 1\n",
      "# K = L - I + 2 = 3\n",
      "D += A1 * (-1.)\n",
      "Q += F1 * (-1.)\n",
      "# I = 2\n",
      "# K = 2\n",
      "D += A2 * 4.0\n",
      "Q += F2 * 4.0\n",
      "# L = 3\n",
      "A3 = -D / V\n",
      "print(A3)\n",
      "# IF (L.EQ.2) GO TO 2\n",
      "# L1 = (L - 2) / 2 = 0\n",
      "# L2 = L1 + 1 = 1\n",
      "# IF(L2.LT.2) GO TO 5\n",
      "# 5\n",
      "# IF(2*L1.EQ.L-2) GO TO 2\n",
      "A2 += A3 * A2\n",
      "print(A2)\n",
      "V += A3 * D\n",
      "F3 = (0. - Q) / V\n",
      "# L3 = L - 1 = 2\n",
      "# J = 1\n",
      "# K = L - J + 1 = 3\n",
      "F1 += F3 * A3\n",
      "# J = 2\n",
      "# K = 2\n",
      "F2 += F3 * A2\n",
      "print(\"F1: {0}\".format(F1))\n",
      "print(\"F2: {0}\".format(F2))\n",
      "print(\"F3: {0}\".format(F3))\n",
      "D = 0.0\n",
      "Q = 0.0\n",
      "# I = 1\n",
      "# K = L - I + 2 = 4\n",
      "D += A1 * (-4.0)\n",
      "Q += F1 * (-4.0)\n",
      "# I = 2\n",
      "# K = L - I + 2 = 3\n",
      "D += A2 * (-1.0)\n",
      "Q += F2 * (-1.0)\n",
      "# I = 3\n",
      "# K = L - I + 2 = 2\n",
      "D += A3 * (4.0)\n",
      "Q += F3 * (4.0)\n",
      "# L = 4\n",
      "A4 = -D / V\n",
      "# L1 = (L - 2) / 2 = 1\n",
      "# L2 = L1 + 1 = 2\n",
      "# IF(L2.LT.2) GO TO 5\n",
      "# J = 2\n",
      "HOLD = A2\n",
      "# K = L - J + 1 = 3\n",
      "A2 += A4 * A3\n",
      "A3 += A4 * HOLD\n",
      "# IF(2*L1.EQ.L-2) GO TO 2\n",
      "# 2\n",
      "V += A4 * D\n",
      "F4 = (0. - Q) / V\n",
      "# L3 = L - 1 = 3\n",
      "# J = 1\n",
      "# K = L - J + 1 = 4\n",
      "F1 += F4 * A4\n",
      "# J = 2\n",
      "# K = L - J + 1 = 3\n",
      "F2 += F4 * A3\n",
      "# J = 3\n",
      "# K = L - J + 1 = 2\n",
      "F3 += F4 * A2\n",
      "print(\"F1: {0}\".format(F1))\n",
      "print(\"F2: {0}\".format(F2))\n",
      "print(\"F3: {0}\".format(F3))\n",
      "print(\"F4: {0}\".format(F4))\n",
      "D = 0.0\n",
      "Q = 0.0\n",
      "# I = 1\n",
      "# K = L - I + 2 = 5\n",
      "D += A1 * (-4.0)\n",
      "Q += F1 * (-4.0)\n",
      "# I = 2\n",
      "# K = L - I + 2 = 4\n",
      "D += A2 * (-4.0)\n",
      "Q += F2 * (-4.0)\n",
      "# I = 3\n",
      "# K = L - I + 2 = 3\n",
      "D += A3 * (-1.0)\n",
      "Q += F3 * (-1.0)\n",
      "# I = 4\n",
      "# K = L - I + 2 = 2\n",
      "D += A4 * (4.0)\n",
      "Q += F4 * (4.0)\n",
      "# L = 5\n",
      "A5 = -D / V\n",
      "# L1 = (L - 2) / 2 = 1\n",
      "# L2 = L1 + 1 = 2\n",
      "# J = 2\n",
      "HOLD = A2\n",
      "# K = L - J + 1 = 4\n",
      "A2 += A5 * A4\n",
      "A4 += A5 * HOLD\n",
      "A3 += A5 * A3\n",
      "V += A5 * D\n",
      "F5 = (0. - Q) / V\n",
      "# L3 = L - 1 = 4\n",
      "# J = 1\n",
      "# K = L - J + 1 = 5\n",
      "F1 += F5 * A5\n",
      "# J = 2\n",
      "# K = L - J + 1 = 4\n",
      "F2 += F5 * A4\n",
      "# J = 3\n",
      "# K = L - J + 1 = 3\n",
      "F3 += F5 * A3\n",
      "# J = 4\n",
      "# K = L - J + 1 = 2\n",
      "F4 += F5 * A2\n",
      "print(\"F1: {0}\".format(F1))\n",
      "print(\"F2: {0}\".format(F2))\n",
      "print(\"F3: {0}\".format(F3))\n",
      "print(\"F4: {0}\".format(F4))\n",
      "print(\"F5: {0}\".format(F5))\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.1489899  -0.05656566  0.02727273  0.03232323  0.02676768]\n",
        "[ 1.         -0.37966102  0.18305085  0.21694915  0.17966102]\n",
        "0.1\n",
        "0.4\n",
        "-0.4\n",
        "F1: 0.119047619048\n",
        "F2: -0.047619047619\n",
        "0.309523809524\n",
        "-0.52380952381\n",
        "F1: 0.131661442006\n",
        "F2: -0.0689655172414\n",
        "F3: 0.0407523510972\n",
        "F1: 0.14418079096\n",
        "F2: -0.0623728813559\n",
        "F3: 0.0223728813559\n",
        "F4: 0.0424858757062\n",
        "F1: 0.14898989899\n",
        "F2: -0.0565656565657\n",
        "F3: 0.0272727272727\n",
        "F4: 0.0323232323232\n",
        "F5: 0.0267676767677\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "LR = 5\n",
      "R = numpy.array([10.0, 4.0, -1.0, -4.0, -4.0])\n",
      "Q = numpy.empty((LR * LR, ))\n",
      "SPACE = numpy.empty((2 * LR, ))\n",
      "(Q, SPACE) = ER.INVTOP(LR, R, Q, SPACE)\n",
      "print(Q)\n",
      "print(SPACE)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.1489899  -0.05656566  0.02727273  0.03232323  0.02676768 -0.05656566\n",
        "  0.16565657 -0.07272727  0.01010101  0.03232323  0.02727273 -0.07272727\n",
        "  0.16363636 -0.07272727  0.02727273  0.03232323  0.01010101 -0.07272727\n",
        "  0.16565657 -0.05656566  0.02676768  0.03232323  0.02727273 -0.05656566\n",
        "  0.1489899 ]\n",
        "[ 0.          0.          0.          0.          1.          1.\n",
        " -0.37966102  0.18305085  0.21694915  0.17966102]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "LR = 5\n",
      "R = numpy.array([10.0, 4.0, 1.0, 4.0, 4.0])\n",
      "Q = numpy.empty((LR * LR, ))\n",
      "SPACE = numpy.empty((2 * LR, ))\n",
      "(Q, SPACE) = ER.INVTOP(LR, R, Q, SPACE)\n",
      "print(numpy.round(Q, 5))\n",
      "print(numpy.round(SPACE, 5))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.15323 -0.06452  0.03763 -0.06452 -0.01344 -0.06452  0.17921 -0.08602\n",
        "  0.0681  -0.06452  0.03763 -0.08602  0.16129 -0.08602  0.03763 -0.06452\n",
        "  0.0681  -0.08602  0.17921 -0.06452 -0.01344 -0.06452  0.03763 -0.06452\n",
        "  0.15323]\n",
        "[ 0.       0.       0.       0.       1.       1.      -0.42105  0.24561\n",
        " -0.42105 -0.08772]\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "LR = 5\n",
      "R = numpy.array([16.72605, 0, 10.75437, 0., 2.44141])\n",
      "Q = numpy.empty((LR * LR, ))\n",
      "SPACE = numpy.empty((2 * LR, ))\n",
      "(Q, SPACE) = ER.INVTOP(LR, R, Q, SPACE)\n",
      "print(numpy.round(Q, 5))\n",
      "print(numpy.round(SPACE, 5))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.12867  0.      -0.12045  0.       0.05867  0.       0.10192  0.\n",
        " -0.06553  0.      -0.12045  0.       0.21468  0.      -0.12045  0.\n",
        " -0.06553  0.       0.10192  0.       0.05867  0.      -0.12045  0.\n",
        "  0.12867]\n",
        "[ 0.       0.       0.       0.       1.       1.       0.      -0.93613\n",
        "  0.       0.45594]\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> TRIANG </h2>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "N = 2\n",
      "TOP = numpy.array([64, -48, -48, 37])\n",
      "S = numpy.empty((N * N, ))\n",
      "SPACE = numpy.empty((N, ))\n",
      "(S, SPACE) = ER.TRIANG(N, TOP, S, SPACE)\n",
      "print(S, SPACE)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(array([ 8.,  0., -6.,  1.]), array([  0.,  36.]))\n"
       ]
      }
     ],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "N = 2\n",
      "TOP = numpy.array([64, -48, -48, 37], \"complex\")\n",
      "S = numpy.empty((N * N, ), \"complex\")\n",
      "SPACE = numpy.empty((N, ), \"complex\")\n",
      "(S, SPACE) = ER.TRIANG_CPLX(N, TOP, S, SPACE)\n",
      "print(S, SPACE)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(array([ 8.+0.j,  0.+0.j, -6.+0.j,  1.+0.j]), array([  0.+0.j,  36.+0.j]))\n"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = numpy.array([[64, -48], [-48, 37]])\n",
      "print(A)\n",
      "print(numpy.linalg.cholesky(A))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 64 -48]\n",
        " [-48  37]]\n",
        "[[ 8.  0.]\n",
        " [-6.  1.]]\n"
       ]
      }
     ],
     "prompt_number": 43
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "A = numpy.matrix([[1., 0, 0], [2, 3, 0], [4, 5, 6]])\n",
      "B = A.T * A\n",
      "TOP = numpy.array(B).flatten(0)\n",
      "print(B, TOP)\n",
      "N = 3\n",
      "# TOP = numpy.array([1, 1, 1, 4, 2, 1., 6, 5, 3.])\n",
      "S = numpy.empty((N * N, ))\n",
      "SPACE = numpy.empty((N, ))\n",
      "(S, SPACE) = ER.TRIANG(N, TOP, S, SPACE)\n",
      "print(S.reshape((N, N), order=\"F\").T)\n",
      "print(SPACE)\n",
      "print(numpy.linalg.cholesky(B))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(matrix([[ 21.,  26.,  24.],\n",
        "        [ 26.,  34.,  30.],\n",
        "        [ 24.,  30.,  36.]]), array([ 21.,  26.,  24.,  26.,  34.,  30.,  24.,  30.,  36.]))\n",
        "[[ 4.58257569  0.          0.        ]\n",
        " [ 5.67366515  1.34518542  0.        ]\n",
        " [ 5.23722937  0.2123977   2.91998558]]\n",
        "[  0.          32.19047619  27.47368421]\n",
        "[[ 4.58257569  0.          0.        ]\n",
        " [ 5.67366515  1.34518542  0.        ]\n",
        " [ 5.23722937  0.2123977   2.91998558]]\n"
       ]
      }
     ],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "S.reshape((3, 3), order=\"F\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 62,
       "text": [
        "array([[ 4.58257569,  5.67366515,  5.23722937],\n",
        "       [ 0.        ,  1.34518542,  0.2123977 ],\n",
        "       [ 0.        ,  0.        ,  2.91998558]])"
       ]
      }
     ],
     "prompt_number": 62
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "A = numpy.matrix([[1+1j, 0, 0], [2, 3, 0], [4, 5, 6]])\n",
      "B = A.T * A\n",
      "TOP = numpy.array(B).flatten(0)\n",
      "print(B, TOP)\n",
      "N = 3\n",
      "# TOP = numpy.array([1, 1, 1, 4, 2, 1., 6, 5, 3.])\n",
      "S = numpy.empty((N * N, ), \"complex\")\n",
      "SPACE = numpy.empty((N, ), \"complex\")\n",
      "(S, SPACE) = ER.TRIANG_CPLX(N, TOP, S, SPACE)\n",
      "print(S, SPACE)\n",
      "print(numpy.linalg.cholesky(B))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(matrix([[ 20.+2.j,  26.+0.j,  24.+0.j],\n",
        "        [ 26.+0.j,  34.+0.j,  30.+0.j],\n",
        "        [ 24.+0.j,  30.+0.j,  36.+0.j]]), array([ 20.+2.j,  26.+0.j,  24.+0.j,  26.+0.j,  34.+0.j,  30.+0.j,\n",
        "        24.+0.j,  30.+0.j,  36.+0.j]))\n",
        "(array([ 4.47770875+0.2233285j ,  0.00000000+0.j        ,\n",
        "        0.00000000+0.j        ,  5.79213275-0.28888622j,\n",
        "        1.40064771+1.19463825j,  0.00000000+0.j        ,\n",
        "        5.34658408-0.2666642j ,  0.72065036+1.59083002j,\n",
        "        3.08297902+0.09059755j]), array([  0.00000000+0.j        ,  33.46534653-3.34653465j,\n",
        "        26.50344828-0.55862069j]))\n",
        "[[  4.47213595e+00+0.j   0.00000000e+00+0.j   0.00000000e+00+0.j]\n",
        " [  5.81377674e+00+0.j   4.47213595e-01+0.j   0.00000000e+00+0.j]\n",
        " [  5.36656315e+00+0.j  -2.68328157e+00+0.j   3.33200094e-07+0.j]]\n"
       ]
      }
     ],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ER.SPUR(N, A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 69,
       "text": [
        "(10+1j)"
       ]
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "(1.02 // 2) * 2 + 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "1.0"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imp.reload(ER)\n",
      "A = numpy.matrix(\"[1, 2, 3, 4; 11, 12, 13, 14; 21, 22, 23, 24]\")\n",
      "A = numpy.array(A).flatten(\"F\")\n",
      "print(A)\n",
      "M, N = 3, 4\n",
      "AT = ER.MATRAN(M, N, A)\n",
      "print(AT)\n",
      "print(AT.reshape((4, 3), order=\"F\"))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 1 11 21  2 12 22  3 13 23  4 14 24]\n",
        "[ 1  2  3  4 11 12 13 14 21 22 23 24]\n",
        "[[ 1 11 21]\n",
        " [ 2 12 22]\n",
        " [ 3 13 23]\n",
        " [ 4 14 24]]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = numpy.zeros((10, 10))\n",
      "%time B = A.T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 9 \u00b5s, sys: 4 \u00b5s, total: 13 \u00b5s\n",
        "Wall time: 13.8 \u00b5s\n"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = A.flatten()\n",
      "%time AT = ER.MATRAN(10, 10, C)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 330 \u00b5s, sys: 124 \u00b5s, total: 454 \u00b5s\n",
        "Wall time: 386 \u00b5s\n"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}