{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# エンコーダを学習する\n",
"\n",
"__Contents:__\n",
"\n",
"- 視覚刺激から特徴量を造る\n",
"- モデルを初期化してアルゴリズムを走らせる\n",
"- 学習機の出来を評価する\n",
"\n",
"___\n",
"\n",
"ここでの目的は、これまでに習得してきた技法や知見を融合させ、正常に機能する学習機を造ることである。その応用先は視覚刺激から脳活動へのエンコーディングである。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"作業の流れは下記の通りになる：\n",
"\n",
" 視覚刺激の読み込み --> ガボールフィルタで特徴量を作って保存 --> スパースな線形回帰\n",
"\n",
"これらのタスクを一つずつこなしていく。やり方を明確に伝えるべく、非常に単純なプロトタイプを作っておく。その後の課題では、このプロトタイプをベースにして、エンコーダとしてちゃんと機能するように改善してもらうことになる。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## 視覚刺激から特徴量を造る"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import math\n",
"import tables\n",
"\n",
"import gaborfil\n",
"import models\n",
"import dataclass\n",
"import algorithms"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"data/vim-2/stimulus_ds.h5 (File) 'vim-2: stimulus'\n",
"Last modif.: 'Tue Mar 27 21:14:47 2018'\n",
"Object Tree: \n",
"/ (RootGroup) 'vim-2: stimulus'\n",
"/test (Array(64, 64, 3, 8100)) 'Testing data'\n",
"/train (Array(64, 64, 3, 108000)) 'Training data'\n",
"\n"
]
}
],
"source": [
"# Establish connection with the file objects.\n",
"h5_X = tables.open_file(\"data/vim-2/stimulus_ds.h5\", mode=\"r\")\n",
"print(h5_X)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Set up the parameters that specify the first filter bank.\n",
"PIX_W = 64\n",
"PIX_H = 64\n",
"max_cycles = 32 # the maximum cycles per image.\n",
"myparas = {\"freqs\": max_cycles/max(PIX_W,PIX_H),\n",
" \"dir\": 0,\n",
" \"amp\": 0.1,\n",
" \"sdev\": max(PIX_W,PIX_H)/20,\n",
" \"phase\": 0}\n",
"mygrid_h = 4\n",
"mygrid_w = 4"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Images processed so far: 0\n",
"Images processed so far: 2160\n",
"Images processed so far: 4320\n",
"Images processed so far: 6480\n",
"Images processed so far: 8640\n",
"Images processed so far: 10800\n",
"Images processed so far: 12960\n",
"Images processed so far: 15120\n",
"Images processed so far: 17280\n",
"Images processed so far: 19440\n",
"Images processed so far: 21600\n",
"Images processed so far: 23760\n",
"Images processed so far: 25920\n",
"Images processed so far: 28080\n",
"Images processed so far: 30240\n",
"Images processed so far: 32400\n",
"Images processed so far: 34560\n",
"Images processed so far: 36720\n",
"Images processed so far: 38880\n",
"Images processed so far: 41040\n",
"Images processed so far: 43200\n",
"Images processed so far: 45360\n",
"Images processed so far: 47520\n",
"Images processed so far: 49680\n",
"Images processed so far: 51840\n",
"Images processed so far: 54000\n",
"Images processed so far: 56160\n",
"Images processed so far: 58320\n",
"Images processed so far: 60480\n",
"Images processed so far: 62640\n",
"Images processed so far: 64800\n",
"Images processed so far: 66960\n",
"Images processed so far: 69120\n",
"Images processed so far: 71280\n",
"Images processed so far: 73440\n",
"Images processed so far: 75600\n",
"Images processed so far: 77760\n",
"Images processed so far: 79920\n",
"Images processed so far: 82080\n",
"Images processed so far: 84240\n",
"Images processed so far: 86400\n",
"Images processed so far: 88560\n",
"Images processed so far: 90720\n",
"Images processed so far: 92880\n",
"Images processed so far: 95040\n",
"Images processed so far: 97200\n",
"Images processed so far: 99360\n",
"Images processed so far: 101520\n",
"Images processed so far: 103680\n",
"Images processed so far: 105840\n",
"(108000, 16)\n"
]
}
],
"source": [
"# Construct features using the specified filter bank (TRAINING).\n",
"X_tr = gaborfil.G2_getfeatures(ims=h5_X.root.train.read(),\n",
" fil_paras=myparas,\n",
" gridshape=(mygrid_h, mygrid_w),\n",
" mode=\"reflect\", cval=0, verbose=True)\n",
"print(X_tr.shape)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Images processed so far: 0\n",
"Images processed so far: 162\n",
"Images processed so far: 324\n",
"Images processed so far: 486\n",
"Images processed so far: 648\n",
"Images processed so far: 810\n",
"Images processed so far: 972\n",
"Images processed so far: 1134\n",
"Images processed so far: 1296\n",
"Images processed so far: 1458\n",
"Images processed so far: 1620\n",
"Images processed so far: 1782\n",
"Images processed so far: 1944\n",
"Images processed so far: 2106\n",
"Images processed so far: 2268\n",
"Images processed so far: 2430\n",
"Images processed so far: 2592\n",
"Images processed so far: 2754\n",
"Images processed so far: 2916\n",
"Images processed so far: 3078\n",
"Images processed so far: 3240\n",
"Images processed so far: 3402\n",
"Images processed so far: 3564\n",
"Images processed so far: 3726\n",
"Images processed so far: 3888\n",
"Images processed so far: 4050\n",
"Images processed so far: 4212\n",
"Images processed so far: 4374\n",
"Images processed so far: 4536\n",
"Images processed so far: 4698\n",
"Images processed so far: 4860\n",
"Images processed so far: 5022\n",
"Images processed so far: 5184\n",
"Images processed so far: 5346\n",
"Images processed so far: 5508\n",
"Images processed so far: 5670\n",
"Images processed so far: 5832\n",
"Images processed so far: 5994\n",
"Images processed so far: 6156\n",
"Images processed so far: 6318\n",
"Images processed so far: 6480\n",
"Images processed so far: 6642\n",
"Images processed so far: 6804\n",
"Images processed so far: 6966\n",
"Images processed so far: 7128\n",
"Images processed so far: 7290\n",
"Images processed so far: 7452\n",
"Images processed so far: 7614\n",
"Images processed so far: 7776\n",
"Images processed so far: 7938\n",
"(8100, 16)\n"
]
}
],
"source": [
"# Construct features using the specified filter bank (TESTING).\n",
"X_te = gaborfil.G2_getfeatures(ims=h5_X.root.test.read(),\n",
" fil_paras=myparas,\n",
" gridshape=(mygrid_h, mygrid_w),\n",
" mode=\"reflect\", cval=0, verbose=True)\n",
"print(X_te.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"これで出来上がった特徴量は、ある一つだけのガボールフィルタの下でできたものである。空間周波数に着目すると、上記ではその値をゼロとしているが、ほかの値もまんべんなく反映すべきであろう。実際、Nishimoto *et al.* (2011)で試される手法のうち、\"static model\"が我々が実装しているものにきわめて近い。ただし、彼らのモデルでは、フィルタの向きが\"0, 45, 90 and 135 degrees\"（appendixより引用）とある。0度の場合は済んでいるので、残りを用意しておこう。"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Adding features using dir = 45.0 degrees\n",
"Adding features using dir = 90.0 degrees\n",
"Adding features using dir = 135.0 degrees\n"
]
}
],
"source": [
"todo_dir = math.pi * np.array([1,2,3]) / 4\n",
"\n",
"for mydir in todo_dir:\n",
" print(\"Adding features using dir =\", mydir*(360/(2*math.pi)), \"degrees\")\n",
" myparas[\"dir\"] = mydir\n",
" \n",
" tmp_X = gaborfil.G2_getfeatures(ims=h5_X.root.train.read(),\n",
" fil_paras=myparas,\n",
" gridshape=(mygrid_h, mygrid_w),\n",
" mode=\"reflect\", cval=0, verbose=False)\n",
" X_tr = np.concatenate((X_tr, tmp_X), axis=1)\n",
" \n",
" tmp_X = gaborfil.G2_getfeatures(ims=h5_X.root.test.read(),\n",
" fil_paras=myparas,\n",
" gridshape=(mygrid_h, mygrid_w),\n",
" mode=\"reflect\", cval=0, verbose=False)\n",
" X_te = np.concatenate((X_te, tmp_X), axis=1)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"# All finished with this file, so close it.\n",
"h5_X.close()\n",
"print(h5_X)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"次に時間軸におけるダウンサンプリングが必要である。応答信号は一秒に一度のペースで観測しているのに対して、視覚刺激の入力は一秒に15枚もある。15フレームおきにこの時系列を分割して、平均を取るという形で望ましいデータ数に縮約する。"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"7200\n",
"540\n"
]
}
],
"source": [
"framerate = 15\n",
"n_tr = X_tr.shape[0]//framerate\n",
"n_te = X_te.shape[0]//framerate\n",
"print(n_tr)\n",
"print(n_te)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"X_tr shape after down-sampling: (7200, 64)\n",
"X_te shape after down-sampling: (540, 64)\n"
]
}
],
"source": [
"# Training data.\n",
"tmp_X = np.zeros((n_tr,X_tr.shape[1]), dtype=X_tr.dtype)\n",
"idx = np.arange(framerate)\n",
"for i in range(n_tr):\n",
" tmp_X[i,:] = np.mean(X_tr[idx,:], axis=0)\n",
" idx += framerate\n",
"X_tr = tmp_X\n",
"print(\"X_tr shape after down-sampling:\", X_tr.shape)\n",
"\n",
"# Testing data.\n",
"tmp_X = np.zeros((n_te,X_te.shape[1]), dtype=X_te.dtype)\n",
"idx = np.arange(framerate)\n",
"for i in range(n_te):\n",
" tmp_X[i,:] = np.mean(X_te[idx,:], axis=0)\n",
" idx += framerate\n",
"X_te = tmp_X\n",
"print(\"X_te shape after down-sampling:\", X_te.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"また、Nishimoto et al. (2011)にしたがって、特徴量のZスコアを用いることにする。記号で表わすと：\n",
"\n",
"\\begin{align}\n",
"z = \\frac{x - \\bar{x} }{\\sqrt{\\widehat{v}}}\n",
"\\end{align}\n",
"\n",
"ここで$\\bar{x}$が経験平均、$\\widehat{v}$が経験分散である。"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# Z-scores\n",
"Z_tr = X_tr - np.mean(X_tr, axis=0)\n",
"Z_tr = Z_tr / np.std(Z_tr, axis=0)\n",
"#print(\"Mean =\", np.mean(Z_tr, axis=0), \"StdDev =\", np.std(Z_tr, axis=0))\n",
"Z_te = X_te - np.mean(X_te, axis=0)\n",
"Z_te = Z_te / np.std(Z_te, axis=0)\n",
"#print(\"Mean =\", np.mean(Z_te, axis=0), \"StdDev =\", np.std(Z_te, axis=0))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"さらに、標準偏差の3倍を超えた値を外れ値と見なして、切断する。学習課題の公平性を保つため、訓練データと検証データを別々に扱う。"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# Truncation of outliers.\n",
"thres = 3\n",
"for j in range(X_tr.shape[1]):\n",
" stdval = np.std(Z_tr[:,j])\n",
" Z_tr[:,j] = np.clip(Z_tr[:,j], a_min=(-thres*stdval), a_max=thres*stdval)\n",
" stdval = np.std(Z_te[:,j])\n",
" Z_te[:,j] = np.clip(Z_te[:,j], a_min=(-thres*stdval), a_max=thres*stdval)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"新しい階層型データファイル`features.h5`を作って、学習と検証用の特徴量を格納しておく。"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"data/vim-2/features.h5 (File) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
"Last modif.: 'Sat Apr 7 17:14:23 2018'\n",
"Object Tree: \n",
"/ (RootGroup) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
"\n"
]
}
],
"source": [
"# Open file connection, writing new file to disk.\n",
"myh5 = tables.open_file(\"data/vim-2/features.h5\",\n",
" mode=\"w\",\n",
" title=\"Features from vim-2 stimulus, via 2D Gabor filter bank\")\n",
"print(myh5)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"data/vim-2/features.h5 (File) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
"Last modif.: 'Sat Apr 7 17:14:34 2018'\n",
"Object Tree: \n",
"/ (RootGroup) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
"/train (Array(7200, 64)) 'Training data'\n",
"\n",
"data/vim-2/features.h5 (File) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
"Last modif.: 'Sat Apr 7 17:14:34 2018'\n",
"Object Tree: \n",
"/ (RootGroup) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
"/test (Array(540, 64)) 'Testing data'\n",
"/train (Array(7200, 64)) 'Training data'\n",
"\n"
]
}
],
"source": [
"# Add arrays.\n",
"myh5.create_array(where=myh5.root, name=\"train\", obj=Z_tr, title=\"Training data\")\n",
"print(myh5)\n",
"myh5.create_array(where=myh5.root, name=\"test\", obj=Z_te, title=\"Testing data\")\n",
"print(myh5)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"# Close the file connection.\n",
"myh5.close()\n",
"print(myh5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"___"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## モデルを初期化してアルゴリズムを走らせる"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"特徴量の準備が整っているなら、Jupyterのカーネルをリセットし、ここから新たな作業を始める。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"前回で作った`Algo_LASSO_CD`をここで本領発揮してもらう。前回のテストとまったく同様に、多数の$\\lambda$候補を用意し、warm startを生かしながら学習させていく。パフォーマンス評価指標として、Nishimoto et al. (2011)より："
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> *\"Prediction accuracy was defined as the correlation between predicted and observed BOLD signals. The averaged accuracy across subjects and voxels in early visual areas (V1, V2, V3, V3A, and V3B) was 0.24, 0.39, and 0.40 for the static, nondirectional, and directional encoding models, respectively.\"*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"教育目的なので我々の学習機があらゆる意味で彼らのものに劣っているのだが、めざす基準としてはこの数値には意義がある。特に相関係数0.24は、2次元のガボールフィルタを使ったときに出た数字なので、より現実的な基準と見てもよかろう。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"算出方法は次の通りである。訓練データを使って、学習アルゴリズムを一通り走らせると、$\\widehat{w}$に対応する`w_est`が定まる。線形モデルを使っているので、新しい入力$X$から$\\widehat{y} = X\\widehat{w}$を出して、$\\widehat{y} \\approx y$で近似してみる。この$X$と$y$が`X_te`と`y_te`に対応する。予測信号と本当の信号の相関が強いほどよいということなので、下記の通りに相関係数を出す："
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\\begin{align}\n",
"\\text{corr}\\,(\\widehat{y},y) = \\frac{\\text{cov}\\,(\\widehat{y},y)}{\\sqrt{\\text{var}\\,(\\widehat{y})\\text{var}\\,(y)}},\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"これは`scipy.stats.pearsonr`で計算できる。また、予測・真の平均的な2乗誤差（root mean squared error; RMSE）を求めることもできる："
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\\begin{align}\n",
"\\text{RMSE}\\,(\\widehat{y},y) = \\left( \\frac{1}{m} \\sum_{i=1}^{m} (\\widehat{y}_{i}-y_{i})^2 \\right)^{1/2},\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"これはモデルオブジェクトの`mod.eval`というメソッドで実装している。$m$は検証データのサンプル数（ここでは$m=540$）。"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import math\n",
"import tables\n",
"\n",
"import gaborfil\n",
"import dataclass"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"data/vim-2/features.h5 (File) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
"Last modif.: 'Sat Apr 7 17:14:42 2018'\n",
"Object Tree: \n",
"/ (RootGroup) 'Features from vim-2 stimulus, via 2D Gabor filter bank'\n",
"/test (Array(540, 64)) 'Testing data'\n",
"/train (Array(7200, 64)) 'Training data'\n",
"\n",
"data/vim-2/response.h5 (File) 'vim-2: BOLD responses'\n",
"Last modif.: 'Mon Apr 9 09:14:08 2018'\n",
"Object Tree: \n",
"/ (RootGroup) 'vim-2: BOLD responses'\n",
"/sub1 (Group) 'Data for subject 1'\n",
"/sub2 (Group) 'Data for subject 2'\n",
"/sub3 (Group) 'Data for subject 3'\n",
"/sub3/idx (Group) 'ROI-specific voxel indices'\n",
"/sub3/idx/v1lh (Array(653,)) ''\n",
"/sub3/idx/v1rh (Array(713,)) ''\n",
"/sub3/idx/v2lh (Array(735,)) ''\n",
"/sub3/idx/v2rh (Array(642,)) ''\n",
"/sub3/idx/v3alh (Array(164,)) ''\n",
"/sub3/idx/v3arh (Array(118,)) ''\n",
"/sub3/idx/v3blh (Array(88,)) ''\n",
"/sub3/idx/v3brh (Array(138,)) ''\n",
"/sub3/idx/v3lh (Array(504,)) ''\n",
"/sub3/idx/v3rh (Array(627,)) ''\n",
"/sub3/resp (Group) 'Response arrays'\n",
"/sub3/resp/test (Array(4381, 540)) 'Testing data'\n",
"/sub3/resp/train (Array(4381, 7200)) 'Training data'\n",
"/sub2/idx (Group) 'ROI-specific voxel indices'\n",
"/sub2/idx/v1lh (Array(470,)) ''\n",
"/sub2/idx/v1rh (Array(573,)) ''\n",
"/sub2/idx/v2lh (Array(733,)) ''\n",
"/sub2/idx/v2rh (Array(926,)) ''\n",
"/sub2/idx/v3alh (Array(135,)) ''\n",
"/sub2/idx/v3arh (Array(202,)) ''\n",
"/sub2/idx/v3blh (Array(83,)) ''\n",
"/sub2/idx/v3brh (Array(140,)) ''\n",
"/sub2/idx/v3lh (Array(714,)) ''\n",
"/sub2/idx/v3rh (Array(646,)) ''\n",
"/sub2/resp (Group) 'Response arrays'\n",
"/sub2/resp/test (Array(4622, 540)) 'Testing data'\n",
"/sub2/resp/train (Array(4622, 7200)) 'Training data'\n",
"/sub1/idx (Group) 'ROI-specific voxel indices'\n",
"/sub1/idx/v1lh (Array(490,)) ''\n",
"/sub1/idx/v1rh (Array(504,)) ''\n",
"/sub1/idx/v2lh (Array(715,)) ''\n",
"/sub1/idx/v2rh (Array(762,)) ''\n",
"/sub1/idx/v3alh (Array(92,)) ''\n",
"/sub1/idx/v3arh (Array(160,)) ''\n",
"/sub1/idx/v3blh (Array(104,)) ''\n",
"/sub1/idx/v3brh (Array(152,)) ''\n",
"/sub1/idx/v3lh (Array(581,)) ''\n",
"/sub1/idx/v3rh (Array(560,)) ''\n",
"/sub1/resp (Group) 'Response arrays'\n",
"/sub1/resp/test (Array(4120, 540)) 'Testing data'\n",
"/sub1/resp/train (Array(4120, 7200)) 'Training data'\n",
"\n"
]
}
],
"source": [
"# Open file connections with data to be used in learning and evaluation.\n",
"h5_X = tables.open_file(\"data/vim-2/features.h5\", mode=\"r\")\n",
"print(h5_X)\n",
"h5_y = tables.open_file(\"data/vim-2/response.h5\", mode=\"r\")\n",
"print(h5_y)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# Subject ID.\n",
"subid = 1\n",
"y_node = h5_y.get_node(h5_y.root, \"sub\"+str(subid))\n",
"\n",
"# Number of voxels.\n",
"num_voxels = y_node.resp.train.nrows\n",
"\n",
"# Set up the model and data objects.\n",
"mod = models.LinearL1()\n",
"data = dataclass.DataSet()\n",
"data.X_tr = h5_X.root.train.read()\n",
"data.X_te = h5_X.root.test.read()\n",
"\n",
"# Basic info.\n",
"n = data.X_tr.shape[0]\n",
"d = data.X_tr.shape[1]\n",
"\n",
"# Dictionaries of performance over all voxels.\n",
"dict_corr_tr = {}\n",
"dict_corr_te = {}\n",
"dict_l0norm = {}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"いよいよ主要の計算だけが待っている。相当な時間はかかるが、全ボクセル分の学習を下記の通りに行なうことができる。"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Voxel: 1 of 4120\n",
"Voxel: 101 of 4120\n",
"Voxel: 201 of 4120\n",
"Voxel: 301 of 4120\n",
"Voxel: 401 of 4120\n",
"Voxel: 501 of 4120\n",
"Voxel: 601 of 4120\n",
"Voxel: 701 of 4120\n",
"Voxel: 801 of 4120\n",
"Voxel: 901 of 4120\n",
"Voxel: 1001 of 4120\n",
"Voxel: 1101 of 4120\n",
"Voxel: 1201 of 4120\n",
"Voxel: 1301 of 4120\n",
"Voxel: 1401 of 4120\n",
"Voxel: 1501 of 4120\n",
"Voxel: 1601 of 4120\n",
"Voxel: 1701 of 4120\n",
"Voxel: 1801 of 4120\n",
"Voxel: 1901 of 4120\n",
"Voxel: 2001 of 4120\n",
"Voxel: 2101 of 4120\n",
"Voxel: 2201 of 4120\n",
"Voxel: 2301 of 4120\n",
"Voxel: 2401 of 4120\n",
"Voxel: 2501 of 4120\n",
"Voxel: 2601 of 4120\n",
"Voxel: 2701 of 4120\n",
"Voxel: 2801 of 4120\n",
"Voxel: 2901 of 4120\n",
"Voxel: 3001 of 4120\n",
"Voxel: 3101 of 4120\n",
"Voxel: 3201 of 4120\n",
"Voxel: 3301 of 4120\n",
"Voxel: 3401 of 4120\n",
"Voxel: 3501 of 4120\n",
"Voxel: 3601 of 4120\n",
"Voxel: 3701 of 4120\n",
"Voxel: 3801 of 4120\n",
"Voxel: 3901 of 4120\n",
"Voxel: 4001 of 4120\n",
"Voxel: 4101 of 4120\n"
]
}
],
"source": [
"for voxidx in range(num_voxels):\n",
" \n",
" # Set up the responses.\n",
" data.y_tr = np.transpose(np.take(a=y_node.resp.train.read(),\n",
" indices=[voxidx],\n",
" axis=0))\n",
" data.y_te = np.transpose(np.take(a=y_node.resp.test.read(),\n",
" indices=[voxidx],\n",
" axis=0))\n",
"\n",
" # Set up for a loop over trials and lambda values.\n",
" todo_lambda = np.logspace(start=math.log10(1/n), stop=math.log10(2.5), num=50)\n",
" num_loops = 15\n",
" t_max = num_loops * d\n",
" \n",
" # Storage for performance metrics.\n",
" corr_tr = np.zeros(todo_lambda.size, dtype=np.float32)\n",
" corr_te = np.zeros(todo_lambda.size, dtype=np.float32)\n",
" l0norm = np.zeros(todo_lambda.size, dtype=np.uint32)\n",
" \n",
" # Initialize and run learning algorithm.\n",
" w_init = 1*np.random.uniform(size=(d,1))\n",
"\n",
" for l in range(todo_lambda.size):\n",
"\n",
" lamval = todo_lambda[l]\n",
" \n",
" if (voxidx % 100 == 0) and (l == 0):\n",
" print(\"Voxel:\", voxidx+1, \"of\", num_voxels)\n",
" #print(\"Lambda value =\", lamval, \"(\", l, \"of\", todo_lambda.size, \")\")\n",
"\n",
" # Use warm starts when available.\n",
" if l > 0:\n",
" w_init = al.w\n",
" \n",
" al = algorithms.Algo_CDL1(w_init=w_init, t_max=t_max, lamreg=lamval)\n",
"\n",
" # Iterate the learning algorithm.\n",
" for onestep in al:\n",
" al.update(model=mod, data=data)\n",
"\n",
" # Record performance based on final output.\n",
" corr_tr[l] = gaborfil.corr(w=al.w, X=data.X_tr, y=data.y_tr)\n",
" corr_te[l] = gaborfil.corr(w=al.w, X=data.X_te, y=data.y_te)\n",
" l0norm[l] = np.nonzero(al.w)[0].size\n",
" \n",
" # Save the performance for this voxel.\n",
" dict_corr_tr[voxidx] = corr_tr\n",
" dict_corr_te[voxidx] = corr_te\n",
" dict_l0norm[voxidx] = l0norm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"実績をディスクに書き込んで保存する。"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"import pickle\n",
"\n",
"# Method name\n",
"mthname = \"DefaultMth\"\n",
"\n",
"# Lambda values used.\n",
"fname = \"results/\"+mthname+\"sub\"+str(subid)+\".lam\"\n",
"with open(fname, mode=\"bw\") as fbin:\n",
" pickle.dump(todo_lambda, fbin)\n",
"\n",
"# Correlation over lambda values.\n",
"fname = \"results/\"+mthname+\"sub\"+str(subid)+\".corrtr\"\n",
"with open(fname, mode=\"bw\") as fbin:\n",
" pickle.dump(dict_corr_tr, fbin)\n",
"fname = \"results/\"+mthname+\"sub\"+str(subid)+\".corrte\"\n",
"with open(fname, mode=\"bw\") as fbin:\n",
" pickle.dump(dict_corr_te, fbin)\n",
"\n",
"# Sparsity over lambda values.\n",
"fname = \"results/\"+mthname+\"sub\"+str(subid)+\".l0norm\"\n",
"with open(fname, mode=\"bw\") as fbin:\n",
" pickle.dump(dict_l0norm, fbin)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"\n",
"## 学習機の出来を評価する\n",
"\n",
"学習が終わって、あとは成績を集計したり、可視化したりするだけである。\n",
"Jupyterのカーネルをリセットし、ここから新たな作業を始める。"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Preparation.\n",
"import math\n",
"import numpy as np\n",
"import pickle\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import tables"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Method name\n",
"mthname = \"DefaultMth\"\n",
"\n",
"# Subject ID.\n",
"subid = 1\n",
"\n",
"# Lambda values used.\n",
"fname = \"results/\"+mthname+\"sub\"+str(subid)+\".lam\"\n",
"with open(fname, mode=\"br\") as fbin:\n",
" todo_lambda = pickle.load(fbin)\n",
"\n",
"# Correlation over lambda values, on the training and test data.\n",
"fname = \"results/\"+mthname+\"sub\"+str(subid)+\".corrtr\"\n",
"with open(fname, mode=\"br\") as fbin:\n",
" dict_corr_tr = pickle.load(fbin)\n",
"fname = \"results/\"+mthname+\"sub\"+str(subid)+\".corrte\"\n",
"with open(fname, mode=\"br\") as fbin:\n",
" dict_corr_te = pickle.load(fbin)\n",
"\n",
"# Sparsity over lambda values.\n",
"fname = \"results/\"+mthname+\"sub\"+str(subid)+\".l0norm\"\n",
"with open(fname, mode=\"br\") as fbin:\n",
" dict_l0norm = pickle.load(fbin)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"まず、パフォーマンス指標をボクセルごとに見てみよう。"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Results: subject 1 , voxel id 3\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAzwAAAG9CAYAAADQhXO6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl4lNXZx/HvnRDIyr6ooIKCW1ERAXdwF7TiUndt1aqgrbZvUau0al1bbKtSFTdcW+v2ulSqaFFfqRsiYKmCqICiIAqRBEhCWBLO+8d5RsYwgUmYmeeZye9zXblmMvPMzJ0oM/k955z7mHMOERERERGRXJQXdgEiIiIiIiLposAjIiIiIiI5S4FHRERERERylgKPiIiIiIjkLAUeERERERHJWQo8IiIiIiKSsxR4JOeY2QIzO7yZjz3IzD5JdU3pYmYHmNlcM6s2s+PNrJuZvWFmVWZ2i5n9xszuT+J57jGzqzNRs4iINM7MZpvZwRl4nWa/72/J56xIGBR4JOXM7Awzmx78Ef61mb1kZgeGXVciZubMrHfse+fcm865ncOsqYmuB+50zpU65/4BjAC+Bdo65y51zv3eOXf+5p7EOXehc+6GLS3GzA42s0Vb+jwiIplmZgea2TtmtsLMKszsbTMbmOk6nHM/cM5NDmq61sweTdPrpOR9H8DMfmVm3wS/uwfNrE0qnlckVRR4JKXMbBQwFvg90A3YDrgLOK4Zz9UqmdtauO2B2Q2+/8hpR2ERkaSZWVvgBeAOoCPQHbgOWJOG18qpzzEzOwq4EjgM6AnsgP/dpft1zcz0d6wkRf+jSMqYWTv8iMPPnXPPOudqnHPrnHP/dM5dHhzTxszGmtni4Gts7ExQbHTAzK4ws2+AhxLdFhz7QzObaWbLgzNyezRS0yAzmxIc97WZ3WlmrYP73ggO+28wGnVqwxEKM9vVzCYHj59tZsPj7nvYzMaZ2YvBFLKpZrbjJn4/sbOHy81soZmdE/u9mdlfzazczL4ws6vi38TN7KdmNsfMKs3sX2a2fXD7fPwHyz+D+h8HzgZ+HXx/eMOzg5uo4WEzuzHuuEZ/v8FUhsvM7IPgbN6TZlZoZiXAS8A2wetXm9k2jf0+REQiZCcA59zjzrl651ytc26Sc+4DADM7JxjxuSN43/vYzA6LPdjMzg3ep6vM7DMzGxl3X6LPts5m9kLwHlthZm/G3veD99jDzWwo8Bvg1OD99L9mdrKZzYgv3MwuNbN/NPyBzOw0M5ve4LZfmdmE4Pp37/tm1iGopzz4rHnBzHok+bs7G3jAOTfbOVcJ3ACc09jBm/vsNLP9zWxa8HueZmb7x9032cxuMrO3gVXADsFtNwafVdVm9k8z62RmfzezlcFz9EzyZ5EcpcAjqbQfUAg8t4ljfgvsC/QD9gQGAVfF3b8V/uza9vjpWRvdZmb9gQeBkUAn4F5ggiUeQq8HfgV0Duo7DPgZgHNucHDMnsGUsCfjH2hmBcA/gUlAV+AS4O9mFj/l7XT8mawOwDzgpkQ/tJlthw8DdwBdgp9/ZnD3HUA7fHgZAvwEODd43PH4D7wTg8e9CTwe1L8j8CVwbFD/6cDfgT8G37/ahBrij0vm93sKMBToBewBnOOcqwGGAYuD1y91zi1O9PsQEYmYT4F6M3vEzIaZWYcEx+wDfIb/PPkd8KyZdQzuWwr8EGiLf/++LXgvjWn42XYpsAj/XtwN/z7/vZF559zL+NkSTwbvp3sCE4BeZrZr3KFnAX9LUO8EYGcz6xN32xnAYwmOzcOfUNwePzOjFrgzwXGJ/AD4b9z3/wW6mVmnTTwm4Wdn8Pt8Ebgd//lzK/Big+f6Mf53WAZ8Edx2WnB7d2BHYErw83QE5uD/e0kLpsAjqdQJ+NY5V7eJY84ErnfOLXXOlePf8H4cd/964HfOuTXOudpGbrsAuNc5NzU4E/cIftrBvg1fzDk3wzn3rnOuzjm3AP/H+5Akf559gVJgjHNurXPu//BTHk6PO+ZZ59x7wc/8d3yIaOznfjU4e7jOObfMOTfTzPKBU4HRzrmqoMZb4n4nI4E/OOfmBK/xe6BfbJSniRLWkOC4ZH6/tzvnFjvnKvChsLGfW0Qk8pxzK4ED8aFjPFBuZhPMrFvcYUuBscH755PAJ8AxweNfdM7Nd96/8SfKDop7bMPPsXXA1sD2wfO9mcxUZOfcGuBJfMjBzH6An0b2QoJjVwHPE3xmBcFnF3wQanjsMufcM865Vc65KnwASfazshRYEfd97HrZJh7T2GfnMcBc59zfgs/tx4GPgWPjHvtwMJpU55xbF9z2UPD7X4E/sTffOfdq8Pz/C+yV5M8iOUqBR1JpGdDZNj0/eRs2nJEhuB4/7ancObe6wWMa3rY9cGkwFWC5mS0Htm3wPACY2U7B0Pw3ZrYSHxg6J/nzbAMsdM6tb1Bv97jvv4m7vgr/xp/ItsD8BLd3Blqz8e8k9hrbA3+J+zkrAGtQQ7Iaq6GhZH6/yf7cIiJZITixdI5zrgfQF/+eNzbukK8ahJLvPr+CUaF3g+lpy4Gj+f5nTcPPsT/hRzYmBVPgrmxCqY8AZ5iZ4U+OPRUEoUQeY8NJujOAfwRB6HvMrNjM7jU/rXol8AbQPjgptznV+JGtmNj1KvOdQmNTnO+JO6axz5CGfyPAxp+7CxPUsCTuem2C7/UZ1cIp8EgqTQFWA8dv4pjF+D+oY7YLbotJdIar4W0LgZucc+3jvoqDM0EN3Y0/O9THOdcWP23ANvNzxNe6rX1/UeR2wFdJPr5hzYnW93yLP9PX8HcSe42FwMgGP2uRc+6dFNaQ6Lhkf78NqVmCiGQ959zHwMP44BPTPQgZMdsBi4Ppvs8Afwa6OefaAxP5/mdNw+lqVc530twBP3oxKn5NUGOPCx77LrAWP4J0Bomns8VMwp+I7IcPPomms4GfYrczsE/wWRmb8p3M5+Vs/BT1mD2BJcGo0e/jpjhfmMRzNfwbATb+3NXnjDSZAo+kTDCUfA0wzvyeMMVmVhCc+fpjcNjjwFVm1sXMOgfHN7Xl5njgQjPbx7wSMzvGzBINn5cBK4FqM9sFuKjB/Uvwa2cSmQrU4JsAFJjfF+FY4Ikm1gt+yP5wMzvFzFoFCyr7OefqgaeAm8ysLJiqNooNv5N7gNHBtIVYg4OTm/H6jdaQ4Lim/H4bWgJ0Mt/AQkQkK5jZLuYX//cIvt8WHxDejTusK/CL4PPgZGBXfLBpDbQByoE6MxsGHLmZ1/uhmfUOAtRK/HrT+gSHLgF62sbdyP6KX2NT55x7q7HXCaZ0PY0fUeoIvNLIoWX4kZDlwTqapqx5+StwnpntFqx9ugofFptjIrCT+e0tWpnZqcBuJJiyJ9IUCjySUs65W/F/sF+Ff/NfCFwMxDrI3AhMBz4APgTeD25rymtMx68zuROoxE8LOKeRwy/DnwGrwv8h/2SD+68FHgmmbp3S4HXWAsPxC/G/xbfX/klw5q9JnHNf4qc4XIqfljaTDWfELsEHq8+At/Bn4B4MHvcccDPwRDDNYFZQT5Ntpob445ry+2342I/xofaz4HeqLm0ikg2q8E0JpppZDT7ozMK/X8ZMBfrgPw9uAk4KRjGqgF/gT15V4j9zNlon00Af4FX8dLApwF0u2Hungf8NLpeZ2ftxt/8NP/q0qdGdmMeAw4H/dY2vsR0LFOF/tneBl5N4XuC75gp/BF7HTz/7gmY2CXDOLcM3f7gUP03+18APnXPfNuf5RGIsiTVyIiIiIi2W+Rb+5zvnIrGJtpkV4Zso9HfOzQ27HpGo0wiPiIiISHa5CJimsCOSnJza7VdEREQkl5nZAnwzgU01CBKROJrSJiIiIiIiOUtT2kREREREJGclNaXNzIYCfwHygfudc2Ma3D8Y3+FjD+A059zTcfdtB9yP37jQAUcHu8kn1LlzZ9ezZ8+m/RQiIpJyM2bM+NY51yXsOqJIn1UiIuFL9nNqs4En2GV3HHAEsAiYZmYTnHMfxR32Jb5t7WUJnuKv+E0MXzGzUmB9gmO+07NnT6ZPn765skREJM3MrOGO5xLQZ5WISPiS/ZxKZoRnEDDPOfdZ8MRPAMcB3wWe2IiNmX0vzJjZbkAr59wrwXHVyRQlIiIiIiKSCsms4emO3zwyZlFwWzJ2wu/a+6yZ/cfM/hSMGImIiIiIiKRdMoHHEtyWbGu3VsBB+KluA4EdSLBju5mNMLPpZja9vLw8yacWERERERHZtGSmtC3CNxyI6QEsTvL5FwH/iZsO9w9gX+CB+IOcc/cB9wEMGDBAfbJFJHTr1q1j0aJFrF69OuxS0q6wsJAePXpQUFAQdikiIiIpl0zgmQb0MbNewFfAacAZST7/NKCDmXVxzpUDhwJa5Skikbdo0SLKysro2bMnZokGunODc45ly5axaNEievXqFXY5IiIiKbfZKW3OuTrgYuBfwBzgKefcbDO73syGA5jZQDNbBJwM3Gtms4PH1uOns71mZh/ip8eNT8+PIiKSOqtXr6ZTp045HXYAzIxOnTq1iJEsERFpmZLah8c5NxGY2OC2a+KuT8NPdUv02Ffw+/OIiGSVXA87MS3l5xQRkZYpmaYFIiIiIiIiWUmBR0QkgpYvX85dd93V5McdffTRLF++PA0ViYiIZCcFHhGRCGos8NTX12/ycRMnTqR9+/bpKktERCTrJLWGR0REMuvKK69k/vz59OvXj4KCAkpLS9l6662ZOXMmH330EccffzwLFy5k9erV/PKXv2TEiBEA9OzZk+nTp1NdXc2wYcM48MADeeedd+jevTvPP/88RUVFIf9kIiIimaXAIyKyOf/zPzBzZmqfs18/GDu20bvHjBnDrFmzmDlzJpMnT+aYY45h1qxZ37WOfvDBB+nYsSO1tbUMHDiQH/3oR3Tq1Ol7zzF37lwef/xxxo8fzymnnMIzzzzDWWedldqfQ0REJOIUeEREssCgQYO+t0/O7bffznPPPQfAwoULmTt37kaBp1evXvTr1w+AvffemwULFmSsXhERkahQ4BER2ZxNjMRkSklJyXfXJ0+ezKuvvsqUKVMoLi7m4IMPTriPTps2bb67np+fT21tbUZqFRERiRI1LRARiaCysjKqqqoS3rdixQo6dOhAcXExH3/8Me+++26GqxMREckeGuERSacvv4Yly6BvbygqDLsaySKdOnXigAMOoG/fvhQVFdGtW7fv7hs6dCj33HMPe+yxBzvvvDP77rtviJVKU1RVQefOzXvsCSfA/fdDaWlqaxIRyXUKPCLpsvAb+Pwrf/3DudBvF2hdEG5NklUee+yxhLe3adOGl156KeF9sXU6nTt3ZtasWd/dftlll6W8Pmm61q1h1KimP66yEsaPhzlz4PnnoWfPlJcmIpKzFHhE0uGrJfDZIujSAbbpCh9+CrPmwZ47QX5+2NWJSEjatIE//KF5jz3hBDj1VBg4EJ55BgYPTm1tIiK5Smt4RFLt63KYtxA6tYddekH7Mth1B6iqgTmfg3NhVygiWeioo+C996BTJzjsMLjvvrArEhHJDgo8Iqm0ZBl8+gV0aAu77QB5wT+xzh2g93awbDnM+1KhR0SaZaedYOpUOOIIGDkSLr4Y1q0LuyoRkWhT4BFJlfIK+PhzP6Lzg94bwk5M967QoxssLvfre0REmqFdO/jnP+Hyy2HcOD/yU1MTdlUiItGlwCOSCt8u99PV2pb6jmz5jfzT2qEHdO3omxksWZbZGkUkZ+Tnwx//CI88Aq+/DtdfH3ZFIiLRpcAjsqXW1cGcz6C0CHbvvemmBGawc09oVwafLIDKlZmqUkRy0E9+Aj/9Kdx6K8yeHXY1IiLRpMAjsqUqV8D69bDjdtAqicaHeXnQd0coagMffQar16S/Rsk6y5cv56677mrWY8eOHcuqVatSXJFE1c03Q9u28LOfaXmgiEgiCjwiW2rZCihoBW1Lkn9Mq1Z+nY9b70PP+vXpq0+ykgKPJKtzZxgzBt54Ax59NOxqRESiR/vwiGwJ56BiJXRs66erNUVxIezcCz6aD/MXQp/t01OjZKUrr7yS+fPn069fP4444gi6du3KU089xZo1azjhhBO47rrrqKmp4ZRTTmHRokXU19dz9dVXs2TJEhYvXswhhxxC586def3118P+USQDzjsPHngALrsMfvhD6NBh08evWOGbH4iItAQKPCJbYmUN1NX5PXeao0sH37lt0RLf8KBbp9TWJ6kx70uoTvGISWmxb1XeiDFjxjBr1ixmzpzJpEmTePrpp3nvvfdwzjF8+HDeeOMNysvL2WabbXjxxRcBWLFiBe3atePWW2/l9ddfp3PnzqmtWSIrLw/uvhsGDICrrvLd2xJxzq/3+fWvYfx4v/5HRCTXaUqbyJaoWO4vO7Zt/nP06u7DzqdfQE1tauqSnDJp0iQmTZrEXnvtRf/+/fn444+ZO3cuu+++O6+++ipXXHEFb775Ju10yr5F22svvy/P3XfD9Okb3796NZx9th8FatUKxo7Vmh8RaRk0wiOyJZatgHalyTUraExent+kdMZHMHs+9N8VWm2i05tk3iZGYjLBOcfo0aMZOXLkRvfNmDGDiRMnMnr0aI488kiuueaaECqUqLj+enjqKbjoInj33Q1NIxcvhhNOgPfe88d06+Y3Lp0yBfbfP9yaRUTSTSM8Is21eq0fkemYgrPqbVr70FO7Gj5doNOuQllZGVVVVQAcddRRPPjgg1RXVwPw1VdfsXTpUhYvXkxxcTFnnXUWl112Ge+///5Gj5WWpV07P2Vt+nS47z5/27RpMHCgb1v97LNw9dVwxhlQVgb33BNuvSIimaARHpHmqljhL5u7fqeh9m399LbPv4K2S/3aHmmxOnXqxAEHHEDfvn0ZNmwYZ5xxBvvttx8ApaWlPProo8ybN4/LL7+cvLw8CgoKuPvuuwEYMWIEw4YNY+utt1bTghbotNPg/vth9Gior/dT2LbaCt55B/bYwx9TWgo//rFvdHDbbdBJywdFJIeZi9iZ5AEDBrjpiSYfi0TNrLlQXQv77N70Dm2Ncc5Pa6tYAQN/AEWFqXleabI5c+aw6667hl1GxiT6ec1shnNuQEglRVrUP6s+/tiHm3XrYPBgePpp6NLl+8d88AHsuSfccguMGhVOnSIiWyLZzylNaRNpjvXrobIKOrVLXdgB/1w79PDBZ3l16p5XRFqUXXbxzQt++1t45ZWNww74QLT//n5aW8TOfYqIpJQCj0hzLK/yoScV63caKmrjVxpX1aT+uUWkxTjvPLjxRmjduvFjLrwQ5s4FzXwUkVymwCPSHMtW+O5q7begHXVjzKCsWIEnAqI25TddWsrPCWBm7c3saTP72MzmmNl+ZtbRzF4xs7nB5Wa27cwdJ50EHTuqeYGI5DYFHpGmcs7vv9O+DPLT9E+orMR3gFu/Pj3PL5tVWFjIsmXLcj4MOOdYtmwZhYUtZr3YX4CXnXO7AHsCc4Argdecc32A14LvW4SiIjjnHHjuOViyJOxqRETSQ13aRJpq1WrfknrbrdL3GmUlPlhVr/KbkkrG9ejRg0WLFlFeXh52KWlXWFhIjx49wi4j7cysLTAYOAfAObcWWGtmxwEHB4c9AkwGrsh8heEYMcK3sn7wQd/ZTUQk1yjwiDRVrB11OtbvxJQV+8sqBZ6wFBQU0KtXr7DLkNTaASgHHjKzPYEZwC+Bbs65rwGcc1+bWddEDzazEcAIgO22C3cz2lTaeWc49FC/b8+vf71hs1IRkVyhKW0iTbVsBZQUQWGb9L1Gm9ZQ0ErreERSqxXQH7jbObcXUEMTpq855+5zzg1wzg3okqjtWRa78EJYsAAmTQq7EhGR1FPgEWmKujpYWZ3e0R0IGheUKPCIpNYiYJFzbmrw/dP4ALTEzLYGCC6XhlRfaI47Drp1U/MCEclNCjwiTVGx0q+t6ZTmwAM+8KxaDXX16X8tkRbAOfcNsNDMdg5uOgz4CJgAnB3cdjbwfAjlhap1a9/G+oUXYOHCsKsREUktBR6RpqhYAa3yM7OupqzEX1ZrlEckhS4B/m5mHwD9gN8DY4AjzGwucETwfYtzwQX+fM7994ddiYhIaqlpgUiynPOBp2M7P+Us3drGNS5Ix34/Ii2Qc24mMCDBXYdlupao6dkThg3b0LygpCTsikREUiOpER4zG2pmn5jZPDPbaIGnmQ02s/fNrM7MTkpwf1sz+8rM7kxF0SKhqKqBdXXpX78TU1AAha21jkdEMuY3v4FvvoEbbwy7EhGR1Nls4DGzfGAcMAzYDTjdzHZrcNiX+H0NHmvkaW4A/t38MkUiYFkG2lE3pMYFIpJBBxwA554Lf/4zfPRR2NWIiKRGMiM8g4B5zrnPgk3angCOiz/AObfAOfcBsNG28Ga2N9ANULNLyW6VK6FtiW8XnSllJX6T07XrMveaItKi3XwzlJXBz3/uZ/KKiGS7ZAJPdyC+Z8ui4LbNMrM84Bbg8qaXJhIh6+r8SEuHDK+lid+AVEQkA7p0gTFjYPJkeKyxeRsiIlkkmcCTaHV2sud8fgZMdM5tssmlmY0ws+lmNr28vDzJpxbJoOUr/WWHDE5nAygNVg1rWpuIZND558OgQXDppbB8edjViIhsmWQCzyJg27jvewCLk3z+/YCLzWwB8GfgJ2a2UbvPXN69WnJE5UrIz/dT2jKpVT4UFyrwiEhG5eXB3XdDeTlcdVXY1YiIbJlkAs80oI+Z9TKz1sBp+E3aNss5d6ZzbjvnXE/gMuCvzrmNuryJRJpzfsPR9mWZaUfdUKxxgSbTi0gG9e/v1/HcdRfMmBF2NSIizbfZwOOcqwMuBv4FzAGecs7NNrPrzWw4gJkNNLNFwMnAvWY2O51Fi2RU7RpYsxY6hrQXTlmJX0O0Zm04ry8iLdYNN0DXrnDRRVBfH3Y1IiLNk1S7KefcRGBig9uuibs+DT/VbVPP8TDwcJMrFAlbZWz9TliBJ65xQWGbcGoQkRapXTu49VY480wYPx4uvDDsikREmi6pjUdFWrTKlX4D0LDCRmmxn0qndTwiEoLTT4dDDoHRo2Hp0rCrERFpOgUekU1Zv953aOvQNpz1O+BXD5cUKfCISCjMYNw4qKmB446Db74JuyIRkaZR4BHZlKoaqF+f+XbUDZWV+CltalwgIiHYdVd4/HH44AMYOBDefz/sikREkqfAI7IpsfU77cvCraOsxK8Yrl0Tbh0i0mL96Efw9tt+0PnAA+HJJ8OuSEQkOQo8IptSudKHjYKk+nukz3eNCzStTUTC068fTJvmW1afdprfo2f9+rCrEhHZNAUekcbU1cHKmvC6s8UrKfKnVRV4RCRkXbvCa6/BeefBTTfBiSdCVVXYVYmINE6BR6Qxy4NP8CgEHjPfrU2BR0QioE0b36b69tvhhRdgv/3gs8/CrkpEJDEFHpHGVKyE/DxoWxJ2JV5ZCVSv0vwREYkEM7jkEnj5ZVi82DczeP31sKsSEdmYAo9IYypX+mYFeRH5Z9K2GNY7qFkddiUiIt85/HCYOhW6dYMjjoC77gq7IhGR74vIX3IiEVO7GlavCb8ddbyyYKRJ09pEJGL69IEpU2DoUPj5z+Gii2Dt2rCrEhHxFHhEEom1o47C+p2YwjbQKl+BR0QiqV07eP55uOIKuOceOPJIKC8PuyoREQi5165IRFWuhDatoahN2JVsEGtcUFMbdiUiIgnl58OYMbD77r6L28CBcOONUFi48bGFhX4KXJsIvc2KSG5S4BFpyDmorIIuHXzIiJJW+bCuLuwqREQ26cwzYaed4Pjj4cc/bvy4/feHZ5/1639ERNJFgUekoaoaqK+P1nS2GDN1aRORrDBwIHzyCSxYkPj+99+HCy+EAQP8VLj+/TNanoi0IAo8Ig1VRHD9Tkxenh+BEhHJAqWl0Ldv4vv69vVT3447Dg48EB5+GE45JaPliUgLoaYFIg1VroSyYiiI4PkAM9+aWkQkB+y1F0yb5kd3Tj0Vrr5ag9giknoKPCLx6uphZXU0R3fABx6N8IhIDunWDV57zTc5uPFGOPFEqKoKuyoRySUKPCLxVgSfsu0jGnjyFHhEJPe0aQPjx8Ptt8MLL/hmBp99FnZVIpIrFHhE4i2v8qMo7UrDriQxTWkTkRxlBpdcAi+/DF99BYMGweTJYVclIrlAgUck3ooqaFvimwNEkaa0iUiOO/xwmDoVunb1+/TcfXfYFYlItovoX3UiIairh6pV0L4s7EoaF5vSptAjIjmsTx+YMgWOOgp+9jO46CJYty7sqkQkWynwiMTE1u+0i3DgiW2EqsAjIjmuXTu/P88VV8A99/jRnm+/DbsqEclGCjwiMSuqfaBoWxJ2JY2LTbVT4BGRFiA/H8aMgUcfhXff9ZuZzp0bdlUikm0UeERilldBWYn/hI2q2AiPGheISAty5pnw5puwbBn8/vdhVyMi2UaBRwSC9Ts10V6/A5rSJiIt1sCBcOyx8OKLUF8fdjUikk0UeETAbzYK0Q88eRrhEZGWa/hwKC/3XdxERJLVKuwCslZtLSxYAKtW+a/a2g3X16zxu6gVF0NRkb+MfbVr57eVjvK0qZYotv9OlNfvgEZ4RKRFGzoUWrWCCRP85qQiIslQ4GmO+fPhsMPgiy+a9/j8fNhqK+je3X9ts42/3G472H57/7XNNgpFmbSiCsqKo/87j43wuPXh1iEiEoJ27eDgg33gGTMm7GpEJFso8DTV3LlwyCGwejU89BB07Pj9EZyiIj+6s3Zt4tGfigq/hXTs65NP4P/+D1as+P7r5OdDjx4bAlDPnv6rVy9/2aMHFBSE8AvIQfXB/js9uoVdyeZZMAtVU9pEpIUaPhx+8Qv/cdynT9jViEg2UOBpik8+8WFn3TofUvbYI3XPXV0NCxf6UaMvv/SXseuTJ/twtD7urH5eng89u+0Ggwb51ZwDB/rpctI0K6r9FLGor98BTWkTkRbv2GN94PnnP2HUqLCrEZFsoMCTrI8+gkMP9X9ovv469O2b2ucvLYVdd/Vfiaxb5wPRggUbvj7/HP77X7jxxg1haNttffAZNMhPcB4wwI86SeNiG462LQ23jmSoaYGItHA9e/rzjRMmKPCISHIUeJIxa5Zfs5OX58OZQAmNAAAgAElEQVROY6EknQoKYIcd/FdDNTXw/vswbdqGr2ef3fC4vfeGAw7Y8NW1a2Zrj7rl1X7/nVYRX78DGuEREcFPa/vDH/y+PJ06hV2NiESd2lJvzgcf+GlsrVr5qWVhhJ3NKSmBgw7yp7oefxzmzfN9OydMgF/9yq8HuuMOOPFEP+Vtjz3g5pth0aKwKw9ffZbsvxOjwCMiwvDh/u37pZfCrkREsoECz6bMn+/DTmEh/PvfsPPOYVeUvM6d/UTnm2+Gt96ClSvh7bf996WlcOWVvivcoYf65gsrV4ZdcThW1vjw0C5LAs93U9rUpU1EWq6994att/bn9URENkeBZ1NefNF3VXv5ZejdO+xqtkybNn5Nz69/De+849vb/O53fl3QT3/qR35OO82vAl27NuxqM2d5sH6nXRas34ENXdo0wiMiLVhenj+n9/LLfus7EZFNUeDZlIoKf7nLLuHWkQ69e/vA8+mn8O67cP758Nprfp7AVlvByJF+VCvXRxJi++9kw/odiNuHR4FHRFq24cOhqsp/VImIbEpSgcfMhprZJ2Y2z8yuTHD/YDN738zqzOykuNv7mdkUM5ttZh+Y2ampLD7tKiv9LmdR34xyS5jBPvv4NT6LF/tRraOPhr//3e/utt12cPnl8OGHYVeaevXr/ZS2bJnOBhvW8KhLm0izmNkCM/vQzGaa2fTgto5m9oqZzQ0uO4Rdp2zeoYf67e80rU1ENmezgcfM8oFxwDBgN+B0M9utwWFfAucAjzW4fRXwE+fcD4ChwFgza7+lRWdMRYXfWLSlKCjwYefRR2HJEt8AYa+9YOxY3+jgpz/1t+eKlVm0/06MRnhEUuEQ51w/59yA4Psrgdecc32A14LvJeKKiuDII33g0VuiiGxKMiM8g4B5zrnPnHNrgSeA4+IPcM4tcM59AKxvcPunzrm5wfXFwFKgS0oqz4SWFnjilZRsWNPz9dd+lOfRR2GnneDWW/2+QNluRZat3wGN8Iikx3HAI8H1R4DjQ6xFmmD4cL8U9b//DbsSEYmyZAJPd2Bh3PeLgtuaxMwGAa2B+U19bGhacuCJ17kz/PGPfj+iAw6ASy/1Iz6TJoVd2ZZZXgWlxb7leLZQW2qRLeWASWY2w8xGBLd1c859DRBcJtyszMxGmNl0M5teXl6eoXJlU445xr8talqbiGxKMoHHEtzWpL+2zGxr4G/Auc65jVbBR/ZDRIHn+3baCSZOhBdegLo6OOooOP54WLAg7Mqabn0Wrt8BtaUW2XIHOOf646dp/9zMBif7QOfcfc65Ac65AV26ZM9khVzWtSvst58Cj4hsWjKBZxGwbdz3PYDFyb6AmbUFXgSucs69m+iYyH6IKPAkdswxfrRnzBjf2W3PPeHJJ8Ouqmli++9k0/od0AiPyBYKplfjnFsKPIeftr0kODEXO0G3NLwKpamGD4cZM7SXtog0LpnAMw3oY2a9zKw1cBqQ1LmU4PjngL865/63+WWGYP16BZ5NadMGrrjCd2/bbTe/3ueCC6CmJuzKkpNt++/E5GkfHpHmMrMSMyuLXQeOBGbhP9PODg47G3g+nAqlOYYP95cvvBBuHSISXZsNPM65OuBi4F/AHOAp59xsM7vezIYDmNlAM1sEnAzca2azg4efAgwGzglagM40s35p+UlSrarKhx4Fnk3r2RPeeANGj4YHHoABA+CDD8KuavNWVEFpERRk0fqdeGpaINIc3YC3zOy/wHvAi865l4ExwBFmNhc4IvhessQuu/it5TStTUQak9Rfe865icDEBrddE3d9Gn6qW8PHPQo8uoU1hiO26agCz+YVFMDvfw+HHQZnnQWDBvlObhddtGEKVpQ4B1WroGsW/rc1818a4RFpMufcZ8CeCW5fBhyW+YokFcz8KM+dd0J5OURpZryIRENSG4+2SJWV/rKD9p9L2mGH+d6ghx4KP/85/OhHfqQsalavgfp6KCsOu5LmyTON8IiIxDn/fD8p40rtoCQiCSjwNEYjPM3TtaufSH3LLX5+wSGHwNKIrf+tXuUvS7M08GiER0Tke3bd1e+Y8OCD8PbbYVcjIlGjwNMYBZ7my8uDUaPg+efho4/83j2ffx52VRtUBYGnpCjcOprLDDbu7i4i0qJdfTVsu62fTV1XF3Y1IhIlCjyNUeDZcscc49tWV1TA/vvDzJlhV+RVr/JhJy9L//fPy9OUNhGRBkpK4PbbffPQO+4IuxoRiZIs/YsvA2KBR2t4tsx++8Fbb/nGBkOGwOTJYVfkA0+2TmcDTWkTEWnEccf5c23XXANffRV2NSISFQo8jamogOJiKCwMu5Lst+uuflJ1jx5w1FHwzDPh1bJmLayry+7Ak6fAIyKSiJkf5amr8zOrRURAgadx2nQ0tbbdFt580+/Tc/LJcP/94dSR7Q0LwH+ia0qbiEhCO+wAv/0tPPUUTJoUdjUiEgUKPI1R4Em9jh3hlVf8KM/IkfDqq5mvIRcCj0Z4REQ26fLLoU8fv0PC6tVhVyMiYVPgaYwCT3oUF/vTbrvuCqecAvPnZ/b1q2uhsA20ys/s66aSmd9wQkREEmrTBsaNg3nz4E9/CrsaEQmbAk9jKivVsCBdysp8y2rwK0wzuTlptjcsADUtEBFJwhFHwKmnwk03Zf7cmohEiwJPYzTCk1477uhHeubMgbPPzsyIRV0drF6jwCMi0kLceiu0bg2jR4ddiYiESYGnMQo86Xf44XDLLfDcc3DDDel/vVxYvwPah0dEJEnbbAMXXgjPPguLF4ddjYiERYEnkdpav8pRgSf9fvlLP8Jz7bXwj3+k97VigacsywOPRnhERJI2YgTU18ODD4ZdiYiERYEnkdimowo86WcG99wDAwfCj38Ms2al77Wqa6F1gf/KZurSJiKStN69/Xqe++7zwUdEWh4FnkQUeDKrsNBPayst9U0MKivT8zq50LAAtA+PiEgTXXghLFwIL70UdiUiEgYFnkQUeDKve3c/yfqLL+Cqq1L//PXroaY2NwJPnoFTW2oRkWQdeyxsvbWfUCAiLY8CTyIKPOHYbz9/Gu6ee1I/ta2m1l/mQuDRCI+ISJMUFMD558PEif68moi0LAo8iSjwhOfaa6FtWxg1KrXrVHKlQxuA5WkNj4hIE51/vj9fNH582JWISKYp8CQSW0OijUczr3Nn+N3v4JVX/Km4VKleBfn5UNg6dc8ZljyN8IiINNV228Exx8D998O6dWFXIyKZpMCTSEUFtGrlF9FL5v3sZ7DTTnDppan7VKpeBaVF/vRetlNbahGRZrnwQliyBJ5/PuxKRCSTFHgSiW06mgt/HGej1q39hqSffAJ3373lz+cc1ORIhzbYEHgUekREmuSoo2D77dW8QKSlUeBJJBZ4JDzHHOM3Trj2Wli2bMuea9VqPwWsrCQlpYUuLwjiCjwiIk2Sn+83In3tNfj007CrEZFMUeBJRIEnfGZw662wYgVcd92WPVcuNSyADSOPCjwiIk3205/6Wev33Rd2JSKSKQo8iSjwREPfvjByJNx1F8yZ0/znqV7lR0WKC1NXW5g0wiMi0mxbbQUnnAAPPQSrV4ddjYhkggJPIgo80XHddb55xKWXNv85qldBSY40LADflhrUqU1EpJkuvNB/1D/9dNiViEgmKPAkosATHV26wDXXwEsv+a+mci7o0JYj09lAU9pERLbQIYf4ZqBqXiDSMijwNLRuHVRVKfBEycUXQ+/efpSnrq5pj129FurqcyvwxKa0aYRHRKRZzPyM6bffhlmzwq5GRNJNgaeh5cv9pTYdjY7WrWHMGL+O57HHmvbYXGtYABrhERFJgbPO8m+nzz4bdiUikm4KPA1VVPhLjfBEy4knQv/+vk312rXJPy4WeEpyKPCoaYGIyBbr2hX23RcmTAi7EhFJNwWehhR4oskMbrwRPv8cHnww+cdVr/Ld2fJz6H/12AjP+vXh1iEikuWGD4cZM2DRorArEZF0yqG/AlNEgSe6hg6F/feHG26A2trkHpNrDQtgQ5c2jfCIiGyR4cP95QsvhFuHiKSXAk9DCjzRZQY33QSLFyfXWmftOv+Va4FHTQtERFJi111hxx01rU0k1ynwNKTAE20HHwyHHw5/+ANUV2/62FxsWABqWiAikiJmfpTntdc2/5EiItlLgaehigr/DtiuXdiVSGNuvBHKy+Evf9n0cVU5GnjUtEBEJGWGD/e9cCZNCrsSEUkXBZ6GKiqgfXvIzw+7EmnMPvvAscfCn/4ElZWNH1dVA0VtoKBV5mrLBNOUNhGRVDngAL8Thaa1ieQuBZ6GKis1nS0b3HADrFgBt9zS+DHVNVBWkrmaMuW7KW3q0iYisqUKCuDoo33jgvr6sKsRkXRIKvCY2VAz+8TM5pnZlQnuH2xm75tZnZmd1OC+s81sbvB1dqoKT5uKCm06mg323BNOOQXGjoWlSze+f+06WLMOynJsOhtAnrq0iYik0vDhsGwZTJkSdiUikg6bDTxmlg+MA4YBuwGnm9luDQ77EjgHeKzBYzsCvwP2AQYBvzOzaKeJigqN8GSL667z7alvvnnj+1bW+MtcHuHRlDYRkZQ46ig/0qNpbSK5KZkRnkHAPOfcZ865tcATwHHxBzjnFjjnPgAazrE5CnjFOVfhnKsEXgGGpqDu9FHgyR677AI/+QmMGwdfffX9+6qDwJNrDQtATQtERFKsXTvfBFSBRyQ3JRN4ugML475fFNyWjC15bDgUeLLL737nJ103XMtTVQMlRbnZfEIjPCIiKTd8OHzyif8SkdySTOCxBLcl+5dWUo81sxFmNt3MppeXlyf51Gmwfr2aFmSbnj3hxBPh4Yf99DbwIx9Vq3Jz/Q5oHx4RkTQ49lh/+c9/hluHiKReMoFnEbBt3Pc9gMVJPn9Sj3XO3eecG+CcG9ClS5cknzoNVq70oUeBJ7uMHOmD6tNP++/XrIV1dbm5fgc2TGlbry5tIiKpsv32vh+OprWJ5J5kAs80oI+Z9TKz1sBpQLJvB/8CjjSzDkGzgiOD26KposJfKvBkl0MOgT594N57/fdVsfU7ORp4NMIjIpIWw4fD22/Dt9+GXYmIpNJmA49zrg64GB9U5gBPOedmm9n1ZjYcwMwGmtki4GTgXjObHTy2ArgBH5qmAdcHt0WTAk92MoMRI/yn1OzZfjqbGZQWhV1ZeijwiIikxfDhfvB84sSwKxGRVEpqHx7n3ETn3E7OuR2dczcFt13jnJsQXJ/mnOvhnCtxznVyzv0g7rEPOud6B18PpefHSJHKSn+pwJN9zjkHWrf2ozxVNT7s5OXovrpm/ktNC0REUqp/f9hmG01rE8k1OfoXYTPFRni08Wj26dwZTjoJ/va3IPDk6HS2GDON8Ig0g5nlm9l/zOyF4PteZjY12Bz7yWDqtrRQeXm+ecHLL8Pq1WFXIyKposATT1PastvIkVDWDurX527Dgpg8BR6RZvolfnp2zM3Abc65PkAlcF4oVUlkDB8ONTUweXLYlYhIqijwxNMIT3Y76CA49HB/PVdbUsdoSptIk5lZD+AY4P7gewMOBYIWjzwCHB9OdRIVhx4KxcWa1iaSSxR44lVUQEkJtGkTdiXSHGZwzHC/H8+8T8OuJr00pU2kOcYCvwZiPd07AcuD5jywmc2xI7NnnKRVYSEMGwZPPbXhPKiIZDcFnngVFZrOlu169oL5c+G++8KuJL3yTPvwiDSBmf0QWOqcmxF/c4JDGz2TEJk94yTtrrkGli+H3/wm7EpEJBUUeOIp8GQ352D1OqAeHn0UqqvDrih9LE8jPCJNcwAw3MwWAE/gp7KNBdqbWavgmKZsrC05bI894Be/8OfOpk4NuxoR2VIKPPEUeLJbTa0f9dhzd6iqgieeCLui9FHTApEmcc6NDrZP6InfQPv/nHNnAq8DJwWHnQ08H1KJEjHXXgtbbw0XXQT19WFXIyJbQoEnngJPdquq8Zf994S+ff2ePLlKTQtEUuUKYJSZzcOv6Xkg5HokItq2hdtug//8B+6+O+xqRGRLKPDEq6xU4MlmVasgPx+KCn2L6unTYcaMzT8uG6lpgUizOecmO+d+GFz/zDk3KNgc+2Tn3Jqw65PoOPlkOOII+O1v4Ztvwq5GRJpLgSfGOY3wZLuqGt+O2gzOOguKinJ3lCdPIzwiIulmBnfe6TchveyysKsRkeZS4ImprYU1a7QHT7Zav96v4YltONq+PZx2Gjz2GKxYEW5t6WAGTl3aRETSbaed4Ior4O9/h9dfD7saEWkOBZ6YWLN9jfBkp+pVfpQuFngAfv5zv132gw+GV1e6qGmBiEjGjB4NvXrBz34Ga9eGXY2INJUCT4wCT3arWuUvy4o33Lb33nDQQfCXv0BdXeLHZSvL05Q2EZEMKSryU9s+/hhuvTXsakSkqRR4YhR4sltVDRS0gjatv3/7qFHwxRfwj3+EU1e6qGmBiEhGHX00nHACXH89LFgQdjUi0hQKPDEKPNmtqsZPZ7MGG6cfeyzssEPunZLTlDYRkYwbO9ZPaXtAzctFsooCT4wCT/aqr4dVq78/nS0mPx/+539gyhR4993M15Yu2odHRCTjttsO9toL3ngj7EpEpCkUeGIUeLLXd+t3ShLff+650K6d30EuV2iER0QkFEOGwNSpvlW1iGQHBZ6Yykpo3RqKE4wSSLRV1fjLxgJPaSlccAE884xfz5MLNMIjIhKKwYP9LhbvvRd2JSKSLAWemNimow3XgEj0Va3yzQpaFzR+zCWX+Ms77shMTemmfXhEREJx0EH+LVjT2kSyhwJPTEWFNh3NVlU1idfvxNtuOzjpJBg/HqqqMlNXOuWpLbWISBg6dIDdd4d//zvsSkQkWQo8MbERHskudXWweg2UNjKdLd6oUbByJTz0UPrrSje1pRYRCc2QIfDOO7BuXdiViEgyFHhiFHiyU3Wtvywt2vyxgwbB/vv7vqL19emtK91iUy8VekREMm7wYFi1CmbMCLsSEUmGAk+MAk92qokFniSbTYwaBZ9/DhMmpK+mTMgLAo+mtYmIZNxBB/lLreMRyQ4KPDEKPNmpZhW0yt90w4J4xx8PPXtm/0akGuEREQlNt26wyy4KPCLZQoEH/LbJ1dUKPNmoutaP7iTbXS8/H37xC3jrLZg2Lb21pVNshEed2kREQjF4MLz5ZvbPkBZpCRR4wO/BAwo82cY5P6WtJIn1O/HOOw/KyuDmm9NTVyZY8E9XU9pEREIxZIjvg/PBB2FXIiKbo8ADCjzZqnYNrF+f/PqdmLZt4Ve/8huRTpmSntrSTVPaRERCNXiwv1R7apHoU+ABv34HFHiyTc0qf9nUER6Ayy+HrbaCSy/NztCgpgUiIqHq0QN22EHreESygQIPbAg82ng0u8RaUjcn8JSWwg03+BGep59ObV2ZoBEeEZHQDR7sA896LacUiTQFHtAIT7aqqYXiQshr5v/G557rt8u+4gpYsya1taVbngKPiEjYBg+GZctgzpywKxGRTVHgAQWebFW9CkqauH4nXn4+/PnPfl+eO+9MXV2ZEBvh0WlFEZHQDBniLzWtTSTaFHjABx4zaNcu7EokWXV1sGYtlDZjOlu8I4+EoUPhxhv9abpsoSltIiKh69ULundX4wKRqFPgAR94OnRo/tQoybwtWb/T0J/+5HuL3nDDlj9XpqhpgYhI6Mz8KM8bb+j8k0iU6S988IFH09myS00QeJrakjqRvn393jzjxsHcuVv+fJkQ24dHn7AiIqEaPBi+/hrmzQu7EhFpjAIPKPBko5pV0CofWhek5vmuvx4KC30Dg2ygER4RkUjQOh6R6FPgAb/xqAJPdqmu9aM7sbUsW2qrrXzYee657PjU0hoeEZFI2Hln6NJF63hEoiypwGNmQ83sEzObZ2ZXJri/jZk9Gdw/1cx6BrcXmNkjZvahmc0xs9GpLT9FNMKTXZzzU9pSsX4n3qhRfvXppZdGv/uZAo+ISCSYbdiPR0SiabOBx8zygXHAMGA34HQz263BYecBlc653sBtwM3B7ScDbZxzuwN7AyNjYShSYk0LJDvUrvGBJBXrd+IVF8NNN8H06TB+fGqfO9Xy1JZaRCQqhgyBL77wXyISPcmM8AwC5jnnPnPOrQWeAI5rcMxxwCPB9aeBw8zMAAeUmFkroAhYC6xMSeWp4hwsX67Ak01qVvnLVI/wAPz4x3D44X60J8orUDXCIyISGYMH+0uN8ohEUzKBpzuwMO77RcFtCY9xztUBK4BO+PBTA3wNfAn82TlXsYU1p1ZtrT9LXlYWdiWSrFS2pG4oLw8eeghat/bhp64u9a+RCnnq0iYiEhW77w7t2yvwiERVMoEn0arwhn9lNXbMIKAe2AboBVxqZjts9AJmI8xsuplNLy8vT6KkFKqu9pclJZl9XWm+mlVQXJi+fZN69IC77oJ334UxY9LzGlvK1KVNRCQq8vLgoIPUuEAkqpL5i3ERsG3c9z2AxY0dE0xfawdUAGcALzvn1jnnlgJvAwMavoBz7j7n3ADn3IAuXbo0/afYEjU1/rK0NLOvK81XXQslKV6/09Dpp8Npp8F118GMGel9rebQlDYRkUgZMsRv5fb112FXIiINJRN4pgF9zKyXmbUGTgMmNDhmAnB2cP0k4P+ccw4/je1Q80qAfYGPU1N6isRGeBR4skNdHaxZC6VpmM7W0Lhx0K0bnHWWn/oYJdqHR0QkUrSORyS6Nht4gjU5FwP/AuYATznnZpvZ9WY2PDjsAaCTmc0DRgGx1tXjgFJgFj44PeSc+yDFP8OWUeDJLt+t30nzCA/4VuUPPQQffwxXbtSNPVwa4RERiZS99vJ/SijwiERPq2QOcs5NBCY2uO2auOur8S2oGz6uOtHtkaLAk11qgsCTiREegCOOgF/8Am6/HX74Q/99FHwXeNSWWkQkClq1ggMOUOARiaI0rfrOIgo82aVmlf9UaV2QudccMwZ22QXOOcfv2RQVeXma0iYiEiGDB8OsWfDtt2FXIiLxFHgUeLJLda0f3bFEjQHTpKgIHn0Uli6Fiy/O3OtujpmmtImIRMiQIf7yrbfCrUNEvk+BR22ps4dzfkpbOvbf2Zy994arroLHH4fXXsv86yeSp8AjIhIlAwZAYaGmtYlEjQKP2lJnj9o1fpPY0gw0LEjkiiugVy+/pmfdunBqiGemKW0iIhHSpg3su6/24xGJGgWe2AhPcUh/REvyalb5yzBGeMCfths7Fj76CO68M5wa4mmER0QkcoYMgZkzYcWKsCsRkRgFnupqP50tT7+KyPuuJXVIgQfg2GNh6FC49lpYsiS8OkAjPCIiETR4sJ+M8PbbYVciIjH6K7+6WtPZskXNKiguDDecmsFf/uI3Ig17bx4ztaUWEYmYffeFggKt4xGJEgUeBZ7sUV0b3vqdeDvtBKNGwcMPw5Qp4dWRl6cpbSIiEVNcDAMHKvCIRIkCjwJPdqirgzVrw53OFu+qq2CbbeCSS6C+PpwaNKVNRCSShgyBadM29EUSkXAp8MTW8Ei0fbd+JwIjPOBD8p/+BDNmwIMPhlOD9uERaRIzKzSz98zsv2Y228yuC27vZWZTzWyumT1pZq3DrlWy2+DB/jzdu++GXYmIgAKPP/2iEZ7oqwkCT2lERngATj8dDjoIRo+GiorMv36eRnhEmmgNcKhzbk+gHzDUzPYFbgZuc871ASqB80KsUXLA/vv7WcdqTy0SDQo8mtKWHapqoKAVtC4Iu5INzOCOO6CyEq65JpzX1wiPSNKcF+xFQEHw5YBDgaeD2x8Bjg+hPMkhbdtC//5axyMSFQo8CjzR5xxUroT2Zf6P/CjZc0+46CK4+26/8UIm5alLm0hTmVm+mc0ElgKvAPOB5c65uuCQRUD3sOqT3DF4sJ/StmZN2JWIiAKPAk/01dTC2nXQoV3YlSR2/fXQqRNccIGftJ0palog0mTOuXrnXD+gBzAI2DXRYYkea2YjzGy6mU0vLy9PZ5mSA4YM8WHnvffCrkREFHgUeKKvcqW/7Ng23Doa07Gjn9o2fbrfoydTTG2pRZrLObccmAzsC7Q3s1bBXT2AxY085j7n3ADn3IAuXbpkplDJWgce6C81rU0kfC078NTVwerVCjxRV7nSbzjaJsKNk045BYYPh6uvhnnzMvOaalog0iRm1sXM2gfXi4DDgTnA68BJwWFnA8+HU6Hkko4dYffd1bhAJApaduCJNchXW+roqq+H5VXRnc4WYwZ33eW3177ggsyMvKhpgUhTbQ28bmYfANOAV5xzLwBXAKPMbB7QCXggxBolhwwZAu+8A+vWhV2JSMumwAMa4YmyFdX+j/qoTmeL1707/PnPMHkyjB+f/tfLU+ARaQrn3AfOub2cc3s45/o6564Pbv/MOTfIOdfbOXeyc07LzCUlBg/2f2r85z9hVyLSsrXswFMddCdV4ImuypV+JKNdlvw3Ov98OOQQuPxyWLQova+lpgUiIpE2eLC/1LQ2kXAp8IACT5RVrPBhJz8/7EqSY+ZHd9at8+2q0zkCoyltIiKR1q0b7LyzGheIhE2BBxR4omrNWli1GjpGfP1OQzvuCDfeCC+8AE8+mb7XyTNYr314RESibMgQePNNvyRVRMKhwAMKPFEVa0fdIQvW7zT0y1/CoEFwySWQrv06LPjnq1EeEZHIGjwYVqyADz8MuxKRlkuBB9SlLaoqVkLrAigpCruSpsvPhwce8J9yF14Iy5al/jXM/KUCj4hIZMXW8Wham0h4WnbgUZe26HLOj/B0aLvhD/ts07cvXHstPPssbLUVHH00/PWvPgSlQl7we1HjAhGRyNp2W+jZE956K+xKRFqulh14NKUtuqpX+Y1hs3E6W7zf/AZmzIBRo+Cjj+Dss/0q1hNOgCee2BC6m0MjPCIiWWHgQP9RICLhUOABBZ4oqghGQbI98AD07w833wyffw5Tpvjube+9B6efDl26wIym5JQAACAASURBVCmnwDPPQG1t0543T4FHRCQbDBgAn30GFRVhVyLSMinwtGoFrVuHXYk0VLkSSov9Gp5cYQb77gu33QYLF/qNGc4911+edJIPP2ecAc8/D6tXJ/d8oE5tIiIRN2CAv9Qoj0g4FHhKS7N3jUiuqquHlTW5MbrTmLw8v5J13Dj46it49VU480yYNAmOPx623homTNj8c4BGeEREIq5/f385fXq4dYi0VAo8paUwax58uzzsaiRmeZX/Iz6XA0+8Vq3gsMPg3nvh66/h5Zehd28ffG65pfFAY2paICKSDdq3hz59FHhEwqLAs3V3WLZ8w5oRCV/lCj960a4Frq0qKICjjvLT3H70I7jsMhg5Etat2/hYNS0QEckaAwYo8IiEpWUHnpoa2HFHf33N2nBrkQ0qV0L7sg1Ttlqi4mJ48knf5W38eBg2DCorv3+MmhaIiGSNAQPgyy9h6dKwKxFpeVrwX5T4EZ5te/rrCjzRULvGf3VsIdPZNiUvD266CR5+2O9Yt99+MH/+hvs1pU1EJGuocYFIeBR4tunuryvwRENlrB11u3DriJKzz/ZNDcrLYZ994J13/O3fTWlTlzYRkajbay//tq1pbSKZp8DTdSt/va4e6uvDrUf8dLY2raGoTdiVRMvgwTB1KpSUwKWX+tvyNMIjIpItyspgl10UeETCoMDTsfOGtSJrEiwMl8xZvx4qq/x0NrUK31jv3nDAAbBsmf9eTQtERLKKGheIhKNlBx4Mior9AnnQtLawrazxo2yazta40lIf1EH78IiIZJkBA2DxYv8lIpmTVOAxs6Fm9omZzTOzKxPc38bMngzun2pmPePu28PMppjZbDP70MwKU1f+FnAOOnXx12ML5BV4wvVtpR+1UMOCxsUHHjUtEBHJKmpcIBKOzQYeM8sHxgHDgN2A081stwaHnQdUOud6A7cBNwePbQU8ClzonPsBcDAQjXlja9fCttv56+0VeELnnN/8tWM7yM8Pu5roigUe5zSlTUQky/Tr5wfnNa1NJLOSGeEZBP/f3t3HyV3W9/5/fZK9S3YgQAgRwq2AlCBSNYB3hVZBaK1BbC14ePykRw+Qtng8x2qLbfV4rPbGPqrWFq147IGHpxUUbY2W4q/FG9TDragIghICSABJAiGwudvd7HX+uGbYyWbDzmxm5zvzndfz8djHd+Y739n97LWTfOc91/W9LtaklNamlEaBq4FzphxzDnBV9fa1wGsiIoDXAnemlH4IkFJ6IqXUGTMDjIzAEUfBxE5YOAT9fQaeIo1sze1/4H5FV9LZKpUccLZtc9ICSeoyCxfCCScYeKR2ayTwLAMerru/rrpv2mNSSuPAZmAx8AIgRcTXIuKOiPiDvS+5RWqBZ2w0f1I+OGDgKdKG6qKaiw08z6lSyduREXt4JKkL1SYu8L9uqX0aCTzTTZc19Z/pno7pA14FXFDdnhsRr9ntB0RcHBG3R8TtGzZsaKCkFhgZgSOOhFqH0+CAs7QVaeNTefKI/r6iK+ls9YFnnuvwSFK3WbEC1q+Hhx+e+VhJrdFI4FkHHFZ3/1Bg6vwizx5TvW5nEfBkdf+3UkobU0pbgeuAl0z9ASmlK1JKK1JKK5YsWdL8bzEbz4zA8w6ZfNM42G8PT1G2bINt2+HA/YuupPPt0sNT/efrkDZJ6hq1iQsc1ia1TyOB5zbg2Ig4KiIGgPOB1VOOWQ1cWL39m8DXU0oJ+BrwoohYWA1CpwM/bk3pe2nL1rwdrPYoDA64+GhRNlaHs3n9zsym7eEx8EhSt3jRi6Cvz8AjtdOM44dSSuMRcSk5vMwH/iGldHdEfAC4PaW0GvgM8NmIWEPu2Tm/+txNEfERcmhKwHUppX+do9+lOTvG8m+/oDpL9uBAdf8oLFxQWFk9aeNTsM/w5N9AezbdNTz28EhS1xgaghNPNPBI7dTQBRMppevIw9Hq972v7vZ24E17eO7/IU9N3VkmEoyPw77VYVS1N9vbDTxttX1HnqHt+YcWXUl3qAWeLVuctECSutSKFXDttbuuMCBp7jS08Gg5zYdHHoZ99sl3a4Fn1IkL2mrjU3nrcLbG1PfwQD5TGngkqausWAGbNsEDDxRdidQbejfw9A/Agw9MvoEc7M/b7U5c0FYbN8HwgsmhhXpuw8N5Wws88wImnKVNkrqJExdI7dWbgWdiAgaH4GcPTL6BnDfPxUfbbXQMNo/Yu9OM3Xp45tnDI0ld5oUvhIEBA4/ULr0ZeLbtyG8UH3sU5s+f3D/k4qNt9URtOJvTUTdswYI8jK2+h8fAI0ldZWAATjrJwCO1S28Gnq3b8nbD47vuHzDwtNWGTTA0mIe0qTHz5uVeyfpreJylTZK6zooV8L3vOSpZaofeDDxbtudPxTdv2nX/0ECerlpzb3wcnnomD2dziprmVCr28EhSl1uxAp5+GtasKboSqfx6M/Bs3Q5PP5V7dOoN9OeFR8ddfHTOPbE5v1F3OFvz6gOPPTyS1JWcuEBqnx4NPNtgw/rJC8BrhuoWH9Xc2vhUDpj7DhddSfeZGnjs4ZGkrrN8eV6E1MAjzb3eCzwp5R6exx7ZPfAMGnjaYudOeHKzw9lma7chbQ4Al6Ru09cHL36xgUdqh94LPNt35NDz8EMGnqJsejpfpelwttlxSJsklcKKFXDHHflzQElzp/cCz5btefvA/bsHnoHq4qMGnrm1YRP0zYdFlZmP1e52CTyuwyNJ3WrFCtiyBe66q+hKpHLrvcBTm5L6vp/uHnjmzcuhx5na5s6O0Rx4lhyQ21vNmzqkzR4eSepKZ50F/f1w1VVFVyKVW++949y6PYea9Y/n9UymGuy3h2curXs890gc9ryiK+leTlogSaWwdCm88Y1w5ZWwbVvR1Ujl1ZuBZ8EQbN26ew8P5Ot4DDxzY2wcHt0ABx0ACwaLrqZ7uQ6PJJXGqlWwaRN84QtFVyKVV28FnpTykLaB+fm+gae9Hnk8T1Zw+MFFV9LdKhUYHYWxMSctkKQud/rpcNxx8Pd/X3QlUnn1VuAZHYOdE1CbCXlPgWfnBIyPt7W00hvfCY+sh8X7wfCCoqvpbrXX7ZYt1SFtTkstSd0qIvfy3HQT/PCHRVcjlVNvBZ4t1QGyE9VJCfYUeMCJC1rtsQ059BzutTt7rXbt2chInvjBHh5J6mpveUtehPRTnyq6EqmceivwbK1OSb2jup028Dg1dctNTOTJCvbbB/Z1Kuq9Vnvdjow4aYEklcABB8B558FnPwvPPFN0NVL59F7g6ZufhwLBDD08Bp6W+fnGPJzQa3daoz7wOGmB1LCIOCwivhER90TE3RHxjur+AyLi3yPivurWVZHVdqtW5f/WP/e5oiuRyqfHAs82WDgEW6uBZ7ppqV18tLVSgod/DvsM5x4e7b2pPTwOaZMaNQ78fkrpeOBlwO9FxHLgMuCGlNKxwA3V+1JbnXoqnHRSnrzAz7Gk1uqxwLMdFi6YnNJ3uh6eZxcfNfC0xPonYfto7t2JmPl4zcwhbdKspJQeSyndUb39DHAPsAw4B6gt/XgV8IZiKlQvq01e8P3vw223FV2NVC69E3jGxvI6MAuHnjvwQHVqaict2Gspwc8ey7OyLV5UdDXlMd2QNkOP1JSIOBJ4MXALsDSl9BjkUAQcVFxl6mUXXJD/i3eKaqm1eifwbKlOVDA8Qw8PuBZPqzzxVO5VO+x59u600tQeHjDwSE2IiArwReC/pZSebuJ5F0fE7RFx+4YNG+auQPWsffbJoefqq/NipJJao3cCT22GtoZ6ePrzMCzfRM5erXdnaBAOOqDoaspllx6e6j9hX6tSQyKinxx2/jGl9KXq7scj4uDq4wcD66d7bkrpipTSipTSiiVLlrSnYPWcSy6BbdvyjG2SWqO3As+8ebn3pvZGcWho+mMHB/JUyjt3trfGMnniKXhmq707c6F+HZ5a2zpxgTSjiAjgM8A9KaWP1D20GriwevtC4Mvtrk2qefGL4ZRTnLxAaqXeCTxjY3kygoj8RrFS2fMb8drU1Nsd1jYr23fATx7Mwweft7joasqnvx8GBx3SJjXvlcD/B7w6In5Q/fo14C+AMyPiPuDM6n2pMKtWwT33wLe/XXQlUjn0FV1A24yNQ3/1192yZfopqWvq1+KpLJz72spkYgJ+vDa/AV9+9OSQK7VWpTI5aQEYeKQGpJS+A+ypy/k17axFei7nnQf//b/DJz8Jp51WdDVS9+udd6NjYzBQDTy1Hp49eTbwOFNb09aug2e2wHFH5uulNDdqgcchbZJUOgsXwoUXwhe/COunvaJMUjN6J/CM1vXwzBh4XHx0VjZsgkfWw7KDYIkTFcyp3Xp4JoqtR5LUUpdckj+rvfLKoiuRul9vBJ6UqkPaqkFmpsAT4eKjzdq2PV+3s88wPP/QoqspP3t4JKnUli+H00+HT30qjxaXNHu9EXh2TuTQ02gPD7gWTzN2TsCP788j45c/3+t22mF4OF+LFk5LLUlltWoVrF0L//EfRVcidbfeeGc6Vr0Wp5nAM2Tgadj9P4ORbfALR+V1dzT3pg5ps4dHkkrn3HNhyZI8RbWk2euRwDOet80EnoGBPGmBn5w/t8efgMc25vV2Fu9XdDW9Y+qQNl+nklQ6g4Pw1rfC6tXwyCNFVyN1rx4LPNVreGaalhpgqD8Pmh138dE9enoEfvoQLKrAUcuKrqa3GHgkqSdcfHFeB/0znym6Eql79VjgafIaHnBY255s2wF3rclTfS8/es+LuGpu7DakzStaJamMnv98OOss+PSnYXy86Gqk7tQbgWe0eg3PQB+MjuYvA8/sjY3DXfflXoUTj80z2qm9KpXcU1ljD48kldaqVbBuHVx3XdGVSN2pocATEWdHxE8iYk1EXDbN44MRcU318Vsi4sgpjx8eESMR8a7WlN2ksfH8Sfi8eZNvEg08szMxAXffn3t4TjgGFi4ouqLeVKnkkDO6I9838EhSaf36r8Mhhzh5gTRbMwaeiJgPXA78KrAceHNELJ9y2NuATSmlY4CPAn855fGPAv+29+XOUm0Nnog8DAgamLSgtvjo2NzW1k1SytfsbH4GjjsS9tun6Ip6V+31u21b3jpLmySVVl8fXHQRXH89PPBA0dVI3aeRHp5TgDUppbUppVHgauCcKcecA1xVvX0t8JqIfFFHRLwBWAvc3ZqSZ2FsfNfrd2DmwBPhWjxTPfRYnpXtiENg6eKiq+lttdfv1q15aw+PJJXaf/kv+a3Jpz9ddCVS92kk8CwDHq67v666b9pjUkrjwGZgcUQMA38I/M+9L3UvzCbwAAz2G3hqfr4RHno0B50jDi66Gk0NPPbwSFKpHXoovP71eba2Ud+aSE1pJPBMN/3W1HdXezrmfwIfTSmNPOcPiLg4Im6PiNs3bNjQQElNGhubDDy1a3hmmpYa7OGp2bApD2Xbbx94wRHOyNYJaoGn9nq2h0eSSm/VKli/Hv7lX4quROoujQSedcBhdfcPBR7d0zER0QcsAp4ETgU+HBEPAv8N+KOIuHTqD0gpXZFSWpFSWrFkyZKmf4kZzbqHpxp4evXN5OgY/Pj+/DW8AE44Ok/8oOLtFnicllqSyu61r4Ujj3TyAqlZjbx7vQ04NiKOiogB4Hxg9ZRjVgMXVm//JvD1lP1SSunIlNKRwMeAP0sp/V2Lam/MxATsnJhcdLTZwDORJqe17hUp5Wt1brsbNj4FRx4CL/6FfNWkOkPt9Vt7PTukTZJKb948uOQS+MY34M47i65G6h4zBp7qNTmXAl8D7gE+n1K6OyI+EBErq4d9hnzNzhrgncBuU1cXZrS6StfALHp4arOQbXiy9XV1qh2jeUHRex+ABYPw0uV5kgJ7djrLs4Hnmbzt1V5ISeoxb30r7L8/vOY1cOONRVcjdYeGPrJPKV0HXDdl3/vqbm8H3jTD93j/LOrbe2PV3pnZDGmrLIR9huGxjbBsabmvXZmYyBMTrF2Xr746+jBYdlC5f+duZg+PJPWkgw6Cm2+GlStz6Ln8crj44qKrkjpb+T+2H6v28EwNPI1MWgBwyBLYuh02P+e8C91rfCc8/HO49S6472c54K04AQ4tecDrdrXX75Yt+e9kD48k9YwXvABuuQXOOCMPcbv00snPdyXtrocCT901PIODk/dnsmR/mD8fHpuD2eOKtGM09+bcfGfeDg3AC4+BF70gD2VTZ1uwYHIh3XkGHknqNYsWwVe/Cu96V+7lOesseOKJoquSOlP5r0Kf2sOzZUvjvTuQw87SA/KwtmPqZnvrVlu2wbqfw+NP5jfJS/aHQ58H+zbRJipeRB7WNjKSbzukTZJ6zvz58Fd/BSeeCBddBCefDKtXwwtfWHRlUmfp8nfvDaj18fbNz9uRkcau36l38BJ4dEOeuezQpa2trx1Sgs3PwMOPw5Ob8wQEBy/Jv4u9Od2rPvA4LbUk9ay3vCUPczv3XHj1q2HdOhgYKLoqqXP0xpC2/r7J61FmE3ienbxgQ3cNHUoJ1j8Jd9wDP/wpPLMlTzH9shfBsYcbdrpdLfDMm9ddr0tJUsu97GVwxRWwYQN861tFVyN1lvL38IyOw0Dd9TqzCTyQJy/4yYN58oLadNWdZGIih7ux8bxu0JZt8Mj6fK3OgkE49ghYuhjmlz/j9gyHtEmS6pxxRr7Ec/VqOPPMoquROkf5A8/YlOtuZht4luwPax7OvTxFBp6dO+HpLXmI2uYR2L4j/447pxnStG8FjjkcFi9yxrUyeraHx0kLJEk57Lz2tTnwfPzjnvqlmt4IPJUFk/dHRuDww5v/PkVMXpBS7qEZ2ZrDzeaRfLv25rayMIea/r7qV//k7cF+WDA09zWqOJUKbNxoD48k6VkrV8KXvwx33gknnVR0NVJn6IHAMwb9dT0ys+3hgdlPXpBSXu9mdCz30KSUF/esPVbbbh+Fbdth247Jbe3xiHwd0aFLcw/TvpXJiRjUmyoVePBB1+GRJD3rda/Lp4XVqw08Uk25A08taNT3xjQ7LXW9+skLlh00fV/x9h2w7vEcXkbHJr8afUMaka+5WTAEByyChUP5a5/hfHG6VFM/pG3CWdokSbB0aZ7AYPVqeO97i65G6gzlDjxTFx2FvevhgcnJC54egUV1PUcTE3mSgAcfBVIOLAP9OawM9MNAX97O74NgMiw9uwUGB/KXg27ViF2mpbaHR5KUrVwJ73kPPPIILFtWdDVS8Xok8FR/zZT2PvDUJi94dMNk4Hl6BH76UJ4ZbfF+eaKAISfA1xzbZZY2e3gkSVkt8Hz1q3DJJUVXIxWv3GOkRquLjtYCz7ZtOfTsTeCpTV6wYVO+xua+h+D79+ZwdcLR8MJjDDtqj0plcmFdJy2QJFUdfzwcfXQe1iap7IGn1sNTW4dnZCRv9ybwQJ68ICW47a7c07PsIDj5hXDg/nv3faVm1F7HEzsd0iZJelZE7uW54YbJtz5SL+uNwFPr4WlV4KkshAP2zdNdv+T4PITNGdPUbrXJN3YaeCRJu1q5EnbsgH//96IrkYrXG4GnFka2bMnbvQ08ACe+AF6yPM+eJhWh9joeH3dImyRpF698Jey/v8PaJCh94BnLYac2nXOth2e201JLnaQ+8NjDI0mq098Pv/ZreeKCnTuLrkYqVskDz/iua/C0akib1AmeDTxjztImSdrNypWwcSPcfHPRlUjF6oHAM2UNHjDwqBxqr+OxJha2lST1jLPOym+DHNamXtcDgcceHpVU7XU8OmrgkSTtZtEi+OVfNvBIBh6pW9UHHictkCRNY+VKuPde+OlPi65EKk55A09KBh6VW+11vGNHfr3byyNJmuL1r8/br3yl2DqkIpU38IxX1yYZqAs8W7bk1bgWLCiuLqlVarMN7thebB2SpI51xBFw0kkOa1NvK2/geXbR0SmTFixcODlNtdTN+vpgaAi2VwOPM7VJkqaxciV85zvwxBNFVyIVo7zv/MfG8nbqkDaHs6lMKhXYvi3fdkibJGkaK1fmz8S+/OWiK5GKUeLAU+vhMfCoxCoV2Fbr4THwSJJ295KXwPHHw6WXwjXXFF2N1H4GHqmbVSqwdWu+bQ+PJGka8+bBN7+Zg8/558Of/ImjoNVbeiDwTLmGx8CjMqlUYOuWfNvAIz2niPiHiFgfEXfV7TsgIv49Iu6rbvcvskZprhx0EHz96/C2t8GHPgRvfCM880zRVUntUd7AMzqeP9KYX/crGnhUNvU9PA5pk2ZyJXD2lH2XATeklI4Fbqjel0ppYAA+/Wn4+Mfhq1+Fl78c1q4tuipp7pU38IyN7TolNeRpqQ08KpPh4fy6Bnt4pBmklG4Enpyy+xzgqurtq4A3tLUoqc0i4O1vh+uvh0cfhZNPhm98o+iqpLlV4sAzZdFRyD08tbVLpDKoVCYX1LWHR5qNpSmlxwCq24MKrkdqizPOgFtvhaVL4cwz4ROfKLoiae6UPPD077rPIW0qm0qlrofHK1CluRQRF0fE7RFx+4YNG4ouR9prxxwDN98MZ58Nv/d7sGoVjI4WXZXUeiUPPNP08Bh4VCaVCoxUrzq1h0eajccj4mCA6nb9ng5MKV2RUlqRUlqxZMmSthUozaV9983r81x2GXzqU7m3xzyvsumdwDM+nlekN/CoTCqVyWl2vIZHmo3VwIXV2xcCLs2onjN/Pvz5n8M//mMe5nbyyXDnnUVXJbVOOQPPzp15gvn6wFMb9mPgUZlUKjnMg4FHmkFEfA64CTguItZFxNuAvwDOjIj7gDOr96We9J/+E9x4Y5736RWvgC99qeiKpNbom/mQLrSnNXjAwKNyqVRgfCzfdkib9JxSSm/ew0OvaWshUgc7+WS4/XY491z4jd+Aj30M3vGOoquS9k5DPTwRcXZE/CQi1kTEbmsURMRgRFxTffyWiDiyuv/MiPheRPyoun11a8vfg9Fa4LGHRyVXqeSP4sAeHklSSxx8MHzzm7ByJbz73XDvvUVXJO2dGQNPRMwHLgd+FVgOvDkilk857G3AppTSMcBHgb+s7t8IvD6ldCJ5bPRnW1X4c6q9Aaxfh6fWw+O01CqTSgV21oa0OUubJKk1hobyIqXDw3kGNz9TUzdrpIfnFGBNSmltSmkUuJq8UFu9+oXbrgVeExGRUvp+SunR6v67gaGIGGxF4c9pbJoeHoe0qYwqlcnXu0PaJEktdNBBeTKDr38dPve5oquRZq+RwLMMeLju/rrqvmmPSSmNA5uBxVOO+Q3g+ymlHbMrtQkGHvWK+mt4/PhNktRiF12Ur+t55zth8+aiq5Fmp5HAE9Psm/rO6jmPiYgTyMPcLpn2B7R6MbexcYjI8yzWGHhURvWztNnDI0lqsfnz4ZOfhPXr4b3vLboaaXYaCTzrgMPq7h8KPLqnYyKiD1gEPFm9fyjwz8BbUkr3T/cDWr6YW20NnqjLYQYelZGTFkiS5thLXwq/+7tw+eVwxx1FVyM1r5HAcxtwbEQcFREDwPnkhdrq1S/c9pvA11NKKSL2A/4VeE9K6butKnpGY2O7DmcDA4/KySFtkqQ2+OAH4cAD4Xd+Jy91KHWTGQNP9ZqcS4GvAfcAn08p3R0RH4iIldXDPgMsjog1wDuB2tTVlwLHAO+NiB9Uvw5q+W8x1ej47oGnNi21s7SpTIaH64a0eQaSJM2N/faDv/5ruPVW+F//q+hqpOY0tPBoSuk64Lop+95Xd3s78KZpnvdB4IN7WWPzxsZhwZRgMzICfX0wMND2cqQ5MzQ0GXTs4ZEkzaELLshh57LL8sKkrbgKQWqHhhYe7Tpj0/Tw/OxnsO++u17XI3W7iDysbWKnkxZIkuZUBHziE/DMM/CHf1h0NVLjyhd4JiZg585dA89DD8E118Bv/VZxdUlzpVLJr3t7eCRJc2z5cvj934f//b/hO98puhqpMeULPNOtwfOhD+WPJf7oj4qpSZpLlUoO+QYeSVIbvPe9cPjheQKD2kShUicrceDpz9sHHsgfQ1x0ERx22J6fJ3WrWuBxSJskqQ2Gh+Fv/gbuugv+9m+LrkaaWYkDT7WH54MfzKtmvec9xdUkzaXa1NT28EiS2uScc+B1r4P/8T9g3bqiq5GeWwkDT7Vvtb8P7r8frroKVq2CZcuKrUuaK5VKDvoGHklSm0TAxz+eV0Z45zuLrkZ6buULPKN1PTx/+qd5aJtTiajMKpUc9F2HR5LURs9/PvzxH8MXvgBf+1rR1Uh7Vr7AUxvS9sBa+Oxn4Xd/Fw4+uNiapLlUqcDYqD08kqS2e/e74dhj4dJLYfv2oquRplfOwNPfl6/dGRyEP/iDoiuS5lalAqM7nLRAktR2g4Nw+eWwZg18+MNFVyNNr5yBZ2In/NM/5Y8bli4tuiJpblUqMGoPjySpGGeeCeedB3/2Z/nyaanTlDDwjMHDP4MFC3I/q1R2tcCzc2fRlUiSetRHPgIDA/mzZj9/U6cpX+AZ2QL3/RTe/nZYsqToaqS5V6nkaXLGx4uuRJLUow45BD7wAbj+evjnfy66GmlX5Qw8W7fAu95VdCVSewwPVwOPPTySpOJceimcdFIeYGMvjzpJuQLPnXfC4BAcdxwsXlx0NVJ71KaldkibJKlAfX1w8cWwdi088EDR1UiTyhV4Dj8C5s2DV76i6Eqk9qlUYNx1eCRJxTv99Lz91reKrUOqV67A0z+Qt/vtW2wdUjvVruEx8EiSCnb88XmQzY03Fl2JNKmv6AJaanAATjga9hkuuhKpfWqBxwHTkqSCzZsHp51m4FFnKVcPT998OHD/HHykXlG7hse8I0nqAKedlq/jWbeu6EqkrFyBR+pFz/bwFF2IJEmT1/HYy6NOYeCRul1t0oIouhBJkuBFL4J993XiAnUOA4/U7YaH85A2E48kqQPMnw+vepU9POocBh6p282fDxEwz8AjSeoMp58O994L69cXXYlk4JHKYV5AzHOmNklSRzjttLy1l0edPefFFwAADFtJREFUwMAjlcG86j9lA48kqQO89KWwcKGBR53BwCOVQV91Sa0JA48kqXj9/fCKVzhxgTqDgUcqg/nz8zZNFFuHJElVp58OP/oRPPlk0ZWo1xl4pDKwh0eS1GFOOy2PtP7ud4uuRL3OwCOVQS3weA2PJKlDnHIKDA46rE3FM/BIZTBgD48kqbMMDcGppzpxgYpn4JHKoH8gb+3hkSR1kNNOgzvugGeeKboS9TIDj1QGgwYeSVLnOf102LkT/u//LboS9TIDj1QGA4N5Oz5ebB2SJNV5+cvzZaZex6MiGXikMhiqBp6t24qtQ5KkOsPDeRFSr+NRkQw8UhnUAs82A48kqbOcfjrcequnKBXHwCOVwdBQ3m7dWmwdkiRNcdppMDYGN99cdCXqVQYeqQwWLMjbbduLrUOSpCle9SqIcFibitNQ4ImIsyPiJxGxJiIum+bxwYi4pvr4LRFxZN1j76nu/0lEnNW60iU9qxZ4tjteQJqNmc5zkmZv0SL4xV904gIVZ8bAExHzgcuBXwWWA2+OiOVTDnsbsCmldAzwUeAvq89dDpwPnACcDXyi+v0ktdLChXm7fUexdUhdqMHznKS9cPrpcNNNMDpadCXqRX0NHHMKsCaltBYgIq4GzgF+XHfMOcD7q7evBf4uIqK6/+qU0g7ggYhYU/1+N7WmfEkADC+E0R0wMgKbny66GnWKyjDM9zOmBjRynpO0F047DT72Mfj2t+HUU4uuRp1iaChPWz7XGvkRy4CH6+6vA6a+VJ89JqU0HhGbgcXV/TdPee6yWVcraXr77gubNsDiQ+AHPy26GnWKA4fhhOOLrqIbNHKek7QXfumX8nU8Z5xRdCXqJDfcAK9+9dz/nEYCT0yzb+py7ns6ppHnEhEXAxcDHH744Q2UJGkXRx4Bd/zQdXi0q+PbcBYpB89V0hw78ED44hfh/vuLrkSd5Jhj2vNzGgk864DD6u4fCjy6h2PWRUQfsAh4ssHnklK6ArgCYMWKFbudZCQ14I0ri65A6laeq6Q2OPfcoitQr2pklrbbgGMj4qiIGCBPQrB6yjGrgQurt38T+HpKKVX3n1+dxe0o4Fjg1taULklSSzRynpMkdakZe3iq1+RcCnwNmA/8Q0rp7oj4AHB7Smk18Bngs9VJCZ4knyyoHvd58oWf48DvpZR2ztHvIklS0/Z0niu4LElSizQ0L0JK6Trguin73ld3ezvwpj0890PAh/aiRkmS5tR05zlJUjk0tPCoJEmSJHUjA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0oqUUtE17CIiNgAPAYuAzXUPHQhsnKMfO/Vnteo5Mx2zp8en2z9130z356q9ZtNWjT6vqPbqtNdWo8+zvZp73nMdszdtNd2+Tm6vRp+zCNgvpbSk6ap6QPVc9RTdf56a6bhuPE9N97Na9RzP6809z/NUc89r1b/F6faXtb0aO0+llDryC7hiyv3b2/WzWvWcmY7Z0+PT7Z+mPWa6PyftNZu26vT26rTXlu3V/vbam7bqtvZq9Dmz/Vv00lcZzlMzHdeN56m5bC/P693RXp322mpFezXTVrbX7l+dPKTtKx3+sxp5zkzH7Onx6fZP3TfT/bky259je7X+ebZXc897rmP2pq2m29fJ7dXoc9r5f3C3KsN5aqbjuvH/kdn+LP/fbf3zbK/mnteqf4vT7e+19tpFxw1p25OIuD2ltKLoOrqF7dU426o5tldzbK/e4d+6ObZXc2yvxtlWzemF9urkHp6prii6gC5jezXOtmqO7dUc26t3+Lduju3VHNurcbZVc0rfXl3TwyNJkiRJzeqmHh5JkiRJaoqBR5IkSVJpGXgkSZIklZaBR5IkSVJplSLwRMRwRHwvIn696Fo6XUQcHxF/HxHXRsTvFF1Pp4uIN0TEpyPiyxHx2qLr6XQR8fyI+ExEXFt0LZ2q+v/VVdXX1QVF16P28VzVGM9TzfE81RzPUzMr43mq0MATEf8QEesj4q4p+8+OiJ9ExJqIuKyBb/WHwOfnpsrO0Yr2Sindk1JaBfwWUOo511vUXv+SUroI+G3gvDkst3Ataq+1KaW3zW2lnafJtnsjcG31dbWy7cWqaZ6rGud5qjmep5rjeWr2ev08VXQPz5XA2fU7ImI+cDnwq8By4M0RsTwiToyIr075OigizgB+DDze7uILcCV72V7V56wEvgPc0N7y2+5KWtBeVX9SfV6ZXUnr2qvXXEmDbQccCjxcPWxnG2vU7F2J56pGXYnnqWZcieepZlyJ56nZupIePk/1FfnDU0o3RsSRU3afAqxJKa0FiIirgXNSSn8O7DYMICJ+BRgm/6G2RcR1KaWJOS28IK1or+r3WQ2sjoh/Bf5p7iouVoteXwH8BfBvKaU75rbiYrXq9dWLmmk7YB35ZPIDiv/QSQ3wXNU4z1PN8TzVHM9Ts9fr56lCA88eLGMyVUJu9FP3dHBK6Y8BIuK3gY1lPIHMoKn2iohfJndVDgLXzWllnamp9gLeDpwBLIqIY1JKfz+XxXWgZl9fi4EPAS+OiPdUTzi9ak9t93Hg7yLidcBXiihMLeG5qnGep5rjeao5nqdmr2fOU50YeGKafWmmJ6WUrmx9KV2hqfZKKX0T+OZcFdMFmm2vj5P/4feqZtvrCWDV3JXTVaZtu5TSFuA/t7sYtZznqsZ5nmqO56nmeJ6avZ45T3ViN9U64LC6+4cCjxZUSzewvZpjezXH9po9267c/Ps2zrZqju3VHNtr9nqm7Tox8NwGHBsRR0XEAHA+sLrgmjqZ7dUc26s5ttfs2Xbl5t+3cbZVc2yv5thes9czbVf0tNSfA24CjouIdRHxtpTSOHAp8DXgHuDzKaW7i6yzU9hezbG9mmN7zZ5tV27+fRtnWzXH9mqO7TV7vd52kdKMQ44lSZIkqSt14pA2SZIkSWoJA48kSZKk0jLwSJIkSSotA48kSZKk0jLwSJIkSSotA48kSZKk0jLwqOtFxMgcfM8HI+LAIn72bOqQJHUuz1NSsQw8kiRJkkrLwKNSiojXR8QtEfH9iPiPiFha3f/+iLgqIv7/6qdSb4yID0fEjyLi+ojor/s2746IW6tfx1Sff1RE3BQRt0XEn9b9vEpE3BARd1S/1znT1PQ7EfHhuvu/HRF/W739LxHxvYi4OyIunua5R0bEXXX33xUR76/ePrpa+/ci4tsR8QvV/W+KiLsi4ocRcePetqkkqXU8T3meUvsYeFRW3wFellJ6MXA18Ad1jx0NvA44B/g/wDdSSicC26r7a55OKZ0C/B3wseq+vwE+mVI6Gfh53bHbgXNTSi8BfgX464iIKTVdC7yx7v55wDXV229NKb0UWAH814hY3MTvegXw9urz3wV8orr/fcBZKaWTgJVNfD9J0tzzPOV5Sm3SV3QB0hw5FLgmIg4GBoAH6h77t5TSWET8CJgPXF/d/yPgyLrjPle3/Wj19iuB36je/izwl9XbAfxZRJwGTADLgKXUnWxSShsiYm1EvAy4DzgO+G714f8aEedWbx8GHAs8MdMvGREV4BXAF+rOW4PV7XeBKyPi88CXZvpekqS28jzleUptYuBRWf0t8JGU0uqI+GXg/XWP7QBIKU1ExFhKKVX3T7Drv4nUwO2aC4AlwEurJ6kHgaFpjrsG+C3gXuCfU0qpWt8ZwMtTSlsj4pvTPHecXXtka4/PA55KKf3i1B+UUloVEaeSPw38QUT8YkppxpOTJKktPE95nlKbOKRNZbUIeKR6+8JZfo/z6rY3VW9/Fzi/evuCKT9vffUk8ivAEXv4nl8C3gC8mclhAouATdWTyC8AL5vmeY8DB0XE4ogYBH4dIKX0NPBARLwJILKTqrePTindklJ6H7CR/ImcJKkzeJ7yPKU2sYdHZbAwItbV3f8I+ZOyL0TEI8DNwFGz+L6DEXEL+YOBN1f3vQP4p4h4B/DFumP/EfhKRNwO/ID8ydhuUkqbIuLHwPKU0q3V3dcDqyLiTuAn1XqnPm8sIj4A3EIe9lD//S8APhkRfwL0k8eC/xD4q4g4ljyM4YbqPklS+3me8jylAsVkL6kkSZIklYtD2iRJkiSVloFHkiRJUmkZeCRJkiSVloFHkiRJUmkZeCRJkiSVloFHkiRJUmkZeCRJkiSV1v8DOIW5DGn2jTgAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Single-voxel performance evaluation.\n",
"voxidx = 3\n",
"\n",
"print(\"Results: subject\", subid, \", voxel id\", voxidx)\n",
"\n",
"myfig = plt.figure(figsize=(14,7))\n",
"ax_corr = myfig.add_subplot(1, 2, 1)\n",
"plt.title(\"Correlation coefficient\")\n",
"plt.xlabel(\"Lambda values\")\n",
"ax_corr.set_xscale('log')\n",
"ax_corr.plot(todo_lambda, dict_corr_tr[voxidx], label=\"train\", color=\"red\")\n",
"ax_corr.plot(todo_lambda, dict_corr_te[voxidx], label=\"test\", color=\"pink\")\n",
"ax_corr.legend(loc=1,ncol=1)\n",
"\n",
"ax_spar = myfig.add_subplot(1, 2, 2)\n",
"plt.title(\"Sparsity via l0-norm\")\n",
"plt.xlabel(\"Lambda values\")\n",
"ax_spar.set_xscale('log')\n",
"ax_spar.plot(todo_lambda, dict_l0norm[voxidx], color=\"blue\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"続いて、全ボクセルまで対象範囲を拡げ、ROIごとの平均的な成績を見ていく。まずは、$\\lambda$の候補として多数の値を試しているので、代表的な成績を選ぶ必要がある。訓練、検証、それぞれにおいて相関係数の絶対値がもっとも高かったものを選出する。"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"num_voxels = len(dict_corr_tr)\n",
"\n",
"best_corr_tr = np.zeros(num_voxels, dtype=np.float32)\n",
"best_corr_te = np.zeros(num_voxels, dtype=np.float32)\n",
"\n",
"for v in range(num_voxels):\n",
" \n",
" # Best absolute correlation value.\n",
" best_corr_tr[v] = np.max(np.abs(dict_corr_tr[v]))\n",
" best_corr_te[v] = np.max(np.abs(dict_corr_te[v]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"次はROIの上でのループである。あらかじめ用意した階層型データはここで役に立つ。"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"dict_roi_corr_tr = {}\n",
"dict_roi_corr_te = {}"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"f = tables.open_file(\"data/vim-2/response.h5\", mode=\"r\")\n",
"tocheck = f.get_node((\"/sub\"+str(subid)), \"idx\")\n",
"for idxnode in tocheck._f_iter_nodes():\n",
" idx = idxnode.read()\n",
" roi_name = idxnode._v_name\n",
" dict_roi_corr_tr[roi_name] = np.mean(best_corr_tr[idx])\n",
" dict_roi_corr_te[roi_name] = np.mean(best_corr_te[idx])\n",
"\n",
"f.close()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0IAAAGrCAYAAADtiicXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xu4ZXV95/n3RyqAXCyIXOQmpQ0OrSmwnVI7KkK3bcRG28wEATUKTlpCMzqxWxJN1IjOaIOBJ6bbRIKxhygxhZSaQbFRYgcUIzpVgkCpKJKSsqq8oFBy0ULLb/+x1sHN8Zw6+5y9z6Xq9349z35qXX5rre/ea+9T67N/a62dqkKSJEmSWvKIxS5AkiRJkhaaQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSQOSrEhSSZbNcfk/SvJX465rrpL89yRn7GD+pUn+nx3Mvy/J4+enOs2nJCcm+fYO5l+c5M0LWdMwkuyR5CtJHjPGdb4syafG3XaIdX0kyUnjWJek8TMISZpWkg1JftwfDN+d5KokR4xpvf9mHDUupqkONKvqHVX17xerpsmq6vlV9dcASc5Mcv0sl9+nqu6Yy7b7kPVg//75YZJrkhwzqc3hSf4myQ+S3J/ki0leMKlNJTlqLjVoelV1dlX936OsY5j3VJJTk/xjkgeSXDvEas8CPlNV3+mX32FYH0ZV/U1V/ca42w7hfODtY1qXpDEzCEmayQurah/gEOC7wH9d5HqGNlWvzlx7ejRn7+zfP4cBm4D3TcxI8qvA9cCDwJOAA4A/BT6Y5JRFqHVaSXZb7Bp2Yj8E3kUXCobxu8AHhl35Uv5MV9UXgUclWbXYtUj6ZQYhSUOpqp8Aa4AnTkzrT2G5MMmdSb7bn2rzyH7eAUk+nuSevjfgs0kekeQDwGOBj/U9BX8w1faSvCjJTUl+lOSbE6eXJDk0yZX9Om9P8qqBZc5LsibJZUl+BJw5zbRHJHlDv94fJPlQf1A+VR2vTPLVJPcmuSPJ7/bT9wb+O3Bo/zzu62s7L8llA8v/uyTr+9fh2iT/fGDehiTnJrk5ydYklyfZc6Z9keRx/foe0Y//VZLvDcy/LMlr++Frk/z7frsXA7/e13rPwCr373v77k3yhST/bGBdD/XG9N/M//l0bXekqn4MfAh48sDk/wjcB/xOVX2nqn5cVX9L9w36RUkyxGtxcpIb+/fJxiTnDcy7OsmrJ7X/cpL/vR8+pu+l+mGS25KcOtDu0iTvSfKJJPcD/2pH2+qXeUWSb/XvqTdnoOdzpvdc/x546TTPccrPUj/vYb1lmaL3JN3pmnf19bxsurZJXpDuM3dPuh6cYwfmHZHuNK/v9/W/e4b31EOq6u+r6kPA5qnmT6r1scA/A77Qj58FvAz4g34bH+unb0jy+iQ3A/cnWTbw+t6b7tS6/21gvQ/ruepft7OTfCNdb/efT7zfZtl2tyQX9a/vPyV5dX751NprgZNneu6SFp5BSNJQkuwFnAbcMDD5AuAJdAe3R9F96//H/bzXAd8GDgQOBv4IqKp6OXAnfU9TVb1zim09DXg/8PvAfsCzgQ397L/t13socArwjiTPGVj8RXSBbT/gb6aZ9n8Bvwmc0K/nbuDPp3nq3wNeADwKeCXwp0meUlX3A88HNvfPY5+qetiBXpIn9PW+tn8dPkEXAHcfaHYqcBLwOOBY4Mxp6nhIVf0T8CPgX/STjgfuyy9C1rOB6yYt81XgbODzfa37Dcx+CfBWYH/gdnZ8Ks9s2j4kXXB8Sb/MhOcCH66qn09q/iG6sPyEIVZ9P/AKun17MvAfkvxmP++D/TYnangicCRwVV/PNX2bg/p2f5HkSQPrfmn//Pal67madlv9uv+C7qD9EGA53edhwg7fc1V1bFV9cJrnOOVnaYjXBuAxdD1thwFnAJck+V8mN0ryFOC/0fXGPBr4S+DKdF927AZ8HPgWsKJf1+oZ3lNztRK4o6p+BlBVl9B9Zt/Zb+OFA21fQrcf9uvbf5Pus7Cc7j16WZJDdrCtFwBPBY6j+xw+bw5tX0X3d+DJwFPo9vFkX+2Xk7TEGIQkzeTv+m96f0R34PonAP03oq8C/mNV/bCq7gXeAZzeL/dTugPCI6vqp1X12aoa9uDtd4D/VlXXVNXPq2pTVX0t3fVJzwJeX1U/qaqbgL8CXj6w7Oer6u/65X48zbTfBd5YVd+uqm3AecApmeIUm6q6qqq+WZ3rgE/RHWwN4zTgqv55/BS4EHgk8IyBNv+lqjZX1Q+Bj/HwHpMduQ44Ib+4oHxNP/44utD25SHXA/CRqvpifzD5NzPUMJu2AOf275976fbd4L46ANgyxTJbBubvUFVdW1W39Pv2ZrrgeUI/+6PAk5Mc2Y+/rK9/G92B7Yaq+n+r6mdV9SXgw3ThesL/V1Wf69f9kxm2dQrwsaq6vqoepPtCYPD9PvR7bgqjfJYA3lxV2/r371V0B/KTvQr4y6r6QlVt768r2wb8S+BpdOHt96vq/v61mNW1ZrOwH917ZRj/pao2TnzOq+qK/rP086q6HPgGXe3TOb+q7qmqO4F/YMfv5enangr8Wb9f72bq0//u7Z+XpCXGICRpJr/Zf9O7B/Bq4Lr+4PtAYC9gXX8qzT3A1f106ALT7cCn0p1S9oZZbPMIum93JzsUmAhdE77Fw7953zjFcpOnHQl8dKDurwLb6b5tf5gkz09yQ39K0j3Av2WIA/SBer81MdL3fGycVO93BoYfAPYZct3XASfS9f58hu70mxP6x2en6GXZkdnUMGXb/vSriVMELx5oc2H//lkB/BgY7I24i+4Af7JDBubvUJKnJ/mH/pStrXQ9FAcA9O+Tq/hFOD+dX/QSHgk8feI90O/bl9H1oEx42PtmR9ui29cPta+qB4AfDCw+9HtuCqN8lu7uey8nfKuvdbIjgddNej2O6NseAXxropdmnt1N1wM3jMn75xUDp/bdA/waO/6sjvy+Z9J+n1xTb19gytMGJS0ug5CkofTfEn+E7uDtWXQHqT8GnlRV+/WP5f2F8VTVvVX1uqp6PPBC4D8NnMI207fZG+muE5hsM/CrSQYPlB5LdxH+Q6VOVf4U63/+QN37VdWeVTW4HpLsQddLcCFwcH9A/wlg4tqVmZ7HZroDzIn1he6gctO0SwzvOrqeqRP74euBZ9IFoeumWWY2vQizUt3d8iZOETx7ivl3Ar8H/Fn668iAvwd+K/31LgNOpdtHXx9i0x8ErgSOqKrldNesDF5b9LfAS5L8Ol1v3D/00zcC1016D+xTVf9hsOxZbGsLcPhEw/45Pnpg2aHec1OZ4bP0AN0XEhMm33J6//40wAmPZeprdTYCb59U317VXbO1EXjsNL1X435P3Qw8ftK2ptvGQ9P7Xr/30n1Z8+j+s3orD38vzIeH7Xe6z/dk/5zZ9dBKWiAGIUlDSedFdNeGfLXvcXgv3TUzB/VtDkvyvH74BUmO6g/+f0QXoLb3q/susKPfpnkf8Mokz0l3kflhSY6pqo3APwL/Ocme6S7m/h1+8S3/sC4G3j5xylSSA/vnNtnudD1h3wd+luT5wOBtdb8LPDrJ8mm28yHg5P55/ArdtR7b+ucwo/6i6xOnmldV36ALor9Nd6vhH/X1/BbTB6HvAodPukZpwVTVNXQH4Wf1k/6U7jS+9yV5TL9PXwK8ke40rGEOsvel6yX8SX9t2eQbDnyCLoy+Dbh8oKfs48ATkrw8ya/0j6cOXGc1222tAV6Y5Bn96/tWHn4QvsP3XLqL/8+caqMzfJZuAl6a7qL9k/jFqXqD3ppk9yTH050SeMUUbd4LnN33eiXJ3uluDrEv8EW6A/7z++l7Jnlmv9yM76m+tj2BZcAj+uV/Zaq2VfVtfvmUtpn+XgDsTReMvt9v85V0PULz7UPA7/V/o/YDXj9FmxPobqwiaYkxCEmayceS3Ed3APZ24IyqWt/Pez3dKTs3pLsj29/zi1Ofju7H7wM+D/xFVV3bz/vPwJv6U1jOnbzB6m45+0q6A+WtdAf2Ez0rL6E7zWoz3TUgb+kPsGfjz+i+2f9UknvpbgDx9CnquJfuIvcP0Z2y89J+uYn5X6Prcbijfy6HTlr+Nrqg8l/petBeSHeTiAdnKjDJ4XSv3S07aHYd8IO+t2ViPMCN07T/H8B64DtJZjztbJ78Cd0dwPaoqh/Q9S7uCXyF7lSy/wS8vL/GYxjnAG/r9+Mf0+2rh/TX43wE+Dd0PToT0++lC7Wn072XvkN384895rKt/jPxGmA1XWi4l+5GG9v6JtO+5/oQ8WgefiOSQTv6LP0e3ftq4tS+v5u07Hfo3rub6b4wOLt/3z5MVa2lu07o3X372+lv3FFV2/ttHEV3o5Nv013/BsO9p15OF9rfQ9eL+WO64DWdv+Th15K9D3hi/xmb/Pwm6v8KcBHd6/NdupsufG4H2xiX99JdN3gz3efuE8DP6INqkqcC9/d/0yQtMZnd9ZaSpIWQ5LfpTjv8w8WuRbOXZB+6cHJ0dXf521HbZwH/Z1W9ZEftxi3J+4Hbq+ptC7ndmfSnpN4IPKeqprqZxpLV9xpfXFUTPX8fBt5XVZ9Y3MokTcUgJEnSGCR5IfBpul65i+h6fJ4yyzu8LYj+GpzPAe+uqqF/vFQP118L9q/oeoUOprum8Iaqeu2iFiZpKJ4aJ0nSeLyI7hS0zXSns52+FENQ7zt0PVYfXuxCdnKhux7sbrperK/yi99Sk7TE2SMkSZIkqTn2CEmSJElqzjC/aL1kHHDAAbVixYrFLkOSJEm7qB/8oPst5Ec/+tEztNRStW7duruq6sCZ2u1UQWjFihWsXbt2scuQJEnSLurSSy8F4Mwzz1zUOjR3Sb41TDtPjZMkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKas2yxC5iNWzZtZcUbrlrsMiRJkjSkDeefvNglSFOyR0iSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmjPWIJTkyCTrktyUZH2Ss3fQdkOSA8a5fUmSJEkaxrIxr28L8Iyq2pZkH+DWJFdW1ebBRkl2G/N2JUmSJGloc+4RSnJBknMGxs8DXlNV2/pJewyuP8l9Sd6W5AvAr/eTX5PkS0luSXLMXGuRJEmSpNkY5dS41cBpA+OnAlckOSLJzcBG4IKB3qC9gVur6ulVdX0/7a6qegrwHuDcqTaS5Kwka5Os3f7A1hHKlSRJkqTOnINQVd0IHJTk0CTHAXdX1Z1VtbGqjgWOAs5IcnC/yHbgw5NW85H+33XAimm2c0lVraqqVbvttXyu5UqSJEnSQ0a9WcIa4BS6nqHVgzP6nqD1wPH9pJ9U1fZJy0+cRred8V+vJEmSJElTGjUIrQZOpwtDa5IcnuSRAEn2B54J3DbiNiRJkiRprEbqhamq9Un2BTZV1ZYkzwUuSlJAgAur6pZxFCpJkiRJ4zLy6WhVtXJg+Brg2Gna7TNpfMXA8FrgxFFrkSRJkqRhjPUHVSVJkiRpZ2AQkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnNG/kHVhbTysOWsPf/kxS5DkiRJ0k7OHiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpObsVLfPvmXTVla84arFLkOSJEmLYIM/o6IxskdIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkpoz9iCU5Mgk65LclGR9krMH5m1IckCSFUluHfe2JUmSJGkYy+ZhnVuAZ1TVtiT7ALcmubKqNs/DtiRJkiRp1kbqEUpyQZJzBsbPA15TVdv6SXvsYBu7JXlv32v0qSSPHKUWSZIkSRrWqKfGrQZOGxg/FbgiyRFJbgY2AhdM0xt0NPDnVfUk4B7gt6baQJKzkqxNsnb7A1tHLFeSJEmSRgxCVXUjcFCSQ5McB9xdVXdW1caqOhY4CjgjycFTLP5PVXVTP7wOWDHNNi6pqlVVtWq3vZaPUq4kSZIkAeO5WcIa4BS6nqHVgzP6nqD1wPFTLLdtYHg783O9kiRJkiT9knEEodXA6XRhaE2Swyeu90myP/BM4LYxbEeSJEmSxmLkXpiqWp9kX2BTVW1J8lzgoiQFBLiwqm4ZdTuSJEmSNC5jOR2tqlYODF8DHDtNuxX94F3Arw1Mv3AcdUiSJEnSMMb+g6qSJEmStNQZhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktScsfyg6kJZedhy1p5/8mKXIUmSJGknZ4+QJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzdqrbZ9+yaSsr3nDVYpchSZKkOdjgz6BoCbFHSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaM+sglOTIJOuS3JRkfZKzB+bdN+Q6Lk1yymy3LUmSJEnjsGwOy2wBnlFV25LsA9ya5Mqq2jzMwknmsk1JkiRJGpsd9ggluSDJOQPj5wGvqapt/aQ9Jq8jyUVJvpTk00kO7Kddm+QdSa4Dfq9v+uwk/5jkDnuHJEmSJC2kmU6NWw2cNjB+KnBFkiOS3AxsBC4Y6A3aG/hSVT0FuA54y8Cy+1XVCVV1UT9+CPAs4AXA+dMVkOSsJGuTrN3+wNahn5gkSZIkTWeHQaiqbgQOSnJokuOAu6vqzqraWFXHAkcBZyQ5uF/k58Dl/fBldEFnwuU83N9V1c+r6ivAwUyjqi6pqlVVtWq3vZbP4qlJkiRJ0tSGuVnCGuAUup6h1YMz+p6g9cDx0yxbA8P3T5q3bWA4Q9QhSZIkSWMxTBBaDZxOF4bWJDk8ySMBkuwPPBO4bWB9E9f7vBS4frzlSpIkSdLoZryDW1WtT7IvsKmqtiR5LnBRkqLrybmwqm7pm98PPCnJOmArD7++SJIkSZKWhKFuZV1VKweGrwGOnabdPv3gmydNP3HS+JnTLCdJkiRJ827WP6gqSZIkSTs7g5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaM9QPqi4VKw9bztrzT17sMiRJkiTt5OwRkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqTqpqsWsY2h6HHF2HnPGuxS5DkiRJu6iTdv8aAFc/eMyU8zf4Uy5LXpJ1VbVqpnb2CEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkpoz5yCU5Mgk65LclGR9krMH5t03zTKXJjmlH96Q5IC5bl+SJEmS5mrZCMtuAZ5RVduS7APcmuTKqto8ptokSZIkaV4M1SOU5IIk5wyMnwe8pqq29ZP2mLyuJBcl+VKSTyc5cJpVv6Zvc0uSY+ZQvyRJkiTN2rCnxq0GThsYPxW4IskRSW4GNgIXDPQG7Q18qaqeAlwHvGWa9d7Vt3kPcO5UDZKclWRtkrXbH9g6ZLmSJEmSNL2hglBV3QgclOTQJMcBd1fVnVW1saqOBY4CzkhycL/Iz4HL++HLgGdNs+qP9P+uA1ZMs+1LqmpVVa3aba/lw5QrSZIkSTs0m5slrAFOoesZWj04o+8JWg8cP82yNc30iVPrtjPa9UqSJEmSNLTZBKHVwOl0YWhNksOTPBIgyf7AM4HbBtZ7Sj/8UuD68ZQrSZIkSaMbuhemqtYn2RfYVFVbkjwXuChJAQEurKpb+ub3A09Ksg7YysOvL5IkSZKkRTWr09GqauXA8DXAsdO026cffPOk6WcODK8YGF4LnDibWiRJkiRprub8g6qSJEmStLMyCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5s/pB1cW28rDlrD3/5MUuQ5IkSbuoSy/9PgAXn+kx567OHiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOakqha7hqHtccjRdcgZ71rsMiRJkrSLOmn3rwFw9YPHALDBn27Z6SRZV1WrZmpnj5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNWdOQSjJkUnWJbkpyfokZ89hHdcmWTWX7UuSJEnSKJbNcbktwDOqaluSfYBbk1xZVZuHWTjJbnPcriRJkiSNbMYeoSQXJDlnYPw84DVVta2ftMfgepK8J8navqforQPTNyT54yTXAy/uJ784yReTfD3J8eN4QpIkSZI0k2FOjVsNnDYwfipwRZIjktwMbAQuGOgNemNVrQKOBU5IcuzAsj+pqmdV1ep+fFlVPQ14LfCWqTae5Kw+WK3d/sDWWTw1SZIkSZrajEGoqm4EDkpyaJLjgLur6s6q2lhVxwJHAWckObhf5NQkXwJuBJ4EPHFgdZdPWv1H+n/XASum2f4lVbWqqlbtttfyoZ+YJEmSJE1n2GuE1gCnAI+h6yF6SFVtTrIeOD7JOuBc4KlVdXeSS4E9B5rfP2m9E6fXbZ9FLZIkSZI0kmHvGrcaOJ0uDK1JcniSRwIk2R94JnAb8Ci6sLO17yF6/vhLliRJkqTRDNULU1Xrk+wLbKqqLUmeC1yUpIAAF1bVLQBJbgTWA3cAn5unuiVJkiRpzoY+Ha2qVg4MX0N3M4Sp2p05zfQVk8ZPHBi+i2muEZIkSZKkcZvTD6pKkiRJ0s7MICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmDP2DqkvBysOWs/b8kxe7DEmSJO2iLr30+wBcfKbHnLs6e4QkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkpqTqlrsGoa2xyFH1yFnvGuxy5AkSdIu6qTdvwbA1Q8es8iV7Bw2LMGftkmyrqpWzdTOHiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJas6cg1CSI5OsS3JTkvVJzh5imWuTrOqH75vrtiVJkiRpFMtGWHYL8Iyq2pZkH+DWJFdW1eYx1SZJkiRJ82KoHqEkFyQ5Z2D8POA1VbWtn7TH4LqSvCfJ2r6n6K07WO/bk3w5yQ1JDp7bU5AkSZKk2Rn21LjVwGkD46cCVyQ5IsnNwEbggoHeoDdW1SrgWOCEJMdOsc69gRuq6jjgM8CrptpwkrP6ULV2+wNbhyxXkiRJkqY3VBCqqhuBg5IcmuQ44O6qurOqNlbVscBRwBkDvTqnJvkScCPwJOCJU6z2QeDj/fA6YMU0276kqlZV1ard9lo+9BOTJEmSpOnM5hqhNcApwGPoeogeUlWbk6wHjk+yDjgXeGpV3Z3kUmDPKdb306qqfnj7LGuRJEmSpDmbzV3jVgOn04WhNUkOT/JIgCT7A88EbgMeBdwPbO17iJ4/3pIlSZIkaTRD98JU1fok+wKbqmpLkucCFyUpIMCFVXULQJIbgfXAHcDn5qFuSZIkSZqzWZ2OVlUrB4avobsZwlTtzpxm+okDw/sMDK+hO/VOkiRJkubdnH9QVZIkSZJ2VgYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNWdWP6i62FYetpy155+82GVIkiRpF3Xppd8H4OIzPebc1dkjJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUnJ3q9tm3bNrKijdctdhlSJIkaQFs8GdTNI/sEZIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5ow1CCV5cpLPJ1mf5OYkp+2g7bVJVo1z+5IkSZI0jGVjXt8DwCuq6htJDgXWJflkVd0z2CjJbmPeriRJkiQNbc49QkkuSHLOwPh5wAur6hsAVbUZ+B5wYD9/Q5I/TnI98OJ+sRcn+WKSryc5fq61SJIkSdJsjHJq3Gpg8NS3U4ErJkaSPA3YHfjmQJufVNWzqmp1P76sqp4GvBZ4y1QbSXJWkrVJ1m5/YOsI5UqSJElSZ85BqKpuBA5KcmiS44C7q+pOgCSHAB8AXllVPx9Y7PJJq/lI/+86YMU027mkqlZV1ard9lo+13IlSZIk6SGjXiO0BjgFeAxdDxFJHgVcBbypqm6Y1P7+SePb+n+3j6EWSZIkSRrKqOFjNfBe4ADghCS7Ax8F3l9VV+xwSUmSJElaJCPdPruq1gP7ApuqagvddULPBs5MclP/ePIY6pQkSZKksRn5dLSqWjkwfBlw2TTtVkwaP3Fg+C6muUZIkiRJksZtrD+oKkmSJEk7A4OQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktScZYtdwGysPGw5a88/ebHLkCRJkrSTs0dIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5O9Xts2/ZtJUVb7hqscuQJEnSErfBn1zRDOwRkiRJktQcg5BNDeQTAAAMoElEQVQkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmjD0IJXlyks8nWZ/k5iSnDcy7Nsmqfvi+cW9bkiRJkoaxbB7W+QDwiqr6RpJDgXVJPllV98zDtiRJkiRp1kbqEUpyQZJzBsbPA15YVd8AqKrNwPeAA6dZ/u1JvpzkhiQHj1KLJEmSJA1r1FPjVgOnDYyfClwxMZLkacDuwDenWHZv4IaqOg74DPCqqTaQ5Kwka5Os3f7A1hHLlSRJkqQRg1BV3QgclOTQJMcBd1fVnQBJDgE+ALyyqn4+xeIPAh/vh9cBK6bZxiVVtaqqVu221/JRypUkSZIkYDzXCK0BTgEeQ9dDRJJHAVcBb6qqG6ZZ7qdVVf3w9jHVIkmSJEkzGkf4WA28FzgAOCHJ7sBHgfdX1RU7XFKSJEmSFsHIt8+uqvXAvsCmqtpCd53Qs4Ezk9zUP5486nYkSZIkaVzGcjpaVa0cGL4MuGyadicODO8zMLyG7hQ7SZIkSZp3Y/9BVUmSJEla6gxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJas5YflB1oaw8bDlrzz95scuQJEmStJOzR0iSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTmpqsWuYWh7HHJ0HXLGuxa7DEmSJO2iTtr9awBc/eAxY1/3Bn8GZkEkWVdVq2ZqZ4+QJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5hiEJEmSJDVn7EEoydVJ7kny8RnaXZtk1bi3L0mSJEkzmY8eoT8BXr6jBkl2m4ftSpIkSdJQ5hyEklyQ5JyB8fOSvK6qPg3cO0X7DUn+OMn1wIv7yS9O8sUkX09y/FxrkSRJkqTZGKVHaDVw2sD4qcAVMyzzk6p6VlWt7seXVdXTgNcCb5lqgSRnJVmbZO32B7aOUK4kSZIkdeYchKrqRuCgJIcmOQ64u6runGGxyyeNf6T/dx2wYprtXFJVq6pq1W57LZ9ruZIkSZL0kGUjLr8GOAV4DF0P0UzunzS+rf93+xhqkSRJkqShjBo+VgPvBQ4AThi9HEmSJEmafyPdNa6q1gP7ApuqagtAks/SXSv0nCTfTvK80cuUJEmSpPEZ+XS0qlo5aXzKu79V1YpJ4ycODN/FNNcISZIkSdK4zcfvCEmSJEnSkmYQkiRJktQcg5AkSZKk5hiEJEmSJDXHICRJkiSpOQYhSZIkSc0xCEmSJElqjkFIkiRJUnNG/kHVhbTysOWsPf/kxS5DkiRJu6hLL/0+ABef6THnrs4eIUmSJEnNMQhJkiRJao5BSJIkSVJzDEKSJEmSmmMQkiRJktQcg5AkSZKk5qSqFruGoe1xyNF1yBnvWuwyJEmStIs6afevAXD1g8csciXjt6GRn6FJsq6qVs3Uzh4hSZIkSc0xCEmSJElqjkFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWrOvAShJFcnuSfJxydNvzbJqn74vvnYtiRJkiTNZL56hP4EePk8rVuSJEmSRjJSEEpyQZJzBsbPS/K6qvo0cO8Qy789yZeT3JDk4FFqkSRJkqRhjdojtBo4bWD8VOCKIZfdG7ihqo4DPgO8aqpGSc5KsjbJ2u0PbB2pWEmSJEmCEYNQVd0IHJTk0CTHAXdX1Z1DLv4gMHEN0TpgxTTbuKSqVlXVqt32Wj5KuZIkSZIEwLIxrGMNcArwGLoeomH9tKqqH94+plokSZIkaUbjCB+rgfcCBwAnjGF9kiRJkjSvRr5rXFWtB/YFNlXVFoAkn6W7Vug5Sb6d5HmjbkeSJEmSxmUsp6NV1cpJ48dP0+7EgeF9BobX0J1iJ0mSJEnzbr5+R0iSJEmSliyDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUnGWLXcBsrDxsOWvPP3mxy5AkSdIu6tJLvw/AxWd6zLmrs0dIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5BiFJkiRJzTEISZIkSWqOQUiSJElScwxCkiRJkppjEJIkSZLUHIOQJEmSpOYYhCRJkiQ1xyAkSZIkqTkGIUmSJEnNMQhJkiRJak6qarFrGFqSe4HbFrsOzegA4K7FLkI75D5a+txHOwf309LnPto5uJ+Wvp1pHx1ZVQfO1GjZQlQyRrdV1arFLkI7lmSt+2lpcx8tfe6jnYP7aelzH+0c3E9L3664jzw1TpIkSVJzDEKSJEmSmrOzBaFLFrsADcX9tPS5j5Y+99HOwf209LmPdg7up6Vvl9tHO9XNEiRJkiRpHHa2HiFJkiRJGplBSJIkSVJzlkwQSnJSktuS3J7kDVPM3yPJ5f38LyRZMTDvD/vptyV53kLW3ZK57qMkz02yLskt/b//eqFrb8kon6V+/mOT3Jfk3IWquTUj/r07Nsnnk6zvP1N7LmTtLRnhb96vJPnrfv98NckfLnTtrRhiHz07yZeS/CzJKZPmnZHkG/3jjIWrui1z3UdJnjzwt+7mJKctbOVtGeWz1M9/VJJNSd69MBWPSVUt+gPYDfgm8Hhgd+DLwBMntTkHuLgfPh24vB9+Yt9+D+Bx/Xp2W+zntKs9RtxH/wI4tB/+NWDTYj+fXfUxyn4amP9h4Arg3MV+PrviY8TP0jLgZuC4fvzR/r1bkvvppcDqfngvYAOwYrGf0672GHIfrQCOBd4PnDIw/VeBO/p/9++H91/s57SrPUbcR08Aju6HDwW2APst9nPaFR+j7KeB+X8GfBB492I/n9k8lkqP0NOA26vqjqp6EFgNvGhSmxcBf90PrwGekyT99NVVta2q/gm4vV+fxmvO+6iqbqyqzf309cCeSfZYkKrbM8pniSS/SXdAsH6B6m3RKPvoN4Cbq+rLAFX1g6ravkB1t2aU/VTA3kmWAY8EHgR+tDBlN2XGfVRVG6rqZuDnk5Z9HnBNVf2wqu4GrgFOWoiiGzPnfVRVX6+qb/TDm4HvAQcuTNnNGeWzRJL/FTgY+NRCFDtOSyUIHQZsHBj/dj9tyjZV9TNgK923ocMsq9GNso8G/RZwY1Vtm6c6Wzfn/ZRkb+D1wFsXoM6WjfJZegJQST7Zn6LwBwtQb6tG2U9rgPvpvsG+E7iwqn443wU3aJT//z12WBhjeZ2TPI2up+KbY6pLDzfn/ZTkEcBFwO/PQ13zbtliF9DLFNMm39d7ujbDLKvRjbKPupnJk4AL6L7V1vwYZT+9FfjTqrqv7yDS/BhlHy0DngU8FXgA+HSSdVX16fGWKEbbT08DttOdzrM/8Nkkf19Vd4y3xOaN8v+/xw4LY+TXOckhwAeAM6rql3ojNBaj7KdzgE9U1cad8dhhqfQIfRs4YmD8cGDzdG360w2WAz8cclmNbpR9RJLDgY8Cr6gqv9GZP6Psp6cD70yyAXgt8EdJXj3fBTdo1L9311XVXVX1APAJ4CnzXnGbRtlPLwWurqqfVtX3gM8Bq+a94vaM8v+/xw4LY6TXOcmjgKuAN1XVDWOuTb8wyn76deDV/bHDhcArkpw/3vLmz1IJQv8/cHSSxyXZne6i0ysntbkSmLiryynA/6ju6qwrgdP7u/c8Djga+OIC1d2SOe+jJPvR/SH7w6r63IJV3KY576eqOr6qVlTVCuBdwDuqaue6+8vOYZS/d58Ejk2yV3/gfQLwlQWquzWj7Kc7gX+dzt7AvwS+tkB1t2SYfTSdTwK/kWT/JPvTnanwyXmqs2Vz3kd9+48C76+qK+axRo2wn6rqZVX12P7Y4Vy6/fVLd51bshb7bg0TD+DfAl+nO//zjf20twH/rh/ek+5OVrfTBZ3HDyz7xn6524DnL/Zz2VUfc91HwJvozpe/aeBx0GI/n131McpnaWAd5+Fd45bkPgJ+m+5mFrcC71zs57IrP0b4m7dPP309XVD9/cV+LrvqY4h99FS6b7vvB34ArB9Y9v/o993twCsX+7nsqo+57qP+b91PJx07PHmxn8+u+hjlszSwjjPZye4al75wSZIkSWrGUjk1TpIkSZIWjEFIkiRJUnMMQpIkSZKaYxCSJEmS1ByDkCRJkqTmGIQkSZIkNccgJEmSJKk5/xNarj7m2bdAUQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Training\n",
"xvals = list(dict_roi_corr_tr.keys())\n",
"yvals = list(dict_roi_corr_tr.values())\n",
"myfig = plt.figure(figsize=(14,7))\n",
"plt.barh(range(len(dict_roi_corr_tr)), yvals)\n",
"plt.yticks(range(len(dict_roi_corr_tr)), xvals)\n",
"plt.title(\"Best correlation, within-ROI average; subject \"+ str(subid)+\" (training)\" )\n",
"plt.axvline(x=np.mean(np.array(yvals)), color=\"gray\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0cAAAGrCAYAAAALo+xTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X+cXXV95/HXW1L5lQCpBDQFiRZdV0qg7ojrDwy71lpbXfvYRkCtgttKWR6ydbe0ttUq+litsbBdtz90se5SxTZIxC7+WJW6JYordScEgalFFCiBRAUMMRANGj77xzkjh3EmuXPvnbmZ5PV8PO5jzu/zOfneOznv+Z5zbqoKSZIkSdrfPWbUBUiSJEnS3sBwJEmSJEkYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkvZKSVYkqSSL+lz/95P8xbDr6leS/53krN3MvzTJf97N/AeSPHluqtNMkhyY5B+SPH4e9vWqJJ8d0rauTPILw9iWpP2L4UjSXi/JHUm+154gb03yySTHDmm7PzeMGkcpyWlJ7upOq6p3VtWvj6qmqarqxVX1lwBJzk5y7SzXX1xVt/Wz7zZ4PdS+f76T5OokT5uyzDFJPpzkviQPJvlykpdMWaaSHN9PDQvYOcDnq+qbsOcQ26vpwn9Vfbiqfn7QbbfeBbxjSNuStB8xHElaKF5aVYuBJwDfAv5kxPX0bLren357hNS3d7fvn58C7gY+MDkjyU8C1wIPAScARwJ/DPxVktUjqHVGSQ6Y513+BvChed7nwKrqy8BhScZGXYukhcVwJGlBqarvA+uAp09Oay/9uSjJnUm+leR9SQ5u5x2Z5BNJ7m97Db6Q5DFJPgQ8Efh426PwO9PtL8nLktyQ5LtJvjF5qU6S5Umuarf59SSv66xzYZJ1SS5L8l3g7BmmPSbJ77bbvS/JR9oT9enqeG2SrybZnuS2JL/RTj8U+N/A8vY4HmhruzDJZZ31/02Sifbf4Zok/7wz744kFyS5Mcm2JJcnOWhPbZHkSe32HtOO/0WSb3fmX5bkDe3wNUl+vd3v+4Bnt7Xe39nk0rZXcHuSv0/y051t/ajXpu29+LOZlt2dqvoe8BHg5M7k/wg8APxaVX2zqr5XVX9N0/NwcZL08G/xS0k2tu+TTUku7Mz7dJLXT1n+K0n+bTv8tLY36ztJbklyeme5S5O8N8mnkjwI/Kvd7atd5zVJ/ql9T/1BOj2ks3zPPRH4aeDv2/FzgFcBv9O23cfb6cuTfDTJPUluT/IfOts4Jcl4W+u3kvyXdtbn25/3t9t6dqb0KLZtfm6SW9P0GP/ZZFskOSDJxUnubff5+vz4ZajXAL+0+5aTpEczHElaUJIcApwBXNeZvAZ4Ks0J7/E0vQNvaef9FnAXsAw4Gvh9oKrq1cCdtD1SVfXuafZ1CvBB4LeBI4DnA3e0s/+63e5yYDXwziQv6Kz+MpoQdwTw4Rmm/Qfgl4FV7Xa2An82w6F/G3gJcBjwWuCPkzyjqh4EXgxsbo9jcVVtnnIcT23rfUP77/ApmlD42M5ipwO/ADwJWAmcPUMdP1JVtwPfBX62nXQq8EAneD0fWD9lna8C5wJfams9ojP7FcDbgKXA19n9ZVGzWfZH2jD5inadSS8EPlpVD09Z/CM0AfqpPWz6QeA1NG37S8C/T/LL7by/avc5WcPTgeOAT7b1XN0uc1S73J8nOaGz7Ve2x7eEpodrxn212/5zmhDzBOBwms/DpNm8504EbquqHwJU1SU079t3t2330jYYfxz4SrufFwBvSPKidhvvAd5TVYfRBK2PtNOf3/48ot3Wl2ao4SXAM4GTaN6jk9t9Hc37/mTgGe0xTfXVdj1J6pnhSNJC8TdtL8N3aU5m/wig/Uvy64D/WFXfqartwDuBM9v1fkBzknhcVf2gqr5QVdXjPn8N+B9VdXVVPVxVd1fVP6a53+l5wBur6vtVdQPwF8CrO+t+qar+pl3vezNM+w3gTVV1V1XtBC4EVmeaS+6q6pNV9Y1qrAc+SxNGenEG8Mn2OH4AXAQcDDyns8x/q6rNVfUdmpPdk6fZznTWA6vyyA3769rxJ9EEua/0uB2AK6vqy+3J+If3UMNslgW4oH3/bKdpu25bHQlsmWadLZ35u1VV11TVTW3b3kgTRle1sz8GnJzkuHb8VW39O2lO/u+oqv9ZVT+squuBj9IE7kn/q6q+2G77+3vY12rg41V1bVU9RPNHgu77vef3HE342r6HQ38msKyq3l5VD7X3hb2fR3/+jk9yZFU9UFXXzbil6b2rqu6vqjuBv+ORdj6dJnTdVVVbae4xmmp7ewyS1DPDkaSF4pfbXoYDgdcD69sT8mXAIcCG9hKv+4FPt9OhCVFfBz6b5nK0353FPo8FvjHN9OXAZBCb9E88+i/0m6ZZb+q044CPder+KrCLpofrUZK8OMl17aVX9wO/SA8n7Z16/2lypO0h2TSl3m92hncAi3vc9nrgNJqegM/TXMq0qn19YZremN2ZTQ3TLpvmKX2Tlxe+r7PMRe37ZwXwPeCfdebdSxOgp3pCZ/5uJXlWkr9rLy3bRtM7diRA+z75JI8EhjN5pDfxOOBZk++Btm1fBXSfDveo983u9kXT1j9avqp2APd1Vu/5PUfTq7RkD4d+HM0lnd36f7+zvV+j6Xn7xyT/L1MectGDmd4TjzpOpv+8LQHun2a6JM3IcCRpQamqXVV1Jc0J3fNoTly/B5xQVUe0r8Pbm++pqu1V9VtV9WTgpcB/6lz+tqcepE00lwJNtRn4ySTdE8cn0tzo/6NSpyt/mu2/uFP3EVV1UFV1t0OSA2l6Ey4Cjm5P8j8FTN4Ls6fj2ExzEju5vdAEv7tnXKN362l6sE5rh68FnksTjtbPsE6vPXez1j6lb/LywnOnmX8n8JvAe9Lelwb8LfAr7SViXafTtNHXetj1XwFXAcdW1eE091V171X6a+AVSZ5N02v3d+30TcD6Ke+BxVX177tlz2JfW4BjJhdsj/FxnXV7es+1bgSePKVXabr38O1Ttrekqn4RoKpurapX0FwyuAZY115KOOh74FHHSfN+nuqfM7ueS0kyHElaWNJ4Gc29Jl9teybeT3MPzlHtMj81ec9DkpckOb4NBN+lCVW72s19C9jdd+d8AHhtkhe0N7L/VJKnVdUm4P8Cf5jkoCQraf5C/uHdbGs67wPeMXm5VZJl7bFN9ViaHrN7gB8meTHQfeTxt4DHJTl8hv18BPil9jh+guY+rJ3tMexRe6P7adPNq6pbacLpr9I88vm7bT2/wszh6FvAMVPueZo3VXU1TWA8p530xzSXAH4gyePbNn0F8Cbgt3u8DHMJTW/i99t71V45Zf6naALq24HLOz1qnwCemuTVSX6ifT2zc9/WbPe1Dnhpkue0/75v49Ehrdf3HFV1F3ArcEpn8tTPzJeB7yZ5Y5KD2wcl/EySZ7bb/9Uky9rjnezF2UXzXn6Y3X/+ducjwG+2n8kjgDdOs8wqmoeVSFLPDEeSFoqPJ3mAJuC8AzirqibaeW+kuXTuujRPgvtbHrls6int+APAl4A/r6pr2nl/CLy5vRzogqk7rOZxwK+lOXneRnOyP9kD8wqaS7Q209xT8tb2pHs23kPTA/DZJNtpHjLxrGnq2E5zI/1HaC51emW73uT8f6TpmbitPZblU9a/hSa8/AlNT9tLaR5E8dCeCkxyDM2/3U27WWw9cF/bKzM5HmDjDMv/H2AC+GaSPV6yNkf+iOapawdW1X00vZAHAf9AcxnafwJeXVWX97i984C3t+34Fh558AAA7f09VwI/R9PzMzl9O03QPZPmvfRNmh6WA/vZV/uZOB9YS9O7sp3mYR4720V6es91/HcefX/WB4Cnt++zv6mqXTTvp5OB22neX39B8yAIaB7yMdF+dt8DnNneN7WD5nP8xXZb/3I3NUzn/TT33d1I8z77FPBD2j98tOHswfYzLEk9S29/EJMk7Y+S/CrNJYu/N+paNHtJFtP02DylmqcLznb9A2nCxwuqarqHVuwV2t7U91XVZI/YR4EPVNWnRluZpIXGcCRJ0j4kyUuBz9H03l1M0zP0jB4vD1wQ2nup/hVN79HRNPfkXVdVbxhpYZIWPC+rkyRp3/Iymkv0NtNcVnrmvhSMWqG5n2orTc/WV3nku80kqW/2HEmSJEkS9hxJkiRJEgDTfSP2XuvII4+sFStWjLoMSZIkzbP77mu+z/hxj3vcHpbU/m7Dhg33VtWyPS/54xZUOFqxYgXj4+OjLkOSJEnz7NJLLwXg7LPPHmkd2vsl+ad+1/WyOkmSJEnCcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEgCLRl3ArGzfAevHR12FJEmS+rVqbNQVSDOy50iSJEmSMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEjDkcJTkuCQbktyQZCLJubtZ9o4kRw5z/5IkSZLUr0VD3t4W4DlVtTPJYuDmJFdV1ebuQkkOGPJ+JUmSJGkgffccJVmT5LzO+IXA+VW1s510YHf7SR5I8vYkfw88u518fpLrk9yU5Gn91iJJkiRJgxrksrq1wBmd8dOBK5Icm+RGYBOwptNrdChwc1U9q6qubafdW1XPAN4LXDDdTpKck2Q8yfg927YOUK4kSZIkzazvcFRVG4GjkixPchKwtarurKpNVbUSOB44K8nR7Sq7gI9O2cyV7c8NwIoZ9nNJVY1V1diyw5f2W64kSZIk7dagD2RYB6ym6UFa253R9hhNAKe2k75fVbumrD95Cd4uhn//kyRJkiT1bNBwtBY4kyYgrUtyTJKDAZIsBZ4L3DLgPiRJkiRpzg3UW1NVE0mWAHdX1ZYkLwQuTlJAgIuq6qZhFCpJkiRJc2ngS9mq6sTO8NXAyhmWWzxlfEVneBw4bdBaJEmSJKlfQ/0SWEmSJElaqAxHkiRJkoThSJIkSZIAw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRIwhC+BnVdLDoFVY6OuQpIkSdI+yJ4jSZIkScJwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBCy0R3lv3wHrx0ddhSRJ0v7Br1DRfsaeI0mSJEnCcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJwByEoyTHJdmQ5IYkE0nO7cy7I8mRSVYkuXnY+5YkSZKkfi2ag21uAZ5TVTuTLAZuTnJVVW2eg31JkiRJ0lAM1HOUZE2S8zrjFwLnV9XOdtKBu9nHAUne3/YufTbJwYPUIkmSJEmDGPSyurXAGZ3x04Erkhyb5EZgE7Bmhl6jpwB/VlUnAPcDvzLdDpKck2Q8yfg927YOWK4kSZIkTW+gcFRVG4GjkixPchKwtarurKpNVbUSOB44K8nR06x+e1Xd0A5vAFbMsI9LqmqsqsaWHb50kHIlSZIkaUbDeCDDOmA1TQ/S2u6MtsdoAjh1mvV2doZ3MTf3P0mSJElST4YRjtYCZ9IEpHVJjpm8fyjJUuC5wC1D2I8kSZIkzZmBe2uqaiLJEuDuqtqS5IXAxUkKCHBRVd006H4kSZIkaS4N5VK2qjqxM3w1sHKG5Va0g/cCP9OZftEw6pAkSZKkfg39S2AlSZIkaSEyHEmSJEkShiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJwJC+BHbeLDkEVo2NugpJkiRJ+yB7jiRJkiQJw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRKw0B7lvX0HrB8fdRWSJEmab/dvb37uTeeCfsXMPseeI0mSJEnCcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJQB/hKMlxSTYkuSHJRJJzO/Me6HEblyZZPdt9S5IkSdJcWdTHOluA51TVziSLgZuTXFVVm3tZOUk/+5QkSZKkObXbnqMka5Kc1xm/EDi/qna2kw6cuo0kFye5Psnnkixrp12T5J1J1gO/2S76/CT/N8lt9iJJkiRJGrU9XVa3FjijM346cEWSY5PcCGwC1nR6jQ4Frq+qZwDrgbd21j2iqlZV1cXt+BOA5wEvAd41UwFJzkkynmT8nm1bez4wSZIkSZqN3YajqtoIHJVkeZKTgK1VdWdVbaqqlcDxwFlJjm5XeRi4vB2+jCb8TLqcR/ubqnq4qv4BOJoZVNUlVTVWVWPLDl86i0OTJEmSpN718kCGdcBqmh6ktd0ZbY/RBHDqDOtWZ/jBKfN2dobTQx2SJEmSNGd6CUdrgTNpAtK6JMckORggyVLgucAtne1N3j/0SuDa4ZYrSZIkSXNjj0+Oq6qJJEuAu6tqS5IXAhcnKZoen4uq6qZ28QeBE5JsALbx6PuVJEmSJGmv1dNjtavqxM7w1cDKGZZb3A7+wZTpp00ZP3uG9SRJkiRpJGb9JbCSJEmStC8yHEmSJEkShiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJQI9fArvXWHIIrBobdRWSJEmab7ff3Pz0XFBzyJ4jSZIkScJwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBCy0R3lv3wHrx0ddhSRJ0uz5CGppr2fPkSRJkiRhOJIkSZIkwHAkSZIkSYDhSJIkSZIAw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRJgOJIkSZIkwHAkSZIkScAA4SjJcUk2JLkhyUSSczvzHphhnUuTrG6H70hyZL/7lyRJkqRhWjTAuluA51TVziSLgZuTXFVVm4dUmyRJkiTNm556jpKsSXJeZ/xC4Pyq2tlOOnDqtpJcnOT6JJ9LsmyGTZ/fLnNTkqf1Ub8kSZIkDUWvl9WtBc7ojJ8OXJHk2CQ3ApuANZ1eo0OB66vqGcB64K0zbPfedpn3AhdMt0CSc5KMJxm/Z9vWHsuVJEmSpNnpKRxV1UbgqCTLk5wEbK2qO6tqU1WtBI4HzkpydLvKw8Dl7fBlwPNm2PSV7c8NwIoZ9n1JVY1V1diyw5f2Uq4kSZIkzdpsHsiwDlhN04O0tjuj7TGaAE6dYd2aYfrkZXm7GOz+J0mSJEkayGzC0VrgTJqAtC7JMUkOBkiyFHgucEtnu6vb4VcC1w6nXEmSJEmaGz331lTVRJIlwN1VtSXJC4GLkxQQ4KKquqld/EHghCQbgG08+n4lSZIkSdrrzOpStqo6sTN8NbByhuUWt4N/MGX62Z3hFZ3hceC02dQiSZIkScPU95fASpIkSdK+xHAkSZIkSRiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAGz/BLYkVtyCKwaG3UVkiRJkvZB9hxJkiRJEoYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRJgOJIkSZIkYKE9ynv7Dlg/PuoqJEmSNEx+VYv2EvYcSZIkSRKGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkC+gxHSY5LsiHJDUkmkpzbxzauSTLWz/4lSZIkadgW9bneFuA5VbUzyWLg5iRXVdXmXlZOckCf+5UkSZKkObHHnqMka5Kc1xm/EDi/qna2kw7sbifJe5OMtz1Kb+tMvyPJW5JcC7y8nfzyJF9O8rUkpw7jgCRJkiSpH71cVrcWOKMzfjpwRZJjk9wIbALWdHqN3lRVY8BKYFWSlZ11v19Vz6uqte34oqo6BXgD8Nbpdp7knDZsjd+zbessDk2SJEmSerfHcFRVG4GjkixPchKwtarurKpNVbUSOB44K8nR7SqnJ7ke2AicADy9s7nLp2z+yvbnBmDFDPu/pKrGqmps2eFLez4wSZIkSZqNXu85WgesBh5P05P0I1W1OckEcGqSDcAFwDOramuSS4GDOos/OGW7k5fm7ZpFLZIkSZI0dL0+rW4tcCZNQFqX5JgkBwMkWQo8F7gFOIwmAG1re5JePPySJUmSJGn4euqtqaqJJEuAu6tqS5IXAhcnKSDARVV1E0CSjcAEcBvwxTmqW5IkSZKGqudL2arqxM7w1TQPXJhuubNnmL5iyvhpneF7meGeI0mSJEmaD319CawkSZIk7WsMR5IkSZKE4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSMIsvgd0rLDkEVo2NugpJkiRJ+yB7jiRJkiQJw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRKw0B7lvX0HrB8fdRWSJEmaLb+ORQuAPUeSJEmShOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkoABwlGS45JsSHJDkokk5/awzjVJxtrhB/rdtyRJkiQN26IB1t0CPKeqdiZZDNyc5Kqq2jyk2iRJkiRp3vTUc5RkTZLzOuMXAudX1c520oHdbSV5b5LxtkfpbbvZ7juSfCXJdUmO7u8QJEmSJGlwvV5WtxY4ozN+OnBFkmOT3AhsAtZ0eo3eVFVjwEpgVZKV02zzUOC6qjoJ+Dzwuul2nOScNmiN37Nta4/lSpIkSdLs9BSOqmojcFSS5UlOArZW1Z1VtamqVgLHA2d1en9OT3I9sBE4AXj6NJt9CPhEO7wBWDHDvi+pqrGqGlt2+NKeD0ySJEmSZmM29xytA1YDj6fpSfqRqtqcZAI4NckG4ALgmVW1NcmlwEHTbO8HVVXt8K5Z1iJJkiRJQzWbp9WtBc6kCUjrkhyT5GCAJEuB5wK3AIcBDwLb2p6kFw+3ZEmSJEkavp57a6pqIskS4O6q2pLkhcDFSQoIcFFV3QSQZCMwAdwGfHEO6pYkSZKkoZrVpWxVdWJn+GqaBy5Mt9zZM0w/rTO8uDO8juayPUmSJEkaib6/BFaSJEmS9iWGI0mSJEnCcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJmOWXwI7ckkNg1dioq5AkSZK0D7LnSJIkSZIwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAEL7VHe23fA+vFRVyFJkiTwK1a0z7HnSJIkSZIwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSMORwlOTkJF9KMpHkxiRn7GbZa5KMDXP/kiRJktSvRUPe3g7gNVV1a5LlwIYkn6mq+7sLJTlgyPuVJEmSpIH03XOUZE2S8zrjFwIvrapbAapqM/BtYFk7/44kb0lyLfDydrWXJ/lykq8lObXfWiRJkiRpUINcVrcW6F42dzpwxeRIklOAxwLf6Czz/ap6XlWtbccXVdUpwBuAt063kyTnJBlPMn7Ptq0DlCtJkiRJM+s7HFXVRuCoJMuTnARsrao7AZI8AfgQ8Nqqeriz2uVTNnNl+3MDsGKG/VxSVWNVNbbs8KX9litJkiRJuzXoPUfrgNXA42l6kkhyGPBJ4M1Vdd2U5R+cMr6z/blrCLVIkiRJUt8GDSRrgfcDRwKrkjwW+Bjwwaq6YrdrSpIkSdJeZKBHeVfVBLAEuLuqttDcd/R84OwkN7Svk4dQpyRJkiTNqYEvZauqEzvDlwGXzbDciinjp3WG72WGe44kSZIkaT4M9UtgJUmSJGmhMhxJkiRJEoYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRJgOJIkSZIkwHAkSZIkSYDhSJIkSZIAWDTqAmZlySGwamzUVUiSJEnaB9lzJEmSJEkYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkoCF9ijv7Ttg/fioq5AkSdIg/GoW7aXsOZIkSZIkDEeSJEmSBBiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAMORJEmSJAGGI0mSJEkCDEeSJEmSBMxBOEpycpIvJZlIcmOSMzrzrkky1g4/MOx9S5IkSVK/Fs3BNncAr6mqW5MsBzYk+UxV3T8H+5IkSZKkoRio5yjJmiTndcYvBF5aVbcCVNVm4NvAshnWf0eSryS5LsnRg9QiSZIkSYMY9LK6tcAZnfHTgSsmR5KcAjwW+MY06x4KXFdVJwGfB1433Q6SnJNkPMn4Pdu2DliuJEmSJE1voHBUVRuBo5IsT3ISsLWq7gRI8gTgQ8Brq+rhaVZ/CPhEO7wBWDHDPi6pqrGqGlt2+NJBypUkSZKkGQ3jnqN1wGrg8TQ9SSQ5DPgk8Oaqum6G9X5QVdUO7xpSLZIkSZLUl2EEkrXA+4EjgVVJHgt8DPhgVV2x2zUlSZIkaS8x8KO8q2oCWALcXVVbaO47ej5wdpIb2tfJg+5HkiRJkubSUC5lq6oTO8OXAZfNsNxpneHFneF1NJfnSZIkSdJIDP1LYCVJkiRpITIcSZIkSRKGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEnAkL4Edt4sOQRWjY26CkmSJEn7IHuOJEmSJAnDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJErDQHuW9fQesHx91FZIkSZpv929vfu6v54J+nc28sOdIkiRJkjAcSZIkSRJgOJIkSZIkwHAkSZIkSYDhSJIkSZIAw5EkSZIkAYYjSZIkSQIMR5IkSZIEGI4kSZIkCTAcSZIkSRIwB+EoyaeT3J9FpY3sAAAJtklEQVTkE3tY7pokY8PevyRJkiT1Yy56jv4IePXuFkhywBzsV5IkSZL61nc4SrImyXmd8QuT/FZVfQ7YPs3ydyR5S5JrgZe3k1+e5MtJvpbk1H5rkSRJkqRBDdJztBY4ozN+OnDFHtb5flU9r6rWtuOLquoU4A3AW6dbIck5ScaTjN+zbesA5UqSJEnSzPoOR1W1ETgqyfIkJwFbq+rOPax2+ZTxK9ufG4AVM+znkqoaq6qxZYcv7bdcSZIkSdqtRQOuvw5YDTyepidpTx6cMr6z/blrCLVIkiRJUt8GDSRrgfcDRwKrBi9HkiRJkkZjoKfVVdUEsAS4u6q2ACT5As29Ry9IcleSFw1epiRJkiTNrYEvZauqE6eMT/vUuapaMWX8tM7wvcxwz5EkSZIkzYe5+J4jSZIkSVpwDEeSJEmShOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEjCEL4GdV0sOgVVjo65CkiRJ8+32m5ufngtqDtlzJEmSJEkYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkoCF9ijv7Ttg/fioq5AkSdJ8u39789Nzwb3LPvZodXuOJEmSJAnDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQBcxSOknw6yf1JPjFl+jVJxtrhB+Zi35IkSZLUj7nqOfoj4NVztG1JkiRJGrqBwlGSNUnO64xfmOS3qupzwPYe1n9Hkq8kuS7J0YPUIkmSJEmDGLTnaC1wRmf8dOCKHtc9FLiuqk4CPg+8brqFkpyTZDzJ+D3btg5UrCRJkiTNZKBwVFUbgaOSLE9yErC1qu7scfWHgMl7kjYAK2bYxyVVNVZVY8sOXzpIuZIkSZI0o0VD2MY6YDXweJqepF79oKqqHd41pFokSZIkqS/DCCRrgfcDRwKrhrA9SZIkSZp3Az+trqomgCXA3VW1BSDJF2juPXpBkruSvGjQ/UiSJEnSXBrKpWxVdeKU8VNnWO60zvDizvA6msvzJEmSJGkk5up7jiRJkiRpQTEcSZIkSRKGI0mSJEkCDEeSJEmSBBiOJEmSJAkwHEmSJEkSYDiSJEmSJMBwJEmSJEmA4UiSJEmSAFg06gJmZckhsGps1FVIkiRpvt1+c/PTc0HNIXuOJEmSJAnDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQBhiNJkiRJAgxHkiRJkgQYjiRJkiQJMBxJkiRJEmA4kiRJkiTAcCRJkiRJAKSqRl1Dz5JsB24ZdR3qy5HAvaMuQn2x7RY222/hsu0WLttu4bLtFq5u2x1XVcv62cii4dUzL26pqrFRF6HZSzJu2y1Mtt3CZvstXLbdwmXbLVy23cI1rLbzsjpJkiRJwnAkSZIkScDCC0eXjLoA9c22W7hsu4XN9lu4bLuFy7ZbuGy7hWsobbegHsggSZIkSXNlofUcSZIkSdKcMBxJkiRJEntROEryC0luSfL1JL87zfwDk1zezv/7JCs6836vnX5LkhfNZ93qv+2SvDDJhiQ3tT//9XzXvr8b5HPXzn9ikgeSXDBfNasx4O/MlUm+lGSi/fwdNJ+17+8G+J35E0n+sm2zryb5vfmuXT213/OTXJ/kh0lWT5l3VpJb29dZ81e1oP+2S3Jy53fmjUnOmN/KNcjnrp1/WJK7k/zpHndWVSN/AQcA3wCeDDwW+Arw9CnLnAe8rx0+E7i8HX56u/yBwJPa7Rww6mPaX14Dtt3PAsvb4Z8B7h718exPr0HarjP/o8AVwAWjPp796TXg524RcCNwUjv+OH9nLpi2eyWwth0+BLgDWDHqY9qfXj223wpgJfBBYHVn+k8Ct7U/l7bDS0d9TPvLa8C2eyrwlHZ4ObAFOGLUx7S/vAZpu8789wB/Bfzpnva3t/QcnQJ8vapuq6qHgLXAy6Ys8zLgL9vhdcALkqSdvraqdlbV7cDX2+1pfvTddlW1sao2t9MngIOSHDgvVQsG+9yR5Jdp/nOfmKd69YhB2u7ngRur6isAVXVfVe2ap7o1WNsVcGiSRcDBwEPAd+enbLX22H5VdUdV3Qg8PGXdFwFXV9V3qmorcDXwC/NRtIAB2q6qvlZVt7bDm4FvA8vmp2wx2OeOJP8COBr4bC8721vC0U8Bmzrjd7XTpl2mqn4IbKP5i2cv62ruDNJ2Xb8CbKyqnXNUp35c322X5FDgjcDb5qFO/bhBPndPBSrJZ9pLEH5nHurVIwZpu3XAgzR/tb4TuKiqvjPXBetRBjnn8HxltIby75/kFJrei28MqS7tWd9tl+QxwMXAb/e6s0WzKm3uZJppU58xPtMyvayruTNI2zUzkxOANTR/0db8GaTt3gb8cVU90HYkaX4N0naLgOcBzwR2AJ9LsqGqPjfcEjWDQdruFGAXzWU9S4EvJPnbqrptuCVqNwY55/B8ZbQG/vdP8gTgQ8BZVfVjPRSaM4O03XnAp6pqU6/nK3tLz9FdwLGd8WOAzTMt015ScDjwnR7X1dwZpO1IcgzwMeA1VeVfYebXIG33LODdSe4A3gD8fpLXz3XB+pFBf2eur6p7q2oH8CngGXNesSYN0navBD5dVT+oqm8DXwTG5rxidQ1yzuH5ymgN9O+f5DDgk8Cbq+q6Idem3Ruk7Z4NvL49X7kIeE2Sd+1uhb0lHP0/4ClJnpTksTQ3oF41ZZmrgMknu6wG/k81d1hdBZzZPt3nScBTgC/PU90aoO2SHEHzi+b3quqL81axJvXddlV1alWtqKoVwH8F3llVe34CjIZlkN+ZnwFWJjmkPfFeBfzDPNWtwdruTuBfp3Eo8C+Bf5ynutXopf1m8hng55MsTbKU5mqJz8xRnfpxfbddu/zHgA9W1RVzWKOm13fbVdWrquqJ7fnKBTRt+GNPu5u60l7xAn4R+BrNNZxvaqe9Hfg37fBBNE/F+jpN+HlyZ903tevdArx41Meyv736bTvgzTTXz9/QeR016uPZn16DfO4627gQn1a3oNoO+FWaB2ncDLx71Meyv70G+J25uJ0+QRNof3vUx7I/vnpov2fS/KX7QeA+YKKz7r9r2/XrwGtHfSz726vftmt/Z/5gyvnKyaM+nv3pNcjnrrONs+nhaXVpF5YkSZKk/dreclmdJEmSJI2U4UiSJEmSMBxJkiRJEmA4kiRJkiTAcCRJkiRJgOFIkiRJkgDDkSRJkiQB8P8BNMk/WF38GsUAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Testing\n",
"xvals = list(dict_roi_corr_te.keys())\n",
"yvals = list(dict_roi_corr_te.values())\n",
"myfig = plt.figure(figsize=(14,7))\n",
"plt.barh(range(len(dict_roi_corr_te)), yvals, color=\"pink\")\n",
"plt.yticks(range(len(dict_roi_corr_te)), xvals)\n",
"plt.title(\"Best correlation, within-ROI average (testing)\")\n",
"plt.axvline(x=np.mean(np.array(yvals)), color=\"gray\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"上記で作ったエンコーダには、改善する余地が山ほどある。というのも、特徴量が非常に少ないこと、学習アルゴリズムの設定が恣意的だったことなどを踏まえると、パラメータの設定をしっかりと詰めていくと大幅に改善できるはずである。注目すべき点としては、以下のものは用意に考えられる：\n",
"\n",
"- 初期値`w_init`の決め方。\n",
"\n",
"- $\\lambda$の候補の範囲と密度。\n",
"\n",
"- 学習機の反復回数の制限`t_max`。\n",
"\n",
"- フィルタバンクを定めるあらゆるパラメータ(`myparas`の中身)が大事だが、特に重要なのは`freqs`、`dir`、`sdev`であろう。\n",
"\n",
"- フィルタバンクの豊富さ。グリッドの解像度の高低、多様な空間周波数、向きなどをまんべんなく含むことが必須。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__種々の練習課題:__\n",
"\n",
"0. ここでは、全ボクセルではなく、Nishimoto et al.のいう「early visual areas」に焦点を当てている（両半球においてV1、V2、V3、V3A、V3Bのこと）。訓練データを使い、これらの全ボクセル用のエンコーダを学習させ、検証データにおいて、訓練データでもっとも良かった$\\lambda$のときの学習結果の成績を算出すること。ROIごとに、これらの成績（相関係数等）のボクセルの上での平均と分散を求めること（上ではすでにほぼ出来ている）。\n",
"\n",
"0. 先の課題を、全被験者分のデータを使って行なうこと。被験者の間、同じモデルとアルゴリズムを使って、パフォーマンスに差異があるか。被験者の上でも平均的な成績を求めて、引用している論文の成績と勝負できるか。\n",
"\n",
"0. 手法の更なる強化を測るため、時間の遅延を利用して特徴量を作ることができる。つまり、前の時点の特徴量ベクトルを、現時点の特徴量ベクトルに連結させるという方法。現時点の特徴量を$x_{i}$とすると、$x_{i-1}$が前の時点に相当するので、連結した結果が$\\widetilde{x}_{i} = (x_{i},x_{i-1})$になって、この$\\widetilde{x}_{i}$を学習に使うことになる。遅延が$k$時点なら、$\\widetilde{x}_{i} = (x_{i},x_{i-1},\\ldots,x_{i-k})$となる。数ステップを試して、パフォーマンスがもっとも良い遅延はいくらか。\n",
"\n",
"0. 学習結果の良し悪しがROIによってどの程度変わるか。特段良い・悪い成績が見られたROIがあれば、それはどれか。\n",
"\n",
"0. モデルは大事だが、学習アルゴリズムもきわめて大きな役割を果たす。ベストな設定を探るべく、反復回数などの終了条件や$\\lambda$の候補範囲、その他のアルゴリズムの工夫を調べること。どのようなアルゴリズム設定がもっとも良い成績につながったか。独自の改造を行った場合、それも併せて説明すること。\n",
"\n",
"0. これまでは同一被験者を前提として、「汎化能力」を評価してきた。当然ながら、被験者間の汎化能力も注目に値する。つまり、2名の被験者のデータを足し合わせて、大きなデータセットにまとめ、訓練データとして使って学習してから、残り1名の被験者のデータを検証データとする。同じモデルとアルゴリズムでも十分なのか。それとも、新たな工夫が必要だったか。また、このときの汎化能力の解釈が、被験者ごとに学習する場合と比べて、どう違うか。"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## 参考文献：\n",
"\n",
" - Nishimoto, Shinji, et al. \"Reconstructing visual experiences from brain activity evoked by natural movies.\" Current Biology 21.19 (2011): 1641-1646.\n",
" - Description of dataset vim-2 (visual imaging 2), at CRCNS - Collaborative Research in Computational Neuroscience. https://crcns.org/data-sets/vc/vim-2/about-vim-2"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}