Upgrade to Pro — share decks privately, control downloads, hide ads and more …

First try for CAS, SymPy with codegen

95963a44712b973e2fd7766baa419478?s=47 dv
June 05, 2016

First try for CAS, SymPy with codegen

95963a44712b973e2fd7766baa419478?s=128

dv

June 05, 2016
Tweet

Transcript

  1. First try for CAS, SymPy with codegen 許邱翔 Chiu-Hsiang Hsu

    wdv4758h@gmail.com
  2. Who am I ? • 我朋友都叫我 "dv" • Pythonista •

    FOSS 熱愛 • 快畢業的大學生 • 即將進入四個月國軍 Online • 出來後準備找工作
  3. 先備知識 • Python 基本能力 • 一般 Python 套件安裝方式

  4. Talk 目標: 會用 SymPy 協助處理你的數學問題!

  5. 調查

  6. 接下來會有 Bob 和 Alice 兩位同學一起參與

  7. Bob:「等等,CAS 到底是什麼?」

  8. Bob:「難道是 嗎?」

  9. Computer Algebra System

  10. 符號運算

  11. CAS (Computer Algebra System) 一些知名的軟體: • Maxima (GPL) • SageMath

    (GPL) • SymPy (BSD License) • Mathematica (Proprietary) • Magma (Proprietary) • Maple (Proprietary) • MATLAB - Symbolic Math Toolbox (Proprietary)
  12. Bob:「這聽起來好數學喔, 我快睡著了」

  13. 別擔心, 以下都是站在外行人的角度來看

  14. Alice:「那什麼是 SymPy?」

  15. 什麼是 SymPy?

  16. 什麼是 SymPy? • Python 的 Library

  17. 什麼是 SymPy? • Python 的 Library • Open Source (BSD

    License)
  18. 什麼是 SymPy? • Python 的 Library • Open Source (BSD

    License) • 針對數學符號做處理 (Symbolic Computation)
  19. 什麼是 SymPy? • Python 的 Library • Open Source (BSD

    License) • 針對數學符號做處理 (Symbolic Computation) • 安裝方便 ( 依賴的軟體相對少 )
  20. 基本 SymPy 使用

  21. 如何安裝 SymPy?

  22. $ pip install sympy

  23. 為了讓沒安裝的人可以嘗試 ……

  24. SymPy Live

  25. SymPy Live live.sympy.org

  26. Try Jupyter!

  27. Try Jupyter! try.jupyter.org 1. 選 "New" 然後 "Python 3" 2.

    輸入 "import sympy" 然後按 "Shift + Enter" 3. 現在你可以使用任何 SymPy 的功能了
  28. SageMathCloud

  29. SageMathCloud cloud.sagemath.com • 註冊 / 登入 • 建立新專案 • 在新專案中建立

    Jupyter Notebook • 選擇 kernel
  30. 說文解字

  31. SymPy

  32. Python

  33. Symbolic Python

  34. 符號!

  35. 什麼是 Symbolic? Symbolic • π • e Numeric • 3.141592653589793

    • 2.718281828459045
  36. 建立 Symbol 3

  37. x = Symbol ( "x" ) from sympy import Symbol

  38. x, y, z = symbols ( "x, y, z" )

    from sympy import symbols
  39. var ( "x, y, z" ) from sympy import var

  40. 建立 Expression 2

  41. x, y = symbols ( "x, y" ) expr =

    x**2 + 2*x*y + y**2
  42. 各式數學操作 from sympy import sin, cos, ...

  43. expr = sympify ( "x**2 + 2*x*y + y**2" )

    from sympy import sympify
  44. 帶入求值 2

  45. expr = sympify ( "x**2 + 2*x*y + y**2" )

    result = expr.subs ( {"x": 1, "y": 2} )
  46. expr = sympify ( "x**2 + 2*x*y + y**2" )

    result = expr.evalf ( subs = {"x": 1, "y": 2} )
  47. 方程式化簡 1

  48. expr = sympify ( "sin(x)**2 + cos(x)**2" ) new_expr =

    expr.simplify ()
  49. 方程式展開 1

  50. expr = sympify ( "(x + y)**2" ) new_expr =

    expr.expand ()
  51. 解方程式 1

  52. expr = sympify ( "x**4 - 1" ) result =

    solve ( expr ) [-1, 1, -I, I] from sympy import solve
  53. expr1 = sympify ( "x + y - 10" )

    expr2 = sympify ( "x - y + 6" ) result = solve ( [expr1, expr2] ) {y: 8, x: 2}
  54. 等等, 講這麼多我當初到底是 什麼需求來用 SymPy 的?

  55. 讓電腦可以計算三角函數!

  56. But how ?

  57. 逼近!

  58. Chebyshev Approximation

  59. Chebyshev Approximation 那是什麼鬼?我用的到嗎?

  60. Chebyshev Approximation 那是什麼鬼?我用的到嗎? (你可能已經用過了)

  61. Polynomial Approximation

  62. 目標:sin ( x ) 多項式逼近:a n xn + a n-1

    xn-1 + … + a 1 x1 + a 0
  63. 取點

  64. None
  65. None
  66. None
  67. None
  68. None
  69. N 個點 可以解一元 N-1 次方程式 ( N-1 degree )

  70. None
  71. None
  72. static const double half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */

    S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ FreeBSD - msun (k_sin.c) 方程式係數
  73. GNU C Library (s_sin.c) static const double sn3 = -1.66666666666664880952546298448555E-01,

    sn5 = 8.33333214285722277379541354343671E-03, cs2 = 4.99999999999999999999950396842453E-01, cs4 = -4.16666666666664434524222570944589E-02, cs6 = 1.38888874007937613028114285595617E-03; 方程式係數
  74. musl (__sin.c) static const double S1 = -1.66666666666666324348e-01, /* 0xBFC55555,

    0x55555549 */ S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ 方程式係數
  75. 自己用 SymPy 寫 Chebyshev Approximation

  76. None
  77. 轉成可執行 Function 9

  78. lambdify ( args, expr, module=None ) from sympy import lambdify

  79. lambdify ( x, x**x )

  80. lambdify ( (x, y) , x**x + y**y )

  81. lambdify ( x, x**x, "math" ) <function <lambda>>

  82. lambdify ( x, x**x, "numpy" ) <function numpy.<lambda>>

  83. lambdify ( x, x**x, "mpmath" ) <function <lambda>>

  84. lambdify ( x, x**x, "numexpr" ) <function numexpr.<lambda>>

  85. lambdify ( x, x**x, "sympy" ) <function <lambda>>

  86. # NumPy universal function ufuncify ( args, expr, backend="numpy" )

    from sympy.utilities.autowrap import ufuncify
  87. # Python & NumPy's C API # and C compiler

    ufuncify ( x, x**x ) <ufunc 'wrapper_module_0'>
  88. # Cython ufuncify ( x, x**x, backend="cython") <function wrapper_module_0.autofunc_c>

  89. # F2Py & Fortran compiler ufuncify ( x, x**x, backend="f2py")

    <fortran object>
  90. # Theano theano_function ( args, exprs ) from sympy.printing.theanocode import

    theano_function
  91. theano_function ( [x], [x**x] ) <theano.compile.function_module.Function at 0x7fd9a0bede48>

  92. # broadcast and type theano_function ( [x], [x**x], dims={x: 1},

    dtype={x: 'float64'}) <theano.compile.function_module.Function at 0x7fd9a0bede48>
  93. Library 內建的 Chebyshev Approximation

  94. poly, err = chebyfit ( sin, [0, pi/4], 10, error=True)

    from mpmath import chebyfit, sin, pi poly: [2.53705546510064e-6, 4.83882390985545e-7, -0.000198913241162381, 2.96736201694264e-7, 0.00833322702020108, 2.28887435232138e-8, -0.166666669474202, 1.73888384312788e-10, 0.999999999995841, 1.63478544118895e-14]
  95. poly, err = chebyfit ( sin, [0, pi/4], 10, error=True)

    from mpmath import chebyfit, sin, pi err: 1.94337098561508e-14
  96. 找完逼近的方程式了 現在想要讓執行變快 怎麼辦?

  97. NEED FOR SPEED

  98. Bob:「肯定是 Python 直譯器不夠快, 我們來換一個實作吧」

  99. Python 直譯器實作 • CPython ( 官方 ) • PyPy -使用

    Python 的子集來實作,採用 JIT 與不同的 Garbage Collection 技術 • Pyston -基於 LLVM 實作 JIT • Jython -把 Python code 轉成 Java bytecode 跑在 JVM 上 • IronPython -用 .NET Framework 實作 • …
  100. PyPy pypy.org • 專案狀態 ◦ Python 2.7 有支援 ◦ Python

    3 目前指支援到 Python 3.3 (alpha) • 效能比較可以看 speed.pypy.org
  101. Pyston github.com/dropbox/pyston • 由 Dropbox 於 2013 發起的專案 • 使用

    C++11 撰寫 • 基於 LLVM 實作 JIT • 目前方向著重在 Python 2.7
  102. CPython:「你這樣講就不對了, 怎麼可以怪我呢?」

  103. Alice:「不要急著怪別人」

  104. Alice:「其實你不需要 SymPy 了啊, 為什麼不寫成其他形式?」

  105. 自己重寫一遍好麻煩啊

  106. SymPy codegen

  107. codegen ( (name, expr), lang ) from sympy.utilities.codegen import codegen

  108. 要速度?何不試試編譯式語言?

  109. None
  110. C

  111. codegen ( ("func", x**x), "c", header=False) [("func.c", "..."), ("func.h", "...")]

  112. codegen ( ("func", x**x), "c", header=False) #include "func.h" #include <math.h>

    double func(double x) { double func_result; func_result = pow(x, x); return func_result; }
  113. codegen ( ("func", x**x), "c", header=False) #ifndef PROJECT__FUNC__H #define PROJECT__FUNC__H

    double func(double x); #endif
  114. Fortran

  115. codegen ( ("func", x**x), "f95" ) [("func.f90", "..."), ("func.h", "...")]

  116. 喜歡類似 MATLAB 的語言? Octave

  117. codegen ( ("func", x**x), "octave" ) [("func.m", "...")]

  118. 試試針對科學運算的新語言? Julia

  119. codegen ( ("func", x**x), "julia" ) [("func.jl", "...")]

  120. 藏在 codegen 背後的重要資料-AST ( Expression Tree )

  121. SymPy Expression Tree 把 SymPy 的 Expression Tree 轉成 DOT

    語言描述 (接著可以用 Graphviz 這類工具轉成圖檔) from sympy.printing.dot import dotprint
  122. x+2 SymPy Expression Tree

  123. x**2 + 2*x*y + y**2 SymPy Expression Tree

  124. x / (y + 42) SymPy Expression Tree

  125. 優化! Tree 操作!

  126. SymPy Expression Tree sin(x)**2 + cos(x)**2 .simplify()

  127. SymPy Expression Tree x**x + y**y .subs({x: 4})

  128. Visitor Pattern

  129. _print_<EXPR_CLASS>

  130. _print ( expr ) _print_Rational ( expr ) _print_Number (

    expr ) _print_AtomicExpr ( expr ) _print_Atom ( expr ) _print_Basic ( expr )
  131. SymPy 運作流程圖(總結)

  132. Symbol

  133. Symbol Operation

  134. Symbol Symbol Operation Operation Symbol Operation Operation Symbol Symbol

  135. Symbol Operation Operation Symbol Symbol Symbol Operation Symbol

  136. Symbol Operation Symbol C Fortran Octave Julia DOT LaTeX JavaScript

    NumPy Theano Numexpr Python ...
  137. C C compiler ufuncify Fortran Fortran compiler F2PY Cython Python

    & Numpy C API module.so
  138. SymPy 內其他 Modules

  139. SymPy 內其他 Modules • Calculus • Physics • Lie Algebra

    • Linear Algebra • Number Theory • Differential Geometry • Set Theory • ODE • PDE • Statistics • Logic • …
  140. 其他相關專案

  141. SymPy Gamma

  142. SymPy Gamma gamma.sympy.org

  143. SymEngine github.com/symengine/symengine

  144. SymEngine github.com/symengine/symengine (C++ 撰寫的 SymPy,為了提升效能)

  145. 基於 SymPy 的專案

  146. 基於 SymPy 的專案 • PyDy - Multibody Dynamics with Python

    • SageMath - 有把 SymPy 整合進去 • SfePy - Simple Finite Elements in Python • …
  147. 本投影片將放在 speakerdeck.com/wdv4758h/first-try-for-cas- sympy-with-codegen

  148. Thank you ! ! !

  149. Questions ?