calculate_results_final.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74dd14c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from findiff import FinDiff, coefficients\n",
    "from scipy.linalg import solve, det, expm\n",
    "from math import isclose\n",
    "import pandas as pd\n",
    "from scipy.signal import savgol_filter\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import torch.optim as optim\n",
    "from torch.utils.data import DataLoader, TensorDataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60e2c3c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "savefigs_loc = \"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95885f1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "savefigs_loc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c40ffd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "savefigs_loc_nn = \"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78158e35",
   "metadata": {},
   "outputs": [],
   "source": [
    "######### NEURAL NETWORK MODEL ###########\n",
    "\n",
    "#Data as numpy array\n",
    "# t = np.linspace(0,np.pi,50,endpoint=True)\n",
    "# y = np.sin(t)\n",
    "\n",
    "def nn_model(t, time_series, hidden_width, epochs, title, save_fig=False):\n",
    "    print(torch.cuda.is_available())\n",
    "    print(torch.cuda.device_count())\n",
    "    print(torch.cuda.current_device())\n",
    "    print(torch.cuda.device(0))\n",
    "    print(torch.cuda.get_device_name(0))\n",
    "    \n",
    "    #Our input dimension is 1, and each entry in the tensor above is a datapoint for training. So, the length of \n",
    "    #the tensor/numpy array should be the batch size:\n",
    "    n_input, n_hidden, n_output = 1, hidden_width, 1\n",
    "    batch_size = 10\n",
    "    \n",
    "    #Turn into tensor\n",
    "    t_tensor = torch.from_numpy(np.reshape(t, (-1, n_input))).to(torch.float32)\n",
    "    y_tensor = torch.from_numpy(np.reshape(time_series, (-1, n_output))).to(torch.float32)\n",
    "\n",
    "    dataset = TensorDataset(t_tensor, y_tensor)\n",
    "    dataloader = DataLoader(dataset, batch_size=batch_size, shuffle=True)\n",
    "    \n",
    "    dataset_size = len(dataloader.dataset)\n",
    "    \n",
    "    learning_rate = 0.001\n",
    "    \n",
    "    plt.ylabel('y')\n",
    "    plt.xlabel('t')\n",
    "    plt.title(\"Given Time Series\")\n",
    "    plt.plot(t_tensor, y_tensor, marker='.', linestyle='none')\n",
    "    plt.show()\n",
    "\n",
    "    model = nn.Sequential(nn.Linear(n_input, n_hidden), nn.ReLU(), nn.Linear(n_hidden, n_output))\n",
    "    print(model)\n",
    "\n",
    "    loss_function = nn.MSELoss(reduction=\"sum\")\n",
    "    print(loss_function)\n",
    "    \n",
    "    optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)\n",
    "\n",
    "    losses = [] # For plotting purposes\n",
    "    print(\"---------------------------------------------------------\")\n",
    "    print(\"START TRAINING\")\n",
    "    print(f\"HIDDEN LAYER WIDTH: {hidden_width}\")\n",
    "    for epoch in range(epochs):\n",
    "        #print(f\"Epoch {epoch}\\n------------------------\")\n",
    "        \n",
    "        losses_epoch = []\n",
    "        for id_batch, (t_batch, y_batch) in enumerate(dataloader):\n",
    "            #print(\"id_batch: \", id_batch)\n",
    "            #print(\"t_batch: \\n\", t_batch)\n",
    "            #print(\"y_batch: \\n\", y_batch)\n",
    "            \n",
    "            # Forward pass\n",
    "            y_pred = model(t_batch)\n",
    "            \n",
    "            # Calculate the loss\n",
    "            loss = loss_function(y_pred, y_batch)\n",
    "            losses_epoch.append(loss.item())\n",
    "            \n",
    "            # We set the gradients to zero, as PyTorch accumulates gradients on subsequent backward passes.\n",
    "            optimizer.zero_grad()\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "            \n",
    "            # We plot the prediction every 10 epochs.\n",
    "            '''\n",
    "            if epoch%25 == 0:\n",
    "                loss, current = loss.item(), (id_batch + 1)* len(t_batch)\n",
    "                print(f\"loss: {loss:>7f}\")\n",
    "                \n",
    "                with torch.no_grad():\n",
    "                    plt.ylabel('y')\n",
    "                    plt.xlabel('t')\n",
    "                    plt.title(f\"True Time Series VS. Prediction - Iteration {epoch}\")\n",
    "                    plt.plot(t_tensor, y_tensor, marker='.', linestyle='none')\n",
    "                    plt.plot(t_batch, y_pred, marker='.', linestyle='none', color='orange')\n",
    "                    plt.legend([\"True\", \"NN Prediction\"])\n",
    "                    plt.show()\n",
    "            '''\n",
    "            \n",
    "            #model.zero_grad() \n",
    "            # Calculate the derivatives\n",
    "            #loss.backward()\n",
    "\n",
    "            # Update the weights\n",
    "            #learning_rate = 0.01\n",
    "            #for f in model.parameters():\n",
    "            #    f.data.sub_(f.grad.data * learning_rate)\n",
    "        losses.append(np.mean(losses_epoch))\n",
    "        #print(f\"loss: {np.mean(losses_epoch):>7f}\")\n",
    "    print(\"END TRAINING\")\n",
    "    \n",
    "    y_pred_final = model(t_tensor)\n",
    "    rel_error_final = relative_error(y_pred_final.detach().numpy(), y_tensor.detach().numpy())\n",
    "    \n",
    "    with torch.no_grad():\n",
    "        #plt.ylabel('y')\n",
    "        #plt.xlabel('time')\n",
    "        plt.title(f\"{title}\")\n",
    "        \n",
    "        #plt.plot(t_tensor.detach().numpy(), el_price_2018_data.to_numpy(), color=\"gray\",alpha=0.5, label=\"Original\")\n",
    "        \n",
    "        plt.plot(t_tensor, y_tensor, label=\"True\") #, marker='.', linestyle='none')\n",
    "        plt.plot(t_tensor, y_pred_final, label=f\"NN Pred. Width={hidden_width}\") #, marker='.', linestyle='none', color='orange')\n",
    "        plt.legend()\n",
    "        if save_fig:\n",
    "            plt.savefig(savefigs_loc_nn+\"_TRAINING_\"+title+f\"_hw_{hidden_width}__epochs_{epochs}__rel_error_{rel_error_final:.4}_.png\")\n",
    "        plt.show()\n",
    "    \n",
    "    plt.plot(losses)\n",
    "    plt.ylabel('loss')\n",
    "    plt.xlabel('epoch')\n",
    "    plt.title(title+f\"_width={hidden_width}_\"+\"lr = %f\"%(learning_rate))\n",
    "    if save_fig:\n",
    "        plt.savefig(savefigs_loc_nn+\"_LOSS_\"+title+f\"_hw_{hidden_width}__epochs_{epochs}__rel_error_{rel_error_final:.4}_.png\")\n",
    "    plt.show()\n",
    "    \n",
    "    print(f\"Relative error = {rel_error_final:.4}\")\n",
    "    print(\"--------------------------------------\")\n",
    "    return model, rel_error_final, y_pred_final, y_tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68c98701",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_models(hidden_widths, t, time_series, epochs, title, save_fig=False):\n",
    "    models_dict = {}\n",
    "    rel_err_dict = {}\n",
    "    final_preds_dict = {}\n",
    "    ts_tensor_dict={}\n",
    "    \n",
    "    for hw in hidden_widths:\n",
    "        mod, rel_err, y_pred_final, ts_tensor = nn_model(t, time_series, hw, epochs, title, save_fig)\n",
    "        models_dict[hw]=mod\n",
    "        rel_err_dict[hw]=rel_err\n",
    "        \n",
    "        final_preds_dict[hw] = y_pred_final\n",
    "        ts_tensor_dict[hw] = ts_tensor\n",
    "    \n",
    "    #plt.plot(t, el_price_2018_data.to_numpy(), color=\"gray\",alpha=0.5, label=\"Original\")\n",
    "    plt.plot(t, ts_tensor_dict[hidden_widths[0]].detach().numpy(), label=\"True\")\n",
    "    for hw in hidden_widths:\n",
    "        plt.plot(t, final_preds_dict[hw].detach().numpy(), label=f\"NN pred. Width={hw}\")\n",
    "    plt.legend()\n",
    "    plt.title(title)\n",
    "    if save_fig:\n",
    "        plt.savefig(savefigs_loc_nn+\"_TRAINING_ALL_HWs__\"+title+f\"__epochs_{epochs}_.png\")\n",
    "    plt.show()\n",
    "    \n",
    "    \n",
    "    x = y = list(range(hidden_widths[-1]))\n",
    "    idx = rel_err_dict.keys()\n",
    "    new_y = [rel_err_dict[i] for i in idx]\n",
    "    plt.plot(range(len(idx)), new_y, 'o-')\n",
    "    plt.xticks(range(len(idx)),idx)\n",
    "    plt.xlabel(\"Hidden layer width\")\n",
    "    plt.ylabel(\"Training, Relative error\")\n",
    "    plt.legend()\n",
    "    if save_fig:\n",
    "        plt.savefig(savefigs_loc_nn+\"_REL_ERRORS_vs_HWs__\"+title+f\"__epochs_{epochs}_.png\")\n",
    "    plt.show()\n",
    "    \n",
    "    best_width = min(rel_err_dict, key=rel_err_dict.get)\n",
    "    print(\"Best width: \", best_width)\n",
    "    return models_dict, best_width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91f42abc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_test(model,hidden_width, test_time, test_time_series, title, forward_time=False, save_fig=False):\n",
    "    n_input, n_output = 1, 1\n",
    "    test_time_tensor = torch.from_numpy(np.reshape(test_time, (-1, n_input))).to(torch.float32)\n",
    "    test_time_series_tensor = torch.from_numpy(np.reshape(test_time_series, (-1, n_output))).to(torch.float32)\n",
    "    \n",
    "    pred = model(test_time_tensor)\n",
    "    \n",
    "    plt.title(title)\n",
    "    \n",
    "    if forward_time:\n",
    "        middle_index = test_time.shape[0]//2\n",
    "        forward_rel_err = relative_error(pred.detach().numpy()[middle_index:], test_time_series_tensor.detach().numpy()[middle_index:])\n",
    "        \n",
    "        #plt.plot(test_time_tensor.detach().numpy(), el_price_2018_2019_data.to_numpy()[:len(test_time_tensor.detach().numpy())], color=\"gray\",alpha=0.5, label=\"Original\")\n",
    "        \n",
    "        plt.plot(test_time_tensor.detach().numpy(), test_time_series_tensor.detach().numpy(), label=\"True\")\n",
    "        plt.plot(test_time_tensor.detach().numpy(), pred.detach().numpy(), label=f\"Pred. Width={hidden_width}\")\n",
    "        plt.legend()\n",
    "        if save_fig:\n",
    "            plt.savefig(savefigs_loc_nn+\"_FORWARD_\"+title+f\"_hw_{hidden_width}__forward_rel_error_{forward_rel_err:.4}_.png\")\n",
    "        \n",
    "        print(\"Forward, relative error (on future time points): \", forward_rel_err)\n",
    "    else:\n",
    "        test_rel_err = relative_error(pred.detach().numpy(), test_time_series_tensor.detach().numpy())\n",
    "        \n",
    "        plt.plot(test_time_tensor.detach().numpy(), tesla_2021_data.to_numpy(), color=\"gray\",alpha=0.5, label=\"Original\")\n",
    "        \n",
    "        plt.plot(test_time_tensor.detach().numpy(), test_time_series_tensor.detach().numpy(), label=\"Smoothed\")\n",
    "        plt.plot(test_time_tensor.detach().numpy(), pred.detach().numpy(), label=f\"Pred. Width={hidden_width}\")\n",
    "        plt.legend()\n",
    "        if save_fig:\n",
    "            plt.savefig(savefigs_loc_nn+\"_TESTING_\"+title+f\"_hw_{hidden_width}__forward_rel_error_{test_rel_err:.4}_.png\")\n",
    "\n",
    "        \n",
    "        print(\"Testing (same interval), relative error: \", test_rel_err)\n",
    "    \n",
    "    plt.show()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c4b78af",
   "metadata": {},
   "outputs": [],
   "source": [
    "############ ODE MODEL #################\n",
    "\n",
    "### FOR FINITE DIFFERENCE ###\n",
    "\n",
    "def floor_even(num):\n",
    "    return np.floor(num/2)*2\n",
    "\n",
    "def floor_odd(num):\n",
    "    if num%2==1:\n",
    "        return num\n",
    "    else:\n",
    "        return num-1\n",
    "\n",
    "def highest_possible_N(time_series):\n",
    "    number = 2/3 * len(time_series) -2/3\n",
    "\n",
    "    # Due to floating point errors causing np.floor to return the wrong number,\n",
    "    # (for instance number=65.9999 when it is actually =66), we use isclose \n",
    "    if isclose(number,np.ceil(number), abs_tol=1e-8):\n",
    "        return int(round(number))\n",
    "    else:\n",
    "        return int(np.floor(number))\n",
    "\n",
    "def highest_possible_accuracies(time_series, N):\n",
    "    acc_list = []\n",
    "    for deriv in range(1,N+1):\n",
    "        # Different conditions for odd and even derivatives:\n",
    "        if deriv%2==1:\n",
    "            number = 2/3*(len(time_series))+1-deriv\n",
    "            # Handle floating point error:\n",
    "            if isclose(number,np.ceil(number), abs_tol=1e-8):\n",
    "                actual_number=round(number)\n",
    "                acc_list.append(int(floor_even(actual_number)))\n",
    "            else:\n",
    "                acc_list.append(int(floor_even(number)))\n",
    "        else:\n",
    "            number = 2/3*(len(time_series))+4/3-deriv\n",
    "            # Handle floating point error:\n",
    "            if isclose(number,np.ceil(number), abs_tol=1e-8):\n",
    "                actual_number=round(number)\n",
    "                acc_list.append(int(floor_even(actual_number)))\n",
    "            else:\n",
    "                acc_list.append(int(floor_even(number)))\n",
    "                \n",
    "    return acc_list\n",
    "\n",
    "\n",
    "def points_to_accuracy(time_series, p_to_use, N):\n",
    "    print(\"p_to_use: \", p_to_use)\n",
    "    # n+m points used in a findiff approx (order + accuracy)\n",
    "    # Accuracy has to be an even number.\n",
    "    \n",
    "    accuracy_list = []\n",
    "    \n",
    "    if p_to_use == 'few':\n",
    "        accuracy_list = [2]*N\n",
    "        \n",
    "    if p_to_use == 'all':\n",
    "        accuracy_list = highest_possible_accuracies(time_series, N)\n",
    "        \n",
    "    if '%' in p_to_use:\n",
    "        percent = float(p_to_use.split(\"%\")[0])*0.01\n",
    "        ts_len = len(time_series)\n",
    "        num_points = int(np.round(ts_len*percent))\n",
    "        \n",
    "        for n in range(1,N+1):\n",
    "            if n>num_points-2:\n",
    "                m = 2\n",
    "                accuracy_list.append(m)\n",
    "            elif n%2==1:\n",
    "                m = int(floor_even(num_points-n))\n",
    "                while ((n + m - 1)/2 - 1 + n + m) > ts_len:\n",
    "                    m -= 2\n",
    "                accuracy_list.append(m)\n",
    "            elif n%2==0:     \n",
    "                m = int(floor_even(num_points-n))\n",
    "                while ((n + m - 1)/2 - 1 - 1 + n + m) > ts_len:\n",
    "                    m -= 2\n",
    "                accuracy_list.append(m)    \n",
    "    return accuracy_list\n",
    "\n",
    "def derivative_approximations(h, time_series, N, accuracy_list):\n",
    "    fin_diff = [time_series] # the 0-th derivative is the original time series\n",
    "    #np.set_printoptions(precision=2)\n",
    "    for k in range(1,N+1):\n",
    "        derivative = FinDiff(0, h, k, acc=accuracy_list[k-1])\n",
    "        fin_diff.append(derivative(time_series))\n",
    "        \n",
    "        mat_repr_deriv_approx = derivative.matrix(time_series.shape)\n",
    "        #print(f\"FINITE DIFF. MATRIX REPRESENTATION. deriv={k}, acc={accuracy_list[k-1]}\")\n",
    "        #print(mat_repr_deriv_approx.toarray())\n",
    "        #print(\"\")\n",
    "        \n",
    "        \n",
    "    return np.asarray(fin_diff)\n",
    "\n",
    "### FOR DETERMINING COEFFICIENTS ANALYTICALLY ###\n",
    "\n",
    "def coeff_analytical_matrix_vector(fin_diff, coeff_to_lock):\n",
    "    N = np.shape(fin_diff)[0]-1 # Highest order derivative\n",
    "        \n",
    "    # We will return A and b, as we later wish to solve the matrix equation Ax=b, where the solution x \n",
    "    # will contain all the coefficients except the locked one, i.e. c_k: [a, c_0, ..., c_k-1, c_k+1, ..., c_N]:\n",
    "    A = np.zeros((N+1,N+1)) # LHS MATRIX\n",
    "    b = np.zeros((N+1, 1))  # RHS vector\n",
    "        \n",
    "        \n",
    "    # Two cases. If coeff_to_lock is 'a':\n",
    "    if coeff_to_lock == 'a':\n",
    "        temp_matrix = np.copy(fin_diff)\n",
    "        for m in range(N+1):\n",
    "            b[m,0] = -1*np.sum(temp_matrix[m])\n",
    "            A[m,m] = np.sum(temp_matrix[m]*temp_matrix[m]) \n",
    "            for l in range(m+1, N+1):\n",
    "                A_ml = np.sum(temp_matrix[l]*temp_matrix[m])\n",
    "                A[m,l] = A_ml\n",
    "                A[l,m] = A_ml \n",
    "        \n",
    "        # In this case, the final solution x to Ax=b, will look like x=[c_0, c_1, ..., c_k, ..., c_N]\n",
    "        \n",
    "    # coeff_to_lock is c_k:\n",
    "    else:\n",
    "        coeff_to_lock_index = int(coeff_to_lock.split(\"_\")[1]) # e.g. c_0 corresponds to row 0 in fin_diff\n",
    "        \n",
    "        # Code is nicer if we replace row corresponding to locked \n",
    "        # coefficient with ones\n",
    "        dyk_i = np.copy(fin_diff[coeff_to_lock_index])\n",
    "        temp_matrix = np.delete(fin_diff, coeff_to_lock_index, axis=0)\n",
    "        temp_matrix = np.insert(temp_matrix, 0, np.ones_like(dyk_i), axis=0)\n",
    "\n",
    "        for m in range(N+1):\n",
    "            b[m,0] = -1*np.sum(dyk_i*temp_matrix[m])\n",
    "            A[m,m] = np.sum(temp_matrix[m]*temp_matrix[m]) \n",
    "            for l in range(m+1, N+1):\n",
    "                A_ml = np.sum(temp_matrix[l]*temp_matrix[m])\n",
    "                A[m,l] = A_ml\n",
    "                A[l,m] = A_ml \n",
    "                \n",
    "        # In this case, the final solution x to Ax=b will look like x=[a, c_0, ..., c_k-1, c_k+1, ..., c_N]\n",
    "    \n",
    "    return A, b\n",
    "\n",
    "'''\n",
    "# Prevent singular matrix\n",
    "def check_if_deriv_is_zero(fin_diff, tolerance=1e-6):\n",
    "    # We return the index of the first derivative which is practically zero (we know all the following ones will also be zero).\n",
    "    sum_of_rows = np.sum(abs(fin_diff), axis=1)\n",
    "    print(\"Sum of each row in fin_diff: \\n\", sum_of_rows)\n",
    "    for row_number, row_sum in enumerate(sum_of_rows):\n",
    "        if abs(row_sum) <= tolerance:\n",
    "            return row_number\n",
    "    return 'All derivatives nonzero'\n",
    "\n",
    "def prev_sing_matrix(fin_diff, coefficient_to_lock):\n",
    "    ### PREVENT SINGULAR MATRIX: This can happen if the time series is sampled from a polynomial with degree < N.\n",
    "    # If any of the derivative approximations is equal to zero, this will give us a singular matrix below.\n",
    "    # If the derivs after a certain order are zero, we set the coefficients for these derivs to zero.\n",
    "    # Then we only work with fin_diff up to and including the final nonzero deriv approx.   \n",
    "    print(\"PREVENT SINGULAR MATRIX\")\n",
    "    print(\"Check if any of the derivatives are 0: \\n\")\n",
    "    first_index_deriv_zero = check_if_deriv_is_zero(fin_diff, tolerance=1e-6) # Need a tolerance due to floating point error.\n",
    "    print(\"\")\n",
    "    if first_index_deriv_zero!='All derivatives nonzero':\n",
    "        if first_index_deriv_zero == 0:\n",
    "            print(\"The time series is equal to zero.\")\n",
    "            return\n",
    "        \n",
    "        if coefficient_to_lock != 'a' and int(coefficient_to_lock.split(\"_\")[1])>=first_index_deriv_zero:\n",
    "            print(\"The selected coefficient to lock corresponds to a derivative that is equal to zero.\")\n",
    "            print(f\"Select a coefficient c_k to lock with k <= {first_index_deriv_zero-1}, as all derivatives of order\")\n",
    "            print(\"higher than this are equal to zero.\")\n",
    "            return\n",
    "        \n",
    "        fin_diff = fin_diff[:first_index_deriv_zero]\n",
    "        \n",
    "        print(f\"The highest order derivative to include was given as: {N}.\")\n",
    "        print(f\"However, the derivatives of order >= {first_index_deriv_zero} are equal to zero.\")\n",
    "        print(\"The returned list of coefficients will therefore only include the coefficients up to and including:\") \n",
    "        print(f\"c_{first_index_deriv_zero-1} \\n\")\n",
    "    else:\n",
    "        print(\"findiff does not contain rows that are equal to zero.\")\n",
    "    return fin_diff\n",
    "\n",
    "# Prevent linearly dependent matrix\n",
    "def determine_lin_dep_columns(A, tolerance=1e-6):\n",
    "    low = 0\n",
    "    high = A.shape[1]-1\n",
    "    mid = 0\n",
    "    \n",
    "    # We do a binary search until we get the index of the final column to include, whcih corresponds with the \"high\" variable\n",
    "    if abs(det(A))<tolerance:\n",
    "        while low <= high:\n",
    "            mid = (high+low)//2\n",
    "            \n",
    "            if abs(det(A[:mid, :mid]))<tolerance:\n",
    "                high = mid-1\n",
    "            else:\n",
    "                low = mid+1\n",
    "                \n",
    "        return high\n",
    "    else:\n",
    "        return (\"Matrix is linearly independent\")\n",
    "    \n",
    "def prev_lin_dep_matrix(A, b):\n",
    "    ### PREVENT LINEARLY DEPENDENT MATRIX: This can happen if a row in fin_diff is a multiple of \n",
    "    # an earlier row in fin_diff. I.e. what happens when we have a time series sampled from exp(t),\n",
    "    # then every dervative is the same. We solve this by finding the largest submatrix in A s.t.\n",
    "    # A is linearly independent. We do this by checking if the determinant is smaller or bigger\n",
    "    # than some tolerance.\n",
    "    print(\"CHECK IF MATRIX IS LINEARLY INDEPENDENT\")\n",
    "    final_col_to_include = determine_lin_dep_columns(np.copy(A), tolerance=1e-4)\n",
    "    if final_col_to_include != \"Matrix is linearly independent\":\n",
    "        print(\"The matrix was not linearly independent.\")\n",
    "        print(\"We extract the largest submatrix starting in the upper left corner that's linearly independent.\")\n",
    "        print(f\"The final column from A we include is the column with index: {final_col_to_include-1}.\")\n",
    "        print(f\"This gives the matrix A[:{final_col_to_include},:{final_col_to_include}]:\")\n",
    "        \n",
    "        A = A[:final_col_to_include, :final_col_to_include]\n",
    "        b = b[:final_col_to_include]\n",
    "        \n",
    "        print(\"A: \\n\", A)\n",
    "        print(\"b: \\n\", b)\n",
    "        \n",
    "        print(f\"Since the matrix {final_col_to_include} rows/columns corresponds to having originally\")\n",
    "        print(f\"chosen N={final_col_to_include-1} instead, the returned list of coefficients\")\n",
    "        print(f\"will therefore only include coefficients up to and including: c_{final_col_to_include-1}\")\n",
    "    else:\n",
    "        print(final_col_to_include)\n",
    "    return A, b\n",
    "'''\n",
    "\n",
    "\n",
    "def insert_locked_coeff(solved_coeff, coefficient_to_lock):\n",
    "    if coefficient_to_lock == 'a':\n",
    "        solved_all_coeff = np.insert(solved_coeff, 0, 1, axis=0) \n",
    "        # Here, solved_all_coeff looks like [1, c_0, ..., c_k, ..., c_N]\n",
    "    else: \n",
    "        coeff_to_lock_index = int(coefficient_to_lock.split(\"_\")[1])+1 # +1 since 'a' is at index 0.\n",
    "        solved_all_coeff = np.insert(solved_coeff, coeff_to_lock_index, 1, axis=0)\n",
    "        # Here, if c_k is the locked coeff, solved_all_coeff looks like [a, c_0, ..., c_k-1, 1, c_k+1, ..., c_N].\n",
    "    return solved_all_coeff\n",
    "        \n",
    "# Find the coefficients solving the minimization problem\n",
    "def solved_all_coeffs_analytical(time_points, time_series, N, acc_p_to_use, coefficient_to_lock='c_0', accuracy_list=None):\n",
    "    time_points = np.array(time_points)\n",
    "    time_series = np.array(time_series)\n",
    "    \n",
    "    ### Check if chosen N is possible:\n",
    "    highest_pos_N = highest_possible_N(time_series)\n",
    "    if N<=highest_pos_N:\n",
    "        print(f\"The chosen highest order derivative is N={N}, which is less than or equal\")\n",
    "        print(f\"to the highest possible derivative, N_highest = {highest_pos_N}\")\n",
    "    else:\n",
    "        print(f\"The chosen highest order derivative is N={N}, which is greater than\")\n",
    "        print(f\"the highest possible derivative, N_highest = {highest_pos_N}\")\n",
    "        return\n",
    "    \n",
    "    # Check if coefficient to lock is valid\n",
    "    if coefficient_to_lock != 'a' and coefficient_to_lock[:2]!='c_':\n",
    "        print(\"coefficient_to_lock must be either 'a' or 'c_k', for 0 <= k <= N\")\n",
    "        return\n",
    "    \n",
    "    if coefficient_to_lock[:2]=='c_' and (int(coefficient_to_lock.split(\"_\")[1])>N or int(coefficient_to_lock.split(\"_\")[1])<0):\n",
    "        print(\"c_k: k must be between 0 and N\")\n",
    "        return\n",
    "    \n",
    "    ### Create a list with the highest possible accuracy for each derivative approx:\n",
    "    if accuracy_list == None:\n",
    "        accuracy_list = points_to_accuracy(time_series, acc_p_to_use, N)\n",
    "    \n",
    "    print(\"\\n List containing order of accuracy for the derivative approximations.\")\n",
    "    print(\"Index i corresponds to the accuracy of the derivative of order i+1: \\n\", accuracy_list, \"\\n\")\n",
    "    \n",
    "    ### Create an array with all derivative approximations. \n",
    "    ### The row corresponds to the derivative order, and column corresponds to the time point.\n",
    "    ### Ex: value at (2, 1), is the approximation of the second derivative of y at time t_1, that is y_1''.\n",
    "    fin_diff = derivative_approximations(time_points[1]-time_points[0], time_series, N, accuracy_list)  \n",
    "    print(\"\\n fin_diff: \\n\", fin_diff, \"\\n\")\n",
    "    ###  fin_diff - numpy array with  shape (N+1, n+1): \n",
    "    ###    [ [y_0     ,   y_1      ,   ...   ,    y_n     ],\n",
    "    ###      [y_0'    ,   y_1'     ,   ...   ,    y_n'    ],\n",
    "    ###      [y_0''   ,   y_1''    ,   ...   ,    y_n''   ],\n",
    "    ###                     ...,\n",
    "    ###      [y_0^(N) ,   y_1^(N)  ,   ...   ,    y_n^(N) ]  ] \n",
    "    \n",
    "    \n",
    "    #fin_diff = prev_sing_matrix(fin_diff, coefficient_to_lock)\n",
    "    \n",
    "    print(\"CONSTRUCT MATRIX FOR MATRIX EQUATION\")\n",
    "    A, b = coeff_analytical_matrix_vector(fin_diff, coefficient_to_lock) \n",
    "    print(\"A: \\n\", A)\n",
    "    print(\"b: \\n\", b)\n",
    "    print(\"\")    \n",
    "    #A, b = prev_lin_dep_matrix(A,b)\n",
    "    print(\"det of A: \", det(A))\n",
    "    \n",
    "    solved_coeff = solve(A, b, assume_a='sym') # Column vector. We use the convention [\"a\",\"c_0\", \"c_1\", \"c_2\", ..., \"c_N\"]^T\n",
    "    \n",
    "    # The solution to the matrix equation Ax = b does not contain the locked coefficient, so we have to add it back.\n",
    "    # Since this is the coefficient that's \"divided out\", it has value 1:\n",
    "    solved_all_coeff = insert_locked_coeff(solved_coeff, coefficient_to_lock)\n",
    "    \n",
    "    return solved_all_coeff, fin_diff\n",
    "\n",
    "### FOR CALCULATING THE EXPLICIT SOLUTION\n",
    "\n",
    "# Find the matrix and vector for converting ODE to system of first order equations\n",
    "def coeff_list_to_system_matrix(all_coeffs, tolerance):\n",
    "    # all_coeffs = np.array([[a], [1], [c_1], [c_2], ..., [c_N]]) where c_0 is the locked coeff.\n",
    "    \n",
    "    # As we will divide by c_N, we must ensure that the last coeff in the list is nonzero.\n",
    "    # This is OK, becaue if the array ends with c_{k+1}, ..., c_{N} that are (approx) equal to zero\n",
    "    # we in practice have an ODE of order k, so we can just drop these last terms. \n",
    "    while abs(all_coeffs[-1,0])<tolerance:\n",
    "        all_coeffs = all_coeffs[:-1]\n",
    "    \n",
    "    print(\"remove highest order coeff(s) if below threshold: \\n\", all_coeffs)\n",
    "    \n",
    "    # The order of the ODE determines the shape of the matrix.\n",
    "    ode_order = all_coeffs.shape[0]-2 # -2 because we have the constant term and the y^(0) coeff.\n",
    "    print(\"If any coefficients were removed, the effective ODE order will change. The ODE order is now: \", ode_order)\n",
    "    if ode_order == 0:\n",
    "        C = None\n",
    "        A = None\n",
    "        return C, A\n",
    "        \n",
    "    # Create matrix C\n",
    "    C = -1*np.eye(ode_order, k=1)\n",
    "    bottom_row = all_coeffs[1:-1,0].T/all_coeffs[-1,0]\n",
    "    print(\"bottom_row: \", bottom_row)\n",
    "    C[-1, :] = bottom_row \n",
    "    print(\"C: \\n\", C)\n",
    "    # Create vector A\n",
    "    A = np.zeros((ode_order,1))\n",
    "    A[-1,0]=all_coeffs[0,0]\n",
    "    print(\"A: \\n\", A)\n",
    "    return C, A\n",
    "\n",
    "def initial_cond(fin_diff, percentage_of_interval, where_to_start, N):\n",
    "    num_points = floor_odd(int(np.floor(fin_diff.shape[1]*0.01*float(percentage_of_interval.split('%')[0]))))\n",
    "    print(\"num_points: \", num_points)\n",
    "    t0_index = 0\n",
    "    \n",
    "    if where_to_start == 'start':\n",
    "        t0_index = num_points//2 + 1\n",
    "    \n",
    "    if where_to_start == 'middle':\n",
    "        t0_index = fin_diff.shape[1]//2\n",
    "        \n",
    "    if where_to_start == 'end':\n",
    "        t0_index = fin_diff.shape[1] - (num_points//2 + 1)\n",
    "    \n",
    "    print(\"t0_index: \", t0_index)\n",
    "    \n",
    "    y_initial = [np.mean(fin_diff[0,t0_index-num_points//2:t0_index+num_points//2+1])]\n",
    "    for n in range(1,N):\n",
    "        if n%2==1:\n",
    "            new_num_points = num_points-(2+n)\n",
    "        else:\n",
    "            new_num_points = num_points-(3+n)\n",
    "        \n",
    "        y_initial.append(np.mean(fin_diff[n,t0_index-new_num_points//2:t0_index+new_num_points//2+1]))\n",
    "    \n",
    "    return np.reshape(y_initial,(-1,1)), t0_index\n",
    "\n",
    "# Calculate explicit solution to system of first order equations\n",
    "def Cinverse_negA(C,A):\n",
    "    print(\"CINVERSENEGA\")\n",
    "    print(\"C: \\n\", C)\n",
    "    print(\"A: \\n\", A)\n",
    "    print(\"np.linalg.inv(C): \\n\", np.linalg.inv(C))\n",
    "    return np.matmul(np.linalg.inv(C), -A)\n",
    "\n",
    "def first_order_system_sol(t, t_0, y_0, C, Cinverse_negA):\n",
    "    mat_exp = expm(C*(t_0-t))\n",
    "    #if t<0.02:\n",
    "    #    print(\"C: \", C)\n",
    "    #    print(\"mat_exp: \", mat_exp)\n",
    "    #print(\"Cinverse_negA: \", Cinverse_negA)\n",
    "    return np.matmul(mat_exp, y_0-Cinverse_negA)+Cinverse_negA\n",
    "\n",
    "\n",
    "### Code for the model which puts the above together ###\n",
    "\n",
    "def model(time, time_series, N, acc_points_to_use, init_points_to_use, init_pos, coeff_lock='c_0', acc_list=None):\n",
    "    print(\"\")\n",
    "    print(\"\")\n",
    "    print(\"!!!!!!! START !!!!!!!!!\")\n",
    "    num_points_to_use = [1]*N\n",
    "    initial_time = 't_0'\n",
    "    tolerance = 1e-6\n",
    "    \n",
    "    all_coeffs, fin_diff = solved_all_coeffs_analytical(time, time_series, N, acc_points_to_use, coefficient_to_lock=coeff_lock, accuracy_list=acc_list)\n",
    "    print(\"fin_diff: \\n\", fin_diff)\n",
    "    print(\"!!!!!!!!!!!!!\")\n",
    "    print(\"all_coeffs: \\n\", all_coeffs)\n",
    "\n",
    "    C, A = coeff_list_to_system_matrix(all_coeffs, tolerance)\n",
    "    \n",
    "    try:\n",
    "        print(\"Start of try\")\n",
    "        num_derivs = C.shape[0]+1\n",
    "        print(\"num_derivs: \", num_derivs)\n",
    "        fin_diff = fin_diff[:num_derivs]\n",
    "        print(\"If the effective ODE order is different, we correspondingly remove rows from fin_diff: \\n\", fin_diff)\n",
    "        #y_init = initial_conditions(fin_diff, num_points_to_use, initial_time)\n",
    "        y_init, t0_index = initial_cond(fin_diff, init_points_to_use, init_pos, num_derivs-1)\n",
    "        print(\"initial conditions: \\n\", y_init)\n",
    "        #t_0 = time[int(initial_time.split(\"_\")[1])]\n",
    "        t_0 = time[t0_index]\n",
    "        print(\"init time: \", t_0)\n",
    "\n",
    "        Cinv_negA = Cinverse_negA(C,A)\n",
    "        print(\"Cinv_negA: \\n\", Cinv_negA)\n",
    "        y_sol = []\n",
    "\n",
    "        for t in time:\n",
    "            y_sol.append(first_order_system_sol(t, t_0, y_init, C, Cinv_negA))\n",
    "\n",
    "        y_sol_forward = y_sol.copy()\n",
    "        for t in (time+time[-1]):\n",
    "            y_sol_forward.append(first_order_system_sol(t, t_0, y_init, C, Cinv_negA))\n",
    "        print(\"End of try\")\n",
    "        \n",
    "        return np.array(y_sol), np.array(y_sol_forward), (t0_index, t_0, y_init, C, Cinv_negA), all_coeffs\n",
    "\n",
    "    # If the coefficients in front of all the derivatives are 0:\n",
    "    except:\n",
    "        print(\"In except\")\n",
    "        y_sol = np.ones_like(time)*(-1)*all_coeffs[0,0]    \n",
    "        y_sol = np.reshape(y_sol, (-1,1,1))\n",
    "        y_sol_forward = y_sol.copy()\n",
    "        \n",
    "        return np.array(y_sol), np.array(y_sol_forward), (None, None, None, None, None), all_coeffs\n",
    "\n",
    "\n",
    "### Function that loops over and locks each coefficient, comparing which gives the best result ###\n",
    "\n",
    "def best_c_j_to_lock(time, true_ts, N, acc_points_to_use, init_points_to_use, init_pos):\n",
    "    y_sols_c_j_locked = []\n",
    "    rel_errors = []\n",
    "    model_params_list = []\n",
    "    all_coeffs_list = []\n",
    "    \n",
    "    for j in range(N+1):\n",
    "        y_sol, y_sol_forward, model_params, all_coeffs = model(time, true_ts, N, acc_points_to_use, init_points_to_use, init_pos, coeff_lock=f'c_{j}')\n",
    "        y_sols_c_j_locked.append(y_sol)\n",
    "        rel_errors.append(relative_error(y_sol[:,0,0], true_ts))\n",
    "        model_params_list.append(model_params)\n",
    "        all_coeffs_list.append(all_coeffs)\n",
    "    best_j = rel_errors.index(np.nanmin(np.array(rel_errors)))\n",
    "    print(\"\")\n",
    "    print(\"\")\n",
    "    print(\"Relative errors: \\n\", rel_errors)\n",
    "    print(f\"Best c_j to lock: c_{best_j}\", )\n",
    "    return y_sols_c_j_locked[best_j], rel_errors[best_j], model_params_list[best_j], best_j, all_coeffs_list[best_j]\n",
    "\n",
    "### Function that loops through N (up to and including highest_N), and compares the results ###\n",
    "def best_N(time, true_ts, ts_name, highest_N, plot_best_of_each_N, init_points_to_use, init_pos, save_figure=False, acc_points_to_use='few'):\n",
    "    \n",
    "    best_sols_dict = {}\n",
    "    best_errors_dict = {}\n",
    "    best_model_params_dict = {}\n",
    "    best_locked_j_dict = {}\n",
    "    best_all_coeffs_dict = {}\n",
    "    for N in range(1,highest_N+1):\n",
    "        y_best, best_rel_error, best_model_params, best_locked_j, best_all_coeffs = best_c_j_to_lock(time, true_ts, N, acc_points_to_use,  init_points_to_use, init_pos)\n",
    "        best_sols_dict[N] = y_best\n",
    "        best_errors_dict[N] = best_rel_error\n",
    "        best_model_params_dict[N] = best_model_params\n",
    "        best_locked_j_dict[N] = best_locked_j\n",
    "        best_all_coeffs_dict[N] = best_all_coeffs\n",
    "        \n",
    "    best_sols_and_errors = best_errors_dict.items()\n",
    "    \n",
    "    best_N = min(best_errors_dict, key=best_errors_dict.get)\n",
    "    \n",
    "    best_y = best_sols_dict[best_N]\n",
    "    best_err = best_errors_dict[best_N]\n",
    "    best_model_params = best_model_params_dict[best_N]\n",
    "    \n",
    "    #print(\"best_sols_dict: \\n\", best_sols_dict)\n",
    "    print(\"\")\n",
    "    print(\"\")\n",
    "    print(\"!!\")\n",
    "    print(\"best_locked_j_dict: \", best_locked_j_dict)\n",
    "    print(\"best_errors_dict: \\n\", best_errors_dict)\n",
    "    print(\"best_all_coeffs_dict: \\n\", best_all_coeffs_dict)\n",
    "    print(\"best_N: \", best_N)\n",
    "    print(\"best_locked_c_j: \", best_locked_j_dict[best_N])\n",
    "    if plot_best_of_each_N:\n",
    "        for N, best_y in best_sols_dict.items():\n",
    "            plt.title(ts_name)         \n",
    "            #plt.plot(time, el_price_2019_data[:len(time)], color=\"gray\",alpha=0.5, label=\"Original\")\n",
    "            plt.plot(time, true_ts, label='True')\n",
    "            plt.plot(time, best_y[:,0,0], label=f'Pred. N={N}. Locked c_{best_locked_j_dict[N]}.\\nInit. cond: {init_points_to_use}, {init_pos}')\n",
    "            plt.legend()\n",
    "            \n",
    "            if save_figure:\n",
    "                if N==best_N:\n",
    "                    plt.savefig(savefigs_loc+ts_name+f\"__BEST__N_{N}_best_locked_{best_locked_j_dict[N]}_init_{init_points_to_use.split('%')[0]}_{init_pos}__TRAINING__rel_error_{best_errors_dict[N]}.png\")\n",
    "                else:\n",
    "                    plt.savefig(savefigs_loc+\"not_best_figures\\\\\"+ts_name+f\"__NOT_BEST__N_{N}_best_locked_{best_locked_j_dict[N]}_init_{init_points_to_use.split('%')[0]}_{init_pos}__TRAINING__rel_error_{best_errors_dict[N]}.png\")       \n",
    "            plt.show()\n",
    "            print(\"Training, relative error: \", best_errors_dict[N])\n",
    "            print(\"ODE Coefficients: \\n\", best_all_coeffs_dict[N])\n",
    "            \n",
    "    return best_y, best_err, best_model_params, best_N, best_locked_j_dict[best_N]\n",
    "\n",
    "### Function that calculates and plots function given user-specified coefficients and initial values ###\n",
    "def pre_chosen_coeffs_model(chosen_coeffs, chosen_init_conds, time_interval, t0=0):\n",
    "    C, A = coeff_list_to_system_matrix(chosen_coeffs, 1e-10)\n",
    "    Cinv_negA = Cinverse_negA(C,A)\n",
    "    y_pred_chosen=[]\n",
    "    for t in time_interval:\n",
    "        y_pred_chosen.append(first_order_system_sol(t, t0, chosen_init_conds, C, Cinv_negA))\n",
    "        \n",
    "    plt.plot(time_interval, np.array(y_pred_chosen)[:,0,0])\n",
    "    #plt.savefig(savefigs_loc+f\"pre_chosen_coeffs_sol_{time_interval[-1]}.png\")\n",
    "    plt.show()\n",
    "    return np.array(y_pred_chosen)[:,0,0]\n",
    "\n",
    "### Function that takes in model parameters, and gives the model prediction forward in time ###\n",
    "def model_params_testing(time, best_model_params, true, title, forward_in_time, best_N, best_locked_j, save_figure = False):\n",
    "    sol_list = []\n",
    "\n",
    "    t0_index = best_model_params[0]\n",
    "    t0 = best_model_params[1]\n",
    "    y_init = best_model_params[2]\n",
    "    C = best_model_params[3]\n",
    "    Cinv_negA = best_model_params[4]\n",
    "\n",
    "    for t in time:\n",
    "        sol = first_order_system_sol(t, t0, y_init, C, Cinv_negA)\n",
    "        sol_list.append(sol)\n",
    "        \n",
    "    \n",
    "    rel_error=0 \n",
    "    if forward_in_time:\n",
    "        middle_index = time.shape[0]//2\n",
    "        rel_error = relative_error(np.array(sol_list)[middle_index:,0,0],true[middle_index:])\n",
    "        print(\"Testing (forward in time), Relative error (on future time points): \", rel_error)\n",
    "        \n",
    "        #plt.plot(time, tesla_2021_2022_data[:len(time)], color=\"gray\",alpha=0.5, label=\"Original\")\n",
    "        #plt.plot(time, true, label='Smoothed')\n",
    "        plt.plot(time, np.array(sol_list)[:,0,0], label='Pred. data')\n",
    "        plt.title(title)\n",
    "        plt.legend()\n",
    "        if save_figure:\n",
    "            plt.savefig(savefigs_loc+title+f\"__best_N_{best_N}_best_locked_{best_locked_j}__FIT__rel_error_{rel_error}.png\")\n",
    "    else:\n",
    "        rel_error = relative_error(np.array(sol_list)[:,0,0],true[:])\n",
    "        print(\"Testing (same interval), Relative error: \", rel_error)\n",
    "        \n",
    "        #plt.plot(time, tesla_2021_data, color=\"gray\", alpha=0.5, label=\"Original\")\n",
    "        #plt.plot(time, true, label='Smoothed')\n",
    "        plt.plot(time, np.array(sol_list)[:,0,0], label='Pred. data')\n",
    "        plt.title(title)\n",
    "        plt.legend()\n",
    "        if save_figure:\n",
    "            plt.savefig(savefigs_loc+title+f\"__best_N_{best_N}_best_locked_{best_locked_j}__TESTING__rel_error_{rel_error}.png\")\n",
    "    \n",
    "    plt.show()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8d00ee8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def relative_error(prediction, true):\n",
    "    mse = np.mean(np.power(prediction-true,2))\n",
    "    print(\"mse: \", mse)\n",
    "    true_size = np.mean(np.power(true,2))\n",
    "    print(\"true size: \", true_size)\n",
    "    return mse/true_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "228ca090",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Pure functions ###\n",
    "\n",
    "def exp(t):\n",
    "    return np.exp(t)\n",
    "\n",
    "def cos(t):\n",
    "    return np.cos(t*4*np.pi)\n",
    "\n",
    "def quartic(t):\n",
    "    return 0.5*(np.power(3.5*t-1, 4)-3*np.power(3.5*t-1,3)+3*(3.5*t-1))\n",
    "\n",
    "def composite(t):\n",
    "    return (np.cos(15*t+1/2)+np.log(25*t+1/2))/5\n",
    "\n",
    "def pw_linear(t):\n",
    "    result = []\n",
    "    for x in t:\n",
    "        x = x%1\n",
    "        if x < 0.2:\n",
    "            result.append(0)\n",
    "        elif x < 0.4:\n",
    "            result.append(5 * (x - 0.2))\n",
    "        elif x < 0.6:\n",
    "            result.append(1)\n",
    "        elif x < 0.8:\n",
    "            result.append(4-5*x)\n",
    "        else:\n",
    "            result.append(0)\n",
    "    return np.array(result)\n",
    "\n",
    "\n",
    "### Noisy functions ###\n",
    "\n",
    "def exp_noisy(t):\n",
    "    num=len(t)\n",
    "    std = 0.01\n",
    "    noise = np.random.normal(0,std, num)\n",
    "    return np.exp(t)+noise\n",
    "\n",
    "def cos_noisy(t):\n",
    "    num=len(t)\n",
    "    std = 0.01\n",
    "    noise = np.random.normal(0,std, num)\n",
    "    return np.cos(t*4*np.pi)+noise\n",
    "\n",
    "def quartic_noisy(t):\n",
    "    num=len(t)\n",
    "    std = 0.01\n",
    "    noise = np.random.normal(0,std, num)\n",
    "    return 0.5*(np.power(3.5*t-1, 4)-3*np.power(3.5*t-1,3)+3*(3.5*t-1))+noise\n",
    "\n",
    "def composite_noisy(t):\n",
    "    num=len(t)\n",
    "    std = 0.01\n",
    "    noise = np.random.normal(0,std, num)\n",
    "    return (np.cos(15*t+1/2)+np.log(25*t+1/2))/5 + noise\n",
    "\n",
    "def pw_linear_noisy(t):\n",
    "    num=len(t)\n",
    "    std = 0.01\n",
    "    noise = np.random.normal(0,std, num)\n",
    "    result = []\n",
    "    for x in t:\n",
    "        x = x%1\n",
    "        if x < 0.2:\n",
    "            result.append(0)\n",
    "        elif x < 0.4:\n",
    "            result.append(5 * (x - 0.2))\n",
    "        elif x < 0.6:\n",
    "            result.append(1)\n",
    "        elif x < 0.8:\n",
    "            result.append(4-5*x)\n",
    "        else:\n",
    "            result.append(0)\n",
    "    return np.array(result)+noise\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7ed5f4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## training interval\n",
    "\n",
    "time = np.linspace(0,1,1001)\n",
    "h=time[1]-time[0]\n",
    "\n",
    "ts_exp = exp(time)\n",
    "ts_cos = cos(time)\n",
    "ts_quartic = quartic(time)\n",
    "ts_composite = composite(time)\n",
    "ts_pw_linear = pw_linear(time)\n",
    "\n",
    "# Training\n",
    "time_training = time[::2]\n",
    "\n",
    "ts_exp_training = ts_exp[::2]\n",
    "ts_cos_training = ts_cos[::2]\n",
    "ts_quartic_training = ts_quartic[::2]\n",
    "ts_composite_training = ts_composite[::2]\n",
    "ts_pw_linear_training = ts_pw_linear[::2]\n",
    "\n",
    "# Testing\n",
    "time_testing = time[1::2]\n",
    "\n",
    "ts_exp_testing = ts_exp[1::2]\n",
    "ts_cos_testing = ts_cos[1::2]\n",
    "ts_quartic_testing = ts_quartic[1::2]\n",
    "ts_composite_testing = ts_composite[1::2]\n",
    "ts_pw_linear_testing = ts_pw_linear[1::2]\n",
    "\n",
    "\n",
    "# Pre-chosen_coefficients\n",
    "chosen_coeffs = np.reshape([0,20,5,1],(-1,1))\n",
    "chosen_inits = np.reshape([1,5],(-1,1))\n",
    "ts_pre_chosen_coeffs = pre_chosen_coeffs_model(chosen_coeffs, chosen_inits, time)\n",
    "\n",
    "ts_pre_chosen_coeffs_training = ts_pre_chosen_coeffs[::2]\n",
    "ts_pre_chosen_coeffs_testing = ts_pre_chosen_coeffs[1::2]\n",
    "\n",
    "\n",
    "# training\n",
    "\n",
    "\n",
    "# Forward in time. Interval is now twice the length of the original interval (ends at 2 instead of 1)\n",
    "\n",
    "time_forward = list(time).copy()\n",
    "time_forward.extend(list((time+time[-1])[1:]))\n",
    "time_forward = np.array(time_forward)\n",
    "\n",
    "# testing\n",
    "time_forward_testing = time_forward[1::2]\n",
    "\n",
    "ts_exp_forward = exp(time_forward_testing)\n",
    "ts_cos_forward = cos(time_forward_testing)\n",
    "ts_quartic_forward = quartic(time_forward_testing)\n",
    "ts_composite_forward = composite(time_forward_testing)\n",
    "ts_pw_linear_forward = pw_linear(time_forward_testing)\n",
    "ts_pre_chosen_coeffs_forward = pre_chosen_coeffs_model(chosen_coeffs, chosen_inits, time_forward_testing)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Noisy time series:\n",
    "\n",
    "ts_exp_noisy = exp_noisy(time)\n",
    "ts_cos_noisy = cos_noisy(time)\n",
    "ts_quartic_noisy = quartic_noisy(time)\n",
    "ts_composite_noisy = composite_noisy(time)\n",
    "ts_pw_linear_noisy = pw_linear_noisy(time)\n",
    "\n",
    "num=len(ts_pre_chosen_coeffs)\n",
    "std = 0.01\n",
    "noise = np.random.normal(0,std, num)\n",
    "ts_pre_chosen_coeffs_noisy = ts_pre_chosen_coeffs + noise\n",
    "\n",
    "# Training\n",
    "ts_exp_noisy_training = ts_exp_noisy[::2]\n",
    "ts_cos_noisy_training = ts_cos_noisy[::2]\n",
    "ts_quartic_noisy_training = ts_quartic_noisy[::2]\n",
    "ts_composite_noisy_training = ts_composite_noisy[::2]\n",
    "ts_pw_linear_noisy_training = ts_pw_linear_noisy[::2]\n",
    "ts_pre_chosen_coeffs_noisy_training = ts_pre_chosen_coeffs_noisy[::2]\n",
    "\n",
    "# smoothed training\n",
    "ts_quartic_noisy_training_smoothed = savgol_filter(ts_quartic_training, int(np.floor(len(ts_quartic_noisy_training)*0.1)), 3)\n",
    "ts_pw_linear_noisy_training_smoothed = savgol_filter(ts_pw_linear_noisy_training, int(np.floor(len(ts_pw_linear_noisy_training)*0.1)), 3)\n",
    "ts_pre_chosen_coeffs_noisy_training_smoothed = savgol_filter(ts_pre_chosen_coeffs_noisy_training, int(np.floor(len(ts_pre_chosen_coeffs_noisy_training)*0.1)), 3)\n",
    "\n",
    "\n",
    "# Testing\n",
    "ts_exp_noisy_testing = ts_exp_noisy[1::2]\n",
    "ts_cos_noisy_testing = ts_cos_noisy[1::2]\n",
    "ts_quartic_noisy_testing = ts_quartic_noisy[1::2]\n",
    "ts_composite_noisy_testing = ts_composite_noisy[1::2]\n",
    "ts_pw_linear_noisy_testing = ts_pw_linear_noisy[1::2]\n",
    "ts_pre_chosen_coeffs_noisy_testing = ts_pre_chosen_coeffs_noisy[1::2]\n",
    "\n",
    "#smoothed testing\n",
    "ts_quartic_noisy_testing_smoothed = savgol_filter(ts_quartic_noisy_testing, int(np.floor(len(ts_quartic_noisy_testing)*0.1)), 3)\n",
    "ts_pw_linear_noisy_testing_smoothed = savgol_filter(ts_pw_linear_noisy_testing, int(np.floor(len(ts_pw_linear_noisy_testing)*0.1)), 3)\n",
    "ts_pre_chosen_coeffs_noisy_testing_smoothed = savgol_filter(ts_pre_chosen_coeffs_noisy_testing, int(np.floor(len(ts_pre_chosen_coeffs_noisy_testing)*0.1)), 3)\n",
    "\n",
    "\n",
    "# Forward in time. Interval is now twice the length of the original interval (ends at 2 instead of 1)\n",
    "\n",
    "# testing\n",
    "ts_exp_noisy_forward = exp_noisy(time_forward_testing)\n",
    "ts_cos_noisy_forward = cos_noisy(time_forward_testing)\n",
    "ts_quartic_noisy_forward = quartic_noisy(time_forward_testing)\n",
    "ts_composite_noisy_forward = composite_noisy(time_forward_testing)\n",
    "ts_pw_linear_noisy_forward = pw_linear_noisy(time_forward_testing)\n",
    "\n",
    "\n",
    "#smoothed forward\n",
    "ts_quartic_noisy_forward_smoothed = savgol_filter(ts_quartic_noisy_forward, int(np.floor(len(ts_quartic_noisy_testing)*0.1)), 3)\n",
    "ts_pw_linear_noisy_forward_smoothed = savgol_filter(ts_pw_linear_noisy_forward, int(np.floor(len(ts_pw_linear_noisy_testing)*0.1)), 3)\n",
    "#ts_pre_chosen_coeffs_noisy_forward_smoothed = savgol_filter(ts_pre_chosen_coeffs_noisy_forward, int(np.floor(len(ts_pre_chosen_coeffs_noisy_testing)*0.1)), 3)\n",
    "\n",
    "\n",
    "num2=len(ts_pre_chosen_coeffs_forward)\n",
    "std2 = 0.01\n",
    "noise2 = np.random.normal(0,std, num2)\n",
    "ts_pre_chosen_coeffs_noisy_forward = ts_pre_chosen_coeffs_forward+noise2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e22c0c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "##exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cf4573f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "269aa24c",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(ts_exp_training)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f187d693",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_exp_training, \"ts_exp_training\",1, True, '10%', 'middle', save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_exp_training, \"ts_exp_training\", 1, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a9984d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_exp_testing, \"ts_exp_testing\", False, best_N_middle, best_j_middle)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_exp_forward, \"ts_exp_forward\", True, best_N_end, best_j_end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95decdc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e2bdaae",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_exp_training,250,\"ts_exp_training__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_exp_testing, \"ts_exp_testing__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_exp_forward, \"ts_exp_forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b3c6805",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9454d6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cos # "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fc58855",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_cos_training, \"ts_cos_training\",2, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_cos_training, \"ts_cos_training\", 2, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60ee0d2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_cos_testing, \"ts_cos_testing\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_cos_forward, \"ts_cos_forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f13e4d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5549c289",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_cos_training,250,\"ts_cos_training__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_cos_testing, \"ts_cos_testing__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_cos_forward, \"ts_cos_forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ef2afac",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2669325f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# quartic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf984f18",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_quartic_training, \"ts_quartic_training\",4, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_quartic_training, \"ts_quartic_training\", 4, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "376650a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_quartic_testing, \"ts_quartic_testing\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_quartic_forward, \"ts_quartic_forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "540c909d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdeb59f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_quartic_training,250,\"ts_quartic_training__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_quartic_testing, \"ts_quartic_testing__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_quartic_forward, \"ts_quartic_forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eda7ce4b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a468da1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# composite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f349e9a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_composite_training, \"ts_composite_training\",20, True, '10%', 'middle',save_figure=True)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_composite_training, \"ts_composite_training\", 20, True, '10%', 'end', save_figure=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a556c68",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_composite_testing, \"ts_composite_testing\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_composite_forward, \"ts_composite_forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b850ddee",
   "metadata": {},
   "outputs": [],
   "source": [
    "#nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e37a956",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_composite_training,250,\"ts_composite_training__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_composite_testing, \"ts_composite_testing__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_composite_forward, \"ts_composite_forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86734381",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c95ce2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pw linear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "727f73d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_pw_linear_training, \"ts_pw_linear_training\",20, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_pw_linear_training, \"ts_pw_linear_training\", 20, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a963eed",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_pw_linear_testing, \"ts_pw_linear_testing\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_pw_linear_forward, \"ts_pw_linear_forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1472d843",
   "metadata": {},
   "outputs": [],
   "source": [
    "# nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d814a16",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_pw_linear_training,350,\"ts_pw_linear_training__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_pw_linear_testing, \"ts_pw_linear_testing__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_pw_linear_forward, \"ts_pw_linear_forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc388307",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08fa0ee4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08a34b28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pre_chosen_coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23ae1dae",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_pre_chosen_coeffs_training, \"ts_chosen_coeffs_training\",4, True, '10%', 'middle',save_figure=True)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_pre_chosen_coeffs_training, \"ts_chosen_coeffs_training\", 4, True, '10%', 'end', save_figure=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c1a6b19",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_pre_chosen_coeffs_testing, \"ts_chosen_coeffs_testing\", False, best_N_middle, best_j_middle, save_figure=True)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_pre_chosen_coeffs_forward, \"ts_chosen_coeffs_forward\", True, best_N_end, best_j_end, save_figure=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06cba6a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "735c688c",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_pre_chosen_coeffs_training,350,\"ts_pre_chosen_coeffs_training__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_pre_chosen_coeffs_testing, \"ts_pre_chosen_coeffs_testing__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_pre_chosen_coeffs_forward, \"ts_pre_chosen_coeffs_forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1695de3d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4700be68",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba2da104",
   "metadata": {},
   "outputs": [],
   "source": [
    "####### NOISY ##########"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "184f589f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Quartic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40790bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_quartic_noisy_training, \"ts_quartic_noisy_training\",20, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_quartic_noisy_training, \"ts_quartic_noisy_training\", 20, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b558654",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_quartic_noisy_testing, \"ts_quartic_noisy_testing\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing[:], best_model_params_end, ts_quartic_noisy_forward[:], \"ts_quartic_noisy_forward\", True, best_N_end, best_j_end, save_figure=False)\n",
    "model_params_testing(time_forward_testing[:-350], best_model_params_end, ts_quartic_noisy_forward[:-350], \"ts_quartic_noisy_forward1.3forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "774d7f16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5e38838",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_quartic_noisy_training,500,\"ts_quartic_noisy_training__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_quartic_noisy_testing, \"ts_quartic_noisy_testing__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_quartic_noisy_forward, \"ts_quartic_noisy_forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b15b9539",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79e44159",
   "metadata": {},
   "outputs": [],
   "source": [
    "# quartic smoothed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fc2087f",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_quartic_noisy_training_smoothed, \"ts_quartic_noisy_training_smoothed\",20, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_quartic_noisy_training_smoothed, \"ts_quartic_noisy_training_smoothed\", 20, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce9090cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_quartic_noisy_testing_smoothed, \"ts_quartic_noisy_testing_smoothed\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_quartic_noisy_forward_smoothed, \"ts_quartic_noisy_forward_smoothed\", True, best_N_end, best_j_end, save_figure=False)\n",
    "model_params_testing(time_forward_testing[:-350], best_model_params_end, ts_quartic_noisy_forward_smoothed[:-350], \"ts_quartic_noisy_forward_smoothed_1.3forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41f800f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2b87f65",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_quartic_noisy_training_smoothed,500,\"ts_quartic_noisy_training_smoothed__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_quartic_noisy_testing_smoothed, \"ts_quartic_noisy_testing_smoothed__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_quartic_noisy_forward_smoothed, \"ts_quartic_noisy_forward_smoothed__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96cc647c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f6a132b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80c1473b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pw_linear noisy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7a77bb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_pw_linear_noisy_training, \"ts_pw_linear_noisy_training\",20, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_pw_linear_noisy_training, \"ts_pw_linear_noisy_training\", 20, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36a0f2f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_pw_linear_noisy_testing, \"ts_pw_linear_noisy_testing\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_pw_linear_noisy_forward, \"ts_pw_linear_noisy_forward\", True, best_N_end, best_j_end, save_figure=False)\n",
    "model_params_testing(time_forward_testing[:-350], best_model_params_end, ts_pw_linear_noisy_forward[:-350], \"ts_pw_linear_noisy_forward_1.3forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a30d70d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef2a8a6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_pw_linear_noisy_training,500,\"ts_pw_linear_noisy_training__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_pw_linear_noisy_testing, \"ts_pw_linear_noisy_testing__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_pw_linear_noisy_forward, \"ts_pw_linear_noisy_forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd81f540",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2d09414",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pw_linear noisy smoothed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1029995",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_pw_linear_noisy_training_smoothed, \"ts_pw_linear_noisy_training_smoothed\",20, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_pw_linear_noisy_training_smoothed, \"ts_pw_linear_noisy_training_smoothed\", 20, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5123e1a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_pw_linear_noisy_testing_smoothed, \"ts_pw_linear_noisy_testing_smoothed\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_pw_linear_noisy_forward_smoothed, \"ts_pw_linear_noisy_forward_smoothed\", True, best_N_end, best_j_end, save_figure=False)\n",
    "model_params_testing(time_forward_testing[:-350], best_model_params_end, ts_pw_linear_noisy_forward_smoothed[:-350], \"ts_pw_linear_noisy_forward_1.3forward_smoothed\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a62259f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32fadf4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_pw_linear_noisy_training_smoothed,500,\"ts_pw_linear_noisy_training_smoothed__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_pw_linear_noisy_testing_smoothed, \"ts_pw_linear_noisy_testing_smoothed__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_pw_linear_noisy_forward_smoothed, \"ts_pw_linear_noisy_forward_smoothed__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7aa91c8f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fc24162",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pre-chosen coeffs noisy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9aefcd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_pre_chosen_coeffs_noisy_training, \"ts_pre_chosen_coeffs_noisy_training\",20, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_pre_chosen_coeffs_noisy_training, \"ts_pre_chosen_coeffs_noisy_training\", 20, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81fce938",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_pre_chosen_coeffs_noisy_testing, \"ts_pre_chosen_coeffs_noisy_testing\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_pre_chosen_coeffs_noisy_forward, \"ts_pre_chosen_coeffs_noisy_forward\", True, best_N_end, best_j_end, save_figure=False)\n",
    "model_params_testing(time_forward_testing[:-350], best_model_params_end, ts_pre_chosen_coeffs_noisy_forward[:-350], \"ts_pre_chosen_coeffs_noisy_forward_1.3forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d302e53",
   "metadata": {},
   "outputs": [],
   "source": [
    "# nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b302fa8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_pre_chosen_coeffs_noisy_training,500,\"ts_pre_chosen_coeffs_noisy_training__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_pre_chosen_coeffs_noisy_testing, \"ts_pre_chosen_coeffs_noisy_testing__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_pre_chosen_coeffs_noisy_forward, \"ts_pre_chosen_coeffs_noisy_forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65f4fa86",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28b22664",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pre_chosen coeffs noisy smoothed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2690a92e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_training, ts_pre_chosen_coeffs_noisy_training_smoothed, \"ts_pre_chosen_coeffs_noisy_training_smoothed\",20, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_training, ts_pre_chosen_coeffs_noisy_training_smoothed, \"ts_pre_chosen_coeffs_noisy_training_smoothed\", 20, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7797ee2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_testing, best_model_params_middle, ts_pre_chosen_coeffs_noisy_testing_smoothed, \"ts_pre_chosen_coeffs_noisy_testing_smoothed\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_forward_testing, best_model_params_end, ts_pre_chosen_coeffs_noisy_forward_smoothed, \"ts_pre_chosen_coeffs_noisy_forward_smoothed\", True, best_N_end, best_j_end, save_figure=False)\n",
    "model_params_testing(time_forward_testing[:-350], best_model_params_end, ts_pre_chosen_coeffs_noisy_forward_smoothed[:-350], \"ts_pre_chosen_coeffs_noisy_forward_1.3forward_smoothed\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "443d39cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fa574a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_training, ts_pre_chosen_coeffs_noisy_training_smoothed,500,\"ts_pre_chosen_coeffs_noisy_training_smoothed__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test on same interval\n",
    "nn_test(best_width_model, best_width, time_testing, ts_pre_chosen_coeffs_noisy_testing_smoothed, \"ts_pre_chosen_coeffs_noisy_testing_smoothed__NN\",forward_time=False, save_fig=save)\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_forward_testing, ts_pre_chosen_coeffs_noisy_forward_smoothed, \"ts_pre_chosen_coeffs_noisy_forward_smoothed__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2e43e72",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "355d066f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Real world data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd053a08",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5b3212a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tesla stocks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e195482",
   "metadata": {},
   "outputs": [],
   "source": [
    "tesla_2021_df = pd.read_csv(r\"C:\\Users\\Robin\\OneDrive\\Dokumenter\\Universitet\\Fordypningsprosjekt - TMA4500\\STOCK_US_XNAS_TSLA_2021.csv\")\n",
    "tesla_2021_2022_df = pd.read_csv(r\"C:\\Users\\Robin\\OneDrive\\Dokumenter\\Universitet\\Fordypningsprosjekt - TMA4500\\STOCK_US_XNAS_TSLA_2021_2022.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cae3a5ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "tesla_2021_df = tesla_2021_df[::-1].reset_index(drop=True)\n",
    "tesla_2021_2022_df = tesla_2021_2022_df[::-1].reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71d3e271",
   "metadata": {},
   "outputs": [],
   "source": [
    "tesla_2021_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "979a1215",
   "metadata": {},
   "outputs": [],
   "source": [
    "tesla_2021_data = tesla_2021_df['Open'].dropna()\n",
    "tesla_2021_2022_data = tesla_2021_2022_df['Open'].dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fda7a55a",
   "metadata": {},
   "outputs": [],
   "source": [
    "tesla_2021_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67c9ecb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(tesla_2021_data)\n",
    "plt.title(\"Tesla Stock Price (at open) in 2021\")\n",
    "plt.ylabel(\"Tesla Stock Price in USD\")\n",
    "plt.xlabel(\"Stock market (NASDAQ), days open, starting from Jan. 4th 2021\")\n",
    "#plt.savefig(\"C:\\\\Users\\\\Robin\\\\OneDrive\\\\Dokumenter\\\\Universitet\\\\Fordypningsprosjekt - TMA4500\\\\tesla_stock_2021.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a2799ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(tesla_2021_2022_data)\n",
    "plt.title(\"Tesla Stock Price (at open) in 2021 and 2022\")\n",
    "plt.ylabel(\"Tesla Stock Price in USD\")\n",
    "plt.xlabel(\"Stock market (NASDAQ), days open, starting from Jan. 4th 2021\")\n",
    "#plt.savefig(\"C:\\\\Users\\\\Robin\\\\OneDrive\\\\Dokumenter\\\\Universitet\\\\Fordypningsprosjekt - TMA4500\\\\tesla_stock_2021_2022.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6d08aec",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_tesla = np.linspace(0,1, len(tesla_2021_data))\n",
    "time_tesla_forward = np.linspace(0,2, len(tesla_2021_2022_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05241137",
   "metadata": {},
   "outputs": [],
   "source": [
    "tesla_2021_2022_data = tesla_2021_2022_data/max(tesla_2021_data)\n",
    "tesla_2021_data = tesla_2021_data/max(tesla_2021_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41b6161f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_tesla, tesla_2021_data)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f0dd072",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_tesla_forward, tesla_2021_2022_data)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bede0d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "tesla_2021_data_smoothed = savgol_filter(tesla_2021_data, int(np.floor(len(tesla_2021_data)*0.1)), 3)\n",
    "tesla_2021_2022_data_smoothed = savgol_filter(tesla_2021_2022_data, int(np.floor(len(tesla_2021_data)*0.1)), 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "941c025c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_tesla, tesla_2021_data_smoothed)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a67fc745",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_tesla_forward, tesla_2021_2022_data_smoothed)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b07217fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_tesla, tesla_2021_data, color=\"gray\", alpha=0.5, label=\"Data\")\n",
    "plt.plot(time_tesla, tesla_2021_data_smoothed, label = \"Smoothed\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f0fb38b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_tesla_forward, tesla_2021_2022_data, color=\"gray\", label=\"Data\")\n",
    "plt.plot(time_tesla_forward, tesla_2021_2022_data_smoothed, label = \"Smoothed\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69e6fd80",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_tesla, tesla_2021_data_smoothed, \"tesla_2021_data_smoothed\",20, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_tesla, tesla_2021_data_smoothed, \"tesla_2021_data_smoothed\", 20, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4692c6b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#h=time_tesla[1]-time_tesla[0]\n",
    "#time_tesla_testing = time_tesla+h/2\n",
    "#time_tesla_testing = time_tesla_testing[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "657e830c",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_tesla, best_model_params_middle, tesla_2021_data_smoothed, \"tesla_2021_data_testing\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_tesla_forward, best_model_params_end, tesla_2021_2022_data_smoothed, \"tesla_2021_2022_data\", True, best_N_end, best_j_end, save_figure=False)\n",
    "model_params_testing(time_tesla_forward[:-200], best_model_params_end, tesla_2021_2022_data_smoothed[:-200], \"tesla_2021_2022_datad_1.2forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "588e4e78",
   "metadata": {},
   "outputs": [],
   "source": [
    "# nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60ba1719",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_tesla, tesla_2021_data.to_numpy(),500,\"tesla_2021_data__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_tesla_forward, tesla_2021_2022_data.to_numpy(), \"tesla_2021_2022_data_forward__NN\",forward_time=True, save_fig=save)\n",
    "nn_test(best_width_model, best_width, time_tesla_forward[:-200], tesla_2021_2022_data.to_numpy()[:-200], \"tesla_2021_2022_data_forward_1.2forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ee166be",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_tesla, tesla_2021_data_smoothed,500,\"tesla_2021_data_smoothed__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_tesla_forward, tesla_2021_2022_data_smoothed, \"tesla_2021_2022_data_smoothed_forward__NN\",forward_time=True, save_fig=save)\n",
    "nn_test(best_width_model, best_width, time_tesla_forward[:-200], tesla_2021_2022_data_smoothed[:-200], \"tesla_2021_2022_data_smoothed_forward_1.2forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ade4a35",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac14e8dc",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de7a613b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Electricity price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95d77ddd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adf16eec",
   "metadata": {},
   "outputs": [],
   "source": [
    "el_price = pd.read_csv(r\"C:\\Users\\Robin\\OneDrive\\Dokumenter\\Universitet\\Fordypningsprosjekt - TMA4500\\european_wholesale_electricity_price_data_daily.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac9c4806",
   "metadata": {},
   "outputs": [],
   "source": [
    "el_price = el_price[el_price[\"Country\"] == \"Norway\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "201e1a0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "el_price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07ed0588",
   "metadata": {},
   "outputs": [],
   "source": [
    "el_price_2018 = el_price[-5*365-1:-4*365-1]\n",
    "el_price_2018_2019 = el_price[-5*365:-3*365]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e45e33aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_elprice_2018 = np.arange(0,len(el_price_2018[\"Price (EUR/MWhe)\"]))\n",
    "time_elprice_2018_2019 = np.arange(0,len(el_price_2018_2019[\"Price (EUR/MWhe)\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "702c7d84",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_elprice_2018, el_price_2018[\"Price (EUR/MWhe)\"])\n",
    "plt.title(\"Electricity prices in Norway 2018\")\n",
    "plt.ylabel(\"Price (EUR/MWh)\")\n",
    "plt.xlabel(\"Day\")\n",
    "#plt.savefig(\"C:\\\\Users\\\\Robin\\\\OneDrive\\\\Dokumenter\\\\Universitet\\\\Fordypningsprosjekt - TMA4500\\\\el_prices_2018.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e3f61a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_elprice_2018_2019, el_price_2018_2019[\"Price (EUR/MWhe)\"])\n",
    "plt.title(\"Electricity prices in Norway 2018 and 2019\")\n",
    "plt.ylabel(\"Price (EUR/MWh)\")\n",
    "plt.xlabel(\"Day\")\n",
    "#plt.savefig(\"C:\\\\Users\\\\Robin\\\\OneDrive\\\\Dokumenter\\\\Universitet\\\\Fordypningsprosjekt - TMA4500\\\\el_prices_2018_2019.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "374c99b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "el_price_2018_data = el_price_2018[\"Price (EUR/MWhe)\"]/max(el_price_2018[\"Price (EUR/MWhe)\"])\n",
    "el_price_2018_2019_data = el_price_2018_2019[\"Price (EUR/MWhe)\"]/max(el_price_2018[\"Price (EUR/MWhe)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d11f21c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_elprice_2018, el_price_2018_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d8fd51c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_elprice_2018_2019, el_price_2018_2019_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f75a7eb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "el_price_2018_smoothed = savgol_filter(el_price_2018_data, int(np.floor(len(el_price_2018_data)*0.1)), 3)\n",
    "el_price_2018_2019_smoothed = savgol_filter(el_price_2018_2019_data, int(np.floor(len(el_price_2018_data)*0.1)), 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a2344e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_elprice_2018, el_price_2018_data, color=\"gray\", label=\"Original\")\n",
    "plt.plot(time_elprice_2018, el_price_2018_smoothed, label=\"Smoothed\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2d28b8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_elprice_2018_2019, el_price_2018_2019_data, color=\"gray\", label=\"Original\")\n",
    "plt.plot(time_elprice_2018_2019, el_price_2018_2019_smoothed, label=\"Smoothed\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df470597",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_y_middle, best_err_middle, best_model_params_middle, best_N_middle, best_j_middle = best_N(time_elprice_2018, el_price_2018_smoothed, \"el_price_2018_data\",13, True, '10%', 'middle',save_figure=False)\n",
    "best_y_end, best_err_end, best_model_params_end, best_N_end, best_j_end = best_N(time_elprice_2018, el_price_2018_smoothed, \"el_price_2018_data\", 13, True, '10%', 'end', save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60a78c32",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params_testing(time_elprice_2018, best_model_params_middle, el_price_2018_smoothed, \"el_price_2018_data\", False, best_N_middle, best_j_middle, save_figure=False)\n",
    "model_params_testing(time_elprice_2018_2019, best_model_params_end, el_price_2018_2019_smoothed, \"el_price_2018_2019_data\", True, best_N_end, best_j_end, save_figure=False)\n",
    "model_params_testing(time_elprice_2018_2019[:-183], best_model_params_end, el_price_2018_2019_smoothed[:-183], \"el_price_2018_2019_data_1.5forward\", True, best_N_end, best_j_end, save_figure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "223483af",
   "metadata": {},
   "outputs": [],
   "source": [
    "#nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30bede11",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_elprice_2018, el_price_2018_data.to_numpy(),500,\"el_price_2018_data__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_elprice_2018_2019, el_price_2018_2019_data.to_numpy(), \"el_price_2018_2019_data_forward__NN\",forward_time=True, save_fig=save)\n",
    "nn_test(best_width_model, best_width, time_elprice_2018_2019[:-183], el_price_2018_2019_data.to_numpy()[:-183], \"el_price_2018_2019_data_forward_1.2forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7beb84f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "hw_list = [5,15,50,100]\n",
    "# Train models with widths in hw_list\n",
    "\n",
    "save = False\n",
    "\n",
    "mods, best_width = nn_models(hw_list, time_elprice_2018, el_price_2018_smoothed,500,\"el_price_2018_smoothed__NN\", save_fig=save)\n",
    "\n",
    "best_width_model = mods[best_width]\n",
    "\n",
    "# test forward in time\n",
    "nn_test(best_width_model, best_width, time_elprice_2018_2019, el_price_2018_2019_smoothed, \"el_price_2018_2019_smoothed_forward__NN\",forward_time=True, save_fig=save)\n",
    "nn_test(best_width_model, best_width, time_elprice_2018_2019[:-183], el_price_2018_2019_smoothed[:-183], \"el_price_2018_2019_smoothed_1.5forward__NN\",forward_time=True, save_fig=save)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a46bc956",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}