#!/usr/bin/env python
# coding: utf-8

# # Jupyter 科學計算實驗室
# 
# ![Python](images/python2.png)
# 
# 
# #### 蔡炎龍
# #### 政治大學應用數學系
# 
# ### 2015 台南大學應用數學系

# # 本次演講內容在這裡
# 
# ![演講講義](images/url.png)

# # 1. Python 這個程式語言

# ## 創始人是 Guido van Rossum
# 
# ![gudio](images/gudio.png)
# 
# * 1989 年聖誕節宅在家裡的作品
# * 簡單易學、但功能強大

# ## Python 粉絲團
# 
# * Youtube, Blender (蜘蛛人), NASA, 全美各大學
# * 數據分析 Big 3!

# ## 我也加入粉絲團
# 
# 2000 年的時候, 我特別「贊助」了 Python 程式語言。
# 
# ![psa](images/psa.png)

# ## 2015 年 TIOBE Index
# 
# ![TIOBE](images/TIOBE.png)

# # 2. Python 有很多套件管理...
# 
# 好像你可以自由選, 結果弄得有些麻煩!

# ## 大救星 Anaconda
# 
# 很多套件, 尤其科學計算的都裝好了!
# 
# [Anaconda 下載點](https://www.continuum.io/downloads)

# ## 於是你只要打入神秘的
# 
#     ipython notebook
#     
# 就可以用瀏覽器打開 Jupyter!
# 
# [IPython 的安裝與執行 (影片)](https://youtu.be/WB8hg6UZeY0)

# # 3. 等等, Jupyter?

# ## 原本叫 IPython
# 

# ## 原創者是 Fernado Perez 博士
# 
# [Fernado 在台灣介紹 IPython 的影片](http://youtu.be/ZD5n9s8PGtI)

# ## 你可以一邊寫程式、一邊做筆記, 甚至簡報也不用特別做!!

# ## 因為 IPython 太好用, 其他程式語也也想用這個介面!

# ## 於是這個界面改叫 Jupyter!
# 
# #### 名稱來自 Julia, Python, R; 當然還有木星!

# # 4. It's Magic!
# 
# #### 沒想到 Jupyter 這麼炫

# ## 就這樣
# 
# 就讀入基本的科學計算套件 (包括 `numpy`, `matplotlib`)

# In[1]:


get_ipython().run_line_magic('pylab', 'inline')


# ## 我們來畫串資料
# 
# 以台南每月平均溫度為例。

# In[2]:


temp = [17.6, 18.6, 21.2, 24.5, 27.2, 28.5, 29.2, 
        28.8, 28.1, 26.1, 22.8, 19.1]


# In[3]:


plot(temp)


# ## plot 的基本用法
# 
#     plot(x 的座標串列, y　的座標串列)

# ## 偷懶法
# 
# 我們也可省略 x 座標串列, 此時內設 x = [0, 1, 2, ..., n], n 是 y 的長度。
# 
#     plot(y　的座標串列)

# ## 自動生串列
# 
# 因為 Python 的串列太重要, 我們不時要自動生成一下。

# In[4]:


mylist = range(12)


# In[5]:


mylist


# ## 從 1 開始的串列

# In[6]:


x = range(1, 13)


# In[7]:


x


# ## 再度畫出來

# In[8]:


plot(x, temp)


# ## 長條圖也可以

# In[9]:


bar(x, temp)


# ## 自動生資料
# 
# 我們什麼資料都不餵, 當場生出來!

# In[10]:


plot(randn(100))


# # 5. Python 是真的可愛
# 
# 這次純耍寶

# ![xkcd Python](http://imgs.xkcd.com/comics/python.png)
# 
# from <a href="http://xkcd.com/">xkcd</a>

# ## 進入 xkcd 模式

# In[11]:


plt.xkcd()


# ## 可愛的事發生了

# In[12]:


plot(x, temp)


# ## 長條圖也是

# In[13]:


bar(x, temp)


# #### 這裡只是讓一切回復...

# In[14]:


rcdefaults()
get_ipython().run_line_magic('pylab', 'inline')


# # 6. 做筆記好簡單!
# 
# #### Markdown + LaTeX 都可以!!

# ## Markdown
# 
# 快速打出有格式的文件

# ## 標題和引文
# 
# ### headings
# 
#     # 大標題
#     ## 次標題
#     
# ### blockquotes
# 
#     > 引用文章
#     
# > 引用文章

# ## 條列
# 
# ```list
# * item 1
# * item 2
# * item 3```
# 
# * item 1
# * item 2
# * item 3

# ## 程式碼
# 
#     行內的程式碼, 比如說 `cd` 這個指令。
#     
# 行內的程式碼, 比如說 `cd` 這個指令。
# 
# 空格可以照你的排列。
# 
# ```
#     for i in range(10):
#         print(i)
# ```
# 
# 結果:
# 
#     for i in range(10):
#         print(i)
#         
# 還可以根據程式語言高亮化!!
# 
# ```python
#     for i in range(10):
#         print(i)
# ```

# ## 超連結和圖形
# 
# ### 超連結是這樣
# 
#     [web name](http://foo.com.tw)
# 
# ### 圖形是這樣
# 
#     ![image name](filename)

# ## 基本上可以在網頁呈現的都可以
# 
# ![示範 gif](images/symm.gif)

# ## 內建支援 LaTeX!
# 
# 和我們的展示、隨文用法一樣。比如說
# 
# ```latex
# $$f(x) = \int_a^b f(x) \, dx$$
# ```
# 
# 結果是:
# 
# $$f(x) = \int_a^b f(x) \, dx$$

# ## 離線使用 LaTeX
# 
# LaTeX 事實上是用 MathJax 網路服務來做的, 要離線使用 (一台電腦做一次即可), 可以這樣下指令:
# 
# ```python
# from IPython.external.mathjax import install_mathjax
# install_mathjax()
# ```

# # 7. Hello, World!
# 
# Python 的字串

# ## 我們可以這樣 Hello, World!

# In[15]:


print("hello, world")


# ## 用字串來

# In[16]:


message = "hello, world"
print message


# ## 分解字串
# 
# 把資料做各種拆解是科學計學分析的根本。

# In[17]:


message.split(",")


# ## 位置這樣算
# 
# ![hello world](images/hello.png)

# ## 例如第 0 個位置

# In[18]:


message[0]


# ## 負的也可以

# In[19]:


message[-5:-1]


# ## 練習
# 
# 給一個字串, 把所有字母從後面到前面列出。

# In[20]:


message[::-1]


# # 8. list 是我們分析的基本結構
# 
# 我們可以想成 list 是一串的資料, 在 numpy 叫 array, 在 pandas 叫 series, 雖然嚴格說它們有點不一樣。

# ## 前面用過的生資料方法

# In[21]:


range(10)


# ## Python 的迴圈基本上用 list 達成

# In[22]:


for i in range(10):
    print " "*randint(1,15) + "*"


# ## 用 `for` 迴圈創造 list

# In[23]:


L = []
for i in range(6):
    L.append(2*i)


# In[24]:


L


# ## Python 風一行文!

# In[25]:


S = [2*i for i in range(6)]


# In[26]:


S


# # 9. 可以下系統指令
# 
# 如果你熟悉 Unix Like 指令列的動作, 在 Jupyter 中真的非常方便!

# ## 看目前資料夾中的檔案
# 
# 你會發現指接打就好!!

# In[27]:


ls


# ## 沒有直接支援的, 打 ! 就可以

# In[28]:


get_ipython().system('ifconfig en0')


# ## 更棒的是, Python 有在看...

# In[29]:


out = get_ipython().getoutput('ifconfig en0 | grep inet')


# In[30]:


out


# In[31]:


print "My IP is: " + out[1].split()[1]


# # 10. 存檔和執行
# 
# 我們打這麼高興, 可不可以把我們邊玩邊寫的存成一個 .py 檔?

# In[32]:


def love(n):
    return "love" * n


# In[33]:


love(3)


# In[34]:


print love(5)


# ## 存檔

# In[35]:


get_ipython().run_line_magic('save', 'spam.py 32 34')


# ## 執行

# In[36]:


get_ipython().run_line_magic('run', 'spam.py')


# # 11. Array 導向程式
# 
# 在科學計算、資料分析, array 型的思考方式, 常常讓程式快速方便!

# ## 雖然 array 很像 list

# In[37]:


iamlist = [1,3,-2,5]


# In[38]:


iamlist * 2


# In[39]:


iamarray = array([1,3,-2,5])


# In[40]:


iamarray * 2


# ## 這樣也可以!

# In[41]:


sin(iamarray)


# ## 用 array 來畫圖
# 
# 記得畫圖基本上是:
# 
#     plot(X, Y)

# ## 畫 sinc
# 
# 在 information theory/digital signal processing):
# 
# $$\mbox{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$$

# ## 先生出 x 的座標

# In[42]:


x = linspace(-10, 10, 1000)


# ## Array 式的畫圖!

# In[43]:


plot(x, sinc(x));


# ## 想標出函數值大於零的部份

# In[44]:


y = sinc(x)
plot(x, y);
plot(x[y>0], y[y>0], 'ro');


# ## 改變 array 的 shape

# In[45]:


x.shape = (20,50)


# In[46]:


imshow(x);


# ## 當然還是可以直接套用函數

# In[47]:


imshow(sinc(x));


# ## 延伸學習
# 
# Travis Oliphant (Numpy, SciPy, Anaconda 開發者), 介紹 "array-oriented" 的計算技巧。這段影片是 Travis 在 2012 年, 在第一次 PyCon Taiwan 的中演講。
# 
# [Large-scale array-oriented computing with Python (YouTube)](http://youtu.be/vrPRwUOt-7k)
# 
# [投影片](http://www.slideshare.net/pycontw/largescale-arrayoriented-computing-with-python)

# # 12. 網路資料神速抓
# 
# 我們來看個有名的網路爬蟲例子。

# ## Mac Benchmark...
# 
# [Mac Benchmark](https://browser.primatelabs.com/mac-benchmarks)

# ## requests package
# 
# 有名的 http request 套件。

# In[48]:


import requests


# In[49]:


res = requests.get("https://browser.primatelabs.com/mac-benchmarks")


# ## 看看結果
# 
# 基本上就是整個 html 的檔案。

# In[50]:


res.text


# # 13. 數據分析 pandas
# 
# 著名的數據分析套件

# In[51]:


import pandas as pd


# ## 讀進網頁資訊
# 
# 我們把網頁資訊讀入, pandas 會自己找出一個個表格!

# In[52]:


dfs = pd.read_html(res.text)


# In[53]:


dfs[3]


# ## 只想看機型和分數

# In[54]:


df = dfs[3][["Mac", "Score"]]


# In[55]:


df.head()


# ## 只看前十名成績!

# In[56]:


df.ix[:10].plot(kind='bar')


# ## 延伸學習
# 
# Julia Evans 是加拿大一位專業資料分析師, 她用實際的例子介紹 Pandas (還有 IPython!) 簡短有趣, 非常推薦!
# 
# [Diving into Open Data with IPython & Pandas (YouTube)](http://youtu.be/rEalbu8UGeo)

# # 14. 快速簡單 GUI 界面!!
# 
# 我們來介紹 IPython 2.0 起非常炫的互動功能。

# ## 互動開始!

# In[57]:


from IPython.html.widgets import interact


# ## 隨便來個函數
# 
# 然後就等著互動!!

# In[58]:


def f(x):
    print x


# ## interact!!

# In[59]:


interact(f, x=3);


# 有趣的是, 我們改變 x 的定法, 就可以有不同 GUI 界面! 試試:
# 
# * x=3.
# * x=(1,15)
# * x=u"你好"
# * x=[u"台南", u"台北"]
# * x={u"台南":"01", u"台北":"04"}

# ## 換個函數
# 
# 我們來個「認真點」的函數。

# In[60]:


def move(x):
    print " "*x + "ooo"


# ## 簡簡單單讓它動!

# In[61]:


interact(move, x=(0,40));


# # 15. 認真來數據分析

# ## 基本上就是學函數
# 
# 用數學的角度來說, 資料分析大概可以看成「找函數」的工作。
# 
#     f(x) = y

# ## [例子]
# 
# 輸入一個時間 x, 我就出現那個時間台南大學的溫度。
# 
#     f(x) = 在 x 時台南大學的溫度

# ## [例子]
# 
# 某公司做客戶分析, 把客戶分成三類:
# 
# * 1: 不用花什麼力氣就是我們忠實粉絲
# * 2: 有可能成為客戶但要努力說服
# * 3: 別傻了, 這人不會理我們的
# 
# 我們自然希望把某一個人 x 放進我們的神秘函數, 函數就告訴我們他是哪一類的人:
# 
#     f(x) = x 這個人是哪一類
#     
# 很棒吧? 等等, x 這個人電腦怎麼「輸入」? 所以我們會用某些可量化數據方式表式, 例如:
# 
#     x = {'血型':'O', '星座':'射手座', 'A產品':27, 
#     'B產品':0, ...}
#     
# 耶, 血型星座不是數字啊, 但我們可以用 1 表 A 型, 2 表 B 型, 等等。所以 x 總結可能是一個 list:
# 
#     x = [3, 1, 27, 0, ...]

# ## [例子]
# 
# 於是我們知道, 人可以「輸入」後, 世界又更寛廣了。比方說我們還可以做個「配對函數」, 輸入 x, y 兩個人, 得到:
# 
#     f(x,y) = x 和 y 的素配指數

# ## [型式]
# 
# 基本上我們學函數大約有兩種方式:
# 
# 1. regression (迴歸)
# 2. interpolation (插值)

# # 16. 假數據, 真迴歸

# ## 來學這個函數
# 
# 我們假設真實世界某個現象的「理想函數」長這樣:
# 
# $$f(x) = 1.2x + 0.8$$

# ## 模擬 50 個學習數據
# 
# 就像真實的世界, 我們加入 noise。

# In[62]:


x = linspace(0,5,50)
y = 1.2*x + 0.8 + randn(50)


# ## 看看我們的數據
# 
# 看看我們目前取的點, 附上「完全正確的」函數。

# In[63]:


scatter(x,y)
plot(x,1.2*x + 0.8,'b');


# ## 用 SciPy 來迴歸
# 
# 線性迴歸很多套件會做, 我們這裡介紹的方法好處是你只要預想好這些資料長得很像什麼函數, 我們都可以迴歸! 不一定要線性函數!
# 在我們線性的例子很容易。我們先做一個目標函數
# 
# $$f(x) = ax + b$$
# 
# 再來算 a, b 這樣。

# ## 使用 SciPy
# 
# 這裡用了 SciPy 曲線逼近的函數, 並且設定我們要的函數型式。

# In[64]:


from scipy.optimize import curve_fit


# In[65]:


def f(x, a, b):
    return a*x + b


# ## 真正逼近!
# 
# 我們「學習的資料」就是剛剛的 x, y。傳回的 popt 這裡就是最小平方法下得到的 a, b。

# In[66]:


popt, pcov = curve_fit(f, x, y)


# In[67]:


popt


# In[68]:


a, b = popt


# ## 畫出結果
# 
# 紅色是我們模擬出來的, 藍色是「理想狀況函數」。

# In[69]:


scatter(x,y)
plot(x,1.2*x + 0.8,'b',label="original")
plot(x, a*x + b, 'r',label="fit")
ax = gca()
ax.legend();


# # 17. 插值法

# ## 插值法通常要完全命中
# 
# 我們來準備個要完全命中的函數:
#     
# $$f(x) = \sin(x) + 0.5x$$

# ## 準備訓練數據
# 
# 我們準備 15 筆訓練數據, 但為了畫圖, 再做一個叫 `xx` 的 array, 包含 100 個 x 軸上的點。

# In[70]:


x = linspace(0, 10, 15)
y = sin(x) + 0.5 * x


# In[71]:


xx = linspace(0, 10, 100)


# ## 照例畫個圖

# In[72]:


scatter(x, y, s=50)
plot(xx, sin(xx) + 0.5*xx, 'b');


# ## 記得插值法要全部命中!!
# 
# 但其實全部命中其實不是什麼了不起的事...

# ## 做個函數展示

# In[73]:


def draw(a):
    scatter(x, y, s=50)
    xxx = linspace(0, 10, 1000)
    plot(xxx, 3.2*sin(a*xxx)+3, 'r');


# ## 全部命中...

# In[74]:


interact(draw, a=(1,100));


# ## 通常指定不要太多曲折
# 
# 我們準備再呼叫一下 SciPy 來做 "Cubic Spline" 的插值法!

# ## SciPy 的插值法套件

# In[75]:


from scipy.interpolate import interp1d


# ## 指定要 cubic!

# In[76]:


f = interp1d(x, y, kind='cubic')


# ## 這就學成了

# In[77]:


scatter(x, y, s=50)
plot(xx, sin(xx) + 0.5*xx, 'b') # 真正的
plot(xx, f(xx), 'r');


# # 18. 機器學習之 SVM

# ## 兩種機器學習方式
# 
# 在機器學習的分類 (其實也是一個學函數的過程) 大致有
# 
# * supervised learning (監督式學習)
# * unsupervised learning (非監督式學習)

# ## SVM 支持向量機
# 
# SVM (Support Vector Machines) 是用直線、曲線、平面、曲面等等把各類資料分開的方式。

# ## [例子]
# 
# 第一類:
# (-1, -1), (-2, -1)
# 
# 第二類:
# (1, 1), (2, 1)

# ## 教一下 Python
# 
# 告訴 Python 哪些點是第一類, 哪些是第二類。

# In[78]:


X = array([[-1, -1], [-2, -1], [1, 1], [2, 1]])
Y = array([1, 1, 2, 2])


# ## 畫個圖看看

# In[79]:


scatter(X[:,0], X[:,1],s=50,c=Y);


# ## 動用有名的 SciKit-leran!
# 
# 機器學習的重要套件

# In[80]:


from sklearn.svm import SVC


# ## 打開一個「學習機」
# 
# 打開一個「學習機」, 然後真的餵我們的訓練資料去學習!

# In[81]:


clf = SVC()


# In[82]:


clf.fit(X,Y)


# ## 「預測」沒學過的

# In[83]:


clf.predict([[-0.8,-1]])


# ## 餵一堆資料預測

# In[84]:


gd = array([[i,j] for i in arange(-2.5, 2.5, 0.2) 
 for j in arange(-2.5, 2.5, 0.2)])


# In[85]:


gdc = clf.predict(gd)


# ## 畫出來看看

# In[86]:


scatter(gd[:,0], gd[:,1], s=50, c=gdc);


# ## 生出像樣點的「假資料」
# 
# 剛剛的例子太簡單, 我們照例來生生假的訓練資料。sciket-leran 居然也想到要幫我們做這件事。

# In[87]:


from sklearn.datasets import make_classification


# ## 隨便分
# 
# 在平面上隨便取點分三類。

# In[88]:


X, Y = make_classification(n_features=2, n_redundant=0, n_informative=2, 
                           n_clusters_per_class=1, n_classes=3)


# ## 照例畫個圖

# In[89]:


scatter(X[:,0], X[:,1], s=50, c=Y);


# ## 再度開個學習機來學學

# In[90]:


clf = SVC()


# In[91]:


clf.fit(X,Y)


# ## 放一堆點預測

# In[92]:


gdc = clf.predict(gd)
scatter(gd[:,0], gd[:,1], s=50, c=gdc);


# # 19. 非監督式學習之 K-Mean

# ## 沒有標準答案的學習
# 
# 所謂非監督式學習就是「沒有標準解答」的學習。比如說們處理一些數據, 本來沒有分類, 讓電腦自動分類。
# 
# K-mean 是這種的方法。

# ## 亂數取 50 個點

# In[93]:


X = randn(50,2)


# In[94]:


scatter(X[:,0], X[:,1], s=55);


# ## SciKit-Learn 的 KMean
# 
# 做一個 K-Mean 學習機, 告訴 Python 我們想分三類。

# In[95]:


from sklearn.cluster import KMeans


# In[96]:


clf = KMeans(n_clusters=3)


# ## 開始分類

# In[97]:


clf.fit(X)


# In[98]:


clf.labels_


# ## 畫出結果

# In[99]:


scatter(X[:,0], X[:,1], s=55, c=clf.labels_)


# ## 當然沒在訓練資料中的也可以

# In[100]:


gdc = clf.predict(gd)
scatter(gd[:,0], gd[:,1], s=50, c=gdc);


# # 20. 用數學的方式算數學
# 
# 我們來做所謂的「符號式」運算。

# ## 符號式計算套件 SymPy

# In[101]:


from sympy import *


# ## 讓式子美化的魔術

# In[102]:


init_printing()


# ## 要先設什麼是變數
# 
# 我們要告訴 Python, 例如 x 是一個變數。

# In[103]:


x = symbols("x")


# ## 求反導數

# In[104]:


integrate(1/(x**2-1))


# ## 定積分

# In[105]:


integrate(1/(x**2-1), (x, 2,3))


# ## 秀式子

# In[106]:


Integral(1/(x**2-1), (x,2,3))


# ## 互動一下

# In[107]:


def ntaylor(n):
    expr = sin(x).series(x,0,n).removeO()
    display(expr)
    f = lambdify(x, expr, "numpy")
    a = linspace(-10,10,300)
    plt.plot(a, f(a));
    plt.ylim(-1,1);


# ## 開始!

# In[108]:


interact(ntaylor, n = (1,30));


# ![問題](images/question.png)