The notebook is viewable as such here. The storage isn't quite as efficient as it could be, but I think it's about as good as can be expected with simple indexing. The size of that storage only goes as jmax**4.
Last active
August 29, 2015 14:10
-
-
Save moble/1cef404fc6c7b56eb650 to your computer and use it in GitHub Desktop.
Show how to efficiently store and access Wigner 3j coefficients when one of the m values is <=2
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "metadata": { | |
| "name": "", | |
| "signature": "sha256:2e8a13d9f563d142e2019cdeef729a1269912cabeaa79aa1b4239334669f7c0a" | |
| }, | |
| "nbformat": 3, | |
| "nbformat_minor": 0, | |
| "worksheets": [ | |
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Given that you know that one $m_i$ value (let's say $m_3$) will have magnitude at most 2, which means that $m_2$ must be within 2 of $-m_1$. Also, the triangle inequality helps cut down on the storage, so we can cut it down to scaling as $j_{max}^4$.\n", | |
| "\n", | |
| "Because of symmetries, we are free to stipulate that $j_2 \\leq j_1$. And let's suppose that $j_1$ can go up to $j_{max} = 80$. We can also cut off about half of the values by using symmetry to stipulate that $m_1 \\geq 0$. The naive way to do all this (shown below) includes some invalid indices. But that should be okay, since you'll just never access them." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "from __future__ import division, print_function\n", | |
| "import numpy as np\n", | |
| "from sympy import init_printing, symbols, summation, horner\n", | |
| "from sympy.physics.wigner import wigner_3j\n", | |
| "init_printing()\n", | |
| "\n", | |
| "try:\n", | |
| " from numba import njit\n", | |
| " from numba.utils import IS_PY3\n", | |
| "except ImportError:\n", | |
| " print(\"Couldn't find numba; that `length` function below will be slow...\")\n", | |
| " import sys\n", | |
| " def _identity_decorator_outer(*args, **kwargs):\n", | |
| " def _identity_decorator_inner(fn):\n", | |
| " return fn\n", | |
| " return _identity_decorator_inner\n", | |
| " njit = _identity_decorator_outer\n", | |
| " IS_PY3 = (sys.version_info[:2] >= (3, 0))\n", | |
| "if IS_PY3:\n", | |
| " xrange = range\n", | |
| "else:\n", | |
| " xrange = xrange" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 1 | |
| }, | |
| { | |
| "cell_type": "heading", | |
| "level": 1, | |
| "metadata": {}, | |
| "source": [ | |
| "Implementation" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Construct an array of Wigner 3j values as follows:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "jmax=10" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 2 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "@njit\n", | |
| "def length(jmax):\n", | |
| " i=0\n", | |
| " for j1 in xrange(jmax+1):\n", | |
| " for m1 in xrange(j1+1):\n", | |
| " for j2 in xrange(j1+1):\n", | |
| " for m2 in xrange(-m1-2,-m1+2+1):\n", | |
| " for j3 in xrange((j1-j2),(j1+j2)+1): # Assume j1>=j2\n", | |
| " i += 1\n", | |
| " return i\n", | |
| "print(length(jmax))" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "21780\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 3 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "%%time\n", | |
| "\n", | |
| "w3j = np.empty((length(jmax),), dtype=float)\n", | |
| "\n", | |
| "i=0\n", | |
| "for j1 in range(jmax+1):\n", | |
| " for m1 in range(j1+1):\n", | |
| " for j2 in range(j1+1):\n", | |
| " for m2 in range(-m1-2,-m1+2+1):\n", | |
| " for j3 in range((j1-j2),(j1+j2)+1): # Assume j1>=j2\n", | |
| " for m3 in [-m1-m2]:\n", | |
| " tmp = wigner_3j(j1,j2,j3,m1,m2,m3)\n", | |
| " try:\n", | |
| " w3j[i] = float(tmp.evalf(n=18))\n", | |
| " except AttributeError:\n", | |
| " w3j[i] = float(tmp)\n", | |
| " i += 1\n", | |
| "if i!= w3j.size:\n", | |
| " raise ValueError(\"Something went wrong in assigning the values!\")" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "CPU times: user 14.9 s, sys: 194 ms, total: 15.1 s\n", | |
| "Wall time: 15 s\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 4 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "print('I expect the computation with jmax=80 will take at least {0:.2g} hours.'.format(15*length(80)/(3600*length(10))))" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "I expect the computation with jmax=80 will take at least 10 hours.\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 5 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "You can save that array to a text file with something like" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "np.savetxt('Wigner3j_jmax_80_m3_2.txt', w3j)" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 6 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Then, load that array from the text file into whatever program you're using, and express the 3j values as a function with something like the following code:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "w3j = np.loadtxt('Wigner3j_jmax_80_m3_2.txt')\n", | |
| "\n", | |
| "def wigner3j(j1,j2,j3,m1,m2,m3):\n", | |
| " if abs(m1)>j1 or abs(m2)>j2 or abs(m3)>j3 or abs(j1-j2)>j3 or j1+j2<j3 or m1+m2+m3!=0:\n", | |
| " return 0.0\n", | |
| " if abs(m3)>2: # Re-order so that abs(m3) is <=2, if possible\n", | |
| " if abs(m2)<=2:\n", | |
| " return wigner3j(j1,j3,j2,m1,m3,m2) * (-1)**(j1+j2+j3)\n", | |
| " elif abs(m1)<=2:\n", | |
| " return wigner3j(j3,j2,j1,m3,m2,m1) * (-1)**(j1+j2+j3)\n", | |
| " else:\n", | |
| " raise ValueError(\"wigner3j was not designed to handle indices {0}\".format((j1,j2,j3,m1,m2,m3)))\n", | |
| " if j2>j1: # Re-order so that j1>=j2\n", | |
| " return wigner3j(j2,j1,j3,m2,m1,m3) * (-1)**(j1+j2+j3)\n", | |
| " if m1<0: # Use symmetry to eliminate ~1/2 the storage\n", | |
| " return wigner3j(j1,j2,j3,-m1,-m2,-m3) * (-1)**(j1+j2+j3)\n", | |
| " # The // division by 2 is standard c-style integer division\n", | |
| " return w3j[(j1*(j1*(j1*(5*j1 + 10) + 20*m1 + 5) + 40*m1 - 4) + j2*(20*j2 + 8*m1 + 8*m2 + 20) + 4*j3 + 24*m1 + 4*m2 + 8) // 4]" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [], | |
| "prompt_number": 7 | |
| }, | |
| { | |
| "cell_type": "heading", | |
| "level": 1, | |
| "metadata": {}, | |
| "source": [ | |
| "Explanation of the index" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Here, you can do a simple loop to get the indices. It will include some invalid sets of indices, but that should be fine; you'll just never use them." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "indices = [[j1,j2,j3,m1,m2,m3]\n", | |
| " for j1 in range(jmax+1) for m1 in range(j1+1)\n", | |
| " for j2 in range(j1+1) for m2 in range(-m1-2,-m1+2+1)\n", | |
| " for j3 in range(j1-j2,j1+j2+1) for m3 in [-m1-m2]]\n", | |
| "print(\"Total number of elements:\",len(indices))" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "Total number of elements: 21780\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 8 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The array is created by a nested loop, where the innermost loop goes over values of $j_3$, and the outermost loop goes over values of $j_1$. [Actually, the innermost loop as written is a trivial loop over one value of $m_3$, so we ignore it.] So to index it, you have to start with the number of iterations required to get the innermost loop from the initial value of $j_3$ to the desired value of $j_3$. For the second-to-innermost loop, you have to increase the index by the number of iterations of both loops required to increase the value of $m_2$ from its initial value to the desired value. And so on.\n", | |
| "\n", | |
| "This is achieved by the following sum:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "j1,j2,j3,m1,m2,m3,j_max = symbols('j1,j2,j3,m1,m2,m3,j_{max}', integer=True)\n", | |
| "index_value = (\n", | |
| " summation(1, (j3,j1-j2,j3-1))\n", | |
| " +summation(summation(1, (j3,j1-j2,j1+j2)), (m2,-m1-2,m2-1))\n", | |
| " +summation(summation(summation(1, (j3,j1-j2,j1+j2)), (m2,-m1-2,-m1+2)), (j2,0,j2-1))\n", | |
| " +summation(summation(summation(summation(1, (j3,j1-j2,j1+j2)), (m2,-m1-2,-m1+2)), (j2,0,j1)), (m1,0,m1-1))\n", | |
| " +summation(summation(summation(summation(summation(1, (j3,j1-j2,j1+j2)), (m2,-m1-2,-m1+2)), (j2,0,j1)), (m1,0,j1)), (j1,0,j1-1))\n", | |
| ").expand().simplify()\n", | |
| "\n", | |
| "index_value" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "latex": [ | |
| "$$\\frac{5 j_{1}^{4}}{4} + \\frac{5 j_{1}^{3}}{2} + 5 j_{1}^{2} m_{1} + \\frac{5 j_{1}^{2}}{4} + 10 j_{1} m_{1} - j_{1} + 5 j_{2}^{2} + 2 j_{2} m_{1} + 2 j_{2} m_{2} + 5 j_{2} + j_{3} + 6 m_{1} + m_{2} + 2$$" | |
| ], | |
| "metadata": {}, | |
| "output_type": "pyout", | |
| "png": "iVBORw0KGgoAAAANSUhEUgAAAsoAAAAwBAMAAAD6LkElAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAzXYQMplE74mrIma7\n3VSKKnSYAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJQUlEQVR4Ae1ab4hUVRQ/u/N3Z3Znd0NBqGQ1\nDBWhUQgliN0IoQ+hI6UoKs6HPggRrhkRJTgtCRqJ22aQ2p8h8ItBuwQh1QeXzOhD4RpIJlZLECKU\nisZmFm3n3H/v3vvuffPeSBuOXtiZe88953f+vLv3vfebC6C1Tq2vd/ML39CHCfu7T5YTWrS2+mVP\nOTbCpuYTz9c7epq3bj3L1H5PldfD8eaz7aqnbzRv3XqW7f2eKgOsaT7bXLUw1bx1q1jeVV5ZxVwu\nlGGVWeW2+9LrcKK0AyD1XOJsA1goXUts3XIGu6avUk6fANTMKuemp/twojgMqXnUSdYCWOjqS2ba\nitpzRuo8rWLdrHLb3sUq3+R3PwULsFvB3L6dYZn6MjCr3CEn8Hv1pDaI1VWwUKzEMmhtJVWOX9Zu\n2ainqqqMu/LomD4Tp69gYSs8FcegtXUOf/Qk3t6GxjHL0bKeasdnCwYAViwBOAfHjBldy9dXsOm5\nl074lG4f+QkYrUIGXsL73Kl39bQ7y/QIVsmNwWMLD+kTsfoKFm+iV2JZtLhS5wR8Cv+4kvwGivXu\nSddMDJkfNoZx66lkbkA563w/21JPwbZ6kxn7YZsEvIXNkGTI48aQnwB43EzjdYD9VYB9AOlxc0aO\n8puekF372w9razYxjvDbBNpMmOR6IIPvZm0VWPmO6W87wBa86c2G9Aa8BK72KuCVcDc/rFs/kTTC\nbyKcmVPOAHT1APRPAvxmeq0AnManD9ywS54qn4PVddNGjfywSqX5ToTf5kH/U8tCBY7VAYj0tKp8\nBIovYoVxofuq/Dw9nribH9atn0ga4TcRzgwqH9lLm+tm/LOqXFj/JpY+g+vYV2V+cdyxemHd6gml\ntChuwXYGY7aqzLPoqkRV+dsGqXphG9g1mm7kt5H9/zL/aOF39Ouocro8OmBXWWNE83gJIpoTtg2X\n4YND+L7J2NYIa8a5avMev0c2fdEQK46ORtVqTsNB7lk4YMw3GgjfAIXrnVQuR5VHaydxwtwxAkYU\nPo504YQtnMIqz4WH62iKbGtUQ85Vb26/qWH2vBmJFUcHNKpW92oDL61meoz5BgPpG9VGGE0kq5yq\nK8uOBWPYV1V+iCYCRjRbyw6QxGiatQM2dfK9MhQHIZ8oVLdf5qmjDLlBIwBjEFsHNKqWIWh5cEQu\nmAOd44YHGoR0mYbb9wVhXeoTHflVkmn8TBLF1cEjl34tSx31HbK2YB8oQ+dwwt9P3H6Zp9wEdNB2\n52mxdWDYQgjlwQRF5ztySJdhOX0f+WmEO7JtCpv/7OMzVra7pqe5XP+0rW1YrHL3MGSdvImOo/fd\nfpmntmuNqxxHJ2aVO+WC08MDO2c+GR2f24YsebaKESVRdu/R+SPZoS+pz5rXWkxglfsrkL0q2Fa0\n2Xro/UV9W5dUub36JM5VNIdfnJGeuq75sRw6nogVVSt8cktOCTMRE3TPGlrnydmDi299+MoRbjKy\n8AzLNmBESSHf+Tn03w+4LaTHmYHXWkxglUdrkL7O2VbiTmpf96WuVmlbopFqxLkKTIdfVJOettUY\nFmM3bCyHThCxOI/DdRRVK1gSLiVKWBz6YYL+uyE3GSDoOQdSAxcA43M0GVl4imWLYo0RXYXv51i1\n44rt8FqLCVVlYluJO0lPnsanHCheMZkU4lwlg+Lwi1FIT+eBMbfEboSwHDoqYnkeR+oIqlawJFxK\nQYpDP0zQfw0yOxSCjI9NKamNi/FNm40qy/G7zYkrNCOzFYwoXclyf5Vez99GswlSEdamMdu+OSzQ\njlGjHYOxrfhck0pN0c2wA8HoKWfpD9RmAeNcBabl11wrGTQkLGI3LCwZDeg6KmJ5HkcEBoKqFSwJ\nlxKwOPTDBN0TUJxSCEbOSmrhkm+75Xt7793e21ux5QAv9/Z+19s7m3i4/VViRPmVvMxez5/lVfZa\nBxNY5e4KlPDul8fLQnUtXaH3+Fyf/cS+j2OG/LK1EgAuolARi7EbBpZHR0bMzuNInYCqJRwp5UHi\noR8pQK6x+BdIBFZlOaWkBi4Ai49iNBu7aKZIjNiaUowov5JrAM5C9kbKuK4uawGLVe4ahyJnW1mV\n8TG/q4a/FpSt96LZvMqIZfo110qpAkcZc8tYFRuLOzV1VMTieZLpYPkEAyzYGREuUsLi0A8T4DNG\ncYqOWjlyVlIdF7jvcEEEfniCZ4tuT2uM6AuA//GZng8TVLl9EA0CtjVXA9x21hQmzSprnCursvIr\nz4TxOFcCLOZYgOyGjeXSURGL8zhMR6NqCQcvcB8rAFLCqXmsywS4tWUGQSEYu6SS6rjA42NYxofA\nN2R8wLI1GVH8BRbvW221WoIqw4+wfIAvGdoxRqt02nFO3toxNM7V8quvlcLZPfQjw+Uy/nPjdbCx\nWC6mThCxOI/DdQIGmLMzogoEzI+8csFXsHwyQNCrHEgNXB4fr5/+KfB1keizbE1GNDsB7fgq99ZA\nzCpv2PX9GKyYfwARJdv6DNu7Lr5iVVnjXC2/+lppw/ssVpmwiFWxsVgupk4QcREvCzaeb0DVcnZG\nVIGA+aEfLsivPQABgl7lQKrjCt+Egg0XtmwCnw0Vo8RGLFuuFmZEdY8SygfL5s/gJ63loBkjzYPl\nV18rwhixLFaFY+m5hHTEeRy3jpCeoV15dAzd6GoyZGfO8pyPwwAf6mTLEiZvGqNEgqeFGFyMKL1W\nYNOsceSGJT1Oi0qSgyTIMvIv+hQeOKblV6yVwBPDslgVjhWlI8/juHWYlAGLQz+BWhClM+cwrjRo\n3yl7xneHh/FyMKIB26EheGBRg9OiiuRgNvqIeQhhcr/2mTCOtctgVXQsAnfohM7jhHS4wH/oJxQf\nSyOEy6T08cE9qqt3fIyXYET5StMNrL4HlrQ4LWrpB0PuIRjzHpPKtRJMNsBiik3pxDEKwmjYq7mr\nHMlmua+k4coDa+gkHvjXSmKomTUo1t1VxijcjFK88CJg4wG0ltYy8FZ5W635VCNgmwe9dS1r/iqf\nv4msImBvAvVWNS0NeKvsYpTiphkBGxeilfRWgLfKHkYpVvYRsLHsW0zp4sGDf4gf/azMfIySpeYe\n+mHd+q0v3elO0ccoubXDUg9sWPH2kPztTFOwWc65WEI3bCzTFlR6bfqwKyuTUXJpRMs8sNFGd2bv\nVKCZCvwLSoWAC9rS7V8AAAAASUVORK5CYII=\n", | |
| "prompt_number": 9, | |
| "text": [ | |
| " 4 3 2 \n", | |
| "5\u22c5j\u2081 5\u22c5j\u2081 2 5\u22c5j\u2081 2 \n", | |
| "\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500 + 5\u22c5j\u2081 \u22c5m\u2081 + \u2500\u2500\u2500\u2500\u2500 + 10\u22c5j\u2081\u22c5m\u2081 - j\u2081 + 5\u22c5j\u2082 + 2\u22c5j\u2082\u22c5m\u2081 + 2\u22c5j\u2082\u22c5m\u2082 +\n", | |
| " 4 2 4 \n", | |
| "\n", | |
| " \n", | |
| " \n", | |
| " 5\u22c5j\u2082 + j\u2083 + 6\u22c5m\u2081 + m\u2082 + 2\n", | |
| " " | |
| ] | |
| } | |
| ], | |
| "prompt_number": 9 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "We can check that this actually correctly indexes the `indices` array like this:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "# Insert the result into the `index` function and the `wigner3j` function above\n", | |
| "print(horner(4*index_value))" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "j1*(j1*(j1*(5*j1 + 10) + 20*m1 + 5) + 40*m1 - 4) + j2*(20*j2 + 8*m1 + 8*m2 + 20) + 4*j3 + 24*m1 + 4*m2 + 8\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 10 | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "def index(j1,j2,j3,m1,m2,m3):\n", | |
| " return (j1*(j1*(j1*(5*j1 + 10) + 20*m1 + 5) + 40*m1 - 4) + j2*(20*j2 + 8*m1 + 8*m2 + 20) + 4*j3 + 24*m1 + 4*m2 + 8) // 4 \n", | |
| "\n", | |
| "def isvalid(j1,j2,j3,m1,m2,m3):\n", | |
| " if abs(m1)>j1 or abs(m2)>j2 or abs(m3)>j3 or abs(j1-j2)>j3 or j3>j1+j2:\n", | |
| " return False\n", | |
| " return True\n", | |
| "\n", | |
| "for i in range(len(indices)):\n", | |
| " assert i == index(*indices[i]), '{0}, {1}, {2}'.format(i,indices[i],index(*indices[i]))\n", | |
| "print(\"That seems to have worked\")" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "That seems to have worked\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 11 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "As an additional check, we can calculate the total size of the array, as the index with $j_1 = j_{max}+1$ and everything else as the lowest possible indices for that value of $j_1$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "n_indices = index_value\n", | |
| "n_indices = n_indices.subs({j2:0,j3:j1,m1:0,m2:-2,m3:2}).subs(j1,j_max+1).expand().simplify()\n", | |
| "display(n_indices)\n", | |
| "assert (n_indices.subs(j_max,jmax)) == len(indices), ((n_indices.subs(j_max,jmax)),len(indices))\n", | |
| "print(\"The size {0} seems to be correct.\".format(n_indices.subs(j_max,jmax)))" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "latex": [ | |
| "$$\\frac{5 j_{{max}}^{4}}{4} + \\frac{15 j_{{max}}^{3}}{2} + \\frac{65 j_{{max}}^{2}}{4} + 15 j_{{max}} + 5$$" | |
| ], | |
| "metadata": {}, | |
| "output_type": "display_data", | |
| "png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAAwBAMAAAAsveaxAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAzXYQMplE74mrIma7\n3VSKKnSYAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFmklEQVRoBeWZW2gcVRjH/7szm52dvSTboog3\nloRVUqIMEWkVJUOl6Itm1dY2RMk++NDSh7TYl2qp62rAgKXpWtEQ0AHJgxfY9Y5VyEItPqh0C94e\nioSiiC+6NRI1ta7fzM7O7s71nE3z5HnY+c73/c43/z2Zmf0PATpGoiMODouzSjDUJmbOcOHthQzR\nbzytI5VIlaFnC5G02EArvtJH4TiX8EK8zqEgqYlrHDgXGh7jEQ5EeLYwmpdXudRwwA9wCt+qcjQH\nUitceCC8SdmRJ+hnBYUg4ZuB0JD4CNGpJ4DRTf69peJznTiSGX+et/p043d9yceIawHC36YtizYa\nGaLjc/Thf3MKV+PBLnyGllzJcU1Ja7Ybhb/wuyZIeGh2i3Vy8V8rdAlCBXzZicdzLtB6UvreGeOn\nXVOTrdj1GCbhMauyVZUvWROXYEmjZBvHXjzuQq0jZQkHyopvn27hSS1OF7r3+FwvtYWLN/562hvu\npbLw7k5AKFbpwj37qm8DQ/gnWRXYNgKIxVnND794XzaPmIXTzVH3w/lrp1HOI4KnglfqwhOK/jjO\nRStBuHCxgqPseFA713qihlPwvdGay3ThNL6i50//shH6fAgNBbv1S48N9+nkXYqsQelj+D02hU9p\nAqY173Zm5TKwlKeYEQ/sZwfI+kj0x5dq9oJzrgs/ARzPA8ecVXvmIAmvsOP25YHz6AAiJCiUCySh\nCyc5U3QBXBVM01NlKc+OBze0EREgOQCMLdvyLlNdOH2/c/QQYrgjyCPv5sBdzuefknP4SAPoNB+8\nP5gb/Ay4ffAkpOz2k451uvBFxJ80/FJ8cns2PKJCeGhYRXEnfdhGoiAc5cBtqxmmi7MPE/UY8OEl\nHFaGIXwtrmAoNvCtfW3k5csLkPe8RFdKpIb4dCX6YngOd6vjSlg5kczbccxMqDy4Yz1b4hvIozWs\nYg8iObEGNVTxW5fM4fVb0J+PFLAL90PEdX40mPDNXdbTr5/pUHXkHvkP+nnOiAfwA6LLtI/o1/S8\n+xCVsgqUcC+Sy/QjM0yU38OUDbdbT/dTG1nToVIs/5Wgmy6qhqtyXR3TEu8pGIdqLhU0M7AO5cIZ\nig9hBGVRWMMXt0Gsi61qj7jDerb6dR7vbE46HGppklJl2vN4tRBVx/fL9WMomEtSGTOwDrEsXUdy\nnfb8UwmHxNVKWapJrWqvuH7nty1Zq1v38UJzasf2IaTJryh9pQuTQmnvoGIucioxCn1zuBZvPIt9\nt2Yrd9xUHGqdold8HcJbp7YdPZTYKGvaK24It7yk1a4raO24gXVV3Ca9KnHr5cxZ3XXhQdbTFN7E\nnL1sGau1Le8x7RXXhdPw9ZKmcANDw3sYndBS0u8NGpX6unDDCFEHXy/ZFj6lNc/m+Sml0zccTKfp\nYck01oHrO+7rJZ9Jp79Lp3V3Z2JBilo7HsSZ9V5xXbjlJT2MkLnjJhYkqFclQX2bdau7Lpz+sOcM\n6+lhhEzhTSywv9U6kDSAXnFduGU9PYyQKbyJBcrpVUlg467v2W09PYyQKbzpULv6iy53oYfwxUfJ\nursMDxzY4QKTs8840uQlPYyQKdyxgt6Fqs5kX8WZowtxznj3dJbcceLOO1nKOHDDS3oYof2uLfTk\nzVXPkq0QUxA9YMv5TsNHfMtW0bCeXkbIomyB/E7VlvGcRmuIkYdnH29dz8Ya1pMNbVNSqNqe+Eeh\nFU7hBUbh/qf1qJ5iF04dkvQAYx5xbQOFCxku4dMFZtn07wtsoHAJXMLPc+im96kNFP4jl3D/f6TY\nvlRK3UDhQoFLuP5yzzy2YQOFp+bnXzucYdWSyuFNVhb4ZX7+zxI7zk1Gq8xL6Bd8CzOsg0e4aE64\nv8q6QP7++YkaK2xw/3DRfLB09m9XY+LSJUSva1zCX2gsuLT5H6b+AyIZ6zNXOK+rAAAAAElFTkSu\nQmCC\n", | |
| "text": [ | |
| " 4 3 2 \n", | |
| "5\u22c5j_{max} 15\u22c5j_{max} 65\u22c5j_{max} \n", | |
| "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 + 15\u22c5j_{max} + 5\n", | |
| " 4 2 4 " | |
| ] | |
| }, | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "The size 21780 seems to be correct.\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 12 | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Just for fun, let's actually test that the result is the same as sympy's for every value you'll need up to $j_{1,2}=10$, and for every such value with $j_{1,2}=j_{max}$. (That limited range should be convincing enough, and the rest just take too long.)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "collapsed": false, | |
| "input": [ | |
| "for j1 in range(10) + [jmax]:\n", | |
| " for m1 in range(-j1,j1+1):\n", | |
| " for j2 in range(10) + [jmax]:\n", | |
| " for m2 in range(-j2,j2+1):\n", | |
| " for j3 in range(10) + [jmax]:\n", | |
| " for m3 in range(-j3,j3+1):\n", | |
| " if abs(m3)>2 and abs(m2)>2 and abs(m1)>2:\n", | |
| " continue\n", | |
| " a = wigner_3j(j1,j2,j3,m1,m2,m3)\n", | |
| " try:\n", | |
| " a = float(a.evalf(n=18))\n", | |
| " except AttributeError:\n", | |
| " a = float(a)\n", | |
| " b = wigner3j(j1,j2,j3,m1,m2,m3)\n", | |
| " assert a==b, 'sympy...wigner_3j(*{0})={1}; wigner3j(*{0})={2}'.format([j1,j2,j3,m1,m2,m3], a, b)\n", | |
| "print(\"That seems to have worked\")" | |
| ], | |
| "language": "python", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "stream": "stdout", | |
| "text": [ | |
| "That seems to have worked\n" | |
| ] | |
| } | |
| ], | |
| "prompt_number": 13 | |
| } | |
| ], | |
| "metadata": {} | |
| } | |
| ] | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment