{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Python packages used in this code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import os\n",
    "import pickle\n",
    "import time\n",
    "import sklearn\n",
    "import platform\n",
    "import sys\n",
    "from sklearn.base import BaseEstimator, RegressorMixin\n",
    "from sklearn.model_selection import train_test_split, GridSearchCV\n",
    "from sklearn.linear_model import Ridge, LinearRegression, Lasso, ElasticNet\n",
    "from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score\n",
    "\n",
    "from sklearn.linear_model._base import LinearModel\n",
    "from spmimage.linear_model import admm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Environments\n",
    "\n",
    "--Platform--\n",
    "OS : Windows-10-10.0.19044-SP0\n",
    "--Version--\n",
    "python :  3.9.12 (main, Apr  4 2022, 05:22:27) [MSC v.1916 64 bit (AMD64)]\n",
    "numpy : 1.23.1\n",
    "pandas : 1.4.3\n",
    "sklearn : 1.1.1\n",
    "\"\"\"\n",
    "\n",
    "print('--Platform--')\n",
    "print('OS :', platform.platform())\n",
    "print('--Version--')\n",
    "print('python : ', sys.version)\n",
    "print('numpy :', np.__version__)\n",
    "print('pandas :', pd.__version__)\n",
    "print('sklearn :', sklearn.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the model cladd for the log-difference model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fused ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "class cls_FusedRidge(BaseEstimator, RegressorMixin):\n",
    "    def __init__(self, lambda_scale=1, lambda_trend=1):\n",
    "        \"\"\"\n",
    "        Define the model class for log-difference model\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "            lambda_scale : reguralization parameter for the scale of gamma\n",
    "            lambda_trend : reguralization parameter for the trend of gamma\n",
    "        \"\"\"\n",
    "        self.lambda_scale = lambda_scale\n",
    "        self.lambda_trend = lambda_trend\n",
    "    \n",
    "    def make_D(self, n_features):\n",
    "        trend_matrix = np.eye(n_features - 1, n_features, k=1) - np.eye(n_features - 1, n_features)\n",
    "        trend_matrix[[9, 29, 49, 69, 89, 109, 129, 149, 169]] = np.zeros(n_features)\n",
    "        \n",
    "        if self.lambda_scale == 0:\n",
    "            return self.lambda_trend * trend_matrix\n",
    "        elif self.lambda_trend == 0:\n",
    "            return self.lambda_scale * np.identity(n_features)\n",
    "        else:\n",
    "            generated = np.vstack([self.lambda_scale * np.identity(n_features),\n",
    "                                   self.lambda_trend * trend_matrix])\n",
    "            return generated\n",
    "        \n",
    "    def fit(self, X, y=None):\n",
    "        \"\"\"\n",
    "        Model fitting\n",
    "        \n",
    "        Required grobal variables\n",
    "        -----------------------\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "\n",
    "        \"\"\"\n",
    "        # dataset\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        \n",
    "        # dimension\n",
    "        self.n_sample, self.dim_X = self.X.shape\n",
    "        \n",
    "        # Matrices\n",
    "        self.D = self.make_D(self.dim_X)\n",
    "        self.D2 = np.transpose(self.D).dot(self.D)\n",
    "        \n",
    "        # Fit\n",
    "        self.intercept = np.mean(self.y)\n",
    "        self.theta = np.linalg.pinv(self.X.T.dot(self.X) + self.D2).dot(self.X.T).dot(self.y-self.intercept)\n",
    "\n",
    "        return self\n",
    "    \n",
    "    def predict(self, X):\n",
    "        \"\"\"\n",
    "        Prediction function\n",
    "            \n",
    "        Returns\n",
    "        -------\n",
    "        \"\"\"\n",
    "        return X.dot(self.theta) + self.intercept\n",
    "\n",
    "    def score(self, X, y=None):\n",
    "        \"\"\"\n",
    "        Score function for cross-validation\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "            -\\sum(y-\\hat{y})/n (Consider the minus value because 'GridSearchCV' maximize the score.)\n",
    "        \"\"\"\n",
    "        return -sum((y.values - self.predict(X).values)**2)/self.n_sample\n",
    "    \n",
    "    def get_params(self, deep=True):\n",
    "        \"\"\"\n",
    "        Create parameter dictionary for cross-validation\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "        \"\"\"\n",
    "        return {'lambda_scale' : self.lambda_scale,\n",
    "                'lambda_trend' : self.lambda_trend}\n",
    "    \n",
    "    def set_params(self, **parameters):\n",
    "        \"\"\"\n",
    "        For cross-validation\n",
    "        \"\"\"\n",
    "        for parameter, value in parameters.items():\n",
    "            setattr(self, parameter, value)\n",
    "        return self     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the model class proposed in the paper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fused ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "class cls_InvTrans_FR(BaseEstimator, RegressorMixin):\n",
    "    def __init__(self, lambda1=1, lambda2=1, lambda3_scale=1, lambda3_trend=1):\n",
    "        \"\"\"\n",
    "        Define the model class proposed in the paper\n",
    "            h(x) = \\alpha_1*ys + \\alpha_0 + (\\beta*ys + 1)<\\gamma, x>\n",
    "                x  : discriptor\n",
    "                ys : CP(MD)\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "            lambda1       : regularization parameter for alpha\n",
    "            lambda2       : regularization parameter for beta\n",
    "            lambda3_scale : reguralization parameter for the scale of gamma\n",
    "            lambda3_trend : reguralization parameter for the trend of gamma\n",
    "        \"\"\"\n",
    "        self.lambda1 = lambda1\n",
    "        self.lambda2 = lambda2\n",
    "        self.lambda3_scale = lambda3_scale\n",
    "        self.lambda3_trend = lambda3_trend\n",
    "    \n",
    "    def estimation_alpha(self):\n",
    "        \"\"\"\n",
    "        Optimization with respect to alpha\n",
    "            [\\hat{\\alpha}_1, \\hat{\\alpha}_2]^T \n",
    "                = (Ys^TYs + n\\Lambda_1)^{-1} Ys^T (y + (\\beta*ys+1)*<\\gamma, x>)\n",
    "            \n",
    "                    Ys = | ys_1  1 | \\in R^{n*2}, \\Lambda_1 = | \\lambda1  0 |, ys = | ys_1 |\n",
    "                         | ys_2  1 |                          |        0  0 |       | ys_2 |\n",
    "                              :                                                     |   :  |\n",
    "                         | ys_n  1 |                                                | ys_n |\n",
    "                         \n",
    "        Note that the regularization applies only to \\alpha_1 and not to the intercept \\alpha_1.\n",
    "        \"\"\"\n",
    "        self.alpha = self.InvMat.dot(self.X_source1.T).dot(self.y + (self.Mat2+1)*(self.Mat3))\n",
    "        self.Mat1 = self.X_source1.dot(self.alpha)\n",
    "        self.result_alpha[self.i_count] = self.alpha\n",
    "        return self\n",
    "    \n",
    "    def estimation_beta(self):\n",
    "        \"\"\"\n",
    "        Optimization with respect to beta\n",
    "            \\hat{\\beta}\n",
    "                = -(ys^T Diag(X\\gamma)^2 ys + n\\lambda_2)^{-1} ys^T Diag(X\\gamma) (y - Ys\\alpha + X\\gamma)\n",
    "        \"\"\"\n",
    "        tmp_mat1 = np.linalg.pinv(self.X_source2.T.dot(np.diag(self.Mat3)).dot(np.diag(self.Mat3)).values.dot(self.X_source2.values) + self.n_sample*self.lambda2*np.diag(np.ones(self.dim_X_source2)), hermitian=True)\n",
    "        self.beta = -tmp_mat1.dot(self.X_source2.T).dot(np.diag(self.Mat3)).dot(self.y-self.Mat1+self.Mat3)\n",
    "        self.Mat2 = self.X_source2.dot(self.beta)\n",
    "        self.result_beta[self.i_count] = self.beta\n",
    "        return self\n",
    "    \n",
    "    def estimation_gamma(self):\n",
    "        \"\"\"\n",
    "        Optimization with respect to gamma\n",
    "        \"\"\"\n",
    "        tmp_x = pd.DataFrame(np.diag(self.X_source2.values.reshape(-1)*self.beta+1).dot(self.X), index=self.X.index, columns=self.X.columns)\n",
    "        tmp_y = self.y-self.Mat1\n",
    "        \n",
    "        self.tmp_x = tmp_x\n",
    "        self.tmp_y = tmp_y\n",
    "                \n",
    "        fix_seed(373)\n",
    "        self.gamma = -np.linalg.pinv(tmp_x.T.dot(tmp_x) + self.D2).dot(tmp_x.T).dot(tmp_y)\n",
    "        \n",
    "        self.Mat3 = self.X.dot(self.gamma)\n",
    "        self.result_gamma[self.i_count] = self.gamma\n",
    "        return self\n",
    "    \n",
    "    def make_diff(self, w_new, w_old):\n",
    "        \"\"\"\n",
    "        Function to calculate parameter changes for algorithm convergence determination\n",
    "        We use \\max{|w_new - w_old|}/\\max{|w_old|} for determining the convergence.\n",
    "        This criterion is used in some algorithms in scikit-learn, for example, see https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html\n",
    "        We apply this criterion to each of \\alpha, \\beta and \\gamma separately and use their maximum value for the convergence decision.\n",
    "        \"\"\"\n",
    "        diff1 = np.max(np.abs(w_new-w_old))\n",
    "        diff2 = np.max(np.abs(w_old))\n",
    "        if diff2 < 1e-10: #Avoiding division by zero\n",
    "            out = 0\n",
    "        else:\n",
    "            out = diff1/diff2\n",
    "        return out\n",
    "        \n",
    "    def make_D(self, n_features):\n",
    "        trend_matrix = np.eye(n_features - 1, n_features, k=1) - np.eye(n_features - 1, n_features)\n",
    "        trend_matrix[[9, 29, 49, 69, 89, 109, 129, 149, 169]] = np.zeros(n_features)\n",
    "        \n",
    "        if self.lambda3_scale == 0:\n",
    "            return self.lambda3_trend * trend_matrix\n",
    "        elif self.lambda3_trend == 0:\n",
    "            return self.lambda3_scale * np.identity(n_features)\n",
    "        else:\n",
    "            generated = np.vstack([self.lambda3_scale * np.identity(n_features),\n",
    "                                   self.lambda3_trend * trend_matrix])\n",
    "            return generated\n",
    "        \n",
    "    def fit(self, X, y=None):\n",
    "        \"\"\"\n",
    "        Model fitting\n",
    "        \n",
    "        Required grobal variables\n",
    "        -----------------------\n",
    "            dim_x         : dimension of the discriptor\n",
    "            ini_alpha     : initial value for \\alpha_1\n",
    "            ini_intercept : initial value for \\alpha_0\n",
    "            ini_beta      : initial value for \\beta\n",
    "            ini_gamma     : initial value for \\gamma\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "            i_count       : counter for the iterations\n",
    "            n_loop        : maximum number of iterations\n",
    "            convergence   : flag indicating whether the algorithm has converged before the maximum iteration\n",
    "            error         : flag indication whether the algorithm has terminated with an error\n",
    "            \n",
    "            Input         : X=[discriptors, source features], y=output\n",
    "            X             : descriptors\n",
    "            X_source1     : source features + all-one vector\n",
    "            X_source2     : source features\n",
    "            y             : output\n",
    "            n_sample      : number of sumples\n",
    "            dim_X         : dimension of the discriptor\n",
    "            dim_X_source1 : dimension of the source features + intercept\n",
    "            dim_X_source2 : dimension of the source features\n",
    "\n",
    "            result_alpha  : dataframe to store alpha in all iterations\n",
    "            result_beta   : dataframe to store beta in all iterations\n",
    "            result_gamma  : dataframe to store gamma in all iterations\n",
    "            diff_i        : series to store the difference between \\alpha_0_new and \\alpha_0_old \n",
    "            diff_a        : series to store the difference between \\alpha_1_new and \\alpha_1_old \n",
    "            diff_b        : series to store the difference between \\beta_new and \\beta_old\n",
    "            diff_c        : series to store the difference between \\gamma_new and \\gamma_old \n",
    "            diff          : series to store the difference between AllParams_new and AllParams_old \n",
    "\n",
    "            InvMat        : (Ys^TYs + n\\Lambda_1)^{-1}\n",
    "            Mat1          : Ys\\alpha (updated with every update of \\alpha)\n",
    "            Mat2          : ys\\beta (updated with every update of \\alpha)\n",
    "            Mat3          : X\\gamma (updated with every update of \\alpha)\n",
    "        \"\"\"\n",
    "        # setting\n",
    "        self.i_count = 0        \n",
    "        self.n_loop = 1000\n",
    "        self.convergence = False\n",
    "        self.error = False\n",
    "        \n",
    "        # dataset\n",
    "        self.X = X.iloc[:,:dim_x]\n",
    "        self.X_source1 = X.iloc[:,dim_x:].copy()\n",
    "        self.X_source1['Intercept'] = 1\n",
    "        self.X_source2 = X.iloc[:,dim_x:].copy()\n",
    "        self.y = y\n",
    "        \n",
    "        # dimension\n",
    "        self.n_sample, self.dim_X = self.X.shape\n",
    "        self.dim_X_source1 = self.X_source1.shape[1]\n",
    "        self.dim_X_source2 = self.X_source2.shape[1]\n",
    "        \n",
    "        # for storing the results\n",
    "        self.result_alpha = np.zeros([self.n_loop+1, self.dim_X_source1])\n",
    "        self.result_beta = np.zeros([self.n_loop+1, self.dim_X_source2])\n",
    "        self.result_gamma = np.zeros([self.n_loop+1, self.dim_X])\n",
    "        self.diff_i = np.zeros(self.n_loop+1)\n",
    "        self.diff_a = np.zeros(self.n_loop+1)\n",
    "        self.diff_b = np.zeros(self.n_loop+1)\n",
    "        self.diff_c = np.zeros(self.n_loop+1)\n",
    "        self.diff = np.zeros(self.n_loop+1)\n",
    "        \n",
    "        # initialize the parameters\n",
    "        self.alpha = np.array([ini_alpha, ini_intercept])\n",
    "        self.beta = np.array([ini_beta])\n",
    "        self.gamma = -ini_gamma\n",
    "        self.result_alpha[0] = self.alpha\n",
    "        self.result_beta[0] = self.beta\n",
    "        self.result_gamma[0] = self.gamma\n",
    "        self.diff[0] = np.nan\n",
    "        \n",
    "        # Matrices\n",
    "        self.InvMat = np.linalg.pinv(self.X_source1.T.dot(self.X_source1) + self.n_sample*self.lambda1*np.diag(np.ones(self.dim_X_source1-1).tolist()+[0]), hermitian=True)\n",
    "        self.Mat1 = self.X_source1.dot(self.alpha)\n",
    "        self.Mat2 = self.X_source2.dot(self.beta)\n",
    "        self.Mat3 = self.X.dot(self.gamma)\n",
    "        self.D = self.make_D(self.dim_X)\n",
    "        self.D2 = np.transpose(self.D).dot(self.D)\n",
    "\n",
    "        # try:\n",
    "        # Repeat until convergence\n",
    "        for i_loop in range(self.n_loop):\n",
    "            self.i_count += 1\n",
    "\n",
    "            # Update\n",
    "            self.estimation_alpha()\n",
    "            self.estimation_beta()\n",
    "            self.estimation_gamma()\n",
    "\n",
    "            # Compute changes of parameters\n",
    "            diff_i = self.make_diff(w_new=self.result_alpha[self.i_count][self.dim_X_source1-1], w_old=self.result_alpha[self.i_count-1][self.dim_X_source1-1])\n",
    "            diff_a = self.make_diff(w_new=self.result_alpha[self.i_count][:(self.dim_X_source1-1)], w_old=self.result_alpha[self.i_count-1][:(self.dim_X_source1-1)])\n",
    "            diff_b = self.make_diff(w_new=self.result_beta[self.i_count], w_old=self.result_beta[self.i_count-1])\n",
    "            diff_c = self.make_diff(w_new=self.result_gamma[self.i_count], w_old=self.result_gamma[self.i_count-1])\n",
    "            diff = np.max([diff_i, diff_a, diff_b, diff_c]) # We use the maximum value of {diff_i, diff_a, diff_b, diff_c}.\n",
    "            # Store\n",
    "            self.diff_i[self.i_count] = diff_i\n",
    "            self.diff_a[self.i_count] = diff_a\n",
    "            self.diff_b[self.i_count] = diff_b\n",
    "            self.diff_c[self.i_count] = diff_c\n",
    "            self.diff[self.i_count] = diff\n",
    "\n",
    "            # Check the convergence\n",
    "            if diff < 1e-3:\n",
    "                self.convergence = True\n",
    "                break\n",
    "        # except:\n",
    "            # self.error = True\n",
    "        \n",
    "        # Cut off the unused parts of the dataframes\n",
    "        self.result_alpha = self.result_alpha[:(self.i_count+1),:]\n",
    "        self.result_beta = self.result_beta[:(self.i_count+1),:]\n",
    "        self.result_gamma = self.result_gamma[:(self.i_count+1),:]\n",
    "        self.diff_i = self.diff_i[:(self.i_count+1)]\n",
    "        self.diff_a = self.diff_a[:(self.i_count+1)]\n",
    "        self.diff_b = self.diff_b[:(self.i_count+1)]\n",
    "        self.diff_c = self.diff_c[:(self.i_count+1)]\n",
    "        self.diff = self.diff[:(self.i_count+1)]\n",
    "        \n",
    "        return self\n",
    "    \n",
    "    def predict(self, X):\n",
    "        \"\"\"\n",
    "        Prediction function\n",
    "            h(x) = \\alpha_1*ys + \\alpha_0 + (\\beta*ys + 1)<\\gamma, x>\n",
    "            \n",
    "        Returns\n",
    "        -------\n",
    "            pred1  : \\alpha_1*ys + \\alpha_0\n",
    "            pred2  : \\beta*ys + 1\n",
    "            pred3  : <\\gamma, x>\n",
    "            \n",
    "            y_pred : \\alpha_1*ys + \\alpha_0 + (\\beta*ys + 1)<\\gamma, x>\n",
    "        \"\"\"\n",
    "        # dataset\n",
    "        X_source_pred1 = X.iloc[:,dim_x:].copy()\n",
    "        X_source_pred1['Intercept'] = 1\n",
    "        X_source_pred2 = X.iloc[:,dim_x:].copy()\n",
    "        X_pred = X.iloc[:,:dim_x]\n",
    "        \n",
    "        # Compute each term\n",
    "        self.pred1 = X_source_pred1.dot(self.alpha)\n",
    "        self.pred2 = X_source_pred2.dot(self.beta) + 1\n",
    "        self.pred3 = X_pred.dot(self.gamma)\n",
    "        y_pred = self.pred1 - self.pred2*self.pred3\n",
    "        \n",
    "        return y_pred\n",
    "\n",
    "    def score(self, X, y=None):\n",
    "        \"\"\"\n",
    "        Score function for cross-validation\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "            -\\sum(y-\\hat{y})/n (Consider the minus value because 'GridSearchCV' maximize the score.)\n",
    "        \"\"\"\n",
    "        return -sum((y.values - self.predict(X).values)**2)/self.n_sample\n",
    "    \n",
    "    def get_params(self, deep=True):\n",
    "        \"\"\"\n",
    "        Create parameter dictionary for cross-validation\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "            {'lambda1', 'lambda2', 'lambda3', 'l1_ratio'}\n",
    "        \"\"\"\n",
    "        return {'lambda1' : self.lambda1,\n",
    "                'lambda2' : self.lambda2,\n",
    "                'lambda3_scale' : self.lambda3_scale,\n",
    "                'lambda3_trend' : self.lambda3_trend}\n",
    "    \n",
    "    def set_params(self, **parameters):\n",
    "        \"\"\"\n",
    "        For cross-validation\n",
    "        \"\"\"\n",
    "        for parameter, value in parameters.items():\n",
    "            setattr(self, parameter, value)\n",
    "        return self     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fix_seed function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def fix_seed(seed):\n",
    "    # Numpy\n",
    "    np.random.seed(seed)\n",
    "    # for built-in random\n",
    "    random.seed(seed)\n",
    "    # for hash seed\n",
    "    os.environ[\"PYTHONHASHSEED\"] = str(seed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Search parameters in cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# For difference-learning\n",
    "SearchParams_diff_FR = {\n",
    "    'lambda_scale' : 10**np.linspace(-2, 2, 25),\n",
    "    'lambda_trend' : [50, 100, 150]\n",
    "}\n",
    "\n",
    "# For proposed method\n",
    "SearchParams_InvTrans_FR = {\n",
    "    'lambda1' : [0],\n",
    "    'lambda2' : [1],\n",
    "    'lambda3_scale' : 10**np.linspace(-2, 2, 25),\n",
    "    'lambda3_trend' : [50, 100, 150]\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create output directories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.isdir('../30_Output/10_Model/300_TransferLearning'):\n",
    "    os.makedirs('../30_Output/10_Model/300_TransferLearning')\n",
    "if not os.path.isdir('../30_Output/20_Plot/300_TransferLearning'):\n",
    "    os.makedirs('../30_Output/20_Plot/300_TransferLearning')\n",
    "if not os.path.isdir('../30_Output/30_csv/300_TransferLearning'):\n",
    "    os.makedirs('../30_Output/30_csv/300_TransferLearning')\n",
    "if not os.path.isdir('../30_Output/40_pkl/300_TransferLearning'):\n",
    "    os.makedirs('../30_Output/40_pkl/300_TransferLearning')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main codes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('../30_Output/40_pkl/100_CheckData/100_Data.pkl', 'rb') as f:\n",
    "# with open('../30_Output/40_pkl/100_CheckData/100_Data_masked.pkl', 'rb') as f:\n",
    "    data_list = pickle.load(f)\n",
    "x = data_list['x']\n",
    "y = data_list['y']\n",
    "ys = data_list['ys']\n",
    "PID = data_list['PID']\n",
    "dim_x = x.shape[1]\n",
    "\n",
    "# Scaling parameter\n",
    "#     For stability of the estimation, scaling parameters are calculated using all data.\n",
    "x_Mean = x.mean()\n",
    "x_Std = x.std()\n",
    "y_LogMean = np.log(y).mean()\n",
    "y_LogStd = np.log(y).std()\n",
    "ys_LogMean = np.log(ys).mean()\n",
    "ys_LogStd = np.log(ys).std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User parameter setting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of iterations\n",
    "num_itr = 20\n",
    "# Number of fold of cross-validation\n",
    "n_fold = 5\n",
    "# seed for fix_seed function\n",
    "SEED = 373"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Repeat for different splittings of samples\n",
    "for i_itr in range(num_itr):\n",
    "    print('Itr : ' + str(i_itr))\n",
    "    \n",
    "    # Data splitting\n",
    "    PID_train, PID_test = train_test_split(PID, train_size=60, random_state=int(i_itr))\n",
    "    x_train  = x.loc[PID_train]\n",
    "    y_train  = y.loc[PID_train]\n",
    "    ys_train = ys.loc[PID_train]\n",
    "    x_test   = x.loc[PID_test]\n",
    "    y_test   = y.loc[PID_test]\n",
    "    ys_test  = ys.loc[PID_test]\n",
    "\n",
    "    # Data scaling\n",
    "    x_train_scal  = (x_train - x_Mean)/x_Std\n",
    "    x_test_scal   = (x_test - x_Mean)/x_Std\n",
    "    y_LogMean     = 0\n",
    "    y_LogStd      = 1\n",
    "    y_train_scal  = (np.log(y_train)-y_LogMean)/y_LogStd\n",
    "    ys_train_scal = (np.log(ys_train)-ys_LogMean)/ys_LogStd\n",
    "    ys_test_scal  = (np.log(ys_test)-ys_LogMean)/ys_LogStd\n",
    "\n",
    "    # Model training\n",
    "    ## Simple linear regression\n",
    "    t_tmp = time.time()\n",
    "    model_slr = LinearRegression(fit_intercept=True)\n",
    "    model_slr.fit(ys_train_scal.values.reshape(-1,1), y_train_scal.values)\n",
    "\n",
    "    y_fits_slr_scal = model_slr.predict(ys_train_scal.values.reshape(-1,1))\n",
    "    y_pred_slr_scal = model_slr.predict(ys_test_scal.values.reshape(-1,1))\n",
    "    y_fits_slr = y_fits_slr_scal*y_LogStd + y_LogMean\n",
    "    y_pred_slr = y_pred_slr_scal*y_LogStd + y_LogMean\n",
    "    print('  Simple linear regression has been done.    '+str(time.time()-t_tmp))\n",
    "\n",
    "    ## Learn the difference\n",
    "    y_train_diff      = np.log(y_train) - np.log(ys_train)\n",
    "    y_DiffMean        = 0\n",
    "    y_DiffStd         = 1\n",
    "    y_train_diff_scal = (y_train_diff - y_DiffMean)/y_DiffStd\n",
    "    \n",
    "    gsr_diff = GridSearchCV(\n",
    "        cls_FusedRidge(),\n",
    "        SearchParams_diff_FR,\n",
    "        cv      = n_fold,\n",
    "        n_jobs  = -1,\n",
    "        scoring = 'neg_mean_squared_error',\n",
    "        verbose = False\n",
    "    )\n",
    "    t_tmp = time.time()\n",
    "    fix_seed(SEED)\n",
    "    gsr_diff.fit(x_train_scal, y_train_diff_scal)\n",
    "\n",
    "    model_diff = cls_FusedRidge(\n",
    "        lambda_scale=gsr_diff.best_params_['lambda_scale'],\n",
    "        lambda_trend=gsr_diff.best_params_['lambda_trend']\n",
    "    )\n",
    "\n",
    "    fix_seed(373)\n",
    "    model_diff.fit(x_train_scal, y_train_diff_scal)\n",
    "\n",
    "    y_fits_diff_scal = model_diff.predict(x_train_scal)*y_DiffStd + y_DiffMean\n",
    "    y_pred_diff_scal = model_diff.predict(x_test_scal)*y_DiffStd + y_DiffMean\n",
    "    y_fits_diff      = y_fits_diff_scal+np.log(ys_train)\n",
    "    y_pred_diff      = y_pred_diff_scal+np.log(ys_test)\n",
    "    print('  Learn the difference has been done.    '+str(time.time()-t_tmp))\n",
    "\n",
    "    ## Proposed method\n",
    "    s_train_scal = pd.DataFrame(ys_train_scal)\n",
    "    s_test_scal  = pd.DataFrame(ys_test_scal)\n",
    "    x_train_adds = pd.merge(x_train_scal, s_train_scal, left_index=True, right_index=True)\n",
    "    x_test_adds  = pd.merge(x_test_scal, s_test_scal, left_index=True, right_index=True)\n",
    "    gsr_InvTrans = GridSearchCV(\n",
    "        cls_InvTrans_FR(),\n",
    "        SearchParams_InvTrans_FR,\n",
    "        cv      = n_fold,\n",
    "        n_jobs  = -1,\n",
    "        scoring = 'neg_mean_squared_error',\n",
    "        verbose = False\n",
    "    )\n",
    "    ini_intercept = model_slr.intercept_\n",
    "    ini_alpha     = model_slr.coef_[0]*ys_LogStd\n",
    "    ini_beta      = 0\n",
    "    fix_seed(373)\n",
    "    ini_gamma = model_diff.theta\n",
    "\n",
    "    t_tmp = time.time()\n",
    "    fix_seed(SEED)\n",
    "    gsr_InvTrans.fit(X=x_train_adds, y=y_train_scal)\n",
    "    model_InvTrans = cls_InvTrans_FR(\n",
    "        lambda1       = gsr_InvTrans.best_params_['lambda1'],\n",
    "        lambda2       = gsr_InvTrans.best_params_['lambda2'],\n",
    "        lambda3_scale = gsr_InvTrans.best_params_['lambda3_scale'],    \n",
    "        lambda3_trend = gsr_InvTrans.best_params_['lambda3_trend']\n",
    "    )\n",
    "    fix_seed(SEED)\n",
    "    model_InvTrans.fit(X=x_train_adds, y=y_train_scal)\n",
    "\n",
    "    y_fits_inv_scal = model_InvTrans.predict(x_train_adds)\n",
    "    y_pred_inv_scal = model_InvTrans.predict(x_test_adds)\n",
    "    y_fits_inv      = y_fits_inv_scal*y_LogStd + y_LogMean\n",
    "    y_pred_inv      = y_pred_inv_scal*y_LogStd + y_LogMean\n",
    "    print('  Proposed method has been done.    '+str(time.time()-t_tmp))\n",
    "\n",
    "    # Make figures\n",
    "#     fig = plt.figure(figsize=(15,5))\n",
    "\n",
    "    ## Simple linear regression\n",
    "#     y_prd1 = y_fits_slr\n",
    "#     y_obs1 = np.log(y_train).values\n",
    "#     y_prd2 = y_pred_slr\n",
    "#     y_obs2 = np.log(y_test).values\n",
    "#     ax = fig.add_subplot(1, 3, 1, \n",
    "#                          title='Simple linear regression', \n",
    "#                          xlabel='Prediction', \n",
    "#                          ylabel='Observation')\n",
    "#     ax.scatter(y_prd1, y_obs1, color='steelblue', alpha=0.7)\n",
    "#     ax.scatter(y_prd2, y_obs2, color='darkorange', alpha=1)\n",
    "#     xy_min = min(ax.get_xlim()[0], ax.get_ylim()[0])\n",
    "#     xy_max = max(ax.get_xlim()[1], ax.get_ylim()[1])\n",
    "#     ax.axis('equal')\n",
    "#     ax.axis('square')\n",
    "#     ax.set_xlim([xy_min, xy_max])\n",
    "#     ax.set_ylim([xy_min, xy_max])\n",
    "#     ax.grid(color='gray', linestyle='dotted', linewidth=1, alpha=0.5)\n",
    "#     ax.text(0.6, 0.14, 'Corr : '+str(round(np.corrcoef(y_prd2, y_obs2)[0,1], 4)), size=15, transform=ax.transAxes)\n",
    "#     ax.text(0.6, 0.08, 'MSE : '+str(round(np.mean((y_prd2-y_obs2)**2), 4)), size=15, transform=ax.transAxes)\n",
    "#     ax.text(0.6, 0.02, 'MAE : '+str(round(np.mean(np.abs(y_prd2-y_obs2)), 4)), size=15, transform=ax.transAxes)\n",
    "#     _ = ax.plot([-30000, 30000], [-30000, 30000], color='gray', linewidth=0.5)\n",
    "\n",
    "#     ## log-difference\n",
    "#     y_prd1 = y_fits_diff\n",
    "#     y_obs1 = np.log(y_train)\n",
    "#     y_prd2 = y_pred_diff\n",
    "#     y_obs2 = np.log(y_test)\n",
    "#     ax = fig.add_subplot(1, 3, 2, \n",
    "#                          title='Learning the log-difference', \n",
    "#                          xlabel='Prediction', \n",
    "#                          ylabel='Observation')\n",
    "#     ax.scatter(y_prd1, y_obs1, color='steelblue', alpha=0.7)\n",
    "#     ax.scatter(y_prd2, y_obs2, color='darkorange', alpha=1)\n",
    "#     xy_min = min(ax.get_xlim()[0], ax.get_ylim()[0])\n",
    "#     xy_max = max(ax.get_xlim()[1], ax.get_ylim()[1])\n",
    "#     ax.axis('equal')\n",
    "#     ax.axis('square')\n",
    "#     ax.set_xlim([xy_min, xy_max])\n",
    "#     ax.set_ylim([xy_min, xy_max])\n",
    "#     ax.grid(color='gray', linestyle='dotted', linewidth=1, alpha=0.5)\n",
    "#     ax.text(0.6, 0.14, 'Corr : '+str(round(np.corrcoef(y_prd2, y_obs2)[0,1], 4)), size=15, transform=ax.transAxes)\n",
    "#     ax.text(0.6, 0.08, 'MSE : '+str(round(np.mean((y_prd2-y_obs2)**2), 4)), size=15, transform=ax.transAxes)\n",
    "#     ax.text(0.6, 0.02, 'MAE : '+str(round(np.mean(np.abs(y_prd2-y_obs2)), 4)), size=15, transform=ax.transAxes)\n",
    "#     _ = ax.plot([-30000, 30000], [-30000, 30000], color='gray', linewidth=0.5)\n",
    "\n",
    "#     ## Proposed method\n",
    "#     y_prd1 = y_fits_inv.values\n",
    "#     y_obs1 = np.log(y_train).values\n",
    "#     y_prd2 = y_pred_inv.values\n",
    "#     y_obs2 = np.log(y_test).values\n",
    "#     ax = fig.add_subplot(1, 3, 3, \n",
    "#                          title='Proposed method', \n",
    "#                          xlabel='Prediction', \n",
    "#                          ylabel='Observation')\n",
    "#     ax.scatter(y_prd1, y_obs1, color='steelblue', alpha=0.7)\n",
    "#     ax.scatter(y_prd2, y_obs2, color='darkorange', alpha=1)\n",
    "#     xy_min = min(ax.get_xlim()[0], ax.get_ylim()[0])\n",
    "#     xy_max = max(ax.get_xlim()[1], ax.get_ylim()[1])\n",
    "#     ax.axis('equal')\n",
    "#     ax.axis('square')\n",
    "#     ax.set_xlim([xy_min, xy_max])\n",
    "#     ax.set_ylim([xy_min, xy_max])\n",
    "#     ax.grid(color='gray', linestyle='dotted', linewidth=1, alpha=0.5)\n",
    "#     ax.text(0.6, 0.14, 'Corr : '+str(round(np.corrcoef(y_prd2, y_obs2)[0,1], 4)), size=15, transform=ax.transAxes)\n",
    "#     ax.text(0.6, 0.08, 'MSE : '+str(round(np.mean((y_prd2-y_obs2)**2), 4)), size=15, transform=ax.transAxes)\n",
    "#     ax.text(0.6, 0.02, 'MAE : '+str(round(np.mean(np.abs(y_prd2-y_obs2)), 4)), size=15, transform=ax.transAxes)\n",
    "#     _ = ax.plot([-30000, 30000], [-30000, 30000], color='gray', linewidth=0.5)\n",
    "\n",
    "#     plt.tight_layout()\n",
    "#     plt.savefig('../30_Output/20_Plot/300_TransferLearning/300_Plot_'+str(i_itr)+'.png')\n",
    "#     plt.close()\n",
    "\n",
    "    # Save results\n",
    "    result_list = {\n",
    "        'x_train'     : x_train,\n",
    "        'y_train'     : y_train,\n",
    "        'ys_train'    : ys_train,\n",
    "        'x_test'      : x_test,\n",
    "        'y_test'      : y_test,\n",
    "        'ys_test'     : ys_test,\n",
    "        'PID_train'   : PID_train,\n",
    "        'PID_test'    : PID_test,\n",
    "        'y_fits_slr'  : y_fits_slr,\n",
    "        'y_pred_slr'  : y_pred_slr,\n",
    "        'model_slr'   : model_slr,\n",
    "        'y_fits_diff' : y_fits_diff,\n",
    "        'y_pred_diff' : y_pred_diff,\n",
    "        'model_diff'  : model_diff,\n",
    "        'y_fits_inv'  : y_fits_inv,\n",
    "        'y_pred_inv'  : y_pred_inv,\n",
    "        'model_inv'   : model_InvTrans\n",
    "    }\n",
    "\n",
    "    f = open('../30_Output/40_pkl/300_TransferLearning/300_Results_'+str(i_itr)+'.pkl','wb')\n",
    "    pickle.dump(result_list,f)\n",
    "    f.close()"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
