{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "コード3.1■ソフトマージン線形SVMによるペンギン分類モデル"
      ],
      "metadata": {
        "id": "sgG-26D0cK6q"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 725
        },
        "id": "sinf6RbkcErw",
        "outputId": "95425332-e538-4858-f320-661e49823180"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 812.85x720 with 20 Axes>"
            ],
            "image/png": "[encoded data, 348024 characters]\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "import seaborn as sns\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.svm import SVC\n",
        "\n",
        "# 図の見た目をよくする設定\n",
        "sns.set(style=\"ticks\")\n",
        "# Pandas のDataFrame としてデータ読込み\n",
        "df = sns.load_dataset(\"penguins\")\n",
        "# pairplot の設定\n",
        "sns.pairplot(df, hue='species', markers=['o', 'X', '^'],\n",
        "             palette={'Adelie': 'grey',\n",
        "                      'Chinstrap': 'black',\n",
        "                      'Gentoo': 'black'})\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# 欠損値のある行を削除\n",
        "df = df.copy().dropna()\n",
        "# 種類別に番号を付け別の列に保存\n",
        "df.loc[df[\"species\"] == 'Adelie',[\"species_num\"]] = 1\n",
        "df.loc[df[\"species\"] == 'Chinstrap',[\"species_num\"]] = 0\n",
        "df.loc[df[\"species\"] == 'Gentoo',[\"species_num\"]] = 0\n",
        "# 分類に使う特徴量を2つ決める\n",
        "x1label = \"bill_depth_mm\"\n",
        "x2label = \"bill_length_mm\"\n",
        "# 分類のラベル\n",
        "y1label = \"species\"\n",
        "y2label = \"species_num\"\n",
        "# 特徴量のベクトルをX に入れる\n",
        "X = df[[x1label, x2label]].values\n",
        "# 種別についているラベルをy に格納\n",
        "y = df[[y1label, y2label]].values\n",
        "# トレーニング用データとテスト用データに分割\n",
        "X_train, X_test, y_train, y_test = \\\n",
        "     train_test_split(X, y, test_size=0.2,\n",
        "     random_state=0, stratify=y)\n",
        "# トレーニング用データとテスト用データのデータ型を設定\n",
        "y_train_num = y_train[:,1].astype('int8')\n",
        "y_train = y_train[:,0]\n",
        "y_test_num = y_test[:,1].astype('int8')\n",
        "y_test = y_test[:,0]\n",
        "# 描画キャンパスを作成\n",
        "plt.figure(figsize=(8, 8))\n",
        "# カウンタを初期化\n",
        "cnt = 0\n",
        "cnt1 = 0\n",
        "cnt2 = 0\n",
        "cnt3 = 0\n",
        "# データの長さを取得\n",
        "datalen = len(y_train)\n",
        "# 種類ごとにマーカーや色を変えトレーニング用データを散布図にプロット\n",
        "for cnt in range(datalen):\n",
        "    if y_train[cnt] == 'Adelie':\n",
        "      if cnt1 == 0:\n",
        "        plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt],\n",
        "                    marker='o', s=50, c=\"grey\", label='Adelie')\n",
        "        cnt1 = 1\n",
        "      else:\n",
        "        plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt],\n",
        "                    marker='o', s=50, c=\"grey\")\n",
        "    elif y_train[cnt] == 'Chinstrap':\n",
        "      if cnt2 == 0:\n",
        "        plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt],\n",
        "                    marker='x', s=50, c=\"black\", label='Chinstrap')\n",
        "        cnt2 = 1\n",
        "      else:\n",
        "        plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt],\n",
        "                    marker='x', s=50, c=\"black\")\n",
        "    elif y_train[cnt] == 'Gentoo':\n",
        "      if cnt3 == 0:\n",
        "        plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt],\n",
        "                    marker='^', s=50, c=\"black\", label='Gentoo')\n",
        "        cnt3 = 1\n",
        "      else:\n",
        "        plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt],\n",
        "                    marker='^', s=50, c=\"black\")\n",
        "\n",
        "# 罫線を表示\n",
        "plt.grid()\n",
        "# 点群が表すペンギンの種類を凡例表示\n",
        "plt.legend()\n",
        "# x 軸のデータの最小値・最大値を取得\n",
        "xlim = plt.xlim()\n",
        "# y 軸のデータの最小値・最大値を取得\n",
        "ylim = plt.ylim()\n",
        "# 座標の最小値・最大値間を30等分する\n",
        "xx = np.linspace(xlim[0], xlim[1], 30)\n",
        "yy = np.linspace(ylim[0], ylim[1], 30)\n",
        "# xy平面をメッシュ状に分割\n",
        "XX, YY = np.meshgrid(xx, yy)\n",
        "# すべての格子点を格納\n",
        "xy = np.vstack([XX.ravel(), YY.ravel()]).T\n",
        "# 分類モデルを初期化\n",
        "clf = SVC(C=0.1, kernel='linear')\n",
        "# トレーニング用データから分類モデルを構築\n",
        "clf.fit(X_train,y_train_num)\n",
        "# 分類境界を表す関数を作成\n",
        "Z = clf.decision_function(xy).reshape(XX.shape)\n",
        "# Z=-1, 0,1の3本の境界線を引く\n",
        "plt.contour(XX, YY, Z, colors='k', levels=[-1, 0, 1], alpha=0.5,\n",
        "            linestyles=['--', '-', '--'])\n",
        "# サポートベクトルをプロット\n",
        "plt.scatter(clf.support_vectors_[:, 0], clf.support_vectors_[:, 1],\n",
        "            s=140, linewidth=1,\n",
        "            facecolors='none', edgecolors='k')\n",
        "# カウンタを初期化\n",
        "cnt = 0\n",
        "datalen = len(y_test)\n",
        "# テスト用データの分類結果を散布図にプロット\n",
        "for cnt in range(datalen):\n",
        "    if y_test[cnt] == 'Adelie':\n",
        "      plt.scatter(X_test[:,0][cnt],X_test[:,1][cnt],\n",
        "                  marker='+', s=50, color =\"grey\")\n",
        "    elif y_test[cnt] == 'Chinstrap':\n",
        "      plt.scatter(X_test[:,0][cnt],X_test[:,1][cnt],\n",
        "                  marker='+', s=50, color =\"black\")\n",
        "    elif y_test[cnt] == 'Gentoo':\n",
        "      plt.scatter(X_test[:,0][cnt],X_test[:,1][cnt],\n",
        "                  marker='+', s=50, color =\"black\")\n",
        "plt.xlabel(x1label, fontsize=15)\n",
        "plt.ylabel(x2label, fontsize=15)\n",
        "plt.figsize=(8,8)\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 505
        },
        "id": "HnbPuI7Xc-Yc",
        "outputId": "0f7eb043-bb23-44b0-af08-c0df04942029"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 576x576 with 1 Axes>"
            ],
            "image/png": "[encoded data, 64264 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "コード3.2■図3.30(b)の混同行列の出力"
      ],
      "metadata": {
        "id": "fRo24ZYLdI-I"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from sklearn.metrics import confusion_matrix\n",
        "\n",
        "y_test_pred = clf.decision_function(X_test)\n",
        "y_test_len = len(y_test_pred)\n",
        "\n",
        "# ラベルの個数を取得\n",
        "y_test_pred_num = np.arange(y_test_len)\n",
        "\n",
        "# 分類モデルが出力するy_test_predの値によって0,1どちらに分類するか決定\n",
        "for i in range(y_test_len):\n",
        "  if y_test_pred[i] > 0:\n",
        "    y_test_pred_num[i] = 1\n",
        "  else:\n",
        "    y_test_pred_num[i] = 0\n",
        "\n",
        "# 混同行列の出力\n",
        "cm = confusion_matrix(y_test_num, y_test_pred_num)\n",
        "cm = pd.DataFrame(data=cm, index=[\"others\", \"Adelie\"],\n",
        "                  columns=[\"others\", \"Adelie\"])\n",
        "sns.heatmap(cm, square=True, cbar=True, annot=True, cmap='Blues')\n",
        "plt.yticks(rotation=0)\n",
        "plt.xlabel(\"Predicted\", fontsize=13, rotation=0)\n",
        "plt.ylabel(\"Actual\", fontsize=13)\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 285
        },
        "id": "XC4noLiEdI4A",
        "outputId": "9ef2f834-0595-414e-c312-2864d5c7ed39"
      },
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ],
            "image/png": "[encoded data, 13404 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "図3.19（図3.31(a)）の出力"
      ],
      "metadata": {
        "id": "lZxxLKjzdIh_"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 欠損値のある行を削除\n",
        "df = df.copy().dropna()\n",
        "# 種類別に番号を付け別の列に保存\n",
        "df.loc[df[\"species\"] == 'Adelie',[\"species_num\"]] = 0\n",
        "df.loc[df[\"species\"] == 'Chinstrap',[\"species_num\"]] = 0\n",
        "df.loc[df[\"species\"] == 'Gentoo',[\"species_num\"]] = 1\n",
        "# 分類に使う特徴量を2つ決める\n",
        "x1label = \"bill_depth_mm\"\n",
        "x2label = \"flipper_length_mm\"\n",
        "# 分類のラベル\n",
        "y1label = \"species\"\n",
        "y2label = \"species_num\"\n",
        "#　特徴量のベクトルをXに入れる\n",
        "X = df[[x1label, x2label]].values\n",
        "# 種別についているラベルをyに格納\n",
        "y = df[[y1label, y2label]].values\n",
        "# トレーニング用データとテスト用データに分割\n",
        "X_train, X_test, y_train, y_test = \\\n",
        "      train_test_split(X, y, test_size=0.2, \n",
        "      random_state=0, stratify=y)\n",
        "# トレーニング用データとテスト用データのデータ型を設定\n",
        "y_train_num = y_train[:,1].astype('int8')\n",
        "y_train = y_train[:,0]\n",
        "y_test_num = y_test[:,1].astype('int8')\n",
        "y_test = y_test[:,0]\n",
        "# 描画キャンパスを作成\n",
        "plt.figure(figsize=(8, 8))\n",
        "# カウンタを初期化\n",
        "cnt = 0\n",
        "cnt1 = 0\n",
        "cnt2 = 0\n",
        "cnt3 = 0\n",
        "# データの長さを取得\n",
        "datalen = len(y_train)\n",
        "# 種類ごとにマーカーや色を変えてトレーニング用データを散布図にプロット\n",
        "for cnt in range(datalen):\n",
        "    if y_train[cnt] == 'Adelie':\n",
        "      if cnt1 == 0:\n",
        "        plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt], \n",
        "                    marker='o', s=50, c=\"grey\", label='Adelie')\n",
        "        cnt1 = 1\n",
        "      else:\n",
        "        plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt], \n",
        "                    marker='o', s=50, c=\"grey\")\n",
        "    elif y_train[cnt] == 'Chinstrap':\n",
        "      if cnt2 == 0:\n",
        "       plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt], \n",
        "                   marker='x', s=50, c=\"black\", label='Chinstrap')\n",
        "       cnt2 = 1\n",
        "      else:\n",
        "        plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt], \n",
        "                    marker='x', s=50, c=\"black\")\n",
        "    elif y_train[cnt] == 'Gentoo':\n",
        "      if cnt3 == 0:\n",
        "       plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt], \n",
        "                   marker='^', s=50, c=\"black\", label='Gentoo')\n",
        "       cnt3 = 1\n",
        "      else:\n",
        "       plt.scatter(X_train[:,0][cnt],X_train[:,1][cnt], \n",
        "                   marker='^', s=50, c=\"black\")\n",
        "\n",
        "# 罫線を表示      \n",
        "plt.grid()\n",
        "# 点群が表すペンギンの種類を凡例表示\n",
        "plt.legend()\n",
        "# x軸のデータの最小値・最大値を取得\n",
        "xlim = plt.xlim()\n",
        "# y軸のデータの最小値・最大値を取得\n",
        "ylim = plt.ylim()\n",
        "# 座標の最小値・最大値間を30等分する\n",
        "xx = np.linspace(xlim[0], xlim[1], 30)\n",
        "yy = np.linspace(ylim[0], ylim[1], 30)\n",
        "# xy平面をメッシュ状に分割\n",
        "XX, YY = np.meshgrid(xx, yy)\n",
        "# すべての格子点を格納\n",
        "xy = np.vstack([XX.ravel(), YY.ravel()]).T\n",
        "# 分類モデルを初期化\n",
        "clf = SVC(C=100000, kernel='linear')\n",
        "# トレーニング用データから分類モデルを構築\n",
        "clf.fit(X_train,y_train_num)\n",
        "# 分類境界を表す関数を作成\n",
        "Z = clf.decision_function(xy).reshape(XX.shape)\n",
        "# Z=-1, 0,1の3本の境界線を引く\n",
        "plt.contour(XX, YY, Z, colors='k', levels=[-1, 0, 1], alpha=0.5,\n",
        "                linestyles=['--', '-', '--'])\n",
        "# サポートベクトルをプロット\n",
        "plt.scatter(clf.support_vectors_[:, 0], clf.support_vectors_[:, 1], \n",
        "                s=140, linewidth=1, \n",
        "                facecolors='none', edgecolors='k')\n",
        "# カウンタを初期化\n",
        "cnt = 0\n",
        "datalen = len(y_test)\n",
        "# テスト用データの分類結果を散布図にプロット\n",
        "for cnt in range(datalen):\n",
        "    if y_test[cnt] == 'Adelie':\n",
        "       plt.scatter(X_test[:,0][cnt],X_test[:,1][cnt], \n",
        "                   marker='+', s=50, color =\"grey\")\n",
        "    elif y_test[cnt] == 'Chinstrap':\n",
        "       plt.scatter(X_test[:,0][cnt],X_test[:,1][cnt], \n",
        "                   marker='+', s=50, color =\"black\")\n",
        "    elif y_test[cnt] == 'Gentoo':\n",
        "       plt.scatter(X_test[:,0][cnt],X_test[:,1][cnt], \n",
        "                   marker='+', s=50, color =\"black\")\n",
        "plt.xlabel(x1label, fontsize=15)\n",
        "plt.ylabel(x2label, fontsize=15)\n",
        "plt.figsize=(8,8)\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 505
        },
        "id": "zqFAsCigdH16",
        "outputId": "ae092394-c0ee-45e2-ce82-5ca964983105"
      },
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 576x576 with 1 Axes>"
            ],
            "image/png": "[encoded data, 68756 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "図3.31(b)の混同行列の出力"
      ],
      "metadata": {
        "id": "qaK9hX3HeHXy"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from sklearn.metrics import confusion_matrix\n",
        "\n",
        "# 混同行列\n",
        "y_test_pred = clf.decision_function(X_test)\n",
        "y_test_len = len(y_test_pred)\n",
        "# ラベルの個数を取得\n",
        "y_test_pred_num = np.arange(y_test_len)\n",
        "# 分類モデルが出力するy_test_predの値によって0,1どちらに分類するか決定\n",
        "for i in range(y_test_len):\n",
        "  if y_test_pred[i] > 0:\n",
        "    y_test_pred_num[i] = 1\n",
        "  else:\n",
        "    y_test_pred_num[i] = 0\n",
        "# 混同行列\n",
        "cm = confusion_matrix(y_test_num, y_test_pred_num)\n",
        "cm = pd.DataFrame(data=cm, index=[\"others\", \"Gentoo\"], \n",
        "                           columns=[\"others\", \"Gentoo\"])\n",
        "sns.heatmap(cm, square=True, cbar=True, annot=True, cmap='Blues')\n",
        "plt.yticks(rotation=0)\n",
        "plt.xlabel(\"Predicted\", fontsize=13, rotation=0)\n",
        "plt.ylabel(\"Actual\", fontsize=13)\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 285
        },
        "id": "hHvcdSaPeHP-",
        "outputId": "c6c53ddd-35e8-4fc9-ced1-cf3d4a793892"
      },
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ],
            "image/png": "[encoded data, 14104 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "図3.8の出力"
      ],
      "metadata": {
        "id": "WFS5hGI7efd0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt \n",
        "import numpy as np\n",
        "\n",
        "\n",
        "def zz1(x,y):\n",
        "  #zz=  10 - 1.5*(x - 2)**2 - 1.5*(y - 2)**2\n",
        "  zz = 3 * x * y\n",
        "  return zz\n",
        "\n",
        "def zz2(x,y):\n",
        "  zz= x**2 + y**2 - 4\n",
        "  return zz\n",
        "\n",
        "  \n",
        "fig = plt.figure(figsize= (11.69,8.27), dpi=100 ) #プロット領域の作成\n",
        "xmin = -3\n",
        "xmax = 3\n",
        "ymin = -3\n",
        "ymax = 3\n",
        "ax = fig.gca(projection='3d') #プロット中の軸の取得。gca は\"Get Current Axes\" の略。\n",
        "\n",
        "x = np.arange(xmin, xmax, 0.05) # x点として[-2, 2]まで0.05刻みでサンプル\n",
        "y = np.arange(ymin, ymax, 0.05)  # y点として[-2, 2]まで0.05刻みでサンプル\n",
        "x, y = np.meshgrid(x, y)  # 上述のサンプリング点(x,y)を使ったメッシュ生成\n",
        "\n",
        "z1 = zz1(x,y)  #zz1(x,y) を計算してzz座標へ格納する。\n",
        "#z2 = zz2(x,y)\n",
        "#ax.plot_surface(x, y, z1, rstride=1, cstride=1, cmap='hsv', linewidth=0.3) # 曲面のプロット。rstrideとcstrideはステップサイズ，cmapは彩色，linewidthは曲面のメッシュの線の太さ，をそれぞれ表す。\n",
        "# ax.plot_surface(x, y, z2, rstride=1, cstride=1, cmap='hsv', linewidth=0.3)\n",
        "ax.plot_wireframe(x, y, z1, color='blue',linewidth=0.2) # ワイヤーフレームのプロット。linewidthは曲面のメッシュの線の太さ，をそれぞれ表す。\n",
        "#ax.plot_wireframe(x, y, z2, color='green',linewidth=0.2) # ワイヤーフレームのプロット。linewidthは曲面のメッシュの線の太さ，をそれぞれ表す。\n",
        "ax.view_init(elev=15, azim=310)  #elevはz軸方向から見た仰角，azimはx, y軸方向の方位角を指定\n",
        "for x1 in np.arange(xmin, xmax, 0.005):\n",
        "  for y1 in np.arange(ymin, ymax, 0.005):\n",
        "    if zz2(x1,y1)**2 < 0.00005 :\n",
        "      z3 = zz1(x1,y1)\n",
        "      ax.scatter(x1, y1, z3, c='r',s=1)\n",
        "# ax.scatter(np.sqrt(2), np.sqrt(2), 6, s=100, c=\"r\", linewidths=2, ec=\"g\")\n",
        "# ax.scatter(-np.sqrt(2), -np.sqrt(2), 6, s=100, c=\"r\", linewidths=2, ec=\"g\")\n",
        "# ax.scatter(-np.sqrt(2), np.sqrt(2), -6, s=100, c=\"r\", linewidths=2, ec=\"y\")\n",
        "# ax.scatter(np.sqrt(2), -np.sqrt(2), -6, s=100, c=\"r\", linewidths=2, ec=\"y\")\n",
        "ax.scatter(np.sqrt(2), np.sqrt(2), 6, s=100, c=\"r\",  linewidths=2)\n",
        "ax.scatter(-np.sqrt(2), -np.sqrt(2), 6, s=100, c=\"r\", linewidths=2)\n",
        "ax.scatter(-np.sqrt(2), np.sqrt(2), -6, s=100, c=\"r\", linewidths=2, marker=\"x\")\n",
        "ax.scatter(np.sqrt(2), -np.sqrt(2), -6, s=100, c=\"r\", linewidths=2, marker=\"x\")\n",
        "ax.set_xlabel('X')\n",
        "ax.set_ylabel('Y')\n",
        "ax.set_zlabel('Z')\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 661
        },
        "id": "r7Y5Y7BGefWp",
        "outputId": "fcd737a2-c433-4f29-b309-8552bf26cae9"
      },
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1169x827 with 1 Axes>"
            ],
            "image/png": "[encoded data, 225568 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "図3.9の出力"
      ],
      "metadata": {
        "id": "3egC3Y79efQD"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt \n",
        "import numpy as np\n",
        "\n",
        "def zz1(x, y):\n",
        "    zz = .5 * x - 1.5\n",
        "    return zz\n",
        "\n",
        "\n",
        "def zz2(x, y):\n",
        "    zz = -np.sqrt(x**2 + y**2)\n",
        "    return zz\n",
        "\n",
        "\n",
        "fig = plt.figure(figsize=(11.69, 8.27), dpi=100)  # プロット領域の作成\n",
        "ax = fig.gca(projection='3d')  # プロット中の軸の取得。gca は\"Get Current Axes\" の略。\n",
        "\n",
        "x = np.arange(-4, 2, 0.05)  # x点として[-2, 2]まで0.05刻みでサンプル\n",
        "y = np.arange(-4, 2, 0.05)  # y点として[-2, 2]まで0.05刻みでサンプル\n",
        "x, y = np.meshgrid(x, y)  # 上述のサンプリング点(x,y)を使ったメッシュ生成\n",
        "\n",
        "z1 = zz1(x, y)  # zz1(x,y) を計算してzz座標へ格納する。\n",
        "z2 = zz2(x, y)\n",
        "# ax.plot_surface(x, y, z1, rstride=1, cstride=1, cmap='hsv', linewidth=0.3)  # 曲面のプロット。rstrideとcstrideはステップサイズ，cmapは彩色，linewidthは曲面のメッシュの線の太さ，をそれぞれ表す。\n",
        "# ax.plot_surface(x, y, z2, rstride=1, cstride=1, cmap='hsv', linewidth=0.3)\n",
        "ax.plot_wireframe(x, y, z1, color='blue', linewidth=0.2)   # ワイヤーフレームのプロット。linewidthは曲面のメッシュの線の太さ，をそれぞれ表す。\n",
        "ax.plot_wireframe(x, y, z2, color='black', linewidth=0.2)   # ワイヤーフレームのプロット。linewidthは曲面のメッシュの線の太さ，をそれぞれ表す。\n",
        "ax.view_init(elev=10, azim=120)  # elevはz軸方向から見た仰角，azimはx, y軸方向の方位角を指定\n",
        "for x1 in np.arange(-4.0, 2.0, 0.005):\n",
        "    for y1 in np.arange(-4.0, 2.0, 0.005):\n",
        "        if (zz1(x1, y1) - zz2(x1, y1))**2 < 0.00005:\n",
        "            z3 = zz2(x1, y1)\n",
        "            ax.scatter(x1, y1, z3, c='k', s=1)\n",
        "plt.plot([0, -3], [0, 0], [0, -3], color=\"red\", linewidth=4, linestyle=\":\")\n",
        "plt.plot([0, 1], [0, 0], [0, -1], color=\"red\", linewidth=4, linestyle=\":\")\n",
        "ax.set_xlabel('X')\n",
        "ax.set_ylabel('Y')\n",
        "ax.set_zlabel('Z')\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 661
        },
        "id": "qWdEhnlrefJ9",
        "outputId": "a64d0467-52be-469d-d7ac-7253c40f7da6"
      },
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1169x827 with 1 Axes>"
            ],
            "image/png": "[encoded data, 316592 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "図3.10の出力"
      ],
      "metadata": {
        "id": "tZfKL7NeefDt"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt \n",
        "import numpy as np\n",
        "\n",
        "def zz1(x,y):\n",
        "  zz= - (x-1)**2 - (y-1)**2 + 2 \n",
        "  return zz\n",
        "\n",
        "def zz2(x,y):\n",
        "  zz= -np.sqrt(x**2 + y**2)\n",
        "  return zz\n",
        "fig = plt.figure(figsize= (11.69,8.27), dpi=100 ) #プロット領域の作成\n",
        "ax = fig.gca(projection='3d') #プロット中の軸の取得。gca は\"Get Current Axes\" の略。\n",
        "plt.tick_params(labelsize=11)\n",
        "x = np.arange(0, 2, 0.01) # x点として[-2, 2]まで0.05刻みでサンプル\n",
        "y = np.arange(0, 2, 0.01)  # y点として[-2, 2]まで0.05刻みでサンプル\n",
        "x, y = np.meshgrid(x, y)  # 上述のサンプリング点(x,y)を使ったメッシュ生成\n",
        "\n",
        "z1 = zz1(x,y)  #zz1(x,y) を計算してzz座標へ格納する。\n",
        "z2 = zz2(x,y)\n",
        "\n",
        "ax.plot_wireframe(x, y, z1, color='blue',linewidth=0.2) # ワイヤーフレームのプロット。linewidthは曲面のメッシュの線の太さ，をそれぞれ表す。\n",
        "\n",
        "ax.view_init(elev=18, azim=210)  #elevはz軸方向から見た仰角，azimはx, y軸方向の方位角を指定\n",
        "\n",
        "ax.set_xlabel('$X_1$', fontsize=13)\n",
        "ax.set_ylabel('$X_2$', fontsize=13)\n",
        "ax.set_zlabel('$X_3$', fontsize=13)\n",
        "ax.set_zticks(np.arange(0,2.1,.5))\n",
        "# plt.tight_layout()\n",
        "plt.xticks(np.arange(0,2.1,.5))\n",
        "plt.yticks(np.arange(0,2.1,.5))\n",
        "\n",
        "plt.xlim(0,2)\n",
        "plt.ylim(0,2)\n",
        "plt.subplots_adjust()\n",
        "plt.show() "
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 661
        },
        "id": "N038i01iee-Q",
        "outputId": "e312342d-f7d7-48ca-83f1-a60cab7d90ec"
      },
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1169x827 with 1 Axes>"
            ],
            "image/png": "[encoded data, 332500 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "図3.14の出力"
      ],
      "metadata": {
        "id": "pPZAjyW_ee1W"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt \n",
        "import numpy as np\n",
        "\n",
        "def zz1(x,y):\n",
        "  #zz=  10 - 1.5*(x - 2)**2 - 1.5*(y - 2)**2\n",
        "  zz = x**2 + y**2 + 6 * x * y\n",
        "  return zz\n",
        "\n",
        "def zz2(x,y):\n",
        "  zz= x**2 + y**2 - 4\n",
        "  return zz\n",
        "\n",
        "  \n",
        "fig = plt.figure(figsize= (11.69,8.27), dpi=100 ) #プロット領域の作成\n",
        "xmin = -2.5\n",
        "xmax = 2.5\n",
        "ymin = -2.5\n",
        "ymax = 2.5\n",
        "ax = fig.gca(projection='3d') #プロット中の軸の取得。gca は\"Get Current Axes\" の略。\n",
        "\n",
        "x = np.arange(xmin, xmax, 0.05) # x点として[-2, 2]まで0.05刻みでサンプル\n",
        "y = np.arange(ymin, ymax, 0.05)  # y点として[-2, 2]まで0.05刻みでサンプル\n",
        "x, y = np.meshgrid(x, y)  # 上述のサンプリング点(x,y)を使ったメッシュ生成\n",
        "\n",
        "z1 = zz1(x,y)  #zz1(x,y) を計算してzz座標へ格納する。\n",
        "z2 = zz2(x,y)\n",
        "ax.plot_wireframe(x, y, z1, color='blue',linewidth=0.2) # ワイヤーフレームのプロット。linewidthは曲面のメッシュの線の太さ，をそれぞれ表す。\n",
        "ax.view_init(elev=40, azim=330)  #elevはz軸方向から見た仰角，azimはx, y軸方向の方位角を指定\n",
        "for x1 in np.arange(xmin, xmax, 0.05):\n",
        "  for y1 in np.arange(ymin, ymax, 0.05):\n",
        "    if zz2(x1,y1) <= 0 :\n",
        "      z3 = zz1(x1,y1)\n",
        "      ax.scatter(x1, y1, z3, c='brown',s=1)\n",
        "ax.scatter(np.sqrt(2), np.sqrt(2), 16, s=120, c=\"r\",  linewidths=2, marker=\"s\", edgecolors=\"#ffffff\")\n",
        "ax.scatter(-np.sqrt(2), -np.sqrt(2), 16, s=120, c=\"r\", linewidths=2, marker=\"s\", edgecolors=\"#ffffff\")\n",
        "ax.scatter(-np.sqrt(2), np.sqrt(2), -8, s=200, c=\"r\", linewidths=2, marker=\"^\", edgecolors=\"#ffffff\")\n",
        "ax.scatter(np.sqrt(2), -np.sqrt(2), -8, s=200, c=\"r\", linewidths=2, marker=\"^\", edgecolors=\"#ffffff\")\n",
        "ax.set_xlabel('X')\n",
        "ax.set_ylabel('Y')\n",
        "ax.set_zlabel('Z')\n",
        "plt.show() "
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 661
        },
        "id": "_qXPAqbveer-",
        "outputId": "de9c80fb-353c-41bd-bf7c-368ac3b6c182"
      },
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1169x827 with 1 Axes>"
            ],
            "image/png": "[encoded data, 399688 characters]\n"
          },
          "metadata": {}
        }
      ]
    }
  ]
}