Tweet
Tweet

Slide 1

Slide 1 text

First try for CAS, SymPy with codegen 許邱翔 Chiu-Hsiang Hsu wdv4758h@gmail.com

Slide 2

Slide 2 text

Who am I ? ● 我朋友都叫我 "dv" ● Pythonista ● FOSS 熱愛 ● 快畢業的大學生 ● 即將進入四個月國軍 Online ● 出來後準備找工作

Slide 3

Slide 3 text

先備知識 ● Python 基本能力 ● 一般 Python 套件安裝方式

Slide 4

Slide 4 text

Talk 目標： 會用 SymPy 協助處理你的數學問題！

Slide 5

Slide 5 text

調查

Slide 6

Slide 6 text

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

Slide 7

Slide 7 text

Bob：「等等，CAS 到底是什麼？」

Slide 8

Slide 8 text

Bob：「難道是 嗎？」

Slide 9

Slide 9 text

Computer Algebra System

Slide 10

Slide 10 text

符號運算

Slide 11

Slide 11 text

CAS (Computer Algebra System) 一些知名的軟體： ● Maxima (GPL) ● SageMath (GPL) ● SymPy (BSD License) ● Mathematica (Proprietary) ● Magma (Proprietary) ● Maple (Proprietary) ● MATLAB - Symbolic Math Toolbox (Proprietary)

Slide 12

Slide 12 text

Bob：「這聽起來好數學喔， 我快睡著了」

Slide 13

Slide 13 text

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

Slide 14

Slide 14 text

Alice：「那什麼是 SymPy？」

Slide 15

Slide 15 text

什麼是 SymPy？

Slide 16

Slide 16 text

什麼是 SymPy？ ● Python 的 Library

Slide 17

Slide 17 text

什麼是 SymPy？ ● Python 的 Library ● Open Source (BSD License)

Slide 18

Slide 18 text

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

Slide 19

Slide 19 text

什麼是 SymPy？ ● Python 的 Library ● Open Source (BSD License) ● 針對數學符號做處理 (Symbolic Computation) ● 安裝方便 ( 依賴的軟體相對少 )

Slide 20

Slide 20 text

基本 SymPy 使用

Slide 21

Slide 21 text

如何安裝 SymPy？

Slide 22

Slide 22 text

$ pip install sympy

Slide 23

Slide 23 text

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

Slide 24

Slide 24 text

SymPy Live

Slide 25

Slide 25 text

SymPy Live live.sympy.org

Slide 26

Slide 26 text

Try Jupyter!

Slide 27

Slide 27 text

Try Jupyter! try.jupyter.org 1. 選 "New" 然後 "Python 3" 2. 輸入 "import sympy" 然後按 "Shift + Enter" 3. 現在你可以使用任何 SymPy 的功能了

Slide 28

Slide 28 text

SageMathCloud

Slide 29

Slide 29 text

SageMathCloud cloud.sagemath.com ● 註冊 / 登入 ● 建立新專案 ● 在新專案中建立 Jupyter Notebook ● 選擇 kernel

Slide 30

Slide 30 text

說文解字

Slide 31

Slide 31 text

SymPy

Slide 32

Slide 32 text

Python

Slide 33

Slide 33 text

Symbolic Python

Slide 34

Slide 34 text

符號！

Slide 35

Slide 35 text

什麼是 Symbolic？ Symbolic ● π ● e Numeric ● 3.141592653589793 ● 2.718281828459045

Slide 36

Slide 36 text

建立 Symbol 3

Slide 37

Slide 37 text

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

Slide 38

Slide 38 text

x, y, z = symbols ( "x, y, z" ) from sympy import symbols

Slide 39

Slide 39 text

var ( "x, y, z" ) from sympy import var

Slide 40

Slide 40 text

建立 Expression 2

Slide 41

Slide 41 text

x, y = symbols ( "x, y" ) expr = x**2 + 2*x*y + y**2

Slide 42

Slide 42 text

各式數學操作 from sympy import sin, cos, ...

Slide 43

Slide 43 text

expr = sympify ( "x**2 + 2*x*y + y**2" ) from sympy import sympify

Slide 44

Slide 44 text

帶入求值 2

Slide 45

Slide 45 text

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

Slide 46

Slide 46 text

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

Slide 47

Slide 47 text

方程式化簡 1

Slide 48

Slide 48 text

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

Slide 49

Slide 49 text

方程式展開 1

Slide 50

Slide 50 text

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

Slide 51

Slide 51 text

解方程式 1

Slide 52

Slide 52 text

expr = sympify ( "x**4 - 1" ) result = solve ( expr ) [-1, 1, -I, I] from sympy import solve

Slide 53

Slide 53 text

expr1 = sympify ( "x + y - 10" ) expr2 = sympify ( "x - y + 6" ) result = solve ( [expr1, expr2] ) {y: 8, x: 2}

Slide 54

Slide 54 text

等等， 講這麼多我當初到底是 什麼需求來用 SymPy 的？

Slide 55

Slide 55 text

讓電腦可以計算三角函數！

Slide 56

Slide 56 text

But how ?

Slide 57

Slide 57 text

逼近！

Slide 58

Slide 58 text

Chebyshev Approximation

Slide 59

Slide 59 text

Chebyshev Approximation 那是什麼鬼？我用的到嗎？

Slide 60

Slide 60 text

Chebyshev Approximation 那是什麼鬼？我用的到嗎？ （你可能已經用過了）

Slide 61

Slide 61 text

Polynomial Approximation

Slide 62

Slide 62 text

目標：sin ( x ) 多項式逼近：a n xn + a n-1 xn-1 + … + a 1 x1 + a 0

Slide 63

Slide 63 text

取點

Slide 64

Slide 64 text

No content

Slide 65

Slide 65 text

No content

Slide 66

Slide 66 text

No content

Slide 67

Slide 67 text

No content

Slide 68

Slide 68 text

No content

Slide 69

Slide 69 text

N 個點 可以解一元 N-1 次方程式 ( N-1 degree )

Slide 70

Slide 70 text

No content

Slide 71

Slide 71 text

No content

Slide 72

Slide 72 text

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) 方程式係數

Slide 73

Slide 73 text

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; 方程式係數

Slide 74

Slide 74 text

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 */ 方程式係數

Slide 75

Slide 75 text

自己用 SymPy 寫 Chebyshev Approximation

Slide 76

Slide 76 text

No content

Slide 77

Slide 77 text

轉成可執行 Function 9

Slide 78

Slide 78 text

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

Slide 79

Slide 79 text

lambdify ( x, x**x )

Slide 80

Slide 80 text

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

Slide 81

Slide 81 text

lambdify ( x, x**x, "math" ) >

Slide 82

Slide 82 text

lambdify ( x, x**x, "numpy" ) >

Slide 83

Slide 83 text

lambdify ( x, x**x, "mpmath" ) >

Slide 84

Slide 84 text

lambdify ( x, x**x, "numexpr" ) >

Slide 85

Slide 85 text

lambdify ( x, x**x, "sympy" ) >

Slide 86

Slide 86 text

# NumPy universal function ufuncify ( args, expr, backend="numpy" ) from sympy.utilities.autowrap import ufuncify

Slide 87

Slide 87 text

# Python & NumPy's C API # and C compiler ufuncify ( x, x**x )

Slide 88

Slide 88 text

# Cython ufuncify ( x, x**x, backend="cython")

Slide 89

Slide 89 text

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

Slide 90

Slide 90 text

# Theano theano_function ( args, exprs ) from sympy.printing.theanocode import theano_function

Slide 91

Slide 91 text

theano_function ( [x], [x**x] )

Slide 92

Slide 92 text

# broadcast and type theano_function ( [x], [x**x], dims={x: 1}, dtype={x: 'float64'})

Slide 93

Slide 93 text

Library 內建的 Chebyshev Approximation

Slide 94

Slide 94 text

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]

Slide 95

Slide 95 text

poly, err = chebyfit ( sin, [0, pi/4], 10, error=True) from mpmath import chebyfit, sin, pi err: 1.94337098561508e-14

Slide 96

Slide 96 text

找完逼近的方程式了 現在想要讓執行變快 怎麼辦？

Slide 97

Slide 97 text

NEED FOR SPEED

Slide 98

Slide 98 text

Bob：「肯定是 Python 直譯器不夠快， 我們來換一個實作吧」

Slide 99

Slide 99 text

Python 直譯器實作 ● CPython ( 官方 ) ● PyPy －使用 Python 的子集來實作，採用 JIT 與不同的 Garbage Collection 技術 ● Pyston －基於 LLVM 實作 JIT ● Jython －把 Python code 轉成 Java bytecode 跑在 JVM 上 ● IronPython －用 .NET Framework 實作 ● …

Slide 100

Slide 100 text

PyPy pypy.org ● 專案狀態 ○ Python 2.7 有支援 ○ Python 3 目前指支援到 Python 3.3 (alpha) ● 效能比較可以看 speed.pypy.org

Slide 101

Slide 101 text

Pyston github.com/dropbox/pyston ● 由 Dropbox 於 2013 發起的專案 ● 使用 C++11 撰寫 ● 基於 LLVM 實作 JIT ● 目前方向著重在 Python 2.7

Slide 102

Slide 102 text

CPython：「你這樣講就不對了， 怎麼可以怪我呢？」

Slide 103

Slide 103 text

Alice：「不要急著怪別人」

Slide 104

Slide 104 text

Alice：「其實你不需要 SymPy 了啊， 為什麼不寫成其他形式？」

Slide 105

Slide 105 text

自己重寫一遍好麻煩啊

Slide 106

Slide 106 text

SymPy codegen

Slide 107

Slide 107 text

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

Slide 108

Slide 108 text

要速度？何不試試編譯式語言？

Slide 109

Slide 109 text

No content

Slide 110

Slide 110 text

C

Slide 111

Slide 111 text

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

Slide 112

Slide 112 text

codegen ( ("func", x**x), "c", header=False) #include "func.h" #include double func(double x) { double func_result; func_result = pow(x, x); return func_result; }

Slide 113

Slide 113 text

codegen ( ("func", x**x), "c", header=False) #ifndef PROJECT__FUNC__H #define PROJECT__FUNC__H double func(double x); #endif

Slide 114

Slide 114 text

Fortran

Slide 115

Slide 115 text

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

Slide 116

Slide 116 text

喜歡類似 MATLAB 的語言？ Octave

Slide 117

Slide 117 text

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

Slide 118

Slide 118 text

試試針對科學運算的新語言？ Julia

Slide 119

Slide 119 text

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

Slide 120

Slide 120 text

藏在 codegen 背後的重要資料－AST ( Expression Tree )

Slide 121

Slide 121 text

SymPy Expression Tree 把 SymPy 的 Expression Tree 轉成 DOT 語言描述 （接著可以用 Graphviz 這類工具轉成圖檔） from sympy.printing.dot import dotprint

Slide 122

Slide 122 text

x+2 SymPy Expression Tree

Slide 123

Slide 123 text

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

Slide 124

Slide 124 text

x / (y + 42) SymPy Expression Tree

Slide 125

Slide 125 text

優化！ Tree 操作！

Slide 126

Slide 126 text

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

Slide 127

Slide 127 text

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

Slide 128

Slide 128 text

Visitor Pattern

Slide 129

Slide 129 text

_print_

Slide 130

Slide 130 text

_print ( expr ) _print_Rational ( expr ) _print_Number ( expr ) _print_AtomicExpr ( expr ) _print_Atom ( expr ) _print_Basic ( expr )

Slide 131

Slide 131 text

SymPy 運作流程圖（總結）

Slide 132

Slide 132 text

Symbol

Slide 133

Slide 133 text

Symbol Operation

Slide 134

Slide 134 text

Symbol Symbol Operation Operation Symbol Operation Operation Symbol Symbol

Slide 135

Slide 135 text

Symbol Operation Operation Symbol Symbol Symbol Operation Symbol

Slide 136

Slide 136 text

Symbol Operation Symbol C Fortran Octave Julia DOT LaTeX JavaScript NumPy Theano Numexpr Python ...

Slide 137

Slide 137 text

C C compiler ufuncify Fortran Fortran compiler F2PY Cython Python & Numpy C API module.so

Slide 138

Slide 138 text

SymPy 內其他 Modules

Slide 139

Slide 139 text

SymPy 內其他 Modules ● Calculus ● Physics ● Lie Algebra ● Linear Algebra ● Number Theory ● Differential Geometry ● Set Theory ● ODE ● PDE ● Statistics ● Logic ● …

Slide 140

Slide 140 text

其他相關專案

Slide 141

Slide 141 text

SymPy Gamma

Slide 142

Slide 142 text

SymPy Gamma gamma.sympy.org

Slide 143

Slide 143 text

SymEngine github.com/symengine/symengine

Slide 144

Slide 144 text

SymEngine github.com/symengine/symengine （C++ 撰寫的 SymPy，為了提升效能）

Slide 145

Slide 145 text

基於 SymPy 的專案

Slide 146

Slide 146 text

基於 SymPy 的專案 ● PyDy - Multibody Dynamics with Python ● SageMath - 有把 SymPy 整合進去 ● SfePy - Simple Finite Elements in Python ● …

Slide 147

Slide 147 text

本投影片將放在 speakerdeck.com/wdv4758h/first-try-for-cas- sympy-with-codegen

Slide 148

Slide 148 text

Thank you ! ! !

Slide 149

Slide 149 text

Questions ?