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

# # Table of Contents
#  <p><div class="lev1"><a href="#NumPy-Tutorial-6-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>NumPy Tutorial 6</a></div><div class="lev2"><a href="#Linear-Algebra（線形代数）-1.1"><span class="toc-item-num">1.1&nbsp;&nbsp;</span>Linear Algebra（線形代数）</a></div><div class="lev3"><a href="#Simple-Array-Operations（シンプルな配列演算）-1.1.1"><span class="toc-item-num">1.1.1&nbsp;&nbsp;</span>Simple Array Operations（シンプルな配列演算）</a></div><div class="lev3"><a href="#行列クラス-1.1.2"><span class="toc-item-num">1.1.2&nbsp;&nbsp;</span>行列クラス</a></div><div class="lev3"><a href="#添字アクセス：行列と2次元配列の比較-1.1.3"><span class="toc-item-num">1.1.3&nbsp;&nbsp;</span>添字アクセス：行列と2次元配列の比較</a></div><div class="lev3"><a href="#Vector-Stacking（ベクトルのスタック）-1.1.4"><span class="toc-item-num">1.1.4&nbsp;&nbsp;</span>Vector Stacking（ベクトルのスタック）</a></div><div class="lev3"><a href="#Histograms-1.1.5"><span class="toc-item-num">1.1.5&nbsp;&nbsp;</span>Histograms</a></div><div class="lev2"><a href="#参考リンク-1.2"><span class="toc-item-num">1.2&nbsp;&nbsp;</span>参考リンク</a></div>

# # NumPy Tutorial 6

# - [NumPy Tutorial メモ1(Basic)](http://nbviewer.jupyter.org/gist/Cartman0/9fa48b89664dc08ef82a55877767b5a0)
# - [NumPy Tutorial メモ2(ShapeManipulation)](http://nbviewer.jupyter.org/gist/Cartman0/bc062d1a162589e40a02c21e0ada7776)
# - [NumPy Tutorial メモ3(Copies and Views)](http://nbviewer.jupyter.org/gist/Cartman0/cd39c1bfc88f5edccd8f4523fbd33b91)
# - [NumPy Tutorial メモ4(Less Basic)](http://nbviewer.jupyter.org/gist/Cartman0/e9f50b0f0dbc28f40938b11088ac211d)
# - [NumPy Tutorial メモ5(Fancy indexing and index tricks)](http://nbviewer.jupyter.org/gist/Cartman0/5d24bdc3234a97bfc7c9fed4355a1580)

# ## Linear Algebra（線形代数）

# Work in progress. 
# Basic linear algebra to be included here.
# （基礎の線形代数がここに含まれる。）

# ### Simple Array Operations（シンプルな配列演算）

# See `linalg.py` in numpy folder for more.

# In[1]:


import numpy as np
a = np.array([
        [1.0, 2.0], 
        [3.0, 4.0]
    ])
print(a)


# In[2]:


a.transpose() # 転置


# In[3]:


a.transpose() == a.T # Tも同じ


# In[4]:


np.linalg.inv(a) # 逆行列


# $$
# A * A^{-1} = I
# $$
# 
# $
# A = \left(
#     \begin{array}{cc}
#       a & b \\
#       c & d 
#     \end{array}
#   \right)
# $ , $
# A^{-1} = \left(
#     \begin{array}{cc}
#       a^{-1} & b^{-1} \\
#       c^{-1} & d^{-1} 
#     \end{array}
#   \right)
# $ として、
# 
# $$
# \left(
#     \begin{array}{cc}
#       a & b \\
#       c & d 
#     \end{array}
#   \right)  
#   \left(
#     \begin{array}{cc}
#       a^{-1} & b^{-1} \\
#       c^{-1} & d^{-1} 
#     \end{array}
#   \right)
#   = 
#   \left(
#     \begin{array}{cc}
#       1 & 0 \\
#       0 & 1 
#     \end{array}
#   \right)
# $$
# 
# \begin{eqnarray}
# a * a^{-1} + b * c{^-1} = 1\\
# a * b^{-1} + b * d^{-1} = 0\\
# c * a^{-1} + d * c^{-1} = 0\\
# c * b^{-1} + d * d^{-1} = 1\\
# \end{eqnarray}
# 
# $$
# c * a^{-1} + d * c^{-1} = 0 \\
# c^{-1} = - \frac{c*a^{-1}}{d}
# $$
# 
# \begin{eqnarray}
# a * a^{-1} + b * c{^-1} &=& 1 \\
# a * a^{-1} - b * \frac{c*a^{-1}}{d} &=& 1 \\
# (a - \frac{bc}{d})a^{-1} &=&1\\
# a^{-1} &=& \frac{1}{(a - \frac{bc}{d})}\\
# &=& \frac{d}{(ad - bc)}\\
# \end{eqnarray}
# 
# \begin{eqnarray}
# c^{-1} &=& - \frac{c}{d}*\frac{d}{(ad - bc)}\\
# &=& - \frac{c}{(ad - bc)}\\
# \end{eqnarray}
# 
# \begin{eqnarray}
# c^{-1} &=& - \frac{c}{d}*\frac{d}{(ad - bc)}\\
# &=& - \frac{c}{(ad - bc)}\\
# \end{eqnarray}
# 
# \begin{eqnarray}
# a * b^{-1} + b * d^{-1} &=& 0 \\
# d^{-1} &=& -\frac{a * b^{-1}}{b}
# \end{eqnarray}
# 
# \begin{eqnarray}
# c * b^{-1} + d * d^{-1} &=& 1 \\
# c * b^{-1} + d * -\frac{a * b^{-1}}{b} &=& 1 \\
# (c - d\frac{a}{b})b^{-1} &=& 1\\
# b^{-1} &=& - \frac{1}{d\frac{a}{b} - c}\\
#  &=& - \frac{b}{ad - bc}
# \end{eqnarray}
# 
# \begin{eqnarray}
# d^{-1} &=& -\frac{a * b^{-1}}{b}\\
# &=& -\frac{a}{b} * - \frac{b}{ad - bc}\\
# &=& \frac{a}{ad - bc}\\
# \end{eqnarray}
# 
# $$
# A^{-1} = \frac{1}{ad-bc} \left(
#     \begin{array}{cc}
#       d & -b \\
#       -c & a 
#     \end{array}
#   \right)
# $$

# In[5]:


u = np.eye(2) # unit 2x2 matrix; "eye" represents "I" 単位行列
u


# In[6]:


j = np.array([
        [0.0, -1.0], 
        [1.0, 0.0]
    ])
j


# In[7]:


np.dot(j, j) # matrix product


# In[8]:


j.dot(j)


# In[9]:


np.trace(u) # trace


# In[10]:


y = np.array([
        [5.], 
        [7.]
    ])
y


# In[11]:


np.linalg.solve(a, y) # 2元連立方程式 A^-1 * y


# 固有値の計算
# 
# http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.linalg.eig.html

# In[12]:


eig_val, eig_vec = np.linalg.eig(j)
eig_val, eig_vec


# ### 行列クラス

# [numpy.matrix](http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.matrix.html)

# In[13]:


A = np.matrix('1.0 2.0; 3.0 4.0')
A


# In[14]:


type(A) # クラスを定義しているファイル


# In[15]:


A.T # 転置


# In[16]:


X = np.matrix('5.0 7.0')
X


# In[17]:


Y = X.T
Y


# In[18]:


print(A * Y)  # 行列の乗算


# In[19]:


print(A.I)  # 逆行列


# In[20]:


np.linalg.solve(A, Y)  # 線形方程式を解く


# ### 添字アクセス：行列と2次元配列の比較

# NumPy の配列と行列の間にはいくつかの重要な差異があることに注意。
# 
# NumPy は、
# - N次元配列オブジェクトと
# - ユニバーサル関数オブジェクト
# 
# という2種類の基本的なオブジェクトを用意している。
# その他のオブジェクトはこれら2つの上に構築されている。
# 
# 特に、行列はNumPyのarrayオブジェクトを継承した2次元配列オブジェクトである。
# 配列も行列も、添字は以下の1つもしくは複数の適切な組み合わせで構成されていなければならない：整数スカラー値、省略記号、整数または真偽値のリスト、整数または真偽値のタプル、整数または真偽値の1次元配列。
# 行列を行列のindex として用いることもできるが、普通は配列、リスト、あるいは他の形で事が足りる。
# Pythonでもそうであるように、index は0ベースである。
# 伝統的に我々は2次元配列や行列を行と列の長方形の配列として表す。
# 軸0に沿った動きは行を横切り、軸1に沿った動きは列を横切る。

# In[21]:


A = np.arange(12)
A


# In[22]:


A.shape = (3,4)
A


# In[23]:


M = np.mat(A.copy())
M


# In[24]:


print(type(A), ' ', type(M))


# In[25]:


print(A)


# In[26]:


print(M)


# ここで、いくつか単純なスライスを取ってみる。
# 基本的なスライシングにはスライスオブジェクトまたは整数の組を用いる。
# 
# 例えば、`A[:]` と `M[:]` の評価はPython の index付けに似たものに見えるのだが、重要なのはNumPy配列のスライシングはデータをコピーしないという点である。
# スライシングは同じデータに対し新しいview を提供する。

# In[27]:


print(A[:])


# In[28]:


print(A[:].shape)


# In[29]:


print(M[:]) 


# In[30]:


print(M[:].shape)


# Python のindex付けと異なるが、
# 複数の軸に沿って同時にindex付けする場合にカンマ区切りのindex が使える。

# In[31]:


print(A[:, 1])


# In[32]:


print(A[:, 1].shape)


# In[33]:


print(M[:, 1])


# In[34]:


print(M[:, 1].shape)


# 最後の2つの結果は異なる。
# - 2次元配列に対しコロンを1つ使うと1次元配列ができる
# - 一方、行列に対してだと2次元行列ができる。
# 
# 行列のスライスは常に行列を生み出す。
# 例えばスライス `M[2, :]` は `shape (1, 4)` の行列になるが、
# 対照的に、配列のスライスは常に可能な限り一番低い次元の配列になる。
# 
# 例えば `C` が3次元配列だとすると、
# - `C[..., 1]` は2次元配列に、
# - `C[1, :, 1]` は1次元配列になる。

# 配列の1番目と3番目の列が欲しいとしよう。
# リストを使ってスライスするのが1つの方法である：

# In[35]:


A


# In[36]:


A[:, [1, 3]]


# `take()` メソッドを使うというのもある：

# In[37]:


A.take([1, 3], axis=1)


# 最初の行(第0軸の0番目)を飛ばしたいなら以下のようになる：

# In[38]:


A[1:,].take([1, 3], axis=1)


# あるいは単に `A[1:, [1, 3]]` としてもよい。
# 上のをスライスするもう1つの方法は、`numpy.ix_` の利用がある：

# In[39]:


np.ix_((1, 2), (1, 3))


# In[40]:


A


# In[41]:


A[np.ix_((1, 2), (1, 3))]


# 1行目が1より大きい全ての列を保持したいととする。
# まずは真偽値のindex を作る方法：

# In[42]:


A[0, :] > 1


# しかし行列をindex付けするのはあまり便利ではない。

# In[43]:


M[0, :] > 1


# In[44]:


M[:, M[0, :] > 1]


# 問題としているのは、行列スライスをスライスすることで行列が生成されてしまう点である。
# しかし行列には便利な'A'属性があって、この値はその行列の配列表現になっているので、代わりにこうするだけで良い：

# In[45]:


M[:, M.A[0, :] > 1]


# 行列を条件によって2つの方向にスライスしたいのであれば、
# 戦略をわずかに補正しなければならない。

# In[46]:


A


# In[47]:


A[:, 0] > 2


# In[48]:


A[0, :] > 1


# In[49]:


A[A[:, 0] > 2, A[0, :] > 1]


# の代わりに、外積 `'ix_'` を用いる：

# In[50]:


np.ix_(A[:, 0] > 2, A[0, :] > 1)


# In[51]:


A[np.ix_(A[:, 0] > 2, A[0, :] > 1)]


# In[52]:


M[np.ix_(M.A[:, 0] > 2, M.A[0, :] > 1)]


# ### Vector Stacking（ベクトルのスタック）

# How do we construct a 2D array from a list of equally-sized row vectors? 
# （サイズの等しい行ベクトルのリストから2次元配列を構築するにはどうするか？）
# 
# In MATLAB this is quite easy（MATLABではこれはとても簡単）:
# if `x` and `y` are two vectors of the same length you only need do `m = [x; y]`. 
# （`x` と `y` の2つが同じ長さのベクトルなら、`m=[x;y]` とするだけで良い。）
# 
# In NumPy this works via the functions `column_stack`, 
# `dstack`, `hstack` and `vstack`, 
# depending on the dimension in which the stacking is to be done. 
# （NumPyでは、関数 `column_stack`, `dstack`, `hstack`, `vstack` で同じことができる、
# 関数の選択はスタックが行われる次元に依存する。）
# 
# For example:

# In[53]:


x = np.arange(0, 10, 2) # x=([0,2,4,6,8])
y = np.arange(5)# y=([0,1,2,3,4])


# In[54]:


m = np.vstack([x, y]) # m=([[0,2,4,6,8],
m                    #    [0,1,2,3,4]]), 2 * 5


# In[55]:


m.shape


# In[56]:


xy = np.hstack([x, y])   # xy =([0,2,4,6,8,0,1,2,3,4]) # 1 * 10
xy


# In[57]:


xy.shape


# See alse:
# 
# - [Numpy for Matlab users](https://docs.scipy.org/doc/numpy-dev/user/numpy-for-matlab-users.html)

# ### Histograms

# The NumPy `histogram` function applied to an array returns a pair of vectors（NumPy の`histogram`関数は、2つのベクトルを返す配列に適用される）： 
# the histogram of the array and the vector of bins. 
# （配列のヒストグラムと、瓶(bins)のベクトル、が返ってくる。）
# 
# Beware: `matplotlib` also has a function to build histograms (called `hist`, as in Matlab) that differs from the one in NumPy. 
# （注意：`matplotlib` にもヒストグラムを構築する関数がある
# （Matlabと同様 hist と呼ばれる）が、NumPyのものとは異なる。）
# 
# The main difference is that `pylab.hist` plots the histogram automatically, 
# while `numpy.histogram` only generates the data.
# （主な違いは、`pylab.hist` はヒストグラムを自動的に描画するが、
# `numpy.histogram` はデータを生成するところまでしか行わない点である。

# In[58]:


import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')

# Build a vector of 10000 normal deviates with variance 0.5^2 and mean 2

(mu, sigma) = 2, 0.5
v = np.random.normal(mu, sigma, 10000)
# Plot a normalized histogram with 50 bins
plt.hist(v, bins=50, normed=1)       # matplotlib version (plot)
plt.show()


# In[59]:


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

# Compute the histogram with numpy and then plot it
(n, bins) = np.histogram(v, bins=50, normed=True)  # NumPy version (no plot)
print('n ', n)
print('bins', bins)
print('bins[:-1]', bins[:-1])
plt.plot(.5 * (bins[:-1] + bins[1:]), n)


# ## 参考リンク

# - [QuickStartTutorial](https://docs.scipy.org/doc/numpy-dev/user/quickstart.html)
# - [私訳「暫定的 NumPy チュートリアル」](http://naoyat.hatenablog.jp/entry/2011/12/29/021414)