2020年12月31日木曜日

やり残し事項 2020

GEE コミュニティのチュートリアルを再度実施してみようと手を動かしました。
https://phreeqc.blogspot.com/2020/12/gee-sentinel-1.html

1か月でPart3まで増えていました。仮説検定を用いた判定がメイン。今回は任意の場所でも１～3まで走りました。が、結果はイマイチ。
これまで、IｎSAR、Change Detection ともに、良い結果を得られたことがありません。イマイチなんですよね。やはり、有償ソフトの出番でしょうか。この辺りは来年の短期目標にしましょう。

職種が変わったこともあり、地質関連の作業がなくなりました。帯磁率計などはやり残したままで終わりました。

あらためてやり残し事項を整理しておきましょう。

優先度高：粒子法コードのGPGPU対応
優先度高：SAR・画像の Change Detection
優先度低：PSIｎSAR・・・700万で購入可
優先度低：DtransuのGPU対応・・・100万で購入可
優先度低：地表流＋地下水・・・GSFLOW、HydroGeoSphere
優先度中：OpenMI・・・PHREEQCとDtransuの連成等

中期目標４年目は惨敗。新たな目標を考えましょう。

2020年12月30日水曜日

DataFrames.jl vs. pandas

Python は私のようなノンプロが手を出しやすい簡素な言語です。

一方、ライブラリへの依存が高いため、それらの互換性や version に応じて複数の仮想環境の構築を迫られることがあります。将来的には、作成していたコードがライブラリのVer.UP や開発中止で修正を余儀なくされることもあるでしょう（前者は既に経験済みです）。
個人的には精度や速度が必須となる科学技術計算を Fortran 、データ整理や可視化を Python、といったような使い分けをしています。現状維持で良いともいえるのですが、そろそろ代替え言語を試しておくべきと思い始めました（リスク管理です）。

先日、ある方から Julia を勧められました。Python と同じくらい容易な書き方で、計算が速く、コンパイルができるとのこと。
それは良い。

ということで、試してみました。
少し古いのですが、LTS の v1.0.5 + win10 （i5のラップトップ）です。インストールはプレコンパイル版を入れるだけ。そのままでは Jupyter での利用は不可でしたが、何度か試行し、動くようになりました。
今回の比較は csv 読み込みからの DataFrame 操作。頻繁に使う機能です（Juria と Python の比較ではなく、DataFrames.jl と pandas の比較になりました）。
結果は以下の通り。

csv読み込み	df=DataFrame(CSV.File(open(read,	1170万行	84.8s
	df6=DataFrame(CSV.File(open(read,	6万行	0.50s
	df=pd.read_csv	1170万行	51.4s
	df6=pd.read_csv	6万行	0.14s
カラム名変更	rename!(df,		0.21s
カラム名変更	df.rename(columns= ,inplace=True)		0.02s
inner_join	dfj=innerjoin(df, df6	290万行	15.5s
inner_join	dfj=pd.merge(df, df6	290万行	5.8s
型変換	dfj.＿=Date.(dfj.＿,"yyyy/mm/dd")		86.7s
型変換	dfj['__']=pd.to_datetime(dfj['__'])		1.4s
期間抽出	df0=dfj[dfj[:,"__"] .< _‗, :]	190万行	4.6s
期間抽出	df0=dfj[dfj['＿’] < ‗_]	190万行	3.6s
集約	df0_st=stack(df0,	380万行	5.6s
集約	df0_pv=pd.pivot_table(df0,values=	190万行	1.4s

pandas の圧勝でしたね。強い。

Julia でもパッケージを使うので Python と同じような問題は潜在的に有しているのでしょう。
イタレーションでは Julia の方が速くなるそうですので、扱う作業内容により言語を選択する形になるでしょうか？（今回の結果を見る限り、当面の出番はないかもしれませんが）
https://julialang.org/benchmarks/

引き続きその動向に注目しましょう。

2020年12月29日火曜日

常微分方程式

振動などの物理現象を解く際には、常微分方程式の解法が利用されています。

が、既に学んだことすら記憶にありません。
冬休み中に再読。以下にメモ。学ぶべき時に励んでおけばよかった。

*******************************

未知関数 y(x) と x および導関数 y'(x),y"(x),・・・yn(x) からなる方程式を n階常微分方程式という。
F(x,y(x),y'(x),y"(x),・・・yn(x))=0
このとき、y(x) を常微分方程式の解という。（方程式を解く＝式を求める）

常微分方程式の分類
・線形・非線形
・定数係数・変数係数
・同次・非同次

登坂宣好「微分方程式の解法と応用」より常微分方程式の例
線形‐変数係数‐非同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+x^2y\left(x\right)=e^x\end{align*}
線形‐定数係数‐非同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+3y\left(x\right)=e^x\end{align*}
線形‐変数係数‐同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+x^2y\left(x\right)=0\end{align*}
線形‐定数係数‐同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+3y\left(x\right)=0\end{align*}
非線形‐非同次形
\begin{align*}y\left(x\right)\frac{dy}{dx}\left(x\right)+3y\left(x\right)=x^3\end{align*}

解法（最近、使用した解法のみ記載）
線形-1階-定数係数-非同次形
\begin{align*}\displaystyle\frac{dy}{dx}\left(x\right)=f\left(x\right)\end{align*}
直接積分形：両辺をxで積分（両辺にdxを乗ずる形にして積分)
\begin{align*}f\left(x\right)=x : y\left(x\right)=\int xdx=\frac{1}{2}x^2+c\end{align*}\begin{align*}f\left(x\right)=sinx：y\left(x\right)=\int sinxdx=-cosx+c\end{align*}
線形-1階-定数係数-同次形
\begin{align*}\displaystyle\frac{dy}{dx}\left(x\right)=ay\left(x\right)\end{align*}
指数関数解:
\begin{align*}\frac{d}{dx}e^{\lambda x}=\lambda e^{\lambda x}=ae^{\lambda x}\end{align*}\begin{align*}\left(\lambda-a\right)e^{\lambda x}=0\rightarrow\lambda=a\end{align*}
解として\begin{align*}e^{ax}\end{align*}を有する。定数倍した\begin{align*}Ce^{ax}\end{align*}も解となる。（一般解）

線形-n階-定数係数-同次形
\begin{align*}\frac{d^ny}{dx^n}\left(x\right)+a_1\frac{d^{n-1}y}{dx^{n-1}}\left(x\right)+\ldots+a_ny\left(x\right)=0\end{align*}左辺に\begin{align*}y=C_1y_1+C_2y_2+\ldots+C_ny_n\end{align*}を代入すると\begin{align*}\frac{d^n}{dx^n}\left(C_1y_1+C_2y_2+\ldots\right)+a_1\frac{d^{n-1}}{dx^{n-1}}\left(C_1y_1+C_2y_2+\ldots\right)+a_n\left(C_1y_1+C_2y_2+\ldots\right)\end{align*}\begin{align*}{=C}_1\left({\frac{d^n}{dx^n}y}_1+a_1\frac{d^{n-1}}{dx^{n-1}}y_1+\ldots+a_ny_1\right)+C_2\left({\frac{d^n}{dx^n}y}_2+a_1\frac{d^{n-1}}{dx^{n-1}}y_2+\ldots+a_ny_2\right)+\ldots\end{align*}\begin{align*}{=C}_1・0+C_2・0+\ldots C_n・0=0\end{align*}よって\begin{align*}y=C_1y_1+C_2y_2+\ldots+C_ny_n\end{align*}は一般解となる。
解の形が\begin{align*}y=Ce^{\lambda x}\end{align*}のとき、
\begin{align*}C\ \left(\lambda^n+a_1\lambda^{n-1}+\ldots+a_n\right)e^{\lambda x}=0\end{align*}
C＝0のとき、y＝0は自明な解。\begin{align*}\lambda^n+a_1\lambda^{n-1}+a_2\lambda^{n-2}+\ldots+a_n=0 ：特性方程式\end{align*}のとき、非自明な解\begin{align*}y=C_1e^{α_1ｘ}+C_2e^{α_2ｘ}+⋯+C_ne^{α_nｘ}\end{align*}が存在。

非線形（特別な形式）
変数分離型解法：両辺をg(y)で割ってから積分(g(y)≠0)
\begin{align*}\displaystyle\frac{dy}{dx}\left(x\right)=f\left(x\right)g(y\left(x\right))\end{align*}\begin{align*}\frac{1}{g\left(y\left(x\right)\right)}dy=f\left(x\right)dx\end{align*}\begin{align*}\int\frac{1}{g\left(y\left(x\right)\right)}dy=\int f\left(x\right)dx\end{align*}

2020年12月27日日曜日

Mann-Whitney's U test

Mann-Whitney's U test を実施。

scipy.stats.mannwhitneyu を利用していたのですが、ランキングのつけ方が書かれていません。Σ の計算に n(n+1)/2 が含まれていますので、１位、２位、２位、４位のような同順位のつけ方はまずい。そこは正しく計算してくれているのだろうと思いながら、初めて使うので確認。

ソースでは rankdata が使用されています。

ranked = rankdata(np.concatenate((x, y)))

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rankdata.html?highlight=rank#scipy.stats.rankdata

>>> from scipy.stats import rankdata
>>> rankdata([0, 2, 3, 2])
array([ 1. ,  2.5,  4. ,  2.5])
>>> rankdata([0, 2, 3, 2], method='min')
array([ 1,  2,  4,  2])
>>> rankdata([0, 2, 3, 2], method='max')
array([ 1,  3,  4,  3])
>>> rankdata([0, 2, 3, 2], method='dense')
array([ 1,  2,  3,  2])
>>> rankdata([0, 2, 3, 2], method='ordinal')
array([ 1,  2,  4,  3])

はい、一番上の平均 (default is ‘average’) が利用されていました。OKです。

比較する母集団の平均、標準偏差、サンプル数の違いはどのように影響するのでしょうか？簡単な正規分布のプラス側を取り出してチェックしてみました（プロットは正規化済み）。シード値は設けず、何度かテストを繰り返して両側・有意水準5％にかかわる傾向を見てみました。

x1=np.random.normal(1, 10, 1000)
x2=np.random.normal(1, 10, 1000)
x1=x1[x1 > 0]
x2=x2[x2 > 0]

当然ですが、棄却されない分布です（稀に棄却されるのが気になりますが)。

2群が同じ平均・標準偏差で、一方のサンプルサイズを1000から100000に増やしたケース：検定結果に大きく影響せず（棄却され難い）。

2群が同じ平均・サンプルサイズで、一方の標準偏差を10から20に増やしたケース：検定結果に影響（棄却され易い）。

2群が同じ標準偏差・サンプルサイズで、一方の平均を1から２に増やしたケース：検定結果に大きく影響せず（棄却され難い）。

2群が同じ標準偏差・サンプルサイズで、一方の平均を1から3に増やしたケース：検定結果に影響（棄却され易い）。

2群が同じ標準偏差で、一方の平均を1から3、サンプルサイズを1000から100000に増やしたケース：検定結果に影響（棄却され易い）。

ランダムサンプリングをしているため、棄却・棄却されないという結果は固定になりません。それを繰り返し実行しただけですが、比較的容易に傾向が見えてきました。
平均2、３は見た目が微妙ですね。やや大きいかな？というような違いでも「有意差がないとは言えない」という傾向を出せるようです。人によってはやや大きい、あるいは同じように見える、などと意見が分かれるような分布ですが、ある考え方に従い再現性のある答えを出してくれるのはありがたい。

最近では、土研交流研究員の方が景観に配慮した護岸に関する類似性の評価に、この検定を利用されていたように思います。アイデア次第でいろいろな分野に適用できそうです。

好きではなかった統計的手法ですが、その合理性をもっと利用すべきでしょうとも考えるようになってきました。

2020年12月22日火曜日

shpfile in ZIP

先日、防災科研さんの地すべり地形分布図をプロットしようとして、ちょっとした問題に引っかかりました。

1．ZIP ファイルに５つの SHP ファイルが納められている。
2．ZIP ファイルは１５０個。
3．ZIP ファイルから、目的の SHP のみを解凍せずに取り出したい。
4．取り出した SHP を結合し、１つの SHP として保存したい。

QGIS で試みるも、複数の SHP からその一部を機械的に取り出せません。
ArcGIS Pro 2.5 では ZIP を扱えません。技があるのかとサポートに問い合わせましたが、やはり標準では対応しておらず、開発&有償になるとのこと（余談ですが、サポートさんのレスが以前に比べずいぶん遅くなりました）。

それほど大仰なモノでもないでしょう。
ということで、簡素にできたのがコチラ。

import os
import glob
import pandas as pd
import geopandas as gpd

path = './input/*.zip' #folder 'input'
zfiles = glob.glob(path) #zip files in the folder
gdf_300=gpd.GeoDataFrame()

for zfile in zfiles:
    fname=os.path.splitext(os.path.basename(zfile))[0]
    gdf_temp=gpd.read_file('zip://'+str(zfile)+'!ls300_'+fname+'.shp')
    gdf_300=pd.concat([gdf_300, gdf_temp])

gdf_300=gdf_300.reset_index(drop=True)
gdf_300.to_file('./output/300.shp', encoding='utf-8')

GIS のカジュアルユーザーだからでしょうか、Python との併用が合理的に感じます。

2020年12月17日木曜日

5次メッシュを3次メッシュに集約

AVS30（5次メッシュ）を 3次メッシュ内で平均化したいなあと思案。

全国だと EXCEL では行数が足りないので、正攻法の ArcGISで処理。

時間がかかるので Python でも追っかけ処理。
shpを読んだ後の主な処理はこれだけ。

#AVS30=0(水域)を削除
gdf=gdf[gdf['5th_CODE'] != 0]
#3次メッシュコード作成
gdf['3rd_CODE']=gdf['5th_CODE'].astype(np.int64)//100
#3次メッシュにおける平均を求める
df=gdf.groupby('3rd_CODE').mean()

空間演算が必要ないので、Python の方が早く終わりました。
便利。

2020年12月14日月曜日

scipy.optimize.curve_fit

任意関数の最小二乗法は EXCEL でも Pythonでも可能です。

精度、速度、記載の単純さで scipy.optimize.curve_fit の圧勝です。初期値によっては一発で求まらないこともありますが、 EXCELソルバーに比べると安定かつ高速で解を求めてくれます。
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

変分推論は難しいですね。収束させるのにあたりのつけ方が必要になるようです。で、結局は最小二乗法に戻ってしまったり。慣れが必要なのでしょう。

2020年12月13日日曜日

DataFrame のアダマール積

Pandas の DataFrame でよく使う便利な乗算。
間違えやすいので、( ..)φメモメモ

**************************************

import numpy as np
import pandas as pd

1. DataFrame 作成

df1=pd.DataFrame(np.ones(9).reshape(3, 3),
                 columns=['a','b','c'],
                 index=[1,2,3])
df1

	a	b	c
1	1.0	1.0	1.0
2	1.0	1.0	1.0
3	1.0	1.0	1.0

df2=pd.DataFrame(np.arange(9).reshape(3, 3),
                 columns=['a','b','c'],
                 index=[1,2,3])
df2

	a	b	c
1	0	1	2
2	3	4	5
3	6	7	8

2. アダマール積

#アダマール積
df1*df2

	a	b	c
1	0.0	1.0	2.0
2	3.0	4.0	5.0
3	6.0	7.0	8.0

index, column の名前、型が一致していないとダメ

df2.columns = ['d','e','f']
df1*df2

	a	b	c	d	e	f
1	NaN	NaN	NaN	NaN	NaN	NaN
2	NaN	NaN	NaN	NaN	NaN	NaN
3	NaN	NaN	NaN	NaN	NaN	NaN

df2.columns = ['a','b','c']
df2.index = [4,5,6]
df1*df2

	a	b	c
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN
4	NaN	NaN	NaN
5	NaN	NaN	NaN
6	NaN	NaN	NaN

#indexを文字列として入力
df2.index = ['1','2','3']
df1*df2

	a	b	c
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN

#indexを数値として入力
df2.index = [1,2,3]
df1*df2

	a	b	c
1	0.0	1.0	2.0
2	3.0	4.0	5.0
3	6.0	7.0	8.0

3. ブロードキャストの利用

dataframe * 行、列

#df1の各行に、df2の2行目をかける
df1*df2.iloc[1]

	a	b	c
1	3.0	4.0	5.0
2	3.0	4.0	5.0
3	3.0	4.0	5.0

# df1の各列に、df2のc列をかける・・・よりスマートな方法がありそうですが。
(df1.T*df2['c']).T

	a	b	c
1	2.0	2.0	2.0
2	5.0	5.0	5.0
3	8.0	8.0	8.0

4. スカラー倍

df1*df2.iloc[1,2]

	a	b	c
1	5.0	5.0	5.0
2	5.0	5.0	5.0
3	5.0	5.0	5.0

2020年12月7日月曜日

地下水と MCMC

ようやくMCMC を概ね理解し動かせるようになったので、地下水モデル10章の不確実性に関する部分へ少しだけ立ち戻ることに。
https://phreeqc.blogspot.com/2019/07/4.html

MCMC関連の文献を探してみると、簡単な計算をされている事例がありました。第2著者以降は日本の方ですね。
Julien Boulange, Hirozumi Watanabe, Shinpei Akai (2017) A Markov Chain Monte Carlo technique for parameter estimation and　inference in pesticide fate and transport modeling

気候、水収支、土壌、農薬特性などの40以上のパラメータから、農薬濃度の予測精度に大きな影響を与えることが報告されている以下の4パラメータを選択
農薬溶解速度（kdiss）、水田土壌中の農薬の一次分解速度（kbio）、農薬の脱離速度（kdes）、農薬分配係数（kd）
これらの事前確率分布を作成（最大・最小間で一様分布）
尤度は正規分布(E(θ)は観測値と予測値の平均2乗誤差)
MCMC で事後確率分布を作成

尤度を計算する前に採択されたパラメータでシミュを回す必要があります。大きなモデルを扱う実務向きではないですね。

実務で扱う地下水モデルで不確実性を検討するには、PEST一択なのでしょうか？確かに、市販のパッケージソフトには組み込まれており、使い勝手も良かった印象はあります。加えて、逆解析を古くから扱っているせいか、Bayes も良く出てきます。海外では当たり前なのかな？

もう少し、文献等を読んでみましょう。

2020年12月6日日曜日

change detection

Ohki et al.(2020)
Landslide detection in mountainous forest areas using polarimetry and interferometric coherence
https://earth-planets-space.springeropen.com/articles/10.1186/s40623-020-01191-5

ALOS-2 のデータ（Lバンド、空間解像度約６m)を用いた、土砂移動箇所検出の検討事例です。北海道胆振東部地震（2018年9月）と九州北部豪雨（2017年7月）を対象としています。
提案手法は change detection の範疇ですが、画像データを利用していません。波のデータのみです。そして最後はお決りの機械学習です（といっても決定木です）。

特徴

PolSAR、InSAR、DEM を組み合わせた解析

結果

すべての特徴量（InSAR、PolSAR、DEM）を用いた場合、kappa＞0.6程度の精度。
正確な検出には、少なくとも二重偏波（HHとHV）が必要。特に大きな入射角では四重偏波を推奨。
精度向上には、より高い空間分解能、より高い時間周波数、より多くの方向からの観測が必要である。

留意点

Hサイト（北海道）は、試験地周辺の火山の火砕流堆積物に覆われている。そのため、胆振地震では浅い土砂崩壊が多く、比較的緩やかな斜面(<30°)で発生したと考えられている。Hサイトを用いて構築したモデルは，急峻な地形での多発したFサイトへの適用性が悪い。地質の違いに留意する必要あり。

最近の国交省さんの動向を見ていますと、今後、土砂災害分野における各種衛星データの利用は拡大して行くようですね。関連マニュアルも複数出ていますし、ALOS-3へ向けてのシンポジウムでも以下のような報告がされています。
https://www.pco-prime.com/alos-3_alos-4_sympo/pdf/1505.pdf

機械学習を用いた change detection の精度は高いといえない状況ですが、いずれ向上し、災害時の人手不足問題も（change detection に関しては）回避できるでしょう。
検出アルゴリズムの構築以外は特に難しい作業はないので、備えておきましょう。

2020年12月5日土曜日

スパースモデリング

EXCELをよく使っていた頃、データに多項式をフィッティングする場面で何次式を使うかは感覚的でした。

低次だと合いませんし、高次だとノイズまで合わせてしまいます（いわゆる過学習）。ちょうどよいところを探す方法は、、、古くからありましたね。

LASSO
https://scikit-learn.org/stable/modules/linear_model.html#lasso
\begin{align*}
\min_{w}{ \frac{1}{2n_{\text{samples}}} ||y - X w||_2 ^ 2 + \alpha ||w||_1}\end{align*}

この損失関数はモデルの適合性を表す誤差項と、複雑さを表すペナルティー項で構成されています。αは求めたいパラメータwを制御する上位のパラメータ（ハイパーパラメータ）。αを大きくするほど大きなL1ノルムに対してペナルティーがかかる仕組みです。するとL1ノルムは小さくならざるを得なくなり、スパースな解が導かれる（が適合性は低くなりやすい）という仕組み。逆にαを小さくするほどL1ノルムを大きく取れるようになり、適合性が向上しやすい（が複雑で、過学習に陥りやすい）。αをどのように調整・採用するかがキモです。

その判断基準は以下の通り。

・情報量基準（AIC、BICなど）
・cross‐validation

情報量基準の利用では客観性に優れるものの感覚的にずれていると感じることがあるでしょう。一方、cv は理屈でなく現実的。ベースを統計にしているか、機械学習にしているかなどで好き嫌いが分かれそうです。

スパースを仮定した場合、非ゼロのwを増やすと組み込まれるXも高次まで多く含まれて適合性は良くなる（が複雑になる）といった表現が理想です。が、L0ノルムを扱うと組み合わせ（計算コスト）が膨大になるとのこと。そこでα＊L1ノルムで代用しているそうです。
ちなみに、L1、L2 は機械学習でも多用されていましたね。回帰だと以下の通りです。

・L0正則化・・・非ゼロノルムの個数
・L1正則化・・・L1ノルム利用、LASSO回帰
・L2正則化・・・L2ノルム利用、Ridge回帰
・L1＋L2正則化・・・Elastic Net

scikit-learn を利用すれば、これらのモデリングは容易。ありがたいですね。

==========================================
ラグランジュの未定乗数法を適用した場合の図や説明も web 上には多くあります。そちらの方が直感的で理解し易いでしょう。

2020年12月3日木曜日

GEE＋ Sentinel-1

SAR データを利用した Change Detection を整理。

以前利用していたDocker は－vオプションを使用しなければ動きました。機能が追加されているようです。
https://phreeqc.blogspot.com/2019/04/sar-change-detection-on-google-earth.html

GEE コミュニティのチュートリアルにもありました。Part１は動きましたが、Part２の後半はダメでした。
https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-2

いずれも、GEE＋ Sentinel-1。解析エンジンもデータも海外製です。データは政府事業の Tellus 提供より揃っています。
https://phreeqc.blogspot.com/2019/02/google-earth-engine.html
国内の複数個所で Change Detection を試行しましたが、動作も問題ありません。きちんと結果は出てきます。
扱っているモノがモノなので、まだ精度面では物足りません。が、手軽に加工できるところはありがたいですね。将来的にはアルゴリズムが向上し、精度面でも耐え得るモノになるでしょう。

Tellus でも ALOS2 のデータが公開され始めました。が、未だ直近1年分くらいは公開されていません。過去分も、２～3年前くらいまで。直近のデータがないのは辛い。
しかも、開発環境を割り当ててもらえないと、自由に加工できない状況は継続中。https://phreeqc.blogspot.com/2020/09/tellusar.html

まだまだ海外製には追いつけない感が漂っています。頑張って欲しいですね。

登録: 投稿 (Atom)

	a	b	c
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN
4	NaN	NaN	NaN
5	NaN	NaN	NaN
6	NaN	NaN	NaN

	a	b	c
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN

	a	b	c
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN
4	NaN	NaN	NaN
5	NaN	NaN	NaN
6	NaN	NaN	NaN

	a	b	c
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN

	a	b	c
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN
4	NaN	NaN	NaN
5	NaN	NaN	NaN
6	NaN	NaN	NaN

	a	b	c
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN
1	NaN	NaN	NaN
2	NaN	NaN	NaN
3	NaN	NaN	NaN