Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2016-07-19

MikuHatsune2016-07-19

球面調和関数をrPython でやる

球面調和関数をやりたいがRではいいものがないのでPythonでやった結果をrPython を使ってRでゴリ押し計算するやり方。R を投げ捨ててPython に移行したかったが、結局R でやっちゃう。

rPython について知りたければ下にスクロール。

 

球面上の離散的な分布をなんらかの関数級数展開できれば、関数化できる。球面上の分布fは、球座標表示でf(¥theta, ¥phi)と書ける。これを、ルジャンドル多項式

P_l(x)=¥frac{1}{2^n n!}¥frac{d^l}{dx^l}(x^2-1)^l

ルジャンドル陪関数

P_l^m(x)=(-1)^m(1-x^2)^{¥frac{m}{2}}¥frac{d^m}{dx^m}P_l(x)

を用いて、球面調和関数

Y_l^m(¥theta, ¥phi)=(-1)^{¥frac{m+|m|}{2}}¥sqrt{¥frac{2l+1}{4¥pi}¥frac{(l-|m|)!}{(l+|m|)!}}P_l^{|m|}(¥cos¥theta)exp(im¥phi)

と表すことができる。虚数が入っているので、実数版は

Y_l^m=¥{¥begin{matrix}¥sqrt{2}(-1)^mRe(Y_l^m)&if¥hspace{3}m>0¥¥Y_l^0&if¥hspace{3}m=0¥¥¥sqrt{2}(-1)^{-m}Im(Y_l^{-m})&if¥hspace{3}m<0¥end{matrix}

となる。

l はdegree, m はorder とよく書かれている。-l¥leq m ¥leq lである。

 

Rでは球面調和関数を直接やってくれそうな関数がなさそうなので、組み合わせてやってみる。orthopolynom パッケージはルジャンドル多項式を返す。結果は0 degree からl degree まで出力される。

library(orthopolynom)
l_spharm <- 5 # l; degree
m_spharm <- 0 # m; order
legp <- legendre.polynomials(l_spharm)
[[1]]
1 

[[2]]
x 

[[3]]
-0.5 + 1.5*x^2 

[[4]]
-1.5*x + 2.5*x^3 

[[5]]
0.375 - 3.75*x^2 + 4.375*x^4 

[[6]]
1.875*x - 8.75*x^3 + 7.875*x^5 

legendre.polynomials は直接微分してくれる。複数回の微分(m) はその都度loop を回す。

l_spharm <- 5 # l; degree
m_spharm <- 2 # m; order
legp <- legendre.polynomials(l_spharm)[l_spharm + 1]
for(mm in 1:m_spharm){
  legp <- polynomial.derivatives(legp)
}
[[1]]
-52.5*x + 157.5*x^3 

さてここで、具体的なルジャンドル関数を求めて確かめてみる。

P_2^2(x)=-3(1-x^2)なので、

l_spharm <- 2 # l; degree
m_spharm <- 2 # m; order
legp <- legendre.polynomials(l_spharm)[l_spharm + 1]
for(mm in 1:m_spharm){
  legp <- polynomial.derivatives(legp)
}
P_l_m <- function(x) sum((-1)^m_spharm * (1-x^2)^(m_spharm/2) *polynomial.coefficients(legp)[[1]]* x^(0:(l_spharm-m_spharm)))

x <- seq(-1, 1, length=100)
plot(x, mapply(P_l_m, x))
#lines(x, -3*x*sqrt(1-x^2)) # l <- 2; m <- 1
lines(x, 3*(1-x^2))

確かにあってる。

f:id:MikuHatsune:20160719163025p:image

 

本当なら、球面調和関数同士の積、内積クロネッカーデルタ¥delta_{l,l^{’}}になることを利用して、もとの球面上の関数f(¥theta,¥phi)

f(¥theta,¥phi)=¥sum_{l=0}^¥infty¥sum_{m=-l}^lc_l^mY_l^m(¥theta,¥phi)

c_l^m=¥int_{-{¥pi}}^{¥pi} d¥phi¥int_0^{¥pi} d¥theta ¥sin¥theta f(¥theta, ¥phi)Y_l^m(¥theta, ¥phi)

となる。

が、全然うまくRで再現できない。

上のリンクでわざわざPython コードがあるので、これを再利用することを考える。R からPython を利用するには、rPython パッケージを使う。

# これはPython
import numpy as np
from scipy.special import sph_harm as _sph_harm
def sph_harm(l, m, theta, phi):
    return _sph_harm(m, l, phi, theta)

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib.cm as cm

# R^2(u=θ, v=Φ)
u = np.linspace(0, np.pi, 80)
v = np.linspace(0, 2 * np.pi, 80)

# R^3
x = np.outer(np.sin(u), np.cos(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.cos(u), np.ones(v.size))

fig = plt.figure()
Lmax = 3
for l in range(Lmax + 1):
    for m in range(-l, l+1):
        # 球面調和関数の値を色で表す
        sph = np.abs(sph_harm(l, m, np.array(u, ndmin=2).T, np.array(v)))
# 正規化
#        if sph.max() == sph.min():
#            sph.fill(0.5)
#        else:
#            sph = (sph - sph.min()) / (sph.max() - sph.min())
        colors = cm.jet(sph) #cm.bwr
        ax = fig.add_subplot(Lmax + 1, 2*Lmax+1, l*(2*Lmax+1)+m+Lmax+1, projection="3d")
        ax.set_axis_off()
        surf = ax.plot_surface(x, y, z, facecolors=colors,
                rstride=2, cstride=2, linewidth=10, shade=False)

plt.savefig("sh_c.png", transparent=False)

rPython には4つの関数しかない。しかし、Python ができるのであればこれらの関数だけで十分強力な武器になる。

python.assign はPythonオブジェクトを定義する。

a <- 1:4
python.assign( "a", a )

python.exec はPython 上で演算を実行する。代入もできる。

python.exec( "b = len(a)" )

python.get はPython 上で定義されているオブジェクトを取り出す。ただし、Python でのarray 配列などはR に持ち込んだ時にJSON がウンタラカンタラ言われるので、Python 上でlist にしてR に吐き出すのがいいのかもしれない。

python.get("b")

python.call はPython 上の関数に、R から引数をいれて計算させる。python.assign でいちいちPython になげてやってもいいけれども、面倒だったらpython.exec でPython 上に好きに関数を定義してしまえばR から解決できる。

a <- 1:4
b <- 5:8
python.exec( "def concat(a,b): return a+b" )
python.call( "concat", a, b)

Pythonライブラリを持っていれば、実行できる。

python.exec( "import math" )
python.get( "math.pi" )

さて、RでできないことをPython でゴリ押ししよう。scipy 内に球面調和関数を計算してくれるやつがあるので、これをそのまま読み込む。

リンクではlとm, ¥theta¥phi が逆なのでゴニョゴニョしているが、R で書く行が増えるのでいい感じに帳尻を合わせる。

Python で1行のものはそのままpython.* で挟めばよい。

# R で
# rPython
python.exec("import numpy as np")
python.exec("from scipy.special import sph_harm as _sph_harm")

Python はインデント構造にうるさいので、タブインデントをひたすら横に書いた長い行をいれてもよいが、見栄えを大切にするなら改行コマンドをいれる。

元のコードが

# Python
def sph_harm(l, m, theta, phi):
  sph = np.abs(_sph_harm(m, l, np.array(phi), np.array(theta, ndmin=2).T))
  return [list(d0) for d0 in sph]

ゴリ押しすると

# R
python.exec("def sph_harm(l, m, theta, phi):\
  sph = np.abs(_sph_harm(m, l, np.array(phi), np.array(theta, ndmin=2).T))\
  return [list(d0) for d0 in sph]")

また、上の実部だけ取り出すやつも、元のコードが

# Python
def sph_harm(l, m, theta, phi):
  if m > 0:
      sph = np.sqrt(2) * (-1)**m * _sph_harm( m, l, np.array(phi), np.array(theta, ndmin=2).T).real
  if m == 0:
      sph = _sph_harm(m, l, np.array(phi), np.array(theta, ndmin=2).T).real
  else:
      sph = np.sqrt(2) * (-1)**m * _sph_harm(-m, l, np.array(phi), np.array(theta, ndmin=2).T).imag
  return [list(d0) for d0 in sph]

ゴリ押しすると

# R
python.exec("def sph_harm(l, m, theta, phi):\
  if m > 0:\
      sph = np.sqrt(2) * (-1)**m * _sph_harm( m, l, np.array(phi), np.array(theta, ndmin=2).T).real\
  if m == 0:\
      sph = _sph_harm(m, l, np.array(phi), np.array(theta, ndmin=2).T).real\
  else:\
      sph = np.sqrt(2) * (-1)**m * _sph_harm(-m, l, np.array(phi), np.array(theta, ndmin=2).T).imag\
  return [list(d0) for d0 in sph]")

これで、R上からm, l, ¥theta, ¥phi を与えればPython で球面調和関数を計算して、あたかもR でやってくれたようにできる関数ができたので、確かめる。

# R^2(u=θ, v=Φ)
# パラメータ
u <- seq(0, 180, length=180)*pi/180 # u = np.linspace(0, np.pi, 80)
v <- seq(0, 360, length=360)*pi/180 # v = np.linspace(0, 2 * np.pi, 60)

# プロット用の xyz
x <- c(outer(sin(u), cos(v)))
y <- c(outer(sin(u), sin(v)))
z <- c(outer(cos(u), rep(1, length(v))))
xyz <- cbind(x, y, z)

l <- 6
m <- 6

# ここで rPython の利用
# sph_harm は実数部分版
sph <- python.call("sph_harm", l, m, u, v) # u length list and v length vector
sph <- unlist(sph)
cuti <- 20
cols <- greenred(cuti)[cut(sph, quantile(sph, seq(0, 1, length=cuti + 1)))]

plot3d(xyz, col=cols)
zz <- t(polar2rect(abs(sph), t(expand.grid(v, u)[,2:1])))
plot3d(zz, col=cols)

f:id:MikuHatsune:20160719163024p:image

 

ubuntu14.04, R3.3.0, Python 2.7.1 なので他OS では未確認。

2016-07-14

MikuHatsune2016-07-14

球面上に均等な点を配置したい

球面上に均一な点を配置したい。約N個ならばこんな感じでやれるが、厳密にN 個おきたい。

実装が簡単なこれをやってみる。

 

そもそも、なぜ均一な点を球面上に配置したいかというと、基準点集合として扱えるからである。計算機では3次元データは離散的に持っており、ある球面上の点と別の物体の球面上の点を比較したいときに、両者で揃っていて、なおかつ有限個でありながら無数の点があれば離散でも連続っぽく扱えるため、比較がなんとかできるようになる。だから均一な点をたくさん発生させたい。

 

しかしながら、数学的に厳密な球面上の均一配置は不可能である。正多面体を考えて、その頂点に点を配置させれば、厳密に均一な点が配置できるが、正多面体は四、六、八、十二、二十に限られているので、たくさん均一な点を置きたいときには不便である。

球面上にランダムに点を発生させたときは、球面全体で見れば均一だが、局所的には疎だったり密だったりするのでこれまたよくない。メルカトル図法っぽく格子を作っても、北南極ではグシャっと潰れるのでこれもこの部分の均一性がよろしくない。

というわけでらせんを使うらしい。

球面座標系f(¥theta, ¥phi, r=1)に全部でN個の均一な点を球面上に発生させる。ことを考える。螺旋のスタートとゴールは北極南極ということにする。すなわち、¥theta_{1}=¥pi, ¥theta_{N}=0であり、¥phi_1=0とする。k番目の点にたいして、

h_k=2¥frac{k-1}{N-1}-1

として、

¥theta_{k}=acos{h_k}

¥phi_{k}=¥phi_{k-1}+¥frac{3.6}{¥sqrt{N}}¥frac{1}{sqrt{1-{h_k}^2}}=¥phi_1+¥frac{3.6}{¥sqrt{N}}¥sum¥frac{1}{sqrt{1-{h_k}^2}}

が、均一な点の角度となる。

こうしても、厳密な均一ではない。というのも、北極南極の歪みがあるからである。これを適当に補正するために、¥{2,3,5,6,7¥}¥{N-1,N-2,N-4,N-5,N-6¥}の均一点ベクトルの平均をとって適当に変えてやる。

すると、緑の点が黄色になって、極での均一さがちょっとマシになる。

geometry パッケージのconvhulln を使えば、凸包をして簡単に球状オブジェクトができる。

f:id:MikuHatsune:20160713165743p:image

f:id:MikuHatsune:20160713165742p:image

f:id:MikuHatsune:20160713165741p:image

library(rgl)
library(SphericalCubature)
a <- seq(-90, 90, length=200)*pi/180
b <- seq(-180, 180, length=300)*pi/180
r <- 1
rad <- expand.grid(b, a)

GSS <- function(N=600, r=1){
  hk <- 2*((2:N)-1)/(N-1) - 1
  theta_k <- c(pi, acos(hk))
  phi_k <- 0
  for(k in 2:N){
    phi_k <- c(phi_k, tail(phi_k, 1) + 3.6/sqrt(N)/sqrt(1-hk[k-1]^2))
  }
  phi_k[N] <- phi_k[N-1] + 3.6/sqrt(N)
  # phi_k <- c(0, 3.6/sqrt(N)*mapply(function(z) sum(1/sqrt(1-replace(hk, N-1, 0)[1:z]^2)), 1:(N-1))) # slower
  gss <- t(polar2rect(rep(1, N), rbind(theta_k, phi_k)))
  gss[1,] <- colMeans(gss[c(2, 3, 5, 6, 7),])
  gss[N,] <- colMeans(gss[N-c(1, 2, 4, 5, 6),])
  return(gss)
}

xyz <- t(polar2rect(rep(1, nrow(rad)), t(rad)))
gss <- t(polar2rect(rep(1, N), rbind(theta_k, phi_k)))
cols <- replace(rep("red", N), c(1, N), "green")
plot3d(xyz)
spheres3d(gss, col=cols, radius=0.03)

convhulln(GSS())

2016-07-13

rmarkdwon で出力にGIF を入れたい

plot をfor loop にいれたままうっかりrender すると、出力がすべてプロットされて縦に長いhtml ができてしまう。

せっかくだからGIF 化しよう。

こちらこちらを参考に、fig.show="animate" を指定すればできる。

ただ、ubuntu ではffmpeg が必要だが、簡単にはapt-get できなくなったらしいので、ppa でパクってくる。

sudo add-apt-repository ppa:mc3man/trusty-media
sudo apt-get update
sudo apt-get install ffmpeg

2016-07-11

MikuHatsune2016-07-11

高次元での回転と疎行列

3次元物体を回転させたい。

2次元であれば回転行列は回転角¥thetaを用いてR=¥begin{pmatrix}¥cos¥theta&¥sin¥theta¥¥-¥sin¥theta&¥cos¥theta¥end{pmatrix} と書ける。

高次元での回転ならば、正方対角行列(対角成分は1)Rについて、(i, i), (i, j), (j, i), (j, j) 成分をひらすら上のRで置換していけば、ある軸での回転を表す。

特に3次元では、x軸まわりに¥alpha、y軸まわりに¥beta、z軸まわりに¥gamma 回転することを考えると、v=(x,y,z)を回転させて得られる座標V=(X,Y,Z)

^{t}V=¥begin{pmatrix}¥cos¥gamma&¥sin¥gamma&0¥¥-¥sin¥gamma&¥cos¥gamma&0¥¥0&0&1¥end{pmatrix}¥begin{pmatrix}¥cos¥beta&0&-¥sin¥beta¥¥0&1&0¥¥¥sin¥beta&0&¥cos¥beta¥end{pmatrix}¥begin{pmatrix}1&0&0¥¥0&¥cos¥alpha&¥sin¥alpha¥¥0&-¥sin¥alpha&¥cos¥alpha¥end{pmatrix}{^{t}v}

と表される。

f:id:MikuHatsune:20160711144447p:image

f:id:MikuHatsune:20160711144448p:image

 

これを疎行列でやる。疎行列自体の説明は飛ばして、実際に使ってみる

基本的な疎行列の形は、例えば対角成分だけ何か値があって他は0というような、

S=¥begin{pmatrix}s_1&0&0&0¥¥0&s_2&0&0¥¥0&0&s_3&0¥¥0&0&0&s_4¥end{pmatrix}

というものがまず思いつくやつである。ここで、s_iも行列だと思えば、例えばこれがxyz 座標に相当するとして、s_i=¥begin{pmatrix}x_i&0&0¥¥0&y_i&0¥¥0&0&z_i¥end{pmatrix} と考えれば、

S=¥begin{pmatrix}x_1&0&0&&&&¥¥0&y_1&0&¥dots&&0&¥¥0&0&z_1&&&&¥¥&¥vdots&&¥ddots&&¥vdots&¥¥&&&&x_n&0&0¥¥&0&&¥dots&0&y_n&0¥¥&&&&0&0&z_n¥end{pmatrix}

というような対角ブロック行列ができる。

ブロックに相当する回転行列を持った、これまた疎行列を用意すれば、すべてが行列演算でできる。

しかし、xyz座標をベタに呼び出して計算したほうが速いっぽい…メモリサイズ的には疎行列のほうがよさそうなんだが。

rotate_polar <- function(xyz, alpha=pi/3, beta=pi/3, gamma=pi/3, sparse=FALSE){
  # angle alpha around x axis, beta y axis, and gamma z axis
  X <- matrix(c(cos(gamma), -sin(gamma), 0, sin(gamma), cos(gamma), 0, 0, 0, 1), nr=3)
  Y <- matrix(c(cos(beta), 0, sin(beta), 0, 1, 0, -sin(beta), 0, cos(beta)),     nr=3)
  Z <- matrix(c(1, 0, 0, 0, cos(alpha), -sin(alpha), 0, sin(alpha), cos(alpha)), nr=3)
  if (!sparse){
    res <- t(mapply(function(w) X %*% Y %*% Z %*% xyz[w,], seq(nrow(xyz))))
  } else {
    ij <- seq(nrow(g1)*3)
    M <- sparseMatrix(i=ij, j=ij, x=c(t(g1)))
    zeroX <- sparseMatrix(i=rep(1:(nrow(g1)*3), each=3), j=rep(1:(nrow(g1)*3), 3), x=as.numeric(X))
    zeroY <- sparseMatrix(i=rep(1:(nrow(g1)*3), each=3), j=rep(1:(nrow(g1)*3), 3), x=as.numeric(Y))
    zeroZ <- sparseMatrix(i=rep(1:(nrow(g1)*3), each=3), j=rep(1:(nrow(g1)*3), 3), x=as.numeric(Z))
    res <- matrix(rowSums(zeroX %*% zeroY %*% zeroZ %*% M), nc=3, byrow=TRUE)
  }
  return(res)
}
xyz <- matrix(rnorm(1000*3), nc=3)
# 疎行列にしてみる
system.time(rotate_polar(xyz, sparse=TRUE))
   ユーザ   システム       経過  
     0.257      0.000      0.258 
# 疎行列を使わずベタにfor loop で呼び出す
system.time(rotate_polar(xyz, sparse=FALSE))
   ユーザ   システム       経過  
     0.001      0.001      0.003 

2016-07-08

rmarkdown とknitr による文書作成

markdown 形式のrmarkdown とknitr を使って、授業資料の作成および解析結果の報告、宿題の提出などやっている。

HTML にすることでrgl の3D プロットなども出せるので、3D データを扱っているときにもよさげ。

真面目にやるならpandoc をインストールしてrender して…というので初心者は挫折する可能性があるが、Rstudio ならばknitr ボタンをポチるだけでできる。

 

いくつかハマった箇所。

いい感じの見出しや目次をつけたいとき、一番上に書いておけばいい感じになる。

---
title: "Title"
author: "Author name"
date: "Any date such as created or presented date"
output:
  html_document:
    toc: true
    toc_depth: 4 #デフォルト値は2
    number_section: true #trueでナンバリング
---
Created at `r Sys.time()`.

 

rmarkdownのchunk オプションをまとめて書いておけば一括で変更できる。


```{r echo=FALSE, warning=FALSE, message=FALSE, comment=""}
knitr::opts_chunk$set(echo=TRUE, eval=FALSE, warning=FALSE, message=FALSE, comment="")
```

 

本文自体はmarkdown 記法にしたがってかく。例えば見出しにリンク付きハイライトを書いたとしても、目次では目次番号がHTML 内リンク、見出し文字がリンク先へリンク、となっている。

 

rgl は少し挙動がオカシイときがある。こちらを参考にオプションをセットし

```{r setup, eval=TRUE, echo=FALSE}
library(rgl)
knitr::knit_hooks$set(rgl = hook_webgl)
```

実際にplot3d を出力する部分はrgl=TRUE を設定すればよい。

```{r rgl=TRUE, eval=TRUE}
plot3d()
```

しかし、ubuntu が悪いのかrgl するデータ容量が大きすぎるのが悪いのか、chrome が悪いのか、結構な頻度で

You must enable Javascript to view this page properly.

と言われてhtml 上、表示されないことがある。windows に持っていくと見れることもある。

 

というわけでサンプル。.Rmd として保存して、render するかRstudio でknitr ボタンをポチる。

---
title: "Matrix and Gram-Schmidt orthonormalization process"
author: "@Med_KU"
date: "20160705"
output:
  html_document:
    toc: true
    toc_depth: 4 #デフォルト値は2
    number_section: true #trueでナンバリング
---
Created at `r Sys.time()`.

```{r echo=FALSE, warning=FALSE, message=FALSE, comment=""}
knitr::opts_chunk$set(echo=TRUE, eval=FALSE, warning=FALSE, message=FALSE, comment="")
```
```{r setup, eval=TRUE, echo=FALSE}
library(rgl)
knitr::knit_hooks$set(rgl = hook_webgl)
```
# [vector](https://en.wikipedia.org/wiki/Vector_space)
Vector is mathematically a set of [scalar](https://en.wikipedia.org/wiki/Scalar_(mathematics)). Vector has scalar values and direction.

R keeps data by *vector* structure. This is an essential form. This is defined by `c`.
The data made in R is called *object*.
```{r}
dat <- c(1, 2, 10, 5, 9, 291)
dat
```
You can access to i th element in the *vactor* by `dat[i]`.
```{r}
dat[4]
```
You can modify the *vector object* by adding or reducing its elemetns.
```{r}
c(dat, 4, 1, 5, 12312)
dat[-4]
```

You can sum/diff/multiply/divide the *vector* or each element between *vector*s.
```{r}
# define another vector
dat1 <- c(3, 13, 64, 11, 0, 4)
sum(dat)
dat[3] + dat[1]
dat + dat1

dat[1] - dat[2]
diff(dat)
dat - dat1

dat[1] * dat[5]
prod(dat)
dat * dat1

dat/dat1
```

# [Matrix](https://en.wikipedia.org/wiki/Matrix_(mathematics) "matrix") 

*Matrix* is a 2 dimensional data structure. *Matrix* keeps data as N length vectors by M columns.
This is called N by M (dimension) *matrix*. *Matrix object* is defined by `matrix`.
```{r}
mat <- matrix(dat, nr=2, nc=3)
mat
```
We can access the cell (or element) in the matrix by indicating indices. If you want to take the cell of second row and third column, write `mat[2, 3]`.
```{r eval=TRUE, echo=FALSE}
hoge <- matrix(c("", "", "↓", "→", "→", "i, j"), nr=2, byrow=TRUE)
```
```{r eval=TRUE, echo=TRUE, comment=""}
hoge
```

**Warning!!** Elements in a matrix object should be the same class. In upper case, "", arrows, and 'i, j' are *character*. If there are various type of objects, object type will be forced to transform into *complex* > *character* > *float* > *integer* automatically. To avoid this problem, *data.frame* is sometime used.

Arithmetic operators for *matrix* are two ways. One is an operator for each element. Like arithmetic operators for *vector* in upper cases, simple description of `+`, `-`, `*`, and `/` return results of each element.
That is, there are two matrices $A=\begin{pmatrix}a_{1,1}&a_{1,2}\\a_{2,1}&a_{2,2}\end{pmatrix}$ and $B=\begin{pmatrix}b_{1,1}&b_{1,2}\\b_{2,1}&b_{2,2}\end{pmatrix}$, `+` operator returns $A+B=\begin{pmatrix}a_{1,1}+b_{1,1}&a_{1,2}+b_{1,2}\\a_{2,1}+b_{2,1}&a_{2,2}+b_{2,2}\end{pmatrix}$
```{r echo=FALSE}
mat1 <- matrix(1:4, nr=2)
mat2 <- matrix(5:8, nr=2)
```
```{r}
mat1 + mat2
mat1 * mat2
```

On the other hand, `%*%` operator returns matrix mutiplication $AB=\begin{pmatrix}a_{1,1}b_{1,1}+a_{1,2}b_{2,1}&a_{1,1}b_{1,2}+a_{1,2}b_{2,2}\\a_{2,1}b_{1,1}+a_{2,2}b_{2,1}&a_{2,1}b_{1,2}+a_{2,2}b_{2,2}\end{pmatrix}$.

**Optional**: `%*%` operator for *vector* returns inner product.

# [Gram-Schmidt orthonormalization process](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process)
Gram-Schmidt process is a method to orthonormalize a set of vectors in an inner product space. We generally suppose 3 dimensional space ($\mathbb{R}^3$) but it can be generalized in $\mathbb{R}^n$ space.

You replace a set of basis $\boldsymbol{v}=\{v_1, v_2, v_3\}$ with another set of basis $\boldsymbol{u}=\{u_1, u_2, u_3\}$. During the process, *normalization* and *orthogonalization* can be done.
Given *vector*s $\boldsymbol{x}$ and $\boldsymbol{y}$, these processes are performed as below.

***normalization***: *Normarization* means set length of *vector* (*norm*) 1. Length is difined as $|\boldsymbol{x}|$. That is, *normalization* returns $\boldsymbol{x}=\cfrac{\boldsymbol{x}}{|\boldsymbol{a}|}$.

***orthogonalization***: *Orthoganal* means crossing two *vector*s by 90 degree. Angle between *vector*s $\boldsymbol{x}$ and $\boldsymbol{y}$ is difined as $\theta$, and $\theta$ can be also defined by *inner product*

$cos{\theta}=\cfrac{\boldsymbol{x}\cdot \boldsymbol{y}}{|\boldsymbol{x}||\boldsymbol{y}|}$.

When $\theta=\pi$, $cos{\theta}=0$. That is $\boldsymbol{x}\cdot \boldsymbol{y}=0$.


Let's take a *vector* $\boldsymbol{a}$. This *vector* has three *vector*s `a1 <- c(1, 1, 1)`, `a2 <- c(1, -1, 2)`, and `a3 <- c(-1, 1, 3)` in itself. Data structure looks like *matrix*.
```{r eval=TRUE}
a <- cbind(c(1, 1, 1), c(1, -1, 2), c(-1, 1, 3))
```
```{r eval=TRUE, out.width=480, out.height=480, fig.width=7, fig.height=7}
# install.packages("scatterplot3d")
s <- scatterplot3d::scatterplot3d(t(cbind(0, a)), pch=16, xlab="x", ylab="y", zlab="z")
dt <- seq(0, 1, length=100)
for(i in seq(ncol(a))){
  s$points3d(t(outer(a[, i], dt)), cex=0.01, col=i+1)
}
```

```{r rgl=TRUE, eval=TRUE}
plot3d(t(cbind(0, a)), pch=16, xlab="x", ylab="y", zlab="z")
spheres3d(t(a), radius=0.05)
for(i in seq(ncol(a))){
  lines3d(rbind(0, a[, i]), lwd=2, col=i+1)
}
```

From *vector*s $a$, you want to get Gram-Schmidt *vector*s $u$. Gram-Schmidt process is performed as below.

1. STEP 1: $u_{1}=\cfrac{a_{1}}{|a_{1}|}$

2. STEP 2: $b_{2}=a_{2}-(a_{2}\cdot u_{1})u_{1}$, then $u_{2}=\cfrac{b_{2}}{|b_{2}|}$

3. STEP 3: $b_{3}=a_{3}-(a_{3}\cdot u_{1})u_{1}-(a_{3}\cdot u_{2})u_{2}$, then $u_{3}=\cfrac{b_{3}}{|b_{3}|}$

$b$ is a temporal objects.

Gram-Schmidt process is much major and famous algorithm. Someone has already implemented in R as a function.
[pracma](https://cran.r-project.org/web/packages/pracma/index.html) package has a function `gramSchmidt` to perform it. A precise result is $u_{1}=\cfrac{1}{\sqrt{3}}(1, 1, 1)$, $u_{2}=\cfrac{1}{\sqrt{42}}(1, -5, 4)$, and $u_{3}=\cfrac{1}{\sqrt{14}}(-3, 1, 2)$.
```{r echo=TRUE, eval=TRUE}
# install.packages("pracma")
library(pracma)
(a_g <- pracma::gramSchmidt(a))
```
```{r rgl=TRUE, eval=TRUE}
plot3d(t(cbind(0, a)), pch=16, xlab="x", ylab="y", zlab="z")
spheres3d(t(a), radius=0.05)
for(i in seq(ncol(a))){
  lines3d(rbind(0, a[, i]), lwd=2, col=i+1, lty=3)
}
spheres3d(t(a_g$Q), radius=0.05, col="yellow")
for(i in seq(ncol(a_g$Q))){
  lines3d(rbind(0, a_g$Q[, i]), cex=0.01, lwd=4, col=i+1)
}
```

# Assignment
1. Define any *vector* as you like.

2. Calculate mean of that *vector* by writing down the difinition of mean and using pre-implemented function in R.

3. Define 3 by 3 *matrix* $M$ as you like.

4. Calculate *matrix multiplication* $M^3$.

5. Perform Gram-Schmidt process of $M$. As shown in above, you need 3 steps.

6. Check whether STEP 5 is performed successfully, that is, compute *inner product* of the basis of Gram-Schmidt $M$ and length of it (This procedure needs multiplication and summation. These functions are in R.).