First try for CAS, SymPy with codegen

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. 取點
69. N 個點 可以解一元 N-1 次方程式 ( N-1 degree )
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
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. 要速度？何不試試編譯式語言？
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 ?