sympy-lambdify-vectorization-problem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Purpose\n",
    "\n",
    "Demonstrate an inconsistency in `sympy.lambdify`  (datatype (array or scalar) depends on numerical runtype-data).\n",
    "See sympy issue [#5642](https://github.com/sympy/sympy/issues/5642).\n",
    "\n",
    "\n",
    "# Code\n",
    "\n",
    "First load an extension to easily visualize intermediate results.\n",
    "More Info:\n",
    "http://nbviewer.jupyter.org/github/cknoll/displaytools/blob/master/example1_python3.ipynb\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# this can be commented out savely\n",
    "%load_ext displaytools3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "time.ctime()  := 'Wed Jul 12 11:45:46 2017'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "___\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "sp.__version__  := '1.0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "___\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "M := Matrix([\n",
       "[     x],\n",
       "[a*x**2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "M1 := Matrix([\n",
       "[   x],\n",
       "[x**2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "M0 := Matrix([\n",
       "[x],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "time.ctime() ##:\n",
    "\n",
    "import sympy as sp\n",
    "import numpy as np\n",
    "\n",
    "sp.__version__ ##:\n",
    "\n",
    "a, x = sp.symbols(\"a, x\")\n",
    "\n",
    "# general expression for a vectorfield\n",
    "M = sp.Matrix([x, a*x**2]) ##:\n",
    "\n",
    "# simplification (suppose value of a is determined at runtime)\n",
    "M1 = M.subs(a, 1) ##:\n",
    "M0 = M.subs(a, 0) ##:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# convert Matrix to flat list (funcs will also return a list of results)\n",
    "\n",
    "f1 = sp.lambdify(x, list(M1), modules=\"numpy\")\n",
    "f0 = sp.lambdify(x, list(M0), modules=\"numpy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "a1 := [array([0, 1, 2, 3, 4]), array([ 0,  1,  4,  9, 16])]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "a0 := [array([0, 1, 2, 3, 4]), 0]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---\n"
     ]
    }
   ],
   "source": [
    "xx = np.arange(5)\n",
    "\n",
    "# inconsistency\n",
    "a1 = f1(xx) # list of two arrays ##:\n",
    "a0 = f0(xx) # list of array and scalar ##:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "**convenient casting: (each column is the evaluated vectorfield at x)**\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3,  4],\n",
       "       [ 0,  1,  4,  9, 16]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "___\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "np.array(a1).shape  := (2, 5)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "___\n",
      "**inconvenient casting: -> object array**\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([array([0, 1, 2, 3, 4]), 0], dtype=object)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "___\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "np.array(a0).shape  := (2,)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "___\n"
     ]
    }
   ],
   "source": [
    "print(\"**convenient casting: (each column is the evaluated vectorfield at x)**\")\n",
    "np.array(a1) ##\n",
    "np.array(a1).shape ##:\n",
    "\n",
    "print(\"**inconvenient casting: -> object array**\")\n",
    "\n",
    "np.array(a0) ##\n",
    "np.array(a0).shape ##:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Suggestsion for a workarround:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3, 4],\n",
       "       [0, 0, 0, 0, 0]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(np.broadcast_arrays(*a0))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}