Examples/Example_PDE_system_identification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3bc728b9",
   "metadata": {},
   "source": [
    "## PDE system identification using the SCDT"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "198ef22c",
   "metadata": {},
   "source": [
    "The dynamical system of interest:<br><br>\n",
    "\\begin{equation}\n",
    "    \\xi_0 \\ddot{u}-\\xi_1 u_{xx}-\\xi_2 \\dot{u}_{xx}-\\xi_3 \\ddot{u}_{xx} - \\xi_4 u_{xxxx}+\\beta u_x u_{xx} = 0,\n",
    "\\end{equation}<br>\n",
    "We want to identify the nonlinear parameter $\\beta$ using the sensor measurement $s(t)=\\dot{u}(x_m,t)$ at location $x=x_m$. The sensor measurements were generated by simulating the above PDE using the spectral method. The initial conditions for the displacement and the velocity are given by:\n",
    "\\begin{align}\n",
    "    &u(x,0) = e^{-\\frac{(x-x_0)^2}{2\\sigma^2}} \\nonumber\\\\\n",
    "    &\\dot{u}(x,0) = \\frac{\\nu(x - x_0)}{\\sigma^2}u(x,0),\n",
    "\\end{align}\n",
    "where $\\dot{u}(x,0):=\\frac{\\partial u(x,t)}{\\partial t}|_{t=0}$ and $\\nu$ is the wave speed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0170d895",
   "metadata": {},
   "source": [
    "## Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ce148da9",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Import necessary python packages\n",
    "\n",
    "import numpy as np\n",
    "from numpy import interp\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.linalg import lstsq\n",
    "from scipy import signal\n",
    "import os\n",
    "import h5py\n",
    "from pathlib import Path\n",
    "import random\n",
    "import time\n",
    "\n",
    "import sys\n",
    "sys.path.append('../')\n",
    "from pytranskit.optrans.continuous import SCDT\n",
    "\n",
    "import numpy.linalg as LA\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import confusion_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37dfd00d",
   "metadata": {},
   "source": [
    "### Define the regression model for system identification using the SCDT-NLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9b8a9642",
   "metadata": {},
   "outputs": [],
   "source": [
    "class SCDT_NLS_sysId:\n",
    "    def __init__(self, num_classes):\n",
    "        default = False\n",
    "        self.num_classes = num_classes\n",
    "        self.Nset = []\n",
    "        self.subspaces = []\n",
    "        self.len_subspace = 0\n",
    "        self.k = 1\n",
    "        self.label = []\n",
    "        self.pca_basis = []\n",
    "        self.n_enrV = 1\n",
    "\n",
    "    def fit(self, X, Y, Beta, no_local_enrichment=False):      \n",
    "        Xtrain, Xval, Ytrain, Yval = train_test_split(X, Y, test_size=0.3, random_state=0)\n",
    "        self.bas = []\n",
    "        for class_idx in range(self.num_classes):\n",
    "            # generate the bases vectors\n",
    "            class_data = Xtrain[Ytrain == class_idx]\n",
    "            self.Nset.append(class_data)\n",
    "            self.label.append(class_idx)\n",
    "            bas = []\n",
    "            for j in range(class_data.shape[0]):\n",
    "                flat = np.copy(class_data[j].reshape(1,-1))\n",
    "                u, s, vh = LA.svd(flat,full_matrices=False)\n",
    "                bas.append(vh[:flat.shape[0]])\n",
    "            self.bas.append(bas)\n",
    "            \n",
    "        if Xtrain.shape[0]//self.num_classes < 4:\n",
    "            self.k = 1\n",
    "        else:\n",
    "            smp_class = []\n",
    "            for i in range(len(np.unique(Ytrain))):\n",
    "                smp_class.append(np.count_nonzero(Ytrain == i))\n",
    "            \n",
    "            k_range = range(1,min(min(smp_class),100))\n",
    "            n_range = range(-1,2)  # range(-1,0) means N=-1, i.e., no tuning of N\n",
    "            \n",
    "            self.k, self.n_enrV = self.find_kN(Xval, Yval, k_range, n_range)        \n",
    "        self.Nset = []\n",
    "        self.label = []\n",
    "        self.bas = []\n",
    "        self.Beta = []\n",
    "        for class_idx in range(self.num_classes):\n",
    "            # generate the bases vectors\n",
    "            class_data = X[Y == class_idx]\n",
    "            self.Beta.append(Beta[Y == class_idx])\n",
    "            self.Nset.append(class_data)\n",
    "            self.label.append(class_idx)\n",
    "            bas = []\n",
    "            for j in range(class_data.shape[0]):\n",
    "                if no_local_enrichment:\n",
    "                    flat = np.copy(class_data[j].reshape(1,-1))\n",
    "                else:\n",
    "                    flat = self.add_k_enrich(class_data[j].reshape(1,-1), k=self.n_enrV) # k=0 => translation only\n",
    "                u, s, vh = LA.svd(flat,full_matrices=False)\n",
    "                bas.append(vh[:flat.shape[0]])\n",
    "            self.bas.append(bas)\n",
    "\n",
    "    def predict(self, X, use_gpu=False, k=None, N=None):\n",
    "        if k is not None:\n",
    "            k_opt = k\n",
    "        else:\n",
    "            k_opt = self.k\n",
    "            \n",
    "        if N is not None:\n",
    "            n_opt = N\n",
    "        else:\n",
    "            n_opt = self.n_enrV\n",
    "            \n",
    "        V = np.ones([1, X.shape[1]])        \n",
    "        lmd = 0.001\n",
    "        D = []\n",
    "        Beta = []\n",
    "        \n",
    "        for class_idx in range(self.num_classes):\n",
    "            Xi = self.Nset[class_idx]\n",
    "            Xi_bas = self.bas[class_idx]\n",
    "            Xi_beta = self.Beta[class_idx]\n",
    "            d = np.zeros([X.shape[0],1])\n",
    "            beta_class = np.zeros([X.shape[0],1])\n",
    "            B = []\n",
    "            L_basis = []\n",
    "            for i in range(X.shape[0]):\n",
    "                x = X[i,:]\n",
    "                dist_i = []\n",
    "\n",
    "                for j in range(Xi.shape[0]):\n",
    "                    basj = Xi_bas[j]#[:self.len_subspace,:]\n",
    "                    projR = x @ basj.T  @ basj  # (n_samples, n_features)\n",
    "                    dist_i.append(LA.norm(projR - x))\n",
    "                dist_i = np.stack(dist_i)\n",
    "\n",
    "                indx = dist_i.argsort()[:k_opt]\n",
    "\n",
    "                # calculate est_beta for this class\n",
    "                f_norm = np.sum([1./dist_i[j] for j in indx])\n",
    "                beta_class[i] = np.sum([Xi_beta[j]/dist_i[j] for j in indx])/f_norm\n",
    "\n",
    "                # calculate local subspace distance\n",
    "                Ni = self.add_k_enrich(Xi[indx,:], k=n_opt) # k=0 => translation only\n",
    "\n",
    "                u, s, vh = LA.svd(Ni,full_matrices=False)\n",
    "\n",
    "                cum_s = np.cumsum(s)\n",
    "                cum_s = cum_s/np.max(cum_s)\n",
    "                basis = vh[:Ni.shape[0]]\n",
    "                B.append(basis)\n",
    "                L_basis.append((np.where(cum_s>=0.99)[0])[0]+1)\n",
    "            max_basis = min(L_basis)\n",
    "            Beta.append(np.squeeze(beta_class))\n",
    "            for i in range(X.shape[0]):\n",
    "                x = X[i,:]\n",
    "                basis = B[i][:max_basis,:]\n",
    "\n",
    "                proj = x @ basis.T  # (n_samples, n_basis)\n",
    "                projR = proj @ basis  # (n_samples, n_features)\n",
    "\n",
    "                d[i]=LA.norm(projR - x)                    \n",
    "            D.append(np.squeeze(d))\n",
    "\n",
    "        \n",
    "        D = np.stack(D, axis=0)\n",
    "        preds = np.argmin(D, axis=0)\n",
    "        pred_label = [self.label[i] for i in preds]\n",
    "\n",
    "        return  pred_label\n",
    "\n",
    "    def find_kN(self, X, y, k_range, n_range):\n",
    "        n = X.shape[0]        \n",
    "        max_acc = 0.\n",
    "        score_prev = 0.\n",
    "        k_opt = 1\n",
    "        count = 0\n",
    "        acc_count = 0\n",
    "\n",
    "        ### calculate distances for samples in validation set\n",
    "        indx = []\n",
    "        for i in range(X.shape[0]):\n",
    "            x = np.copy(X[i,:])\n",
    "            indXi = []\n",
    "            for class_idx in range(self.num_classes):\n",
    "                Xi = self.Nset[class_idx]\n",
    "                Xi_bas = self.bas[class_idx]\n",
    "                dist_i = []\n",
    "\n",
    "                for j in range(Xi.shape[0]):\n",
    "                    basj = Xi_bas[j]#[:self.len_subspace,:]\n",
    "                    projR = x @ basj.T  @ basj  # (n_samples, n_features)\n",
    "                    dist_i.append(LA.norm(projR - x))\n",
    "                dist_i = np.stack(dist_i)\n",
    "\n",
    "                indXi.append(dist_i.argsort()[:max(k_range)+1])\n",
    "            indx.append(indXi)\n",
    "\n",
    "        ### tune k using validation set\n",
    "        for k in k_range:\n",
    "            D = []\n",
    "            for class_idx in range(self.num_classes):\n",
    "                Xi = self.Nset[class_idx]\n",
    "                d = np.zeros([X.shape[0],1])\n",
    "                B = []\n",
    "                L_basis = []\n",
    "                for i in range(X.shape[0]):\n",
    "                    x = np.copy(X[i,:])\n",
    "                    ind = indx[i][class_idx]\n",
    "                    Ni = np.copy(Xi[ind[:k],:])\n",
    "                    u, s, vh = LA.svd(Ni,full_matrices=False)\n",
    "                    cum_s = np.cumsum(s)\n",
    "                    cum_s = cum_s/np.max(cum_s)\n",
    "                    basis = vh[:Ni.shape[0]]\n",
    "                    B.append(basis)\n",
    "                    L_basis.append((np.where(cum_s>=0.99)[0])[0]+1)\n",
    "                max_basis = min(L_basis)\n",
    "                for i in range(X.shape[0]):\n",
    "                    x = np.copy(X[i,:])\n",
    "                    basis = B[i][:max_basis,:]\n",
    "                    projR = x @ basis.T @ basis  # (n_samples, n_features)\n",
    "                    d[i]=LA.norm(projR - x)\n",
    "                D.append(np.squeeze(d))\n",
    "            D = np.stack(D, axis=0)\n",
    "            preds = np.argmin(D, axis=0)\n",
    "            pred_label = [self.label[i] for i in preds]\n",
    "            score = (np.sum(pred_label == y))/n\n",
    "            print('Validation accuracy: {} with k = {}'.format(score, k))\n",
    "            if score > score_prev:\n",
    "                count = 0\n",
    "            else:\n",
    "                count = count + 1\n",
    "            if score >= max_acc:\n",
    "                max_acc = score\n",
    "                k_opt = k\n",
    "                acc_count = 0\n",
    "                #count = 0\n",
    "            else:\n",
    "                acc_count = acc_count + 1\n",
    "            \n",
    "            if count == 10 or acc_count == 20:\n",
    "                break\n",
    "            score_prev = score\n",
    "        \n",
    "        ### tune N using validation set\n",
    "        if len(n_range) == 1:\n",
    "            return k_opt, n_range[0]\n",
    "        \n",
    "        n_iter = []\n",
    "        max_acc = 0.\n",
    "        score_prev = 0.\n",
    "        n_opt = 1\n",
    "        count = 0\n",
    "        acc_count = 0\n",
    "        \n",
    "        for n_enr in n_range:  \n",
    "            print('\\nN = {}'.format(n_enr))\n",
    "            n_iter.append(n_enr)\n",
    "    \n",
    "            D = []\n",
    "            for class_idx in range(self.num_classes):\n",
    "                Xi = self.Nset[class_idx]\n",
    "                d = np.zeros([X.shape[0],1])\n",
    "                B = []\n",
    "                L_basis = []\n",
    "                for i in range(X.shape[0]):\n",
    "                    x = np.copy(X[i,:])\n",
    "                    ind = indx[i][class_idx]\n",
    "                    Ni = self.add_k_enrich(Xi[ind[:k_opt],:], k=n_enr) # k=0 => translation only\n",
    "                    u, s, vh = LA.svd(Ni,full_matrices=False)\n",
    "                    cum_s = np.cumsum(s)\n",
    "                    cum_s = cum_s/np.max(cum_s)\n",
    "                    basis = vh[:Ni.shape[0]]\n",
    "                    B.append(basis)\n",
    "                    L_basis.append((np.where(cum_s>=0.99)[0])[0]+1)\n",
    "                max_basis = min(L_basis)\n",
    "                for i in range(X.shape[0]):\n",
    "                    x = X[i,:]\n",
    "                    basis = B[i][:max_basis,:]\n",
    "                    projR = x @ basis.T @ basis  # (n_samples, n_features)\n",
    "                    d[i]=LA.norm(projR - x)\n",
    "                D.append(np.squeeze(d))\n",
    "            D = np.stack(D, axis=0)\n",
    "            preds = np.argmin(D, axis=0)\n",
    "            pred_label = [self.label[i] for i in preds]\n",
    "            score = (np.sum(pred_label == y))/n\n",
    "            print('Validation accuracy: {} with k = {}'.format(score, k_opt))\n",
    "            if score > max_acc:\n",
    "                max_acc = score\n",
    "                acc_count = 0\n",
    "                n_opt = n_enr\n",
    "            else:\n",
    "                acc_count = acc_count + 1\n",
    "            if score > score_prev:\n",
    "                count = 0\n",
    "            else:\n",
    "                count = count + 1\n",
    "            if count == 10 or acc_count == 20:\n",
    "                break\n",
    "            score_prev = score\n",
    "            \n",
    "        return k_opt, n_opt\n",
    "        \n",
    "    def score(self, X, y):\n",
    "#         print('Optimum k: {}'.format(self.k))\n",
    "#         print('Optimum N: {}'.format(self.n_enrV))\n",
    "        n = X.shape[0]\n",
    "        y_pred = self.predict(X)\n",
    "        n_correct = np.sum(y_pred == y)\n",
    "        return n_correct/n, y_pred\n",
    "    \n",
    "    def add_k_enrich(self, scdt_features, k):\n",
    "        # scdt_features: (n_samples, scdt)\n",
    "        if k<0:\n",
    "            return scdt_features\n",
    "        v= np.ones([1, scdt_features.shape[1]]) # add translation\n",
    "        indx = 0\n",
    "        for i in range(-k,k+1):\n",
    "            if i != 0:\n",
    "                vi = scdt_features-np.sin(i*np.pi*scdt_features)/(np.abs(i)*np.pi)\n",
    "                v = np.concatenate((v,vi))            \n",
    "            indx = indx+1\n",
    "        return np.concatenate((scdt_features,v))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21b2fd33",
   "metadata": {},
   "source": [
    "### Read sensor measurements for training and testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "15bd9d2c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "====1D_wave_deform_3class_2k====\n",
      "\n",
      "Number of classes: 3\n",
      "Train set: (6000, 375)\n",
      "Test set: (600, 375)\n"
     ]
    }
   ],
   "source": [
    "datadir = '../../../PDE_nonlinear_systems/Codes/data/PDE' # '/Path/to/data/'\n",
    "\n",
    "## dataset options:\n",
    "# 2 class: 1D_wave_deform_2class_2k\n",
    "# 3 class: 1D_wave_deform_3class_2k\n",
    "# 10 class: 1D_wave_deform_10class_2k\n",
    "\n",
    "dataset = '1D_wave_deform_3class_2k'\n",
    "\n",
    "print('\\n===='+dataset+'====\\n')\n",
    "\n",
    "data_file = os.path.join(datadir, dataset, 'train.hdf5')\n",
    "with h5py.File(data_file, 'r') as f:\n",
    "    x_train = f['x_train'][()]\n",
    "    y_train = f['y_train'][()]\n",
    "    t_train = f['t_train'][()]\n",
    "    p_train = f['p_train'][()]\n",
    "num_classes = len(np.unique(y_train))\n",
    "\n",
    "data_file = os.path.join(datadir, dataset, 'test.hdf5')\n",
    "with h5py.File(data_file, 'r') as f:\n",
    "    x_test = f['x_test'][()]\n",
    "    y_test = f['y_test'][()]\n",
    "    t_test = f['t_test'][()]\n",
    "    p_test = f['p_test'][()]\n",
    "\n",
    "print('Number of classes: {}'.format(num_classes))\n",
    "print('Train set: {}'.format(x_train.shape))\n",
    "print('Test set: {}'.format(x_test.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3731db90",
   "metadata": {},
   "source": [
    "### Define the classes based on $\\beta$ ranges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d382c1f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# class ranges\n",
    "if num_classes==2:\n",
    "    class_int = [[0.0, 0.], [0.01, 0.6]]\n",
    "elif num_classes==3:\n",
    "    class_int = [[0.01, 0.2], [0.21, 0.4], [0.41, 0.6]]\n",
    "elif num_classes==10:\n",
    "    class_int = [[0.0, 0.06], [0.07, 0.12], [0.13, 0.18], [0.19, 0.24], [0.25, 0.30],\n",
    "[0.31, 0.36], [0.37, 0.42], [0.43, 0.48], [0.49, 0.54], [0.55, 0.60]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9001b10e",
   "metadata": {},
   "source": [
    "### Calculate SCDT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "50134943",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loaded from cache file data: x_train (6000, 746) x_test (600, 746)\n"
     ]
    }
   ],
   "source": [
    "N = x_train.shape[1]\n",
    "## define uniform reference signal\n",
    "t0 = np.linspace(0,1,N) # Domain of the reference\n",
    "s0 = np.ones(N)\n",
    "s0 = s0/s0.sum()\n",
    "\n",
    "scdt = SCDT(reference=s0,x0=t0)\n",
    "\n",
    "rm_edge = True\n",
    "cache_file = os.path.join(datadir, dataset, 'scdt_wMass.hdf5')  # 'scdt_wMass.hdf5'\n",
    "\n",
    "# If the SCDTs are already computed, then load from the saved file. \n",
    "# Otherwise, compute SCDTs for all the samples and save them in the\n",
    "# same directory as the dataset\n",
    "\n",
    "if os.path.exists(cache_file):\n",
    "    with h5py.File(cache_file, 'r') as f:\n",
    "        data_train, y_train = f['data_train'][()], f['y_train'][()]\n",
    "        data_test, y_test = f['data_test'][()], f['y_test'][()]\n",
    "        data_train, data_test = data_train.astype(np.float32), data_test.astype(np.float32)\n",
    "        print('loaded from cache file data: x_train {} x_test {}'.format(data_train.shape, data_test.shape))\n",
    "else:\n",
    "    with h5py.File(cache_file, 'w') as f:            \n",
    "        data_train = []\n",
    "        data_test = []\n",
    "        for i in range(x_train.shape[0]):\n",
    "            Ipos, Ineg, Imasspos, Imassneg = scdt.stransform(x_train[i])\n",
    "            if rm_edge:\n",
    "                data_train.append(np.concatenate((Ipos[1:-2],Ineg[1:-2],Imasspos.reshape(1),Imassneg.reshape(1)),axis=0))\n",
    "            else:\n",
    "                data_train.append(np.concatenate((Ipos[:-1],Ineg[:-1],Imasspos.reshape(1),Imassneg.reshape(1)),axis=0))\n",
    "        for i in range(x_test.shape[0]):\n",
    "            Ipos, Ineg, Imasspos, Imassneg = scdt.stransform(x_test[i])\n",
    "            if rm_edge:\n",
    "                data_test.append(np.concatenate((Ipos[1:-2],Ineg[1:-2], Imasspos.reshape(1), Imassneg.reshape(1)),axis=0))\n",
    "            else:\n",
    "                data_test.append(np.concatenate((Ipos[:-1],Ineg[:-1], Imasspos.reshape(1), Imassneg.reshape(1)),axis=0))\n",
    "        data_train = np.stack(data_train)\n",
    "        data_test = np.stack(data_test)            \n",
    "        print('Train set: {}'.format(data_train.shape))\n",
    "        print('Test set: {}'.format(data_test.shape))\n",
    "\n",
    "        data_train, data_test = data_train.astype(np.float32), data_test.astype(np.float32)\n",
    "        f.create_dataset('data_train', data=data_train)\n",
    "        f.create_dataset('y_train', data=y_train)\n",
    "        f.create_dataset('data_test', data=data_test)\n",
    "        f.create_dataset('y_test', data=y_test)\n",
    "        print('saved to {}'.format(cache_file))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e244e0cf",
   "metadata": {},
   "source": [
    "### Train and test the system identification model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "dc8ddd68",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dataset: 1D_wave_deform_3class_2k\n",
      "\n",
      "\n",
      "Train the model...\n",
      "Validation accuracy: 0.9483333333333334 with k = 1\n",
      "Validation accuracy: 0.9522222222222222 with k = 2\n",
      "Validation accuracy: 0.9555555555555556 with k = 3\n",
      "Validation accuracy: 0.9577777777777777 with k = 4\n",
      "Validation accuracy: 0.9566666666666667 with k = 5\n",
      "Validation accuracy: 0.9566666666666667 with k = 6\n",
      "Validation accuracy: 0.9566666666666667 with k = 7\n",
      "Validation accuracy: 0.9583333333333334 with k = 8\n",
      "Validation accuracy: 0.9583333333333334 with k = 9\n",
      "Validation accuracy: 0.9572222222222222 with k = 10\n",
      "Validation accuracy: 0.9544444444444444 with k = 11\n",
      "Validation accuracy: 0.9505555555555556 with k = 12\n",
      "Validation accuracy: 0.9527777777777777 with k = 13\n",
      "Validation accuracy: 0.9516666666666667 with k = 14\n",
      "Validation accuracy: 0.9527777777777777 with k = 15\n",
      "Validation accuracy: 0.9511111111111111 with k = 16\n",
      "Validation accuracy: 0.9488888888888889 with k = 17\n",
      "Validation accuracy: 0.9472222222222222 with k = 18\n",
      "Validation accuracy: 0.9483333333333334 with k = 19\n",
      "Validation accuracy: 0.9444444444444444 with k = 20\n",
      "Validation accuracy: 0.9444444444444444 with k = 21\n",
      "Validation accuracy: 0.9438888888888889 with k = 22\n",
      "Validation accuracy: 0.9438888888888889 with k = 23\n",
      "Validation accuracy: 0.9422222222222222 with k = 24\n",
      "Validation accuracy: 0.9427777777777778 with k = 25\n",
      "Validation accuracy: 0.9422222222222222 with k = 26\n",
      "Validation accuracy: 0.9416666666666667 with k = 27\n",
      "Validation accuracy: 0.9411111111111111 with k = 28\n",
      "Validation accuracy: 0.9416666666666667 with k = 29\n",
      "\n",
      "N = -1\n",
      "Validation accuracy: 0.9583333333333334 with k = 9\n",
      "\n",
      "N = 0\n",
      "Validation accuracy: 0.9477777777777778 with k = 9\n",
      "\n",
      "N = 1\n",
      "Validation accuracy: 0.9466666666666667 with k = 9\n",
      "\n",
      "Test the model...\n",
      "\n",
      "\n",
      "+++++++++++++++++++Test Results+++++++++++++++++++\n",
      "\n",
      "Number of classes: 3\n",
      "\n",
      "Detection Accuracy: 97.17%\n",
      "\n",
      "Confusion matrix:\n",
      "[[197   3   0]\n",
      " [  6 190   4]\n",
      " [  0   4 196]]\n",
      "\n",
      "Estimation Error (MSE): 0.0031443231976238926\n"
     ]
    }
   ],
   "source": [
    "print('Dataset: '+dataset+'\\n')\n",
    "\n",
    "accuracies = []\n",
    "predictions = []\n",
    "\n",
    "\n",
    "cls = SCDT_NLS_sysId(num_classes)\n",
    "print('\\nTrain the model...')\n",
    "cls.fit(data_train, y_train, p_train[:,-1]) # [:,:-2]\n",
    "\n",
    "print('\\nTest the model...')\n",
    "acc, pred = cls.score(data_test[:,:], y_test[:]) # [:,:-2]\n",
    "\n",
    "print('\\n\\n+++++++++++++++++++Test Results+++++++++++++++++++\\n')\n",
    "print('Number of classes: {}\\n'.format(num_classes))\n",
    "print('Detection Accuracy: {:.2f}%'.format(acc*100.))\n",
    "print('\\nConfusion matrix:')\n",
    "print(confusion_matrix(y_test, pred))\n",
    "accuracies.append(acc*100.)\n",
    "predictions.append(pred)\n",
    "\n",
    "se = []\n",
    "for i in range(len(pred)):\n",
    "    beta = p_test[i,-1]\n",
    "    se.append((beta - np.mean(class_int[pred[i]]))**2)\n",
    "\n",
    "MSE_dist = np.mean(se)\n",
    "print('\\nEstimation Error (MSE): {}'.format(MSE_dist))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4489438",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}