8_3_5_deeplearning_gan_sam.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "aoxI3DOK9vm2"
      },
      "source": [
        "# 8.3 SAM（Structural Agnostic Model）による因果探索の実装\n",
        "\n",
        "本ファイルは、8.3節の実装です。\n",
        "\n",
        "7.5節と同じく、「上司向け：部下とのキャリア面談のポイント研修」の疑似データを作成して、SAMによる因果探索を実施します。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "2XdIDbdlejUk"
      },
      "source": [
        "## プログラム実行前の設定など"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 122
        },
        "colab_type": "code",
        "id": "_QZagoIYv44f",
        "outputId": "1c133069-f1aa-4f79-ed78-ab14a644889b"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Looking in links: https://download.pytorch.org/whl/torch_stable.html\n",
            "\u001b[31mERROR: Could not find a version that satisfies the requirement torch==1.4.0+cu92 (from versions: 1.7.1, 1.8.0, 1.8.1, 1.9.0, 1.9.1, 1.10.0, 1.10.1, 1.10.2, 1.11.0, 1.12.0, 1.12.1, 1.13.0, 1.13.1, 2.0.0)\u001b[0m\u001b[31m\n",
            "\u001b[0m\u001b[31mERROR: No matching distribution found for torch==1.4.0+cu92\u001b[0m\u001b[31m\n",
            "\u001b[0m"
          ]
        }
      ],
      "source": [
        "# PyTorchのバージョンを下げる\n",
        "!pip install torch==1.4.0+cu92 torchvision==0.5.0+cu92 -f https://download.pytorch.org/whl/torch_stable.html"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "colab_type": "code",
        "id": "iqh9FyP-wHGa",
        "outputId": "d70b471a-b77c-47b2-af12-990ad611b00d"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "2.0.0\n"
          ]
        }
      ],
      "source": [
        "import torch \n",
        "print(torch.__version__)  # 元は1.5.0+cu101、versionを1.4に下げた"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "XZFKJwcu-_Oj"
      },
      "outputs": [],
      "source": [
        "# 乱数のシードを設定\n",
        "import random\n",
        "import numpy as np\n",
        "\n",
        "np.random.seed(1234)\n",
        "random.seed(1234)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "hx1idArc_F15"
      },
      "outputs": [],
      "source": [
        "# 使用するパッケージ（ライブラリと関数）を定義\n",
        "# 標準正規分布の生成用\n",
        "from numpy.random import *\n",
        "\n",
        "# グラフの描画用\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# その他\n",
        "import pandas as pd\n",
        "\n",
        "# シグモイド関数をimport\n",
        "from scipy.special import expit\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "AWqP6yeQlI_t"
      },
      "source": [
        "## データの作成"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "QBsAEiQ77xww"
      },
      "outputs": [],
      "source": [
        "# データ数\n",
        "num_data = 2000\n",
        "\n",
        "# 部下育成への熱心さ\n",
        "x = np.random.uniform(low=-1, high=1, size=num_data)  # -1から1の一様乱数\n",
        "\n",
        "# 上司が「上司向け：部下とのキャリア面談のポイント研修」に参加したかどうか\n",
        "e_z = randn(num_data)  # ノイズの生成\n",
        "z_prob = expit(-5.0*x+5*e_z)\n",
        "Z = np.array([])\n",
        "\n",
        "# 上司が「上司向け：部下とのキャリア面談のポイント研修」に参加したかどうか\n",
        "for i in range(num_data):\n",
        "    Z_i = np.random.choice(2, size=1, p=[1-z_prob[i], z_prob[i]])[0]\n",
        "    Z = np.append(Z, Z_i)\n",
        "\n",
        "# 介入効果の非線形性：部下育成の熱心さxの値に応じて段階的に変化\n",
        "t = np.zeros(num_data)\n",
        "for i in range(num_data):\n",
        "    if x[i] < 0:\n",
        "        t[i] = 0.5\n",
        "    elif x[i] >= 0 and x[i] < 0.5:\n",
        "        t[i] = 0.7\n",
        "    elif x[i] >= 0.5:\n",
        "        t[i] = 1.0\n",
        "\n",
        "e_y = randn(num_data)\n",
        "Y = 2.0 + t*Z + 0.3*x + 0.1*e_y \n",
        "\n",
        "\n",
        "# 本章からの追加データを生成\n",
        "\n",
        "# Y2：部下当人のチームメンバへの満足度 1から5の5段階\n",
        "Y2 = np.random.choice([1.0, 2.0, 3.0, 4.0, 5.0],\n",
        "                      num_data, p=[0.1, 0.2, 0.3, 0.2, 0.2])\n",
        "\n",
        "# Y3：部下当人の仕事への満足度\n",
        "e_y3 = randn(num_data)\n",
        "Y3 = 3*Y + Y2 + e_y3\n",
        "\n",
        "# Y4：部下当人の仕事のパフォーマンス\n",
        "e_y4 = randn(num_data)\n",
        "Y4 = 3*Y3 + 2*e_y4 + 5\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "BHcdUlW9koTa"
      },
      "source": [
        "## データをまとめた表を作成し、正規化し、可視化する"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 195
        },
        "colab_type": "code",
        "id": "1EMwdGIIIPrK",
        "outputId": "9530a268-91e5-4d68-b8ab-c2d1f62dddcd"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>x</th>\n",
              "      <th>Z</th>\n",
              "      <th>Y</th>\n",
              "      <th>Y2</th>\n",
              "      <th>Y3</th>\n",
              "      <th>Y4</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>-0.616961</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2.286924</td>\n",
              "      <td>2.0</td>\n",
              "      <td>8.732544</td>\n",
              "      <td>30.326507</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>0.244218</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2.864636</td>\n",
              "      <td>3.0</td>\n",
              "      <td>10.743959</td>\n",
              "      <td>37.149014</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>-0.124545</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.198515</td>\n",
              "      <td>3.0</td>\n",
              "      <td>10.569163</td>\n",
              "      <td>38.481185</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>0.570717</td>\n",
              "      <td>1.0</td>\n",
              "      <td>3.230572</td>\n",
              "      <td>3.0</td>\n",
              "      <td>12.312526</td>\n",
              "      <td>43.709229</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>0.559952</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.459267</td>\n",
              "      <td>5.0</td>\n",
              "      <td>12.418739</td>\n",
              "      <td>40.833938</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "          x    Z         Y   Y2         Y3         Y4\n",
              "0 -0.616961  1.0  2.286924  2.0   8.732544  30.326507\n",
              "1  0.244218  1.0  2.864636  3.0  10.743959  37.149014\n",
              "2 -0.124545  0.0  2.198515  3.0  10.569163  38.481185\n",
              "3  0.570717  1.0  3.230572  3.0  12.312526  43.709229\n",
              "4  0.559952  0.0  2.459267  5.0  12.418739  40.833938"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "df = pd.DataFrame({'x': x,\n",
        "                   'Z': Z,\n",
        "                   't': t,\n",
        "                   'Y': Y,\n",
        "                   'Y2': Y2,\n",
        "                   'Y3': Y3,\n",
        "                   'Y4': Y4,\n",
        "                   })\n",
        "\n",
        "del df[\"t\"]  # 変数tは観測できないので削除\n",
        "\n",
        "df.head()  # 先頭を表示\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "1TPIeXDg6QDG"
      },
      "source": [
        "## SAMによる推論を実施"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 386
        },
        "colab_type": "code",
        "id": "edNNPSLY6u6d",
        "outputId": "36f912b0-1dbb-4516-8ea5-d5af39cce918"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: cdt==0.5.18 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (0.5.18)\n",
            "Requirement already satisfied: scikit-learn in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (1.0.2)\n",
            "Requirement already satisfied: numpy in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (1.24.2)\n",
            "Requirement already satisfied: pandas in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (1.4.4)\n",
            "Requirement already satisfied: requests in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (2.28.1)\n",
            "Requirement already satisfied: joblib in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (1.1.0)\n",
            "Requirement already satisfied: statsmodels in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (0.13.5)\n",
            "Requirement already satisfied: skrebate in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (0.62)\n",
            "Requirement already satisfied: GPUtil in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (1.4.0)\n",
            "Requirement already satisfied: tqdm in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (4.64.1)\n",
            "Requirement already satisfied: networkx in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (2.8.4)\n",
            "Requirement already satisfied: scipy in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from cdt==0.5.18) (1.9.1)\n",
            "Requirement already satisfied: python-dateutil>=2.8.1 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from pandas->cdt==0.5.18) (2.8.2)\n",
            "Requirement already satisfied: pytz>=2020.1 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from pandas->cdt==0.5.18) (2022.1)\n",
            "Requirement already satisfied: idna<4,>=2.5 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from requests->cdt==0.5.18) (3.3)\n",
            "Requirement already satisfied: charset-normalizer<3,>=2 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from requests->cdt==0.5.18) (2.0.4)\n",
            "Requirement already satisfied: urllib3<1.27,>=1.21.1 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from requests->cdt==0.5.18) (1.26.11)\n",
            "Requirement already satisfied: certifi>=2017.4.17 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from requests->cdt==0.5.18) (2022.9.24)\n",
            "Requirement already satisfied: threadpoolctl>=2.0.0 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from scikit-learn->cdt==0.5.18) (2.2.0)\n",
            "Requirement already satisfied: patsy>=0.5.2 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from statsmodels->cdt==0.5.18) (0.5.2)\n",
            "Requirement already satisfied: packaging>=21.3 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from statsmodels->cdt==0.5.18) (21.3)\n",
            "Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from packaging>=21.3->statsmodels->cdt==0.5.18) (3.0.9)\n",
            "Requirement already satisfied: six in /Users/sugiyamahiromu/opt/anaconda3/lib/python3.9/site-packages (from patsy>=0.5.2->statsmodels->cdt==0.5.18) (1.16.0)\n"
          ]
        }
      ],
      "source": [
        "!pip install cdt==0.5.18"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "ihTvgRcv1E8s"
      },
      "source": [
        "### SAMの識別器Dの実装"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "sJQ2_9LY8MQ8"
      },
      "outputs": [],
      "source": [
        "# PyTorchから使用するものをimport\n",
        "import torch\n",
        "import torch.nn as nn\n",
        "\n",
        "\n",
        "class SAMDiscriminator(nn.Module):\n",
        "    \"\"\"SAMのDiscriminatorのニューラルネットワーク\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(self, nfeatures, dnh, hlayers):\n",
        "        super(SAMDiscriminator, self).__init__()\n",
        "\n",
        "        # ----------------------------------\n",
        "        # ネットワークの用意\n",
        "        # ----------------------------------\n",
        "        self.nfeatures = nfeatures  # 入力変数の数\n",
        "\n",
        "        layers = []\n",
        "        layers.append(nn.Linear(nfeatures, dnh))\n",
        "        layers.append(nn.BatchNorm1d(dnh))\n",
        "        layers.append(nn.LeakyReLU(.2))\n",
        "\n",
        "        for i in range(hlayers-1):\n",
        "            layers.append(nn.Linear(dnh, dnh))\n",
        "            layers.append(nn.BatchNorm1d(dnh))\n",
        "            layers.append(nn.LeakyReLU(.2))\n",
        "\n",
        "        layers.append(nn.Linear(dnh, 1))  # 最終出力\n",
        "\n",
        "        self.layers = nn.Sequential(*layers)\n",
        "\n",
        "        # ----------------------------------\n",
        "        # maskの用意（対角成分のみ1で、他は0の行列）\n",
        "        # ----------------------------------\n",
        "        mask = torch.eye(nfeatures, nfeatures)  # 変数の数×変数の数の単位行列\n",
        "        self.register_buffer(\"mask\", mask.unsqueeze(0))  # 単位行列maskを保存しておく\n",
        "\n",
        "        # 注意：register_bufferはmodelのパラメータではないが、その後forwardで使う変数を登録するPyTorchのメソッドです\n",
        "        # self.変数名で、以降も使用可能になります\n",
        "        # https://pytorch.org/docs/stable/nn.html?highlight=register_buffer#torch.nn.Module.register_buffer\n",
        "\n",
        "    def forward(self, input, obs_data=None):\n",
        "        \"\"\"　順伝搬の計算\n",
        "        Args:\n",
        "            input (torch.Size([データ数, 観測変数の種類数])): 観測したデータ、もしくは生成されたデータ\n",
        "            obs_data (torch.Size([データ数, 観測変数の種類数])):観測したデータ\n",
        "        Returns:\n",
        "            torch.Tensor: 観測したデータか、それとも生成されたデータかの判定結果\n",
        "        \"\"\"\n",
        "\n",
        "        if obs_data is not None:\n",
        "          # 生成データを識別器に入力する場合\n",
        "            return [self.layers(i) for i in torch.unbind(obs_data.unsqueeze(1) * (1 - self.mask)\n",
        "                                                         + input.unsqueeze(1) * self.mask, 1)]\n",
        "            # 対角成分のみ生成したデータ、その他は観測データに\n",
        "            # データを各変数ごとに、生成したもの、その他観測したもので混ぜて、1変数ずつ生成したものを放り込む\n",
        "            # torch.unbind(x,1)はxの1次元目でテンソルをタプルに展開する\n",
        "            # minibatch数が2000、観測データの変数が6種類の場合、\n",
        "            # [2000,6]→[2000,6,6]→([2000,6], [2000,6], [2000,6], [2000,6], [2000,6], [2000,6])→([2000,1], [2000,1], [2000,1], [2000,1], [2000,1], [2000,1])\n",
        "            # returnは[torch.Size([2000, 1]),torch.Size([2000, 1]),torch.Size([2000, 1], torch.Size([2000, 1]),torch.Size([2000, 1]),torch.Size([2000, 1])]\n",
        "\n",
        "            # 注：生成した変数全種類を用いた判定はしない。\n",
        "            # すなわち、生成した変数1種類と、元の観測データたちをまとめて1つにし、それが観測結果か、生成結果を判定させる\n",
        "\n",
        "        else:\n",
        "            # 観測データを識別器に入力する場合\n",
        "\n",
        "            return self.layers(input)\n",
        "            # returnは[torch.Size([2000, 1])]\n",
        "\n",
        "\n",
        "    def reset_parameters(self):\n",
        "        \"\"\"識別器Dの重みパラメータの初期化を実施\"\"\"\n",
        "        for layer in self.layers:\n",
        "            if hasattr(layer, 'reset_parameters'):\n",
        "                layer.reset_parameters()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "yLyjZsSc1S2i"
      },
      "source": [
        "### SAMの生成器Gの実装"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 72
        },
        "colab_type": "code",
        "id": "pBUh-fKh8X-E",
        "outputId": "b30b006a-9c56-4e93-9185-1d4dc3d2abfb"
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "No GPU automatically detected. Setting SETTINGS.GPU to 0, and SETTINGS.NJOBS to cpu_count.\n"
          ]
        }
      ],
      "source": [
        "from cdt.utils.torch import ChannelBatchNorm1d, MatrixSampler, Linear3D\n",
        "\n",
        "\n",
        "class SAMGenerator(nn.Module):\n",
        "    \"\"\"SAMのGeneratorのニューラルネットワーク\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(self, data_shape, nh):\n",
        "        \"\"\"初期化\"\"\"\n",
        "        super(SAMGenerator, self).__init__()\n",
        "\n",
        "        # ----------------------------------\n",
        "        # 対角成分のみ0で、残りは1のmaskとなる変数skeletonを作成\n",
        "        # ※最後の行は、全部1です\n",
        "        # ----------------------------------\n",
        "        nb_vars = data_shape[1]  # 変数の数\n",
        "        skeleton = 1 - torch.eye(nb_vars + 1, nb_vars)\n",
        "\n",
        "        self.register_buffer('skeleton', skeleton)\n",
        "\n",
        "        # 注意：register_bufferはmodelのパラメータではないが、その後forwardで使う変数を登録するPyTorchのメソッドです\n",
        "        # self.変数名で、以降も使用可能になります\n",
        "        # https://pytorch.org/docs/stable/nn.html?highlight=register_buffer#torch.nn.Module.register_buffer\n",
        "\n",
        "        # ----------------------------------\n",
        "        # ネットワークの用意\n",
        "        # ----------------------------------\n",
        "        # 入力層（SAMの形での全結合層）　\n",
        "        self.input_layer = Linear3D(\n",
        "            (nb_vars, nb_vars + 1, nh))  # nhは中間層のニューロン数\n",
        "        # https://github.com/FenTechSolutions/CausalDiscoveryToolbox/blob/32200779ab9b63762be3a24a2147cff09ba2bb72/cdt/utils/torch.py#L289\n",
        "\n",
        "        # 中間層\n",
        "        layers = []\n",
        "        # 2次元を1次元に変換してバッチノーマライゼーションするモジュール\n",
        "        layers.append(ChannelBatchNorm1d(nb_vars, nh))\n",
        "        layers.append(nn.Tanh())\n",
        "        self.layers = nn.Sequential(*layers)\n",
        "\n",
        "        # ChannelBatchNorm1d\n",
        "        # https://github.com/FenTechSolutions/CausalDiscoveryToolbox/blob/32200779ab9b63762be3a24a2147cff09ba2bb72/cdt/utils/torch.py#L130\n",
        "\n",
        "        # 出力層（再度、SAMの形での全結合層）\n",
        "        self.output_layer = Linear3D((nb_vars, nh, 1))\n",
        "\n",
        "    def forward(self, data, noise, adj_matrix, drawn_neurons=None):\n",
        "        \"\"\"　順伝搬の計算\n",
        "        Args:\n",
        "            data (torch.Tensor): 観測データ\n",
        "            noise (torch.Tensor): データ生成用のノイズ\n",
        "            adj_matrix (torch.Tensor): 因果関係を示す因果構造マトリクスM\n",
        "            drawn_neurons (torch.Tensor): Linear3Dの複雑さを制御する複雑さマトリクスZ\n",
        "        Returns:\n",
        "            torch.Tensor: 生成されたデータ\n",
        "        \"\"\"\n",
        "\n",
        "        # 入力層\n",
        "        x = self.input_layer(data, noise, adj_matrix *\n",
        "                             self.skeleton)  # Linear3D\n",
        "\n",
        "        # 中間層（バッチノーマライゼーションとTanh）\n",
        "        x = self.layers(x)\n",
        "\n",
        "        # 出力層\n",
        "        output = self.output_layer(\n",
        "            x, noise=None, adj_matrix=drawn_neurons)  # Linear3D\n",
        "\n",
        "        return output.squeeze(2)\n",
        "\n",
        "    def reset_parameters(self):\n",
        "        \"\"\"重みパラメータの初期化を実施\"\"\"\n",
        "\n",
        "        self.input_layer.reset_parameters()\n",
        "        self.output_layer.reset_parameters()\n",
        "\n",
        "        for layer in self.layers:\n",
        "            if hasattr(layer, 'reset_parameters'):\n",
        "                layer.reset_parameters()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "2MubteRua0mj"
      },
      "source": [
        "### SAMの誤差関数"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "Hy2GqNNdapc6"
      },
      "outputs": [],
      "source": [
        "# ネットワークを示す因果構造マトリクスMがDAG（有向非循環グラフ）になるように加える損失\n",
        "\n",
        "def notears_constr(adj_m, max_pow=None):\n",
        "    \"\"\"No Tears constraint for binary adjacency matrixes. \n",
        "    Args:\n",
        "        adj_m (array-like): Adjacency matrix of the graph\n",
        "        max_pow (int): maximum value to which the infinite sum is to be computed.\n",
        "           defaults to the shape of the adjacency_matrix\n",
        "    Returns:\n",
        "        np.ndarray or torch.Tensor: Scalar value of the loss with the type\n",
        "            depending on the input.\n",
        "    参考：https://github.com/FenTechSolutions/CausalDiscoveryToolbox/blob/32200779ab9b63762be3a24a2147cff09ba2bb72/cdt/utils/loss.py#L215\n",
        "    \"\"\"\n",
        "    m_exp = [adj_m]\n",
        "    if max_pow is None:\n",
        "        max_pow = adj_m.shape[1]\n",
        "    while(m_exp[-1].sum() > 0 and len(m_exp) < max_pow):\n",
        "        m_exp.append(m_exp[-1] @ adj_m/len(m_exp))\n",
        "\n",
        "    return sum([i.diag().sum() for idx, i in enumerate(m_exp)])\n",
        "    "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "d01nY6IKKmXe"
      },
      "source": [
        "### SAMの学習を実施する関数"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "LdgNruwmJkxj"
      },
      "outputs": [],
      "source": [
        "from sklearn.preprocessing import scale\n",
        "from torch import optim\n",
        "from torch.utils.data import DataLoader\n",
        "from tqdm import tqdm\n",
        "\n",
        "\n",
        "def run_SAM(in_data, lr_gen, lr_disc, lambda1, lambda2, hlayers, nh, dnh, train_epochs, test_epochs, device):\n",
        "    '''SAMの学習を実行する関数'''\n",
        "\n",
        "    # ---------------------------------------------------\n",
        "    # 入力データの前処理\n",
        "    # ---------------------------------------------------\n",
        "    list_nodes = list(in_data.columns)  # 入力データの列名のリスト\n",
        "    data = scale(in_data[list_nodes].values)  # 入力データの正規化\n",
        "    nb_var = len(list_nodes)  # 入力データの数 = d\n",
        "    data = data.astype('float32')  # 入力データをfloat32型に\n",
        "    data = torch.from_numpy(data).to(device)  # 入力データをPyTorchのテンソルに\n",
        "    rows, cols = data.size()  # rowsはデータ数、colsは変数の数\n",
        "\n",
        "    # ---------------------------------------------------\n",
        "    # DataLoaderの作成（バッチサイズは全データ）\n",
        "    # ---------------------------------------------------\n",
        "    batch_size = rows  # 入力データ全てを使用したミニバッチ学習とする\n",
        "    data_iterator = DataLoader(data, batch_size=batch_size,\n",
        "                               shuffle=True, drop_last=True)\n",
        "    # 注意：引数のdrop_lastはdataをbatch_sizeで取り出していったときに最後に余ったものは使用しない設定\n",
        "\n",
        "    # ---------------------------------------------------\n",
        "    # 【Generator】ネットワークの生成とパラメータの初期化\n",
        "    # cols：入力変数の数、nhは中間ニューロンの数、hlayersは中間層の数\n",
        "    # neuron_samplerは、Functional gatesの変数zを学習するネットワーク\n",
        "    # graph_samplerは、Structual gatesの変数aを学習するネットワーク\n",
        "    # ---------------------------------------------------\n",
        "    sam = SAMGenerator((batch_size, cols), nh).to(device)  # 生成器G\n",
        "    graph_sampler = MatrixSampler(nb_var, mask=None, gumbel=False).to(\n",
        "        device)  # 因果構造マトリクスMを作るネットワーク\n",
        "    neuron_sampler = MatrixSampler((nh, nb_var), mask=False, gumbel=True).to(\n",
        "        device)  # 複雑さマトリクスZを作るネットワーク\n",
        "\n",
        "    # 注意：MatrixSamplerはGumbel-Softmaxを使用し、0か1を出力させるニューラルネットワーク\n",
        "    # SAMの著者らの実装モジュール、MatrixSamplerを使用\n",
        "    # https://github.com/FenTechSolutions/CausalDiscoveryToolbox/blob/32200779ab9b63762be3a24a2147cff09ba2bb72/cdt/utils/torch.py#L212\n",
        "\n",
        "    # 重みパラメータの初期化\n",
        "    sam.reset_parameters()\n",
        "    graph_sampler.weights.data.fill_(2)\n",
        "\n",
        "    # ---------------------------------------------------\n",
        "    # 【Discriminator】ネットワークの生成とパラメータの初期化\n",
        "    # cols：入力変数の数、dnhは中間ニューロンの数、hlayersは中間層の数。\n",
        "    # ---------------------------------------------------\n",
        "    discriminator = SAMDiscriminator(cols, dnh, hlayers).to(device)\n",
        "    discriminator.reset_parameters()  # 重みパラメータの初期化\n",
        "\n",
        "    # ---------------------------------------------------\n",
        "    # 最適化の設定\n",
        "    # ---------------------------------------------------\n",
        "    # 生成器\n",
        "\n",
        "    g_optimizer = optim.Adam(sam.parameters(), lr=lr_gen)\n",
        "    graph_optimizer = optim.Adam(graph_sampler.parameters(), lr=lr_gen)\n",
        "    neuron_optimizer = optim.Adam(neuron_sampler.parameters(), lr=lr_gen)\n",
        "\n",
        "    # 識別器\n",
        "    d_optimizer = optim.Adam(discriminator.parameters(), lr=lr_disc)\n",
        "\n",
        "    # 損失関数\n",
        "    criterion = nn.BCEWithLogitsLoss()\n",
        "    # nn.BCEWithLogitsLoss()は、binary cross entropy with Logistic function\n",
        "    # https://pytorch.org/docs/stable/nn.html#bcewithlogitsloss\n",
        "\n",
        "    # 損失関数のDAGに関する制約の設定パラメータ\n",
        "    dagstart = 0.5\n",
        "    dagpenalization_increase = 0.001*10\n",
        "\n",
        "    # ---------------------------------------------------\n",
        "    # forward計算、および損失関数の計算に使用する変数を用意\n",
        "    # ---------------------------------------------------\n",
        "    _true = torch.ones(1).to(device)\n",
        "    _false = torch.zeros(1).to(device)\n",
        "\n",
        "    noise = torch.randn(batch_size, nb_var).to(device)  # 生成器Gで使用する生成ノイズ\n",
        "    noise_row = torch.ones(1, nb_var).to(device)\n",
        "\n",
        "    output = torch.zeros(nb_var, nb_var).to(device)  # 求まった隣接行列\n",
        "    output_loss = torch.zeros(1, 1).to(device)\n",
        "\n",
        "    # ---------------------------------------------------\n",
        "    # forwardの計算で、ネットワークを学習させる\n",
        "    # ---------------------------------------------------\n",
        "    pbar = tqdm(range(train_epochs + test_epochs))  # 進捗（progressive bar）の表示\n",
        "\n",
        "    for epoch in pbar:\n",
        "        for i_batch, batch in enumerate(data_iterator):\n",
        "\n",
        "            # 最適化を初期化\n",
        "            g_optimizer.zero_grad()\n",
        "            graph_optimizer.zero_grad()\n",
        "            neuron_optimizer.zero_grad()\n",
        "            d_optimizer.zero_grad()\n",
        "\n",
        "            # 因果構造マトリクスM（drawn_graph）と複雑さマトリクスZ（drawn_neurons）をMatrixSamplerから取得\n",
        "            drawn_graph = graph_sampler()\n",
        "            drawn_neurons = neuron_sampler()\n",
        "            # (drawn_graph)のサイズは、torch.Size([nb_var, nb_var])。 出力値は0か1\n",
        "            # (drawn_neurons)のサイズは、torch.Size([nh, nb_var])。 出力値は0か1\n",
        "\n",
        "            # ノイズをリセットし、生成器Gで疑似データを生成\n",
        "            noise.normal_()\n",
        "            generated_variables = sam(data=batch, noise=noise,\n",
        "                                      adj_matrix=torch.cat(\n",
        "                                          [drawn_graph, noise_row], 0),\n",
        "                                      drawn_neurons=drawn_neurons)\n",
        "\n",
        "            # 識別器Dで判定\n",
        "            # 観測変数のリスト[]で、各torch.Size([data数, 1])が求まる\n",
        "            disc_vars_d = discriminator(generated_variables.detach(), batch)\n",
        "            # 観測変数のリスト[] で、各torch.Size([data数, 1])が求まる\n",
        "            disc_vars_g = discriminator(generated_variables, batch)\n",
        "            true_vars_disc = discriminator(batch)  # torch.Size([data数, 1])が求まる\n",
        "\n",
        "            # 損失関数の計算（DCGAN）\n",
        "            disc_loss = sum([criterion(gen, _false.expand_as(gen)) for gen in disc_vars_d]) / nb_var \\\n",
        "                + criterion(true_vars_disc, _true.expand_as(true_vars_disc))\n",
        "\n",
        "            gen_loss = sum([criterion(gen,\n",
        "                                      _true.expand_as(gen))\n",
        "                            for gen in disc_vars_g])\n",
        "\n",
        "            # 損失の計算（SAM論文のオリジナルのfgan）\n",
        "            #disc_loss = sum([torch.mean(torch.exp(gen - 1)) for gen in disc_vars_d]) / nb_var - torch.mean(true_vars_disc)\n",
        "            #gen_loss = -sum([torch.mean(torch.exp(gen - 1)) for gen in disc_vars_g])\n",
        "\n",
        "            # 識別器Dのバックプロパゲーションとパラメータの更新\n",
        "            if epoch < train_epochs:\n",
        "                disc_loss.backward()\n",
        "                d_optimizer.step()\n",
        "\n",
        "            # 生成器のGの損失の計算の残り（マトリクスの複雑さとDAGのNO TEAR）\n",
        "            struc_loss = lambda1 / batch_size*drawn_graph.sum()     # Mのloss\n",
        "            func_loss = lambda2 / batch_size*drawn_neurons.sum()   # Aのloss\n",
        "\n",
        "            regul_loss = struc_loss + func_loss\n",
        "\n",
        "            if epoch <= train_epochs * dagstart:\n",
        "                # epochが基準前のときは、DAGになるようにMへのNO TEARSの制限はかけない\n",
        "                loss = gen_loss + regul_loss\n",
        "\n",
        "            else:\n",
        "                # epochが基準後のときは、DAGになるようにNO TEARSの制限をかける\n",
        "                filters = graph_sampler.get_proba()  # マトリクスMの要素を取得（ただし、0,1ではなく、1の確率）\n",
        "                dag_constraint = notears_constr(filters*filters)  # NO TERARの計算\n",
        "\n",
        "                # 徐々に線形にDAGの正則を強くする\n",
        "                loss = gen_loss + regul_loss + \\\n",
        "                    ((epoch - train_epochs * dagstart) *\n",
        "                     dagpenalization_increase) * dag_constraint\n",
        "\n",
        "            if epoch >= train_epochs:\n",
        "                # testのepochの場合、結果を取得\n",
        "                output.add_(filters.data)\n",
        "                output_loss.add_(gen_loss.data)\n",
        "            else:\n",
        "                # trainのepochの場合、生成器Gのバックプロパゲーションと更新\n",
        "                # retain_graph=Trueにすることで、以降3つのstep()が実行できる\n",
        "                loss.backward(retain_graph=True)\n",
        "                g_optimizer.step()\n",
        "                graph_optimizer.step()\n",
        "                neuron_optimizer.step()\n",
        "\n",
        "            # 進捗の表示\n",
        "            if epoch % 50 == 0:\n",
        "                pbar.set_postfix(gen=gen_loss.item()/cols,\n",
        "                                 disc=disc_loss.item(),\n",
        "                                 regul_loss=regul_loss.item(),\n",
        "                                 tot=loss.item())\n",
        "\n",
        "    return output.cpu().numpy()/test_epochs, output_loss.cpu().numpy()/test_epochs/cols  # Mと損失を出力\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "S5SXuXOCUgmg"
      },
      "source": [
        "### GPUの使用可能を確認\n",
        "\n",
        "画面上部のメニュー ランタイム > ランタイムのタイプを変更 で、 ノートブックの設定 を開く\n",
        "\n",
        "ハードウェアアクセラレータに GPU を選択し、 保存 する"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "colab_type": "code",
        "id": "ClTdYzxzXsL2",
        "outputId": "854e7cdb-d51e-4cb4-fbc6-9830415f4c44"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "False"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# GPUの使用確認：True or False\n",
        "torch.cuda.is_available()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "R-FzZ-W3Xseu"
      },
      "source": [
        "### SAMの学習を実施"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 826
        },
        "colab_type": "code",
        "id": "xfqAztolY1fo",
        "outputId": "8489d950-76d7-46fb-c951-6381c4f871d2"
      },
      "outputs": [
        {
          "ename": "AssertionError",
          "evalue": "Torch not compiled with CUDA enabled",
          "output_type": "error",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
            "\u001b[0;32m/var/folders/m2/wc_jy92n5gggmlhsz06ytgs40000gn/T/ipykernel_62156/1278054218.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m     m, loss = run_SAM(in_data=df, lr_gen=0.01*0.5,\n\u001b[0m\u001b[1;32m     10\u001b[0m                       \u001b[0mlr_disc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.01\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m                       \u001b[0;31m#lambda1=0.01, lambda2=1e-05,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/var/folders/m2/wc_jy92n5gggmlhsz06ytgs40000gn/T/ipykernel_62156/99045999.py\u001b[0m in \u001b[0;36mrun_SAM\u001b[0;34m(in_data, lr_gen, lr_disc, lambda1, lambda2, hlayers, nh, dnh, train_epochs, test_epochs, device)\u001b[0m\n\u001b[1;32m     15\u001b[0m     \u001b[0mnb_var\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlist_nodes\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# 入力データの数 = d\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m     \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'float32'\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# 入力データをfloat32型に\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m     \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_numpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdevice\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# 入力データをPyTorchのテンソルに\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     18\u001b[0m     \u001b[0mrows\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcols\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# rowsはデータ数、colsは変数の数\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/torch/cuda/__init__.py\u001b[0m in \u001b[0;36m_lazy_init\u001b[0;34m()\u001b[0m\n\u001b[1;32m    237\u001b[0m                 \"multiprocessing, you must use the 'spawn' start method\")\n\u001b[1;32m    238\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_C\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'_cuda_getDeviceCount'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 239\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mAssertionError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Torch not compiled with CUDA enabled\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    240\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0m_cudart\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    241\u001b[0m             raise AssertionError(\n",
            "\u001b[0;31mAssertionError\u001b[0m: Torch not compiled with CUDA enabled"
          ]
        }
      ],
      "source": [
        "# numpyの出力を小数点2桁に\n",
        "np.set_printoptions(precision=2, floatmode='fixed', suppress=True)\n",
        "\n",
        "# 因果探索の結果を格納するリスト\n",
        "m_list = []\n",
        "loss_list = []\n",
        "\n",
        "for i in range(5):\n",
        "    m, loss = run_SAM(in_data=df, lr_gen=0.01*0.5,\n",
        "                      lr_disc=0.01*0.5*2,\n",
        "                      #lambda1=0.01, lambda2=1e-05,\n",
        "                      lambda1=5.0*20, lambda2=0.005*20,\n",
        "                      hlayers=2,\n",
        "                      nh=200, dnh=200,\n",
        "                      train_epochs=10000,\n",
        "                      test_epochs=1000,\n",
        "                      device='cuda:0')\n",
        "\n",
        "    print(loss)\n",
        "    print(m)\n",
        "\n",
        "    m_list.append(m)\n",
        "    loss_list.append(loss)\n",
        "\n",
        "# ネットワーク構造（5回の平均）\n",
        "print(sum(m_list) / len(m_list))\n",
        "\n",
        "# mはこうなって欲しい\n",
        "#    x Z Y Y2 Y3 Y4\n",
        "# x  0 1 1 0 0 0\n",
        "# Z  0 0 1 0 0 0\n",
        "# Y  0 0 0 0 1 0\n",
        "# Y2 0 0 0 0 1 0\n",
        "# Y3 0 0 0 0 0 1\n",
        "# Y4 0 0 0 0 0 0\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "MGNG7pzi8LI6"
      },
      "source": [
        "以上"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "S9LudNsLxfkd"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "collapsed_sections": [],
      "name": "8_3_5_deeplearning_gan_sam.ipynb",
      "provenance": []
    },
    "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.9.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}