{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Building Neural ODE\n","metadata":{"id":"RFM9rZwdD_XR"}},{"cell_type":"code","source":"!pip install torchdiffeq","metadata":{"id":"CTsd_zKrEl9e","outputId":"eb296e82-fc35-4ad5-f17c-f38c6f86f198","execution":{"iopub.status.busy":"2023-08-27T16:06:31.967590Z","iopub.execute_input":"2023-08-27T16:06:31.968337Z","iopub.status.idle":"2023-08-27T16:06:46.200635Z","shell.execute_reply.started":"2023-08-27T16:06:31.968294Z","shell.execute_reply":"2023-08-27T16:06:46.199490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import math\nimport numpy as np\nimport pandas as pd\nfrom IPython.display import clear_output\nfrom tqdm import tqdm\n\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nsns.color_palette(\"bright\")\nimport matplotlib as mpl\nimport matplotlib.cm as cm\n\nimport torch\nimport torchvision\nfrom torch.utils.data import TensorDataset, DataLoader, random_split\nfrom torch import Tensor\nfrom torch import nn\nfrom torch.nn  import functional as F\nfrom torch.autograd import Variable\nimport torchdiffeq\nfrom torchdiffeq import odeint_adjoint as odeint\n\nuse_cuda = torch.cuda.is_available()\nprint(use_cuda)","metadata":{"id":"YvlFysJcD-iX","outputId":"c7c8a284-cc03-4ab8-de67-5360a5b0ea73","execution":{"iopub.status.busy":"2023-08-27T16:06:46.204328Z","iopub.execute_input":"2023-08-27T16:06:46.204676Z","iopub.status.idle":"2023-08-27T16:06:50.135230Z","shell.execute_reply.started":"2023-08-27T16:06:46.204642Z","shell.execute_reply":"2023-08-27T16:06:50.134221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def ode_solve(z0, t0, t1, f):\n    \"\"\"\n    Simplest Euler ODE initial value solver\n    \"\"\"\n    h_max = 0.05\n    n_steps = math.ceil((abs(t1 - t0)/h_max).max().item())\n\n    h = (t1 - t0)/n_steps\n    t = t0\n    z = z0\n\n    for i_step in range(n_steps):\n        z = z + h * f(z, t)\n        t = t + h\n    return z","metadata":{"id":"FoMVstEbDWW-","execution":{"iopub.status.busy":"2023-08-27T16:06:50.136640Z","iopub.execute_input":"2023-08-27T16:06:50.137564Z","iopub.status.idle":"2023-08-27T16:06:50.144887Z","shell.execute_reply.started":"2023-08-27T16:06:50.137529Z","shell.execute_reply":"2023-08-27T16:06:50.143772Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ODEF(nn.Module):\n    def forward_with_grad(self, z, t, grad_outputs):\n        \"\"\"Compute f and a df/dz, a df/dp, a df/dt\"\"\"\n        batch_size = z.shape[0]\n\n        out = self.forward(z, t)\n\n        a = grad_outputs\n        adfdz, adfdt, *adfdp = torch.autograd.grad(\n            (out,), (z, t) + tuple(self.parameters()), grad_outputs=(a),\n            allow_unused=True, retain_graph=True\n        )\n        # grad method automatically sums gradients for batch items, we have to expand them back\n        if adfdp is not None:\n            adfdp = torch.cat([p_grad.flatten() for p_grad in adfdp]).unsqueeze(0)\n            adfdp = adfdp.expand(batch_size, -1) / batch_size\n        if adfdt is not None:\n            adfdt = adfdt.expand(batch_size, 1) / batch_size\n        return out, adfdz, adfdt, adfdp\n\n    def flatten_parameters(self):\n        p_shapes = []\n        flat_parameters = []\n        for p in self.parameters():\n            p_shapes.append(p.size())\n            flat_parameters.append(p.flatten())\n        return torch.cat(flat_parameters)","metadata":{"id":"AVduwoOfFLjI","execution":{"iopub.status.busy":"2023-08-27T16:06:50.146507Z","iopub.execute_input":"2023-08-27T16:06:50.146937Z","iopub.status.idle":"2023-08-27T16:06:50.159196Z","shell.execute_reply.started":"2023-08-27T16:06:50.146912Z","shell.execute_reply":"2023-08-27T16:06:50.158444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Building the Adjoint Method\nclass ODEAdjoint(torch.autograd.Function):\n    @staticmethod\n    def forward(ctx, z0, t, flat_parameters, func):\n        assert isinstance(func, ODEF)\n        bs, *z_shape = z0.size()\n        time_len = t.size(0)\n\n        with torch.no_grad():\n            z = torch.zeros(time_len, bs, *z_shape).to(z0)\n            z[0] = z0\n            for i_t in range(time_len - 1):\n                z0 = ode_solve(z0, t[i_t], t[i_t+1], func)\n                z[i_t+1] = z0\n\n        ctx.func = func\n        ctx.save_for_backward(t, z.clone(), flat_parameters)\n        return z\n\n    @staticmethod\n    def backward(ctx, dLdz):\n        \"\"\"\n        dLdz shape: time_len, batch_size, *z_shape\n        \"\"\"\n        func = ctx.func\n        t, z, flat_parameters = ctx.saved_tensors\n        time_len, bs, *z_shape = z.size()\n        n_dim = np.prod(z_shape)\n        n_params = flat_parameters.size(0)\n\n        # Dynamics of augmented system to be calculated backwards in time\n        def augmented_dynamics(aug_z_i, t_i):\n            \"\"\"\n            tensors here are temporal slices\n            t_i - is tensor with size: bs, 1\n            aug_z_i - is tensor with size: bs, n_dim*2 + n_params + 1\n            \"\"\"\n            z_i, a = aug_z_i[:, :n_dim], aug_z_i[:, n_dim:2*n_dim]  # ignore parameters and time\n\n            # Unflatten z and a\n            z_i = z_i.view(bs, *z_shape)\n            a = a.view(bs, *z_shape)\n            with torch.set_grad_enabled(True):\n                t_i = t_i.detach().requires_grad_(True)\n                z_i = z_i.detach().requires_grad_(True)\n                func_eval, adfdz, adfdt, adfdp = func.forward_with_grad(z_i, t_i, grad_outputs=a)  # bs, *z_shape\n                adfdz = adfdz.to(z_i) if adfdz is not None else torch.zeros(bs, *z_shape).to(z_i)\n                adfdp = adfdp.to(z_i) if adfdp is not None else torch.zeros(bs, n_params).to(z_i)\n                adfdt = adfdt.to(z_i) if adfdt is not None else torch.zeros(bs, 1).to(z_i)\n\n            # Flatten f and adfdz\n            func_eval = func_eval.view(bs, n_dim)\n            adfdz = adfdz.view(bs, n_dim)\n            return torch.cat((func_eval, -adfdz, -adfdp, -adfdt), dim=1)\n\n        dLdz = dLdz.view(time_len, bs, n_dim)  # flatten dLdz for convenience\n        with torch.no_grad():\n            ## Create placeholders for output gradients\n            # Prev computed backwards adjoints to be adjusted by direct gradients\n            adj_z = torch.zeros(bs, n_dim).to(dLdz)\n            adj_p = torch.zeros(bs, n_params).to(dLdz)\n            # In contrast to z and p we need to return gradients for all times\n            adj_t = torch.zeros(time_len, bs, 1).to(dLdz)\n\n            for i_t in range(time_len-1, 0, -1):\n                z_i = z[i_t]\n                t_i = t[i_t]\n                f_i = func(z_i, t_i).view(bs, n_dim)\n\n                # Compute direct gradients\n                dLdz_i = dLdz[i_t]\n                dLdt_i = torch.bmm(torch.transpose(dLdz_i.unsqueeze(-1), 1, 2), f_i.unsqueeze(-1))[:, 0]\n\n                # Adjusting adjoints with direct gradients\n                adj_z += dLdz_i\n                adj_t[i_t] = adj_t[i_t] - dLdt_i\n\n                # Pack augmented variable\n                aug_z = torch.cat((z_i.view(bs, n_dim), adj_z, torch.zeros(bs, n_params).to(z), adj_t[i_t]), dim=-1)\n\n                # Solve augmented system backwards\n                aug_ans = ode_solve(aug_z, t_i, t[i_t-1], augmented_dynamics) #Using the Euler method to approximate the solution of augmented answer.\n\n                # Unpack solved backwards augmented system\n                adj_z[:] = aug_ans[:, n_dim:2*n_dim]\n                adj_p[:] += aug_ans[:, 2*n_dim:2*n_dim + n_params]\n                adj_t[i_t-1] = aug_ans[:, 2*n_dim + n_params:]\n\n                del aug_z, aug_ans\n\n            ## Adjust 0 time adjoint with direct gradients\n            # Compute direct gradients\n            dLdz_0 = dLdz[0]\n            dLdt_0 = torch.bmm(torch.transpose(dLdz_0.unsqueeze(-1), 1, 2), f_i.unsqueeze(-1))[:, 0]\n\n            # Adjust adjoints\n            adj_z += dLdz_0\n            adj_t[0] = adj_t[0] - dLdt_0\n        return adj_z.view(bs, *z_shape), adj_t, adj_p, None","metadata":{"id":"iLcpkPuMFL_t","execution":{"iopub.status.busy":"2023-08-27T16:06:50.163441Z","iopub.execute_input":"2023-08-27T16:06:50.164139Z","iopub.status.idle":"2023-08-27T16:06:50.185428Z","shell.execute_reply.started":"2023-08-27T16:06:50.164107Z","shell.execute_reply":"2023-08-27T16:06:50.184652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class NeuralODE(nn.Module):\n    def __init__(self, func):\n        super(NeuralODE, self).__init__()\n        assert isinstance(func, ODEF)\n        self.func = func\n\n    def forward(self, z0, t=Tensor([0., 1.]), return_whole_sequence=False):\n        t = t.to(z0)\n        z = ODEAdjoint.apply(z0, t, self.func.flatten_parameters(), self.func)\n        if return_whole_sequence:\n            return z\n        else:\n            return z[-1]","metadata":{"id":"un5-MRlHFTQW","execution":{"iopub.status.busy":"2023-08-27T16:06:50.186899Z","iopub.execute_input":"2023-08-27T16:06:50.187732Z","iopub.status.idle":"2023-08-27T16:06:50.199655Z","shell.execute_reply.started":"2023-08-27T16:06:50.187694Z","shell.execute_reply":"2023-08-27T16:06:50.198712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def norm(dim):\n    return nn.BatchNorm2d(dim)\n\ndef conv3x3(in_feats, out_feats, stride=1):\n    return nn.Conv2d(in_feats, out_feats, kernel_size=3, stride=stride, padding=1, bias=False)\n\ndef add_time(in_tensor, t):\n    bs, c, w, h = in_tensor.shape\n    return torch.cat((in_tensor, t.expand(bs, 1, w, h)), dim=1)","metadata":{"id":"2nqGqls9FYTb","execution":{"iopub.status.busy":"2023-08-27T16:06:50.201182Z","iopub.execute_input":"2023-08-27T16:06:50.201601Z","iopub.status.idle":"2023-08-27T16:06:50.210177Z","shell.execute_reply.started":"2023-08-27T16:06:50.201570Z","shell.execute_reply":"2023-08-27T16:06:50.209257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ConvODEF(ODEF):\n    def __init__(self, dim):\n        super(ConvODEF, self).__init__()\n        self.conv1 = conv3x3(dim + 1, dim)\n        self.norm1 = norm(dim)\n        self.conv2 = conv3x3(dim + 1, dim)\n        self.norm2 = norm(dim)\n\n    def forward(self, x, t):\n        xt = add_time(x, t)\n        h = self.norm1(torch.relu(self.conv1(xt)))\n        ht = add_time(h, t)\n        dxdt = self.norm2(torch.relu(self.conv2(ht)))\n        return dxdt\n\nclass ContinuousNeuralMNISTClassifier(nn.Module):\n    def __init__(self, ode):\n        super(ContinuousNeuralMNISTClassifier, self).__init__()\n        self.downsampling = nn.Sequential(\n            nn.Conv2d(1, 64, 3, 1),\n            norm(64),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(64, 64, 4, 2, 1),\n            norm(64),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(64, 64, 4, 2, 1),\n        )\n        self.feature = ode\n        self.norm = norm(64)\n        self.avg_pool = nn.AdaptiveAvgPool2d((1, 1))\n        self.fc = nn.Linear(64, 10)\n\n    def forward(self, x):\n        x = self.downsampling(x)\n        x = self.feature(x)\n        x = self.norm(x)\n        x = self.avg_pool(x)\n        shape = torch.prod(torch.tensor(x.shape[1:])).item()\n        x = x.view(-1, shape)\n        out = self.fc(x)\n        return out","metadata":{"id":"yow6zYVVFZ-U","execution":{"iopub.status.busy":"2023-08-27T16:06:50.211166Z","iopub.execute_input":"2023-08-27T16:06:50.211475Z","iopub.status.idle":"2023-08-27T16:06:50.226511Z","shell.execute_reply.started":"2023-08-27T16:06:50.211452Z","shell.execute_reply":"2023-08-27T16:06:50.225627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"func = ConvODEF(64)\node = NeuralODE(func)\nmodel = ContinuousNeuralMNISTClassifier(ode)\nif use_cuda:\n    model = model.cuda()","metadata":{"id":"NWHWxwtBFd3Y","execution":{"iopub.status.busy":"2023-08-27T16:06:50.227931Z","iopub.execute_input":"2023-08-27T16:06:50.228309Z","iopub.status.idle":"2023-08-27T16:06:53.068134Z","shell.execute_reply.started":"2023-08-27T16:06:50.228278Z","shell.execute_reply":"2023-08-27T16:06:53.067134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img_std = 0.3081\nimg_mean = 0.1307\n\n\nbatch_size = 32\ntrain_loader = torch.utils.data.DataLoader(\n    torchvision.datasets.MNIST(\"data/mnist\", train=True, download=True,\n                             transform=torchvision.transforms.Compose([\n                                 torchvision.transforms.ToTensor(),\n                                 torchvision.transforms.Normalize((img_mean,), (img_std,))\n                             ])\n    ),\n    batch_size=batch_size, shuffle=True\n)\n\ntest_loader = torch.utils.data.DataLoader(\n    torchvision.datasets.MNIST(\"data/mnist\", train=False, download=True,\n                             transform=torchvision.transforms.Compose([\n                                 torchvision.transforms.ToTensor(),\n                                 torchvision.transforms.Normalize((img_mean,), (img_std,))\n                             ])\n    ),\n    batch_size=128, shuffle=True\n)","metadata":{"id":"DtOzdh7BFf0V","execution":{"iopub.status.busy":"2023-08-27T16:06:53.069730Z","iopub.execute_input":"2023-08-27T16:06:53.070077Z","iopub.status.idle":"2023-08-27T16:06:54.150955Z","shell.execute_reply.started":"2023-08-27T16:06:53.070045Z","shell.execute_reply":"2023-08-27T16:06:54.149995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"optimizer = torch.optim.Adam(model.parameters())","metadata":{"id":"fvICbZxvFl_I","execution":{"iopub.status.busy":"2023-08-27T16:06:54.152534Z","iopub.execute_input":"2023-08-27T16:06:54.152908Z","iopub.status.idle":"2023-08-27T16:06:54.157885Z","shell.execute_reply.started":"2023-08-27T16:06:54.152875Z","shell.execute_reply":"2023-08-27T16:06:54.156916Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train(epoch):\n    num_items = 0\n    train_losses = []\n\n    model.train()\n    criterion = nn.CrossEntropyLoss()\n    print(f\"Training Epoch {epoch}...\")\n    for batch_idx, (data, target) in tqdm(enumerate(train_loader), total=len(train_loader)):\n        if use_cuda:\n            data = data.cuda()\n            target = target.cuda()\n        optimizer.zero_grad()\n        output = model(data)\n        loss = criterion(output, target)\n        loss.backward()\n        optimizer.step()\n\n        train_losses += [loss.item()]\n        num_items += data.shape[0]\n    print('Train loss: {:.5f}'.format(np.mean(train_losses)))\n    return train_losses\n\n\ndef test():\n    accuracy = 0.0\n    num_items = 0\n\n    model.eval()\n    criterion = nn.CrossEntropyLoss()\n    print(f\"Testing...\")\n    with torch.no_grad():\n        for batch_idx, (data, target) in tqdm(enumerate(test_loader),  total=len(test_loader)):\n            if use_cuda:\n                data = data.cuda()\n                target = target.cuda()\n            output = model(data)\n            accuracy += torch.sum(torch.argmax(output, dim=1) == target).item()\n            num_items += data.shape[0]\n    accuracy = accuracy * 100 / num_items\n    print(\"Test Accuracy: {:.3f}%\".format(accuracy))","metadata":{"id":"QyQ84lrfFn_q","execution":{"iopub.status.busy":"2023-08-27T16:06:54.159764Z","iopub.execute_input":"2023-08-27T16:06:54.160487Z","iopub.status.idle":"2023-08-27T16:06:54.172034Z","shell.execute_reply.started":"2023-08-27T16:06:54.160454Z","shell.execute_reply":"2023-08-27T16:06:54.171162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_epochs = 5\ntest()\ntrain_losses = []\nfor epoch in range(1, n_epochs + 1):\n    train_losses += train(epoch)\n    test()","metadata":{"id":"H1lwj2NOFq4K","outputId":"7095874e-b1db-4f75-82ea-00d88efadbd2","execution":{"iopub.status.busy":"2023-08-27T16:06:54.175156Z","iopub.execute_input":"2023-08-27T16:06:54.175472Z","iopub.status.idle":"2023-08-27T16:15:06.049602Z","shell.execute_reply.started":"2023-08-27T16:06:54.175449Z","shell.execute_reply":"2023-08-27T16:15:06.048645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\n\nplt.figure(figsize=(9, 5))\nhistory = pd.DataFrame({\"loss\": train_losses})\nhistory[\"cum_data\"] = history.index * batch_size\nhistory[\"smooth_loss\"] = history.loss.ewm(halflife=10).mean()\nhistory.plot(x=\"cum_data\", y=\"smooth_loss\", figsize=(12, 5), title=\"train error\")","metadata":{"id":"8tNWBoeRFsAh","execution":{"iopub.status.busy":"2023-08-27T16:15:06.055163Z","iopub.execute_input":"2023-08-27T16:15:06.056050Z","iopub.status.idle":"2023-08-27T16:15:06.466998Z","shell.execute_reply.started":"2023-08-27T16:15:06.056022Z","shell.execute_reply":"2023-08-27T16:15:06.466089Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Building the Convolutional Neural Network","metadata":{"id":"FPKPoD10F0-k"}},{"cell_type":"code","source":"def norm(dim):\n    return nn.GroupNorm(min(32, dim), dim)\nclass Flatten(nn.Module):\n\n    def __init__(self):\n        super(Flatten, self).__init__()\n\n    def forward(self, x):\n        shape = torch.prod(torch.tensor(x.shape[1:])).item()\n        return x.view(-1, shape)\n\nclass ResBlock(nn.Module):\n    def __init__(self, in_planes, planes, stride=1):\n        super(ResBlock, self).__init__()\n        self.relu = nn.ReLU(inplace=True)\n        self.conv1 = nn.Conv2d(in_planes, planes, kernel_size=3, stride=stride, padding=1, bias=False)\n        self.norm1 = norm(in_planes)\n        self.conv2 = nn.Conv2d(in_planes, planes, kernel_size=3, stride=stride, padding=1, bias=False)\n        self.norm2 = norm(planes)\n\n    def forward(self, x):\n        out = self.conv1(x)\n        out = self.norm1(out)\n        out = self.relu(out)\n        out = self.conv2(out)\n        out = self.norm2(out)\n        out = self.relu(out)\n        return out + x\n\nclass CNN_MNIST(nn.Module):\n    def __init__(self):\n        super(CNN_MNIST, self).__init__()\n        downsampling_layers = [\n            nn.Conv2d(1, 64, 3, 1),\n            norm(64),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(64, 64, 4, 2, 1),\n            norm(64),\n            nn.ReLU(inplace=True)\n            ]\n        feature_layers = [ResBlock(64, 64) for _ in range(1)]\n        fc_layers = [nn.AdaptiveAvgPool2d((1, 1)), \n            Flatten(), nn.Linear(64, 10)\n            ]\n        self.net = nn.Sequential(*downsampling_layers, *feature_layers, *fc_layers)\n\n    def forward(self, x):\n        return self.net(x)","metadata":{"id":"bInj3F13F0k-","execution":{"iopub.status.busy":"2023-08-27T16:15:06.468453Z","iopub.execute_input":"2023-08-27T16:15:06.469459Z","iopub.status.idle":"2023-08-27T16:15:06.483942Z","shell.execute_reply.started":"2023-08-27T16:15:06.469423Z","shell.execute_reply":"2023-08-27T16:15:06.482945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow.keras.utils import to_categorical\n\n# MNIST dataset\n(X_train, y_train), (X_test, y_test) = tf.keras.datasets.mnist.load_data()","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:15:06.485436Z","iopub.execute_input":"2023-08-27T16:15:06.485783Z","iopub.status.idle":"2023-08-27T16:15:14.690895Z","shell.execute_reply.started":"2023-08-27T16:15:06.485750Z","shell.execute_reply":"2023-08-27T16:15:14.689923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"BATCH_SIZE = 32\n\ntorch_X_train = torch.from_numpy(X_train).type(torch.LongTensor)\ntorch_y_train = torch.from_numpy(y_train).type(torch.LongTensor) # data type is long\n\n# create feature and targets tensor for test set.\ntorch_X_test = torch.from_numpy(X_test).type(torch.LongTensor)\ntorch_y_test = torch.from_numpy(y_test).type(torch.LongTensor) # data type is long\n\n# Pytorch train and test sets\ntrain = torch.utils.data.TensorDataset(torch_X_train,torch_y_train)\ntest = torch.utils.data.TensorDataset(torch_X_test,torch_y_test)\n\n# data loader\ntrain_loaderCNN = torch.utils.data.DataLoader(train, batch_size = BATCH_SIZE, shuffle = False)\ntest_loaderCNN = torch.utils.data.DataLoader(test, batch_size = BATCH_SIZE, shuffle = False)\n","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:15:14.692606Z","iopub.execute_input":"2023-08-27T16:15:14.693612Z","iopub.status.idle":"2023-08-27T16:15:14.962529Z","shell.execute_reply.started":"2023-08-27T16:15:14.693575Z","shell.execute_reply":"2023-08-27T16:15:14.961512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(y_test.shape)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:15:14.964087Z","iopub.execute_input":"2023-08-27T16:15:14.964473Z","iopub.status.idle":"2023-08-27T16:15:14.970968Z","shell.execute_reply.started":"2023-08-27T16:15:14.964440Z","shell.execute_reply":"2023-08-27T16:15:14.969799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def print_num_parameters(model, name=\"model\"):\n  print(\"Number of parameters in %s\" % name,\n        sum(p.numel() for p in model.parameters()))\n\n\nprint_num_parameters(CNN_MNIST(), \"The CNN\")","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:15:14.972329Z","iopub.execute_input":"2023-08-27T16:15:14.973266Z","iopub.status.idle":"2023-08-27T16:15:14.988111Z","shell.execute_reply.started":"2023-08-27T16:15:14.973233Z","shell.execute_reply":"2023-08-27T16:15:14.987052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cnn = CNN_MNIST()\nprint(cnn)\n\nit = iter(train_loader)\nX_batch, y_batch = next(it)\nprint(cnn.forward(X_batch).shape)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:15:14.991247Z","iopub.execute_input":"2023-08-27T16:15:14.991502Z","iopub.status.idle":"2023-08-27T16:15:15.130055Z","shell.execute_reply.started":"2023-08-27T16:15:14.991480Z","shell.execute_reply":"2023-08-27T16:15:15.129068Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_lossCNN = []\ndef fit(model, train_loader):\n    optimizer = torch.optim.Adam(model.parameters())#,lr=0.001, betas=(0.9,0.999))\n    error = nn.CrossEntropyLoss()\n    EPOCHS = 5\n    model.train()\n    for epoch in range(EPOCHS):\n        correct = 0\n        for batch_idx, (X_batch, y_batch) in enumerate(train_loader):\n            var_X_batch = Variable(X_batch).float()\n            var_y_batch = Variable(y_batch)\n            optimizer.zero_grad()\n            output = model(var_X_batch)\n            loss = error(output, var_y_batch)\n            loss.backward()\n            train_lossCNN.append(loss.data)\n            optimizer.step()\n\n            # Total correct predictions\n            predicted = torch.max(output.data, 1)[1] \n            correct += (predicted == var_y_batch).sum()\n            #print(correct)\n            if batch_idx % 50 == 0:\n                print('Epoch : {} [{}/{} ({:.0f}%)]\\tLoss: {:.6f}\\t Accuracy:{:.3f}%'.format(\n                    epoch, batch_idx*len(X_batch), len(train_loader.dataset), 100.*batch_idx / len(train_loader), loss.data, float(correct*100) / float(BATCH_SIZE*(batch_idx+1))))","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:15:15.131690Z","iopub.execute_input":"2023-08-27T16:15:15.132303Z","iopub.status.idle":"2023-08-27T16:15:15.143816Z","shell.execute_reply.started":"2023-08-27T16:15:15.132269Z","shell.execute_reply":"2023-08-27T16:15:15.142850Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def evaluate(model):\n#model = mlp\n    correct = 0 \n    for test_imgs, test_labels in test_loader:\n        #print(test_imgs.shape)\n        test_imgs = Variable(test_imgs).float()\n        output = model(test_imgs)\n        predicted = torch.max(output,1)[1]\n        correct += (predicted == test_labels).sum()\n    print(\"Test accuracy:{:.3f}% \".format( float(correct) / (len(test_loader)*batch_size)))","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:15:15.144768Z","iopub.execute_input":"2023-08-27T16:15:15.145048Z","iopub.status.idle":"2023-08-27T16:15:15.166006Z","shell.execute_reply.started":"2023-08-27T16:15:15.145025Z","shell.execute_reply":"2023-08-27T16:15:15.165020Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"torch_X_train = torch_X_train.view(-1, 1,28,28).float()\ntorch_X_test = torch_X_test.view(-1,1,28,28).float()\nprint(torch_X_train.shape)\nprint(torch_X_test.shape)\n\n# Pytorch train and test sets\ntrain = torch.utils.data.TensorDataset(torch_X_train,torch_y_train)\ntest = torch.utils.data.TensorDataset(torch_X_test,torch_y_test)\n\n# data loader\ntrain_loaderCNN = torch.utils.data.DataLoader(train, batch_size = BATCH_SIZE, shuffle = False)\ntest_loaderCNN = torch.utils.data.DataLoader(test, batch_size = BATCH_SIZE, shuffle = False)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:15:15.167259Z","iopub.execute_input":"2023-08-27T16:15:15.167809Z","iopub.status.idle":"2023-08-27T16:15:15.380403Z","shell.execute_reply.started":"2023-08-27T16:15:15.167777Z","shell.execute_reply":"2023-08-27T16:15:15.379261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fit(cnn, train_loaderCNN)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:15:15.382307Z","iopub.execute_input":"2023-08-27T16:15:15.382686Z","iopub.status.idle":"2023-08-27T16:24:53.461216Z","shell.execute_reply.started":"2023-08-27T16:15:15.382653Z","shell.execute_reply":"2023-08-27T16:24:53.460182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"evaluate(cnn)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:24:53.462477Z","iopub.execute_input":"2023-08-27T16:24:53.462842Z","iopub.status.idle":"2023-08-27T16:25:02.170170Z","shell.execute_reply.started":"2023-08-27T16:24:53.462810Z","shell.execute_reply":"2023-08-27T16:25:02.169014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\n\nplt.figure(figsize=(9, 5))\nhistory = pd.DataFrame({\"loss\": train_lossCNN})\nhistory[\"cum_data\"] = history.index * batch_size\nhistory[\"smooth_loss\"] = history.loss.ewm(halflife=10).mean()\nhistory.plot(x=\"cum_data\", y=\"smooth_loss\", figsize=(12, 5), title=\"train error\")","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:25:02.171794Z","iopub.execute_input":"2023-08-27T16:25:02.172167Z","iopub.status.idle":"2023-08-27T16:25:02.650062Z","shell.execute_reply.started":"2023-08-27T16:25:02.172135Z","shell.execute_reply":"2023-08-27T16:25:02.649148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class AddGaussianNoise(object):\n    def __init__(self, mean=0., std=1.):\n        self.std = std\n        self.mean = mean\n        \n    def __call__(self, tensor):\n        return tensor + torch.randn(tensor.size()) * self.std + self.mean\n    \n    def __repr__(self):\n        return self.__class__.__name__ + '(mean={0}, std={1})'.format(self.mean, self.std)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:25:02.651530Z","iopub.execute_input":"2023-08-27T16:25:02.652574Z","iopub.status.idle":"2023-08-27T16:25:02.659175Z","shell.execute_reply.started":"2023-08-27T16:25:02.652537Z","shell.execute_reply":"2023-08-27T16:25:02.658284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Add Gaussian noise to the data, with sigma = 25","metadata":{}},{"cell_type":"code","source":"path = \"/kaggle\"\nMNIST = torch.utils.data.DataLoader(\n    torchvision.datasets.MNIST(path, train=True, download=True,\n                             transform=torchvision.transforms.Compose([\n                                        torchvision.transforms.ToTensor(),\n                                        torchvision.transforms.Normalize((0.1307,), (0.3081,)),\n                                        torchvision.transforms.RandomApply([AddGaussianNoise(0, 10)], p=0.5)\n                            ])\n    ),\n    batch_size=batch_size, shuffle=True\n)\n\n#evens = list(range)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:25:02.660682Z","iopub.execute_input":"2023-08-27T16:25:02.661040Z","iopub.status.idle":"2023-08-27T16:25:03.499607Z","shell.execute_reply.started":"2023-08-27T16:25:02.661009Z","shell.execute_reply":"2023-08-27T16:25:03.498568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fit(cnn, MNIST)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:25:03.501108Z","iopub.execute_input":"2023-08-27T16:25:03.502006Z","iopub.status.idle":"2023-08-27T16:35:48.479189Z","shell.execute_reply.started":"2023-08-27T16:25:03.501968Z","shell.execute_reply":"2023-08-27T16:35:48.478142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train(epoch, train_loader):\n    num_items = 0\n    train_losses = []\n\n    model.train()\n    criterion = nn.CrossEntropyLoss()\n    print(f\"Training Epoch {epoch}...\")\n    for batch_idx, (data, target) in tqdm(enumerate(train_loader), total=len(train_loader)):\n        if use_cuda:\n            data = data.cuda()\n            target = target.cuda()\n        optimizer.zero_grad()\n        output = model(data)\n        loss = criterion(output, target)\n        loss.backward()\n        optimizer.step()\n\n        train_losses += [loss.item()]\n        num_items += data.shape[0]\n    print('Train loss: {:.5f}'.format(np.mean(train_losses)))\n    return train_losses\n\n\ndef test(test_loader):\n    accuracy = 0.0\n    num_items = 0\n\n    model.eval()\n    criterion = nn.CrossEntropyLoss()\n    print(f\"Testing...\")\n    with torch.no_grad():\n        for batch_idx, (data, target) in tqdm(enumerate(test_loader),  total=len(test_loader)):\n            if use_cuda:\n                data = data.cuda()\n                target = target.cuda()\n            output = model(data)\n            accuracy += torch.sum(torch.argmax(output, dim=1) == target).item()\n            num_items += data.shape[0]\n    accuracy = accuracy * 100 / num_items\n    print(\"Test Accuracy: {:.3f}%\".format(accuracy))","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:35:48.480914Z","iopub.execute_input":"2023-08-27T16:35:48.481315Z","iopub.status.idle":"2023-08-27T16:35:48.491730Z","shell.execute_reply.started":"2023-08-27T16:35:48.481279Z","shell.execute_reply":"2023-08-27T16:35:48.490753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_epochs = 5\ntest(test_loader)\ntrain_losses = []\nfor epoch in range(1, n_epochs + 1):\n    train_losses += train(epoch, MNIST)\n    test(test_loader)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:35:48.493236Z","iopub.execute_input":"2023-08-27T16:35:48.493861Z","iopub.status.idle":"2023-08-27T16:44:31.191748Z","shell.execute_reply.started":"2023-08-27T16:35:48.493826Z","shell.execute_reply":"2023-08-27T16:44:31.190792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# def visualize_model(model,data, typ=\"ode\", sigma = 50.0):\n#     import random\n#     cnt = 0\n#     tr = []; fal = []\n#     for _, dp in list(enumerate(data)):\n#         np_dp = torch.tensor(np.array(dp[0]).reshape(64, 7, 8, 4)).float()\n#         np_dp = np_dp + torch.normal(torch.zeros(np_dp.shape),torch.ones(np_dp.shape) * sigma).float()\n#         # np_dp = np_dp.to(device)\n#         preds = model(np_dp)\n#         preds = torch.argmax(preds,dim=1).item()\n#         if preds == dp[1]:  tr.append((np_dp,dp[1],preds))\n#         else: fal.append((np_dp,dp[1],preds))\n#         cnt+=1\n#         if cnt == 1001: break\n#     print(f\"With sigma = {sigma}, the number of images going to wrong is: {len(fal)} / 2000 (images)\")\n#     random.shuffle(tr)\n#     random.shuffle(fal)\n#     fig, ax = plt.subplots(2,5,sharex=True,sharey=True)\n#     for i in range(5):\n#         ax[0][i].imshow(tr[i][0].cpu().detach().numpy().reshape((28,28)),cmap=\"gray\")\n#         ax[0][i].set_xlabel(f\"truths: {tr[i][1]}\\nPred: {tr[i][2]}\")\n#         ax[1][i].imshow(fal[i][0].cpu().detach().numpy().reshape((28,28)),cmap=\"gray\")\n#         ax[1][i].set_xlabel(f\"truths: {fal[i][1]}\\nPred: {fal[i][2]}\")\n#     #plt.xlabel(f\"{typ} model with sigma {sigma}\")\n#     plt.show()\n\n# sigma_l = [0.0, 50.0, 75.0, 100.0]\n# for sig in sigma_l:\n#     visualize_model(model, MNIST, typ=\"odenet\", sigma=sig)\n#     visualize_model(cnn, MNIST, typ=\"cnnnet\", sigma=sig)","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:44:31.193267Z","iopub.execute_input":"2023-08-27T16:44:31.193846Z","iopub.status.idle":"2023-08-27T16:44:31.200009Z","shell.execute_reply.started":"2023-08-27T16:44:31.193811Z","shell.execute_reply":"2023-08-27T16:44:31.199025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_gaussian_noise(X_train, mean, stddev):\n    ''' \n    INPUT:  (1) 4D numpy array: all raw training image data, of shape \n                (#imgs, #chan, #rows, #cols)\n            (2) float: the mean of the Gaussian to sample noise from\n            (3) float: the standard deviation of the Gaussian to sample\n                noise from. Note that the range of pixel values is\n                0-255; choose the standard deviation appropriately. \n    OUTPUT: (1) 4D numpy array: noisy training data, of shape\n                (#imgs, #chan, #rows, #cols)\n    '''\n    n_imgs = X_train.shape[0]\n    n_chan = X_train.shape[1]\n    n_rows = X_train.shape[2]\n    n_cols = X_train.shape[3]\n    if stddev == 0:\n        noise = np.zeros((n_imgs, n_chan, n_rows, n_cols))\n    else:\n        noise = np.random.normal(mean, stddev/255., \n                                 (n_imgs, n_chan, n_rows, n_cols))\n    noisy_X = X_train + noise\n    clipped_noisy_X = np.clip(noisy_X, 0., 1.)\n    return clipped_noisy_X","metadata":{"execution":{"iopub.status.busy":"2023-08-27T16:44:31.201472Z","iopub.execute_input":"2023-08-27T16:44:31.201807Z","iopub.status.idle":"2023-08-27T16:44:31.214394Z","shell.execute_reply.started":"2023-08-27T16:44:31.201775Z","shell.execute_reply":"2023-08-27T16:44:31.213483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}