A Quick Demo for the MJD Algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "acd8089b",
   "metadata": {},
   "source": [
    "#### A Quick Demo for the Mixture of Joint Distributions (MJD) Algorithm with a One-Hidden-Layer Neural Network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fa4233d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "import torch.nn.functional as F\n",
    "from torchmin import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "22e066da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the one-hidden-layer neural network model\n",
    "class NeuralNet(nn.Module):\n",
    "    def __init__(self, input_size=1000, hidden_size=100, output_size=10):\n",
    "        super(NeuralNet, self).__init__()\n",
    "        self.fc1 = nn.Linear(in_features=input_size, out_features=hidden_size, bias=True)\n",
    "        self.fc2 = nn.Linear(in_features=hidden_size, out_features=output_size, bias=True)\n",
    "    def forward(self, X):\n",
    "        FX = F.relu(self.fc1(X)) # hidden layer activation features\n",
    "        prob = F.softmax(self.fc2(FX), dim=1) # probability output\n",
    "        return FX, prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e778bd4a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# wrap the Mixture of Joint Distributions (MJD) algorithm as a class, following the sklearn style\n",
    "class MJD:\n",
    "    \"\"\"\n",
    "    In the training procedure, the total batch size per iteration is: batch_size * num_class * num_domain\n",
    "    For instance, if there are 5 domains with 68 classes in each domain, \n",
    "    then batch_size=4 means drawing 4 samples from every class in every domain,\n",
    "    resulting in 4*68*5 samples in the total batch size.\n",
    "    \"\"\"\n",
    "    def __init__(self, input_size=1024, hidden_size=512, output_size=68, seed=1000, device=torch.device('cpu'),\n",
    "                 epoch_pretrain=200, epoch=200, lamda=1, epsilon=1e-3, batch_size=4, lr=1e-3, log=False):\n",
    "        args_values = locals()\n",
    "        args_values.pop(\"self\")\n",
    "        for arg,value in args_values.items():\n",
    "            setattr(self, arg, value)\n",
    "            \n",
    "    def fit(self, X_list, y_list, Xt, yt):\n",
    "        class_labels = np.unique(y_list[0])\n",
    "        n, c = len(X_list), len(class_labels) # number of source domains, number of classes\n",
    "        # generate random target labels\n",
    "        yt_rand = np.random.choice(a=class_labels, size=len(Xt), replace=True, p=1.0 * np.ones(c) / c) \n",
    "        X_list.append(Xt), y_list.append(yt_rand)\n",
    "        \n",
    "        # define the neural network instance and the optimizer\n",
    "        torch.manual_seed(self.seed)\n",
    "        net = NeuralNet(input_size=self.input_size, hidden_size=self.hidden_size, output_size=self.output_size).to(self.device)\n",
    "        optimizer = optim.SGD(params=net.parameters(), lr=self.lr, momentum=0.9)\n",
    "        \n",
    "        #=====pretrain the network to estimate the pseudo target labels=======\n",
    "        print('Pretraining...')\n",
    "        for epoch in range(self.epoch_pretrain):\n",
    "            dataset_loaders, l = [], 0\n",
    "            for X, y in zip(X_list, y_list):\n",
    "                for i, counts in zip(*np.unique(y, return_counts=True)):\n",
    "                    dataset = np.hstack((X[y==i], y[y==i][:,None], l * np.ones((counts,1))))\n",
    "                    dataset_loaders.append(torch.utils.data.DataLoader(dataset=torch.tensor(dataset),\n",
    "                                                                       batch_size=self.batch_size, shuffle=True, drop_last=False))\n",
    "                l = l + 1 # source domain labels {0, ..., n-1}, target domain label n\n",
    "            \n",
    "            log_loss, m_log_loss = 0.0, 0.0\n",
    "            for batches in zip(*dataset_loaders):\n",
    "                Xyl = torch.cat(batches, dim=0)\n",
    "                X, y, l = Xyl[:,:-2].to(self.device,torch.float32), Xyl[:,-2].to(self.device,torch.int64), Xyl[:,-1].to(self.device,torch.int64)\n",
    "\n",
    "                FX, prob = net(X)\n",
    "                negative_log = 0.0\n",
    "                # weights for the source domains are identical in the pretraining procedure\n",
    "                for i in range(n):\n",
    "                    negative_log += -1.0 / n * torch.mean(torch.sum(torch.log(prob[l==i]) * F.one_hot(y[l==i], c), dim=1))\n",
    "                loss = negative_log\n",
    "                optimizer.zero_grad()\n",
    "                loss.backward()\n",
    "                optimizer.step()\n",
    "\n",
    "                log_loss += negative_log.item() * len(X[l!=n])\n",
    "                m_log_loss += len(X[l!=n])\n",
    "\n",
    "            with torch.no_grad():\n",
    "                Xt, yt = torch.as_tensor(Xt, dtype=torch.float32, device=self.device), torch.as_tensor(yt, dtype=torch.int64, device=self.device)    \n",
    "                yt_hat = torch.argmax(net(Xt)[1],dim=1)\n",
    "                correct = torch.sum((yt_hat == yt)).item()\n",
    "                m_test = len(yt)\n",
    "            \n",
    "            if True == self.log:\n",
    "                print('epoch ',epoch, ', log loss ',  \"{:.5f}\".format(log_loss / m_log_loss), \n",
    "                      ', test acc. ', \"{:.5f}\".format((correct / m_test) * 100)) \n",
    "        #========================================================\n",
    "    \n",
    "        #=====train the MSDA network====================================\n",
    "        print('Training...')\n",
    "        for epoch in range(self.epoch):\n",
    "            # update the pseudo target labels every epoch\n",
    "            y_list.pop()\n",
    "            y_list.append(yt_hat.cpu().numpy())\n",
    "            dataset_loaders, l = [], 0\n",
    "            for X, y in zip(X_list, y_list):\n",
    "                for i, counts in zip(*np.unique(y, return_counts=True)):\n",
    "                    dataset = np.hstack((X[y==i], y[y==i][:,None], l * np.ones((counts,1))))\n",
    "                    dataset_loaders.append(torch.utils.data.DataLoader(dataset=torch.tensor(dataset),\n",
    "                                                                       batch_size=self.batch_size, shuffle=True, drop_last=False))\n",
    "                l = l + 1 # source domain labels {0, ..., n-1}, target domain label n\n",
    "            \n",
    "            log_loss, m_log_loss = 0.0, 0.0\n",
    "            # n + 1 batches of identical size are drawn from the n + 1 source and target datasets\n",
    "            # each batch contains the same number of samples from each class\n",
    "            for batches in zip(*dataset_loaders):\n",
    "                Xyl = torch.cat(batches, dim=0)\n",
    "                X, y, l = Xyl[:,:-2].to(self.device, torch.float32), Xyl[:,-2].to(self.device, torch.int64), Xyl[:,-1].to(self.device, torch.int64)\n",
    "                \n",
    "                # compute the Gaussian kernel width\n",
    "                pairwise_dist = torch.cdist(X, X, p=2)**2 \n",
    "                sigma = torch.median(pairwise_dist[pairwise_dist!=0]) \n",
    "                \n",
    "                # compute the product kernel matrix\n",
    "                FX, prob = net(X)\n",
    "                FX_norm = torch.sum(FX ** 2, axis = -1)\n",
    "                K = torch.exp(-(FX_norm[:,None] + FX_norm[None,:] - 2 * torch.matmul(FX, FX.t())) / sigma) # feature kernel matrix     \n",
    "                Deltay = torch.as_tensor(y[:,None]==y, dtype=torch.float64, device=FX.device) # label kernel matrix  \n",
    "                P = torch.as_tensor(K * Deltay, dtype=torch.double) # product kernel matrix\n",
    "                \n",
    "                H = 1.0 / len(P[l==n]) * torch.matmul((P[l==n]).t(), P[l==n])\n",
    "                invM = torch.inverse(H + self.epsilon * torch.eye(len(P), device=FX.device))\n",
    "                \n",
    "                # optimize the mixing weights\n",
    "                def ObjAlpha(alpha):\n",
    "                    alpha = torch.softmax(alpha, dim=0)\n",
    "                    negative_log, b = 0.0, 0.0\n",
    "                    for i in range(n):\n",
    "                        negative_log += -alpha[i] * torch.mean(torch.sum(torch.log(prob[l==i]) * F.one_hot(y[l==i], c), dim=1))       \n",
    "                        b += alpha[i] * torch.mean(P[l==i], axis=0)\n",
    "                    theta = torch.matmul(invM, b)\n",
    "                    # the estimated Pearson chi2 divergence as a loss of alpha\n",
    "                    chi2 = 2 * torch.matmul(b, theta) - torch.matmul(theta, torch.matmul(H, theta)) - 1.0\n",
    "                    return negative_log + self.lamda * chi2\n",
    "                alpha0 = 1.0 * torch.ones(n, device=self.device) / n \n",
    "                result = minimize(ObjAlpha, alpha0, method='l-bfgs', max_iter=50)\n",
    "                alpha = torch.softmax(result.x, dim=0)\n",
    "                \n",
    "                negative_log = 0.0\n",
    "                for i in range(n):\n",
    "                    negative_log += -alpha[i] * torch.mean(torch.sum(torch.log(prob[l==i]) * F.one_hot(y[l==i], c), dim=1))\n",
    "                b = 0.0\n",
    "                for i in range(n):     \n",
    "                    b += alpha[i] * torch.mean(P[l==i], axis=0)\n",
    "                theta = torch.matmul(invM, b)\n",
    "                # the estimated Pearson chi2 divergence as a loss of the feature extractor\n",
    "                chi2 = 2 * torch.matmul(b, theta) - torch.matmul(theta, torch.matmul(H, theta)) - 1.0\n",
    "                loss = negative_log + self.lamda * chi2\n",
    "                optimizer.zero_grad()\n",
    "                loss.backward()\n",
    "                optimizer.step()\n",
    "\n",
    "                log_loss += negative_log.item() * len(X[l!=n])\n",
    "                m_log_loss += len(X[l!=n])\n",
    "\n",
    "            with torch.no_grad():\n",
    "                Xt, yt = torch.as_tensor(Xt, dtype=torch.float32, device=self.device), torch.as_tensor(yt, dtype=torch.int64, device=self.device)    \n",
    "                yt_hat = torch.argmax(net(Xt)[1],dim=1)\n",
    "                correct = torch.sum((yt_hat == yt)).item()\n",
    "                m_test = len(yt)\n",
    "                \n",
    "            if True == self.log:\n",
    "                print('epoch ',epoch, ', log loss ',  \"{:.5f}\".format(log_loss / m_log_loss), \n",
    "                      ', chi2 loss ', \"{:.5f}\".format(chi2.item()), \n",
    "                      ', total loss ', \"{:.5f}\".format(log_loss / m_log_loss + self.lamda * chi2.item()),\n",
    "                      ', test acc. ', \"{:.5f}\".format((correct / m_test) * 100)) \n",
    "            #========================================================     \n",
    "        self.net = net # save the network\n",
    "\n",
    "    def score(self, Xt, yt):\n",
    "        with torch.no_grad():\n",
    "            Xt, yt = torch.as_tensor(Xt, dtype=torch.float32,device=self.device), torch.as_tensor(yt, dtype=torch.int64,device=self.device)    \n",
    "            pred = torch.argmax(self.net(Xt)[1],dim=1)\n",
    "            correct = torch.sum((pred == yt)).item()\n",
    "            m_test = len(yt)\n",
    "        return (correct / m_test) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "53dd6a66",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.io as sio\n",
    "import numpy.linalg as la\n",
    "from sklearn.preprocessing import scale,LabelEncoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fb76649f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def readData(tg, domains):\n",
    "    data = sio.loadmat('PIE/' + tg + '.mat')\n",
    "    Xt, yt = data['fea'].astype(np.float64), data['gnd'].ravel()\n",
    "    yt = LabelEncoder().fit(yt).transform(yt).astype(np.float64)\n",
    "    Xt = scale(Xt / Xt.sum(axis=1,keepdims=True))\n",
    "    \n",
    "    Xs_list, ys_list = [], []\n",
    "    for sc in domains:\n",
    "        if sc != tg:\n",
    "            data = sio.loadmat('PIE/' + sc + '.mat')\n",
    "            Xs, ys = data['fea'].astype(np.float64), data['gnd'].ravel()\n",
    "            ys = LabelEncoder().fit(ys).transform(ys).astype(np.float64)\n",
    "            Xs = scale(Xs / Xs.sum(axis=1,keepdims=True))\n",
    "            Xs_list.append(Xs), ys_list.append(ys)        \n",
    "    \n",
    "    return Xs_list, ys_list, Xt, yt\n",
    "\n",
    "domains = ['C05', 'C07', 'C09', 'C27', 'C29']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c67cb105",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pretraining...\n",
      "epoch  0 , log loss  4.22435 , test acc.  3.21128\n",
      "epoch  1 , log loss  4.10045 , test acc.  9.39376\n",
      "epoch  2 , log loss  3.92304 , test acc.  16.53661\n",
      "epoch  3 , log loss  3.72881 , test acc.  22.38896\n",
      "epoch  4 , log loss  3.52655 , test acc.  26.53061\n",
      "epoch  5 , log loss  3.31548 , test acc.  31.66267\n",
      "epoch  6 , log loss  3.09318 , test acc.  37.30492\n",
      "epoch  7 , log loss  2.87622 , test acc.  42.10684\n",
      "epoch  8 , log loss  2.65817 , test acc.  47.98920\n",
      "epoch  9 , log loss  2.43428 , test acc.  52.40096\n",
      "epoch  10 , log loss  2.22172 , test acc.  56.48259\n",
      "epoch  11 , log loss  2.00899 , test acc.  59.84394\n",
      "epoch  12 , log loss  1.80526 , test acc.  63.02521\n",
      "epoch  13 , log loss  1.63748 , test acc.  64.97599\n",
      "epoch  14 , log loss  1.46537 , test acc.  66.62665\n",
      "epoch  15 , log loss  1.32250 , test acc.  67.25690\n",
      "epoch  16 , log loss  1.19497 , test acc.  68.27731\n",
      "epoch  17 , log loss  1.06854 , test acc.  69.11765\n",
      "epoch  18 , log loss  0.96872 , test acc.  69.77791\n",
      "epoch  19 , log loss  0.88122 , test acc.  70.01801\n",
      "epoch  20 , log loss  0.80998 , test acc.  70.55822\n",
      "epoch  21 , log loss  0.74276 , test acc.  71.27851\n",
      "epoch  22 , log loss  0.67740 , test acc.  71.78872\n",
      "epoch  23 , log loss  0.63180 , test acc.  72.02881\n",
      "epoch  24 , log loss  0.58780 , test acc.  72.32893\n",
      "epoch  25 , log loss  0.55118 , test acc.  72.77911\n",
      "epoch  26 , log loss  0.50810 , test acc.  73.28932\n",
      "epoch  27 , log loss  0.48264 , test acc.  73.82953\n",
      "epoch  28 , log loss  0.44990 , test acc.  74.12965\n",
      "epoch  29 , log loss  0.42320 , test acc.  74.27971\n",
      "Training...\n",
      "epoch  0 , log loss  0.40095 , chi2 loss  2.52246 , total loss  25.62558 , test acc.  73.85954\n",
      "epoch  1 , log loss  0.47141 , chi2 loss  0.67007 , total loss  7.17215 , test acc.  72.08884\n",
      "epoch  2 , log loss  0.82032 , chi2 loss  0.77990 , total loss  8.61935 , test acc.  75.51020\n",
      "epoch  3 , log loss  0.71630 , chi2 loss  0.43454 , total loss  5.06171 , test acc.  78.57143\n",
      "epoch  4 , log loss  0.77536 , chi2 loss  0.21958 , total loss  2.97113 , test acc.  79.32173\n",
      "epoch  5 , log loss  0.83904 , chi2 loss  0.11468 , total loss  1.98580 , test acc.  84.81393\n",
      "epoch  6 , log loss  0.83981 , chi2 loss  0.05879 , total loss  1.42767 , test acc.  85.92437\n",
      "epoch  7 , log loss  0.86088 , chi2 loss  0.04810 , total loss  1.34185 , test acc.  87.99520\n",
      "epoch  8 , log loss  0.77841 , chi2 loss  0.03758 , total loss  1.15419 , test acc.  88.62545\n",
      "epoch  9 , log loss  0.74202 , chi2 loss  0.02587 , total loss  1.00072 , test acc.  89.49580\n",
      "epoch  10 , log loss  0.66272 , chi2 loss  0.01920 , total loss  0.85473 , test acc.  90.21609\n",
      "epoch  11 , log loss  0.59824 , chi2 loss  0.02114 , total loss  0.80967 , test acc.  90.57623\n",
      "epoch  12 , log loss  0.57993 , chi2 loss  0.01194 , total loss  0.69934 , test acc.  91.11645\n",
      "epoch  13 , log loss  0.52893 , chi2 loss  0.01176 , total loss  0.64652 , test acc.  91.38655\n",
      "epoch  14 , log loss  0.50896 , chi2 loss  0.00914 , total loss  0.60037 , test acc.  91.65666\n",
      "epoch  15 , log loss  0.48927 , chi2 loss  0.00781 , total loss  0.56739 , test acc.  91.80672\n",
      "epoch  16 , log loss  0.46759 , chi2 loss  0.01037 , total loss  0.57131 , test acc.  91.95678\n",
      "epoch  17 , log loss  0.44162 , chi2 loss  0.00486 , total loss  0.49025 , test acc.  92.13685\n",
      "epoch  18 , log loss  0.42585 , chi2 loss  0.00754 , total loss  0.50130 , test acc.  92.40696\n",
      "epoch  19 , log loss  0.40440 , chi2 loss  0.00405 , total loss  0.44491 , test acc.  92.97719\n",
      "epoch  20 , log loss  0.38737 , chi2 loss  0.01026 , total loss  0.49000 , test acc.  92.70708\n",
      "epoch  21 , log loss  0.39049 , chi2 loss  0.00143 , total loss  0.40482 , test acc.  93.03721\n",
      "epoch  22 , log loss  0.36686 , chi2 loss  0.00529 , total loss  0.41979 , test acc.  93.09724\n",
      "epoch  23 , log loss  0.35201 , chi2 loss  0.00946 , total loss  0.44656 , test acc.  93.33733\n",
      "epoch  24 , log loss  0.34346 , chi2 loss  0.00035 , total loss  0.34697 , test acc.  93.45738\n",
      "epoch  25 , log loss  0.34186 , chi2 loss  0.00113 , total loss  0.35316 , test acc.  93.54742\n",
      "epoch  26 , log loss  0.32995 , chi2 loss  0.00211 , total loss  0.35106 , test acc.  93.54742\n",
      "epoch  27 , log loss  0.31350 , chi2 loss  0.00306 , total loss  0.34409 , test acc.  93.54742\n",
      "epoch  28 , log loss  0.31358 , chi2 loss  -0.00036 , total loss  0.31000 , test acc.  93.66747\n",
      "epoch  29 , log loss  0.29374 , chi2 loss  0.00029 , total loss  0.29665 , test acc.  93.87755\n",
      "epoch  30 , log loss  0.29039 , chi2 loss  0.00409 , total loss  0.33128 , test acc.  93.75750\n",
      "epoch  31 , log loss  0.29073 , chi2 loss  -0.00017 , total loss  0.28904 , test acc.  93.96759\n",
      "epoch  32 , log loss  0.28079 , chi2 loss  -0.00818 , total loss  0.19901 , test acc.  93.96759\n",
      "epoch  33 , log loss  0.27350 , chi2 loss  -0.00187 , total loss  0.25478 , test acc.  93.99760\n",
      "epoch  34 , log loss  0.26590 , chi2 loss  -0.00494 , total loss  0.21654 , test acc.  94.02761\n",
      "epoch  35 , log loss  0.25689 , chi2 loss  -0.00779 , total loss  0.17902 , test acc.  94.08764\n",
      "epoch  36 , log loss  0.24889 , chi2 loss  -0.00713 , total loss  0.17761 , test acc.  94.08764\n",
      "epoch  37 , log loss  0.24757 , chi2 loss  -0.00293 , total loss  0.21827 , test acc.  94.17767\n",
      "epoch  38 , log loss  0.24903 , chi2 loss  -0.00561 , total loss  0.19292 , test acc.  94.20768\n",
      "epoch  39 , log loss  0.24257 , chi2 loss  -0.00876 , total loss  0.15494 , test acc.  94.11765\n",
      "epoch  40 , log loss  0.23506 , chi2 loss  -0.00824 , total loss  0.15262 , test acc.  94.11765\n",
      "epoch  41 , log loss  0.23002 , chi2 loss  -0.00439 , total loss  0.18607 , test acc.  94.26771\n",
      "epoch  42 , log loss  0.22430 , chi2 loss  -0.01065 , total loss  0.11780 , test acc.  94.26771\n",
      "epoch  43 , log loss  0.22407 , chi2 loss  -0.00805 , total loss  0.14360 , test acc.  94.23770\n",
      "epoch  44 , log loss  0.21993 , chi2 loss  -0.01048 , total loss  0.11509 , test acc.  94.29772\n",
      "epoch  45 , log loss  0.21286 , chi2 loss  -0.00813 , total loss  0.13158 , test acc.  94.35774\n",
      "epoch  46 , log loss  0.21392 , chi2 loss  -0.01028 , total loss  0.11109 , test acc.  94.29772\n",
      "epoch  47 , log loss  0.20735 , chi2 loss  -0.00841 , total loss  0.12322 , test acc.  94.38776\n",
      "epoch  48 , log loss  0.20112 , chi2 loss  -0.00841 , total loss  0.11700 , test acc.  94.35774\n",
      "epoch  49 , log loss  0.20470 , chi2 loss  -0.01347 , total loss  0.07004 , test acc.  94.38776\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "94.38775510204081"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DEVICE = torch.device('cpu') # 'cuda:0'\n",
    "Xs_list, ys_list, Xt, yt = readData('C05', domains)\n",
    "instance = MJD(input_size=1024, hidden_size=512, output_size=68, seed=0, device=DEVICE,\n",
    "                         epoch_pretrain=30, epoch=50, lamda=10, epsilon=1e-3, batch_size=3, lr=1e-2, log=True)\n",
    "instance.fit(Xs_list, ys_list, Xt, yt)\n",
    "instance.score(Xt, yt)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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": 5
}