schur_regression.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyO77TLMelJMKiGB+oPcbHUr",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/microprediction/precise/blob/main/schur_regression.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "YjKZY6bcVe7J"
      },
      "outputs": [],
      "source": [
        "\n",
        "!pip install --upgrade setuptools\n",
        "!pip install --upgrade wheel\n",
        "!pip install --upgrade ecos   # <--- Try conda install ecos if this fails\n",
        "!pip install --upgrade osqp   # <-- Can be tricky on some systems see https://github.com/cvxpy/cvxpy/issues/1190#issuecomment-994613793\n",
        "!pip install --upgrade pyportfolioopt # <--- Skip if you don't plan to use it\n",
        "!pip install --upgrade riskparityportfolio\n",
        "!pip install --upgrade scipy\n",
        "!pip install --upgrade precise"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Illustrates how Schur Portfolio Theory can be used for regression"
      ],
      "metadata": {
        "id": "WyX1KwPXZz4e"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install cvxpy"
      ],
      "metadata": {
        "id": "fUyx9pIJpus2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import cvxpy as cp\n",
        "\n",
        "# Define Q (covariance matrix) and b (expected returns)\n",
        "Q = np.array([[0.04, 0.006, 0.004], [0.006, 0.09, 0.006], [0.004, 0.006, 0.025]])\n",
        "b = np.array([0.10, 0.12, 0.15])\n",
        "\n",
        "# Analytical approach\n",
        "def analytical_port_b(Q, b):\n",
        "    Q_inv = np.linalg.inv(Q)\n",
        "    lambda_val = 2 / (b.T @ Q_inv @ b)\n",
        "    w = (lambda_val / 2) * Q_inv @ b\n",
        "    return w\n",
        "\n",
        "def analytical_port_1(Q):\n",
        "    \"\"\"\n",
        "        Special case where b is vector of ones\n",
        "    \"\"\"\n",
        "    n = Q.shape[0]\n",
        "    v = np.ones(shape=(n,))\n",
        "    return analytical_port_b(Q=Q, b=v)\n",
        "\n",
        "\n",
        "def analytical_port_b_alternative(Q, b):\n",
        "    \"\"\"\n",
        "        Use a change of coordinate to reduce to special case analytical_port_1\n",
        "    \"\"\"\n",
        "    Qb = Q / np.outer(b,b)\n",
        "    w = analytical_port_1(Qb)\n",
        "    return w / b\n",
        "\n",
        "\n",
        "# Example covariance matrix Q and vector b\n",
        "Q = np.array([[0.04, 0.006, 0.004], [0.006, 0.09, 0.006], [0.004, 0.006, 0.025]])\n",
        "b = np.array([0.10, 0.12, 0.15])\n",
        "\n",
        "# Calculate weights using both methods\n",
        "w_original = analytical_port_b(Q, b)\n",
        "w_alternative = analytical_port_b_alternative(Q, b)\n",
        "\n",
        "# Comparing the results\n",
        "comparison_result = np.allclose(w_original, w_alternative, atol=1e-6)\n",
        "if not comparison_result:\n",
        "    print({'w_original':w_original,\n",
        "       'w_alternative':w_alternative})\n",
        "else:\n",
        "    print('Hoorah')\n",
        "\n",
        "\n",
        "def cvxpy_port_1(Q):\n",
        "    \"\"\"\n",
        "        Special case where b is vector of ones\n",
        "    \"\"\"\n",
        "    n = Q.shape[0]\n",
        "    v = np.ones(shape=(n,))\n",
        "    return cvxpy_port_b(Q=Q, b=v)\n",
        "\n",
        "\n",
        "\n",
        "# Using cvxpy for comparison\n",
        "def cvxpy_port_b(Q, b):\n",
        "    n = Q.shape[0]\n",
        "    w = cp.Variable(n)\n",
        "    objective = cp.Minimize(cp.quad_form(w, Q))\n",
        "    constraints = [cp.sum(cp.multiply(w, b)) == 1]\n",
        "    problem = cp.Problem(objective, constraints)\n",
        "    problem.solve()\n",
        "    return w.value\n",
        "\n",
        "# Calculate weights using both methods\n",
        "w_analytical = analytical_port_b(Q, b)\n",
        "w_cvxpy = cvxpy_port_b(Q, b)\n",
        "\n",
        "# Output the results and comparison\n",
        "print(\"Analytical weights:\", w_analytical)\n",
        "print(\"CVXPY weights:\", w_cvxpy)\n",
        "print(\"Are the weights close?\", np.allclose(w_analytical, w_cvxpy, atol=1e-6))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "gN5SEGRbp3xP",
        "outputId": "f796d65a-1bf5-46b4-a575-b0472c4ec0c0"
      },
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Hoorah\n",
            "Analytical weights: [1.64284911 0.7614118  4.96230449]\n",
            "CVXPY weights: [1.64284911 0.7614118  4.96230449]\n",
            "Are the weights close? True\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "from precise.skaters.portfoliostatic.allstaticport import LONG_PORT\n",
        "import numpy as np\n",
        "from sklearn.linear_model import LinearRegression\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "import sklearn.metrics as metrics\n",
        "from scipy.optimize import minimize\n",
        "\n",
        "def estimate_quadratic_form_with_b(X, y):\n",
        "    \"\"\"\n",
        "      Robustly estimate the expression (X^T y)^T (X^T X)^{-1} (X^T y) for given X and y.\n",
        "\n",
        "      Parameters:\n",
        "      - X: The matrix X.\n",
        "      - y: The vector y.\n",
        "\n",
        "      Returns:\n",
        "      - The estimate of the expression.\n",
        "    \"\"\"\n",
        "    # Compute b = X^T y\n",
        "    b = X.T @ y\n",
        "\n",
        "    # Compute Q = X^T X and its pseudo-inverse\n",
        "    Q = X.T @ X\n",
        "    Q_pinv = np.linalg.pinv(Q)\n",
        "\n",
        "    # Compute and return the expression\n",
        "    return b.T @ Q_pinv @ b\n",
        "\n",
        "\n",
        "def port_regression_special(X, y, port):\n",
        "    b =  X.T @ y\n",
        "    Q =  X.T @ X\n",
        "    if False:\n",
        "      w = analytical_port_b(Q = Q, b=b)\n",
        "    else:\n",
        "      Qb = Q / (b @ b.T )\n",
        "      w = port(Q) / b\n",
        "\n",
        "    inv_nu = b.T @ np.linalg.pinv(Q) @ b\n",
        "    beta = inv_nu * w\n",
        "    return beta\n",
        "\n",
        "\n",
        "def squared_error(lambda_, X, y, beta):\n",
        "    \"\"\"\n",
        "    Calculate the squared error of the predictions with rescaled coefficients.\n",
        "    \"\"\"\n",
        "    predictions = X @ (beta * lambda_)\n",
        "    return np.sum((y - predictions) ** 2)\n",
        "\n",
        "\n",
        "def find_optimal_lambda(X, y, beta):\n",
        "    \"\"\"\n",
        "    Find the lambda that minimizes the squared error.\n",
        "    \"\"\"\n",
        "    result = minimize(squared_error, x0=1.0, args=(X, y, beta))\n",
        "    return result.x[0]\n",
        "\n",
        "\n",
        "def port_regression(X,y, port, line_search=True):\n",
        "    \"\"\"\n",
        "\n",
        "        Demonstrates how to use any portfolio method for regression\n",
        "\n",
        "    \"\"\"\n",
        "    X = X - np.mean(X, axis=0)\n",
        "    std_devs = np.std(X, axis=0)\n",
        "    std_devs[std_devs ==0] = 1\n",
        "    X_scaled = X / std_devs\n",
        "\n",
        "    beta_scaled = port_regression_special(X_scaled, y, port)\n",
        "\n",
        "    beta = beta_scaled / std_devs\n",
        "\n",
        "    if line_search:\n",
        "       lambda_ = find_optimal_lambda(X_scaled, y, beta)\n",
        "       beta = beta*lambda_\n",
        "\n",
        "    return beta\n",
        "\n",
        "\n",
        "def schur_regression(X,y):\n",
        "    from precise.skaters.portfoliostatic.schurport import schur_unit_weak_s5_g100_port as s5\n",
        "    return port_regression(X=X, y=y, port=s5)\n",
        "\n",
        "def hrp_regression(X,y):\n",
        "    from precise.skaters.portfoliostatic.hrpport import hrp_diag_diag_s5_long_port as h5\n",
        "    return port_regression(X=X, y=y, port=h5)\n",
        "\n",
        "def unit_port_pure(cov):\n",
        "    return unitary_from_pre(np.linalg.inv(cov))\n",
        "\n",
        "def unitary_from_pre(pre, phi=1.0):\n",
        "    \"\"\" Signed min var portfolio summing to unity \"\"\"\n",
        "    # No optimization required\n",
        "    n_dim = np.shape(pre)[1]\n",
        "    wones = np.ones(shape=(n_dim, 1))\n",
        "    from precise.skaters.locationutil.vectorfunctions import normalize\n",
        "    w = normalize(np.squeeze(np.matmul(pre, wones)))\n",
        "    sum_w = np.sum(w)\n",
        "    if abs(sum_w-1)<1e-4:\n",
        "        return np.array(w)\n",
        "    else:\n",
        "        print('Warning. Unitary portfolio failed.')\n",
        "        return np.ones_like(w)/n_dim\n",
        "\n",
        "\n",
        "def unit_regression(X,y):\n",
        "    from precise.skaters.portfoliostatic.unitport import unit_port as up\n",
        "    return port_regression_special(X=X, y=y, port=unit_port_pure)\n",
        "\n",
        "import numpy as np\n",
        "\n",
        "def scaled_regression(X,y):\n",
        "    \"\"\"\n",
        "        OLS with scaling\n",
        "    \"\"\"\n",
        "    X = X - np.mean(X, axis=0)\n",
        "    std_devs = np.std(X, axis=0)\n",
        "    std_devs[std_devs ==0] = 1\n",
        "    X_scaled = X / std_devs\n",
        "    beta_scaled = np.linalg.inv(X_scaled.T @ X_scaled) @ X_scaled.T @ y\n",
        "    beta = beta_scaled / std_devs\n",
        "    return beta\n",
        "\n",
        "def test_scaled_regression():\n",
        "    np.random.seed(0)\n",
        "    X = np.random.randn(1000, 50)\n",
        "    beta_true = np.array([1, 1, 1, 1, -1, 2, 2, 2, 2, 2]*5)\n",
        "    y = X @ beta_true + np.random.randn(1000) * 0.5\n",
        "\n",
        "    beta_custom = scaled_regression(X, y)\n",
        "    beta_port = unit_regression(X, y)\n",
        "    beta_hrp = hrp_regression(X,y)\n",
        "    beta_schur = schur_regression(X,y)\n",
        "\n",
        "    scaler = StandardScaler()\n",
        "    X_scaled = scaler.fit_transform(X)\n",
        "    model = LinearRegression(fit_intercept=False)\n",
        "    model.fit(X_scaled, y)\n",
        "    beta_sklearn = model.coef_ / scaler.scale_\n",
        "\n",
        "    print(\"Custom implementation coefficients:\", beta_custom)\n",
        "    print(\"Port regression special coefficients:\", beta_port)\n",
        "    print(\"scikit-learn implementation coefficients:\", beta_sklearn)\n",
        "    print(\"Schur coefficients:\", beta_schur)\n",
        "    print('HRP coefficients:', beta_hrp)\n",
        "\n",
        "    assert np.allclose(beta_custom, beta_sklearn, atol=1e-6), \"The coefficients do not match closely enough.\"\n",
        "\n",
        "    print(\"Test passed: The coefficients match closely between the custom implementation, port_regression_special, and scikit-learn's.\")\n",
        "\n",
        "test_scaled_regression()"
      ],
      "metadata": {
        "id": "8B31yjVTZkLl",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "9f36141b-e01c-470f-ede9-66afd4719434"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Custom implementation coefficients: [ 1.00394404  1.00101244  0.98604402  0.99411553 -1.00779756  2.02206025\n",
            "  1.99890455  2.01712081  1.99839708  2.02811794  0.98476052  1.01679591\n",
            "  1.02976976  0.9915431  -0.96640799  1.99517801  1.99337859  2.01313501\n",
            "  1.987767    2.0096497   0.98144382  1.00055326  1.00462235  0.98135243\n",
            " -0.97076441  1.99569078  1.97222357  2.02214309  2.01916881  2.0046913\n",
            "  0.98070427  0.99183814  1.01016456  0.98930339 -1.00618418  1.99341331\n",
            "  2.02156476  2.0054314   2.00081339  1.97143532  1.00290907  0.9917496\n",
            "  0.99776909  0.97276666 -1.00624605  2.0166339   2.01275879  1.98304077\n",
            "  2.00271798  2.0072813 ]\n",
            "Port regression special coefficients: [ 4.46338197 19.54082224  2.43600167  0.97207529 -1.95083178  2.50991688\n",
            "  0.72087277  0.9779235   2.12461801  1.80481131  2.28203133  1.78702586\n",
            "  1.76189318  3.14697399 -2.72344694  0.93155611  0.82338879  1.67662111\n",
            "  1.826796    1.2175634   2.01600797  2.03934161  1.94067613  0.88730001\n",
            " -2.77285038  2.85275244  0.61859678  0.97036272  1.69154151  1.38391542\n",
            "  1.0448419   1.66509204  1.21889236  1.70378682 -2.98451065  0.84923986\n",
            "  2.29483783  1.28563161  0.72145235  0.89240931  2.13951342  2.87589934\n",
            "  3.06024907  4.11702605 -2.62540832  2.06415536  1.40117481  1.97646201\n",
            "  1.39261067  2.24386753]\n",
            "scikit-learn implementation coefficients: [ 1.00394404  1.00101244  0.98604402  0.99411553 -1.00779756  2.02206025\n",
            "  1.99890455  2.01712081  1.99839708  2.02811794  0.98476052  1.01679591\n",
            "  1.02976976  0.9915431  -0.96640799  1.99517801  1.99337859  2.01313501\n",
            "  1.987767    2.0096497   0.98144382  1.00055326  1.00462235  0.98135243\n",
            " -0.97076441  1.99569078  1.97222357  2.02214309  2.01916881  2.0046913\n",
            "  0.98070427  0.99183814  1.01016456  0.98930339 -1.00618418  1.99341331\n",
            "  2.02156476  2.0054314   2.00081339  1.97143532  1.00290907  0.9917496\n",
            "  0.99776909  0.97276666 -1.00624605  2.0166339   2.01275879  1.98304077\n",
            "  2.00271798  2.0072813 ]\n",
            "Schur coefficients: [ 1.5724316   3.4613774   1.05503369  0.56880521 -1.10433689  0.78972225\n",
            "  0.34087852  0.2949954   0.22942507  0.56448177  0.813367    0.74572824\n",
            "  0.86451896  1.09108956 -0.71425911  0.5432924   0.31658495  0.60288936\n",
            "  0.6883931   0.41828572  0.79291325  1.0102172   0.81083353  0.55257468\n",
            " -2.93601587  0.98144051  0.3499671   0.31606223  0.51076867  0.48559068\n",
            "  0.31573976  0.90073354  0.49938143  0.73320764 -2.01331845  0.31646474\n",
            "  0.44692304  0.4631825   0.26836566  0.34236191  0.70106617  1.19332079\n",
            "  0.72649329  1.14037177 -1.27290563  0.47593914  0.26566713  0.33867174\n",
            "  0.49284868  0.68593187]\n",
            "HRP coefficients: [ 1.32759671  4.34753561  0.84332339  0.76193132 -0.90503266  0.69012824\n",
            "  0.40914892  0.47291593  0.63387188  0.45705633  1.1588359   0.74984858\n",
            "  0.99657267  1.18711717 -1.96417688  0.44103879  0.57717381  0.50802365\n",
            "  0.91535127  0.46237184  0.84308098  0.87086725  0.90113741  0.55004213\n",
            " -1.26308378  1.03507658  0.2143655   0.83844186  0.39894913  0.81944263\n",
            "  0.57712489  1.32072232  0.4815809   1.11317623 -1.24730975  0.67372384\n",
            "  0.46458964  0.70489529  0.23692178  0.54816551  0.4364378   1.82245703\n",
            "  0.7926926   2.18482109 -0.85577614  0.98550698  0.42667327  0.80393637\n",
            "  0.38962683  0.99493901]\n",
            "Test passed: The coefficients match closely between the custom implementation, port_regression_special, and scikit-learn's.\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Testing out of sample"
      ],
      "metadata": {
        "id": "MUUZBs-ZYKGZ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "import numpy as np\n",
        "from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "from sklearn.metrics import mean_squared_error\n",
        "from scipy.optimize import nnls\n",
        "from sklearn.datasets import make_spd_matrix\n",
        "\n",
        "def nnls_regression(X, y):\n",
        "    beta_nnls, rnorm = nnls(X, y)\n",
        "    return beta_nnls\n",
        "\n",
        "def test_scaled_regression():\n",
        "    np.random.seed(0)\n",
        "    # Generate correlated X\n",
        "    n_samples =  4000\n",
        "    n_features = 20\n",
        "    n_informative = 20\n",
        "    beta_true = np.concatenate([\n",
        "    np.random.normal(loc=0.0, scale=0.1, size=n_features-n_informative),\n",
        "    np.random.normal(loc=1.5, scale=0.1, size=n_informative)\n",
        "    ])\n",
        "\n",
        "    # Create a covariance matrix\n",
        "    n_samples = 400\n",
        "\n",
        "\n",
        "    # Initialize the covariance matrix\n",
        "    cov_matrix = np.zeros((n_features, n_features))\n",
        "\n",
        "    # Fill the covariance matrix with values that decrease with distance\n",
        "    n_blocks = 1  # Define the number of blocks\n",
        "\n",
        "    # Ensure each block has an equal number of features for simplicity\n",
        "    block_size = n_features // n_blocks\n",
        "    standard_deviations = np.linspace(0.1, 0.1, n_features)\n",
        "\n",
        "    # Initialize the covariance matrix\n",
        "    cov_matrix = make_spd_matrix(n_features)\n",
        "    cov_matrix = (cov_matrix + cov_matrix.T) / 2\n",
        "\n",
        "\n",
        "    # Generate X data with hierarchical covariance structure\n",
        "    X = np.random.multivariate_normal(mean=np.zeros(n_features), cov=cov_matrix, size=n_samples)\n",
        "\n",
        "    # Ensure the matrix is symmetric and positive-definite\n",
        "    np.fill_diagonal(cov_matrix, standard_deviations**2)\n",
        "    cov_matrix = (cov_matrix + cov_matrix.T) / 2\n",
        "    X = np.random.multivariate_normal(mean=np.zeros(n_features), cov=cov_matrix, size=n_samples)\n",
        "\n",
        "\n",
        "    # Generate y with some noise\n",
        "    y = X @ beta_true + np.random.randn(n_samples) * 5.\n",
        "\n",
        "    # Split the data into training and testing sets\n",
        "    train_proportion = 0.05\n",
        "    split_idx = int(n_samples * train_proportion)\n",
        "    X_train, X_test = X[:split_idx], X[split_idx:]\n",
        "    y_train, y_test = y[:split_idx], y[split_idx:]\n",
        "\n",
        "    # Fit models on the training set\n",
        "    beta_custom = scaled_regression(X_train, y_train)\n",
        "    beta_port = unit_regression(X_train, y_train)\n",
        "    beta_hrp = hrp_regression(X_train, y_train)\n",
        "    beta_schur = schur_regression(X_train, y_train)\n",
        "    beta_nnls = nnls_regression(X_train, y_train)\n",
        "\n",
        "    # Regularized regressions\n",
        "    ridge_model = Ridge(alpha=1.0, fit_intercept=False)\n",
        "    ridge_model.fit(X_train, y_train)\n",
        "    beta_ridge = ridge_model.coef_\n",
        "\n",
        "    lasso_model = Lasso(alpha=0.1, fit_intercept=False)\n",
        "    lasso_model.fit(X_train, y_train)\n",
        "    beta_lasso = lasso_model.coef_\n",
        "\n",
        "    # Scale X for comparison with scikit-learn\n",
        "    scaler = StandardScaler()\n",
        "    X_train_scaled = scaler.fit_transform(X_train)\n",
        "    model = LinearRegression(fit_intercept=False)\n",
        "    model.fit(X_train_scaled, y_train)\n",
        "    beta_sklearn = model.coef_ / scaler.scale_\n",
        "\n",
        "    # Evaluate model performance on the test set\n",
        "    y_pred_sklearn = X_test @ beta_sklearn\n",
        "    y_pred_custom = X_test @ beta_custom\n",
        "    y_pred_port = X_test @ beta_port\n",
        "    y_pred_hrp = X_test @ beta_hrp\n",
        "    y_pred_schur = X_test @ beta_schur\n",
        "    y_pred_nnls = X_test @ beta_nnls\n",
        "\n",
        "    # Calculate and print out-of-sample mean squared error\n",
        "    mse_zero = mean_squared_error(y_test,np.zeros_like(y_test))\n",
        "    mse_custom = mean_squared_error(y_test, y_pred_custom)\n",
        "    mse_sklearn = mean_squared_error(y_test, y_pred_sklearn)\n",
        "    mse_port = mean_squared_error(y_test, y_pred_port)\n",
        "    mse_hrp = mean_squared_error(y_test, y_pred_hrp)\n",
        "    mse_schur = mean_squared_error(y_test, y_pred_schur)\n",
        "    mse_nnls = mean_squared_error(y_test, y_pred_nnls)\n",
        "    mse_ridge = mean_squared_error(y_test, X_test @ beta_ridge)\n",
        "    mse_lasso = mean_squared_error(y_test, X_test @ beta_lasso)\n",
        "\n",
        "\n",
        "    mse_dict = {\n",
        "    'Sklearn lin reg      ': mse_sklearn,\n",
        "    'Zero prediction      ': mse_zero,\n",
        "    'Ridge                ': mse_ridge,\n",
        "    'Lasso                ': mse_lasso,\n",
        "    'Custom Implementation': mse_custom,\n",
        "    'Port Regression.     ': mse_port,\n",
        "    'HRP Regression.      ': mse_hrp,\n",
        "    'Schur Regression.    ': mse_schur,\n",
        "    'NNLS Regression.     ': mse_nnls\n",
        "    }\n",
        "\n",
        "    # Sort the MSE values from lowest to highest\n",
        "    sorted_mse = sorted(mse_dict.items(), key=lambda item: item[1])\n",
        "\n",
        "    # Print the sorted MSE values\n",
        "    print(\"MSE values sorted from lowest to highest:\")\n",
        "    for name, mse in sorted_mse:\n",
        "        print(f\"{name} MSE: {mse}\")\n",
        "\n",
        "\n",
        "    # Add assertions or comparisons as needed, for example:\n",
        "    if not np.allclose(beta_custom, beta_sklearn, atol=1e-4):\n",
        "         print(\"The coefficients do not match closely enough between sklearn and custom regression\")\n",
        "\n",
        "    print(\"Test passed: The coefficients match closely between the custom implementation and scikit-learn's.\")\n",
        "\n",
        "\n",
        "    import matplotlib.pyplot as plt\n",
        "    import seaborn as sns\n",
        "    plt.figure(figsize=(10, 8))\n",
        "    corr_matrix = np.zeros_like(cov_matrix)\n",
        "\n",
        "    # Convert covariance matrix to correlation matrix\n",
        "    for i in range(n_features):\n",
        "        for j in range(n_features):\n",
        "            corr_matrix[i, j] = cov_matrix[i, j] / (np.sqrt(cov_matrix[i, i] * cov_matrix[j, j]))\n",
        "\n",
        "    sns.heatmap(corr_matrix, cmap='coolwarm', cbar_kws={'label': 'Correlation'})\n",
        "    plt.title('Covariance Matrix')\n",
        "    plt.xlabel('Feature Index')\n",
        "    plt.ylabel('Feature Index')\n",
        "    plt.show()\n",
        "\n",
        "\n",
        "# Call the test function\n",
        "test_scaled_regression()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 961
        },
        "id": "cD3IWsejZlCe",
        "outputId": "4c6605c8-5bea-434a-f71a-9f9c1e56b972"
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "<ipython-input-3-635c1080d43b>:48: RuntimeWarning: covariance is not symmetric positive-semidefinite.\n",
            "  X = np.random.multivariate_normal(mean=np.zeros(n_features), cov=cov_matrix, size=n_samples)\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "MSE values sorted from lowest to highest:\n",
            "Ridge                 MSE: 69.6753566554431\n",
            "NNLS Regression.      MSE: 81.81922397054139\n",
            "Lasso                 MSE: 100.12337380064562\n",
            "Zero prediction       MSE: 113.50160175440578\n",
            "Schur Regression.     MSE: 130.02694734922122\n",
            "HRP Regression.       MSE: 134.14715903964054\n",
            "Sklearn lin reg       MSE: 398.2385706677219\n",
            "Port Regression.      MSE: 2376.887906890707\n",
            "Custom Implementation MSE: 111077.8474894728\n",
            "The coefficients do not match closely enough between sklearn and custom regression\n",
            "Test passed: The coefficients match closely between the custom implementation and scikit-learn's.\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1000x800 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzMAAAK9CAYAAAAQZjKqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB5VUlEQVR4nO3deZyN5f/H8fc92xmGGTtjG8vIkp0WSlORoZJpkUoh5fsllb2ahNS3JopIRVKWSqKvChUhRJEaIUK27PsyY5szy7l/f/Q1vzlhzGFmrrlnXs/H437U3Pd9rvd15swc8znXdV+3Zdu2LQAAAABwGD/THQAAAACAy0ExAwAAAMCRKGYAAAAAOBLFDAAAAABHopgBAAAA4EgUMwAAAAAciWIGAAAAgCNRzAAAAABwJIoZAAAAAI5EMQMAV6hKlSrq2rWr6W7kS3/99Zcsy9LkyZNNdwUAkAdRzADIE7Zt26Z///vfqlatmoKDgxUaGqobbrhBY8aM0dmzZ013r0CwLEuWZenxxx+/4PFBgwaln3PkyBGf2//mm2/04osvXmEvAQD4f5Zt27bpTgAo2L7++mt16NBBLpdLnTt3Vt26dZWcnKzly5frv//9r7p27aoJEyaY7uZFud1u+fn5KTAw0HRXrohlWQoODlZwcLAOHjyooKAgr+PVqlXT/v37lZSUpMOHD6tUqVI+tf/kk0/qnXfekS//7Ni2LbfbrcDAQPn7+/uUBwDI/xiZAWDUjh079MADDygiIkJ//PGHxowZo+7du6tXr1769NNP9ccff+jqq6823c3z2LadPmLkcrkcX8ic06ZNGyUmJurbb7/12v/TTz9px44duuOOO3KlH6mpqUpOTk4vsChkAAAXQjEDwKgRI0bo1KlT+uCDDxQeHn7e8cjISPXu3Tv969TUVL388suqXr26XC6XqlSpoueff15utzv9nDvvvFPVqlW7YF6zZs3UtGnT9K8nTZqkW2+9VWXKlJHL5VKdOnU0bty48x5XpUoV3XnnnZo/f76aNm2qQoUK6b333ks/lvGamWPHjmnAgAGqV6+eihQpotDQULVt21Zr1671anPJkiWyLEszZszQK6+8oooVKyo4OFgtW7bU1q1bz+vDzz//rNtvv13FixdXSEiI6tevrzFjxnids2nTJt13330qUaKEgoOD1bRpU82ePfuC34sLqVChgm666SZNmzbNa/8nn3yievXqqW7duuc9ZtmyZerQoYMqV64sl8ulSpUqqW/fvl7TA7t27ap33nlH0v9PZ7MsS9L/XxfzxhtvaPTo0emv7R9//HHeNTOHDh1S6dKldfPNN3uN8GzdulUhISHq2LFjlp8rAMD5Akx3AEDBNmfOHFWrVk3NmzfP0vmPP/64pkyZovvuu0/9+/fXzz//rLi4OG3cuFFffPGFJKljx47q3LmzfvnlF11zzTXpj925c6dWrlyp119/PX3fuHHjdPXVV+uuu+5SQECA5syZoyeeeEIej0e9evXyyt68ebMefPBB/fvf/1b37t1Vs2bNC/Zx+/bt+vLLL9WhQwdVrVpVBw8e1HvvvaeoqCj98ccfKl++vNf5r732mvz8/DRgwAAlJCRoxIgR6tSpk37++ef0cxYsWKA777xT4eHh6t27t8qVK6eNGzdq7ty56cXehg0bdMMNN6hChQp67rnnFBISohkzZigmJkb//e9/dffdd2fpe/zQQw+pd+/eOnXqlIoUKaLU1FTNnDlT/fr1U1JS0nnnz5w5U2fOnFHPnj1VsmRJrVq1SmPHjtWePXs0c+ZMSdK///1v7du3TwsWLNBHH310wdxJkyYpKSlJ//rXv+RyuVSiRAl5PB6vc8qUKaNx48apQ4cOGjt2rJ5++ml5PB517dpVRYsW1bvvvpul5wgAyCdsADAkISHBlmS3b98+S+evWbPGlmQ//vjjXvsHDBhgS7K///779HZdLpfdv39/r/NGjBhhW5Zl79y5M33fmTNnzsuJjo62q1Wr5rUvIiLClmTPmzfvvPMjIiLsLl26pH+dlJRkp6WleZ2zY8cO2+Vy2S+99FL6vsWLF9uS7Nq1a9tutzt9/5gxY2xJ9u+//27btm2npqbaVatWtSMiIuzjx497tevxeNL/v2XLlna9evXspKQkr+PNmze3a9SocV6//0mS3atXL/vYsWN2UFCQ/dFHH9m2bdtff/21bVmW/ddff9lDhw61JdmHDx9Of9yFvodxcXHnfa979eplX+ifnR07dtiS7NDQUPvQoUMXPDZp0iSv/Q8++KBduHBh+88//7Rff/11W5L95ZdfXvI5AgDyF6aZATAmMTFRklS0aNEsnf/NN99Ikvr16+e1v3///pL+XkhAUvq0rhkzZnhNRfrss890/fXXq3Llyun7ChUqlP7/CQkJOnLkiKKiorR9+3YlJCR45VStWlXR0dGX7KfL5ZKf399vr2lpaTp69KiKFCmimjVravXq1eed/+ijj3pdbN+iRQtJf4/wSNJvv/2mHTt2qE+fPipWrJjXY89N1Tp27Ji+//573X///Tp58qSOHDmiI0eO6OjRo4qOjtaWLVu0d+/eS/ZdkooXL642bdro008/lSRNmzZNzZs3V0RExAXPz/g9PH36tI4cOaLmzZvLtm399ttvWcqUpHvvvVelS5fO0rlvv/22wsLCdN9992nw4MF65JFH1L59+yxnAQDyB4oZAMaEhoZKkk6ePJml83fu3Ck/Pz9FRkZ67S9XrpyKFSumnTt3pu/r2LGjdu/erRUrVkj6e+nn+Pj4866p+PHHH9WqVSuFhISoWLFiKl26tJ5//nlJumAxkxUej0dvvvmmatSoIZfLpVKlSql06dJat27deW1K8iqupL+LCUk6fvx4et8lXfB6lXO2bt0q27Y1ePBglS5d2msbOnSopL+vN8mqhx56SAsWLNCuXbv05Zdf6qGHHrroubt27VLXrl1VokQJFSlSRKVLl1ZUVJSk87+Hmcnq91eSSpQoobfeekvr1q1TWFiY3nrrrSw/FgCQf3DNDABjQkNDVb58ea1fv96nx50bjchMu3btVLhwYc2YMUPNmzfXjBkz5Ofnpw4dOqSfs23bNrVs2VK1atXSqFGjVKlSJQUFBembb77Rm2++ed71GhlHIDLz6quvavDgwerWrZtefvlllShRQn5+furTp895bUq66Epdtg9LGJ9rd8CAARcdPfpnEZiZu+66Sy6XS126dJHb7db9999/wfPS0tJ022236dixY3r22WdVq1YthYSEaO/everatesFn+/FZPX7e878+fMl/V307dmz57xRKwBA/kcxA8CoO++8UxMmTNCKFSvUrFmzTM+NiIiQx+PRli1bVLt27fT9Bw8e1IkTJ7ymQYWEhOjOO+/UzJkzNWrUKH322Wdq0aKF18X3c+bMkdvt1uzZs71GRxYvXnxFz+nzzz/XLbfcog8++MBr/4kTJ3y+N4skVa9eXZK0fv16tWrV6oLnnFu9LTAw8KLn+KJQoUKKiYnRxx9/rLZt216037///rv+/PNPTZkyRZ07d07fv2DBgvPOzUoRmlXz5s3TxIkT9cwzz+iTTz5Rly5d9PPPPysggH/WAKAgYZoZAKOeeeYZhYSE6PHHH9fBgwfPO75t27b05Ydvv/12SdLo0aO9zhk1apQknXcPlI4dO2rfvn2aOHGi1q5de94Us3MjIhlHQBISEjRp0qQrek7+/v7njarMnDkzy9es/FPjxo1VtWpVjR49WidOnPA6di6nTJkyuvnmm/Xee+9p//7957Vx+PBhn3MHDBigoUOHavDgwRc950LfQ9u2z1syWvq7wJR03nPw1YkTJ/T444/r2muv1auvvqqJEydq9erVevXVV6+oXQCA8/ARFgCjqlevrmnTpqljx46qXbu2OnfurLp16yo5OVk//fSTZs6cmX4PlwYNGqhLly6aMGGCTpw4oaioKK1atUpTpkxRTEyMbrnlFq+2b7/9dhUtWlQDBgyQv7+/7r33Xq/jrVu3VlBQkNq1a6d///vfOnXqlN5//32VKVPmggVBVt1555166aWX9Oijj6p58+b6/fff9cknn1z03jeX4ufnp3Hjxqldu3Zq2LChHn30UYWHh2vTpk3asGFD+nSrd955RzfeeKPq1aun7t27q1q1ajp48KBWrFihPXv2nHefm0tp0KCBGjRokOk5tWrVUvXq1TVgwADt3btXoaGh+u9//5t+vU9GTZo0kSQ9/fTTio6Olr+/vx544AGf+iRJvXv31tGjR7Vw4UL5+/urTZs2evzxx/Wf//xH7du3v2SfAQD5iLF11AAggz///NPu3r27XaVKFTsoKMguWrSofcMNN9hjx471Wmo4JSXFHjZsmF21alU7MDDQrlSpkh0bG+t1TkadOnWyJdmtWrW64PHZs2fb9evXt4ODg+0qVarYw4cPtz/88ENbkr1jx4708yIiIuw77rjjgm1caGnm/v372+Hh4XahQoXsG264wV6xYoUdFRVlR0VFpZ93bmnmmTNnerV3seWIly9fbt9222120aJF7ZCQELt+/fr22LFjvc7Ztm2b3blzZ7tcuXJ2YGCgXaFCBfvOO++0P//88wv2PSP9b2nmzFxoaeY//vjDbtWqlV2kSBG7VKlSdvfu3e21a9ee9xxSU1Ptp556yi5durRtWVb6Ms3nnu/rr79+Xt4/vxdfffWVLckeOXKk13mJiYl2RESE3aBBAzs5OfmSzxUAkD9Ytu3DFaYAAAAAkEdwzQwAAAAAR6KYAQAAAOBIFDMAAAAAHIliBgAAAIAjUcwAAAAAcCSKGQAAAACORDEDAAAAwJECTHcgJ5xc9bWx7J+DWxnLbmz9Yix7ubuZsezGoRuNZbtSzxjLPukqaSz7SEopY9m1Tq80lh2YcMhY9srS9xnLvv7ol8ayE8teZSw7JHGfsezkBXONZQe1usNYdsBRc9/zLeM+M5bd03rRWHbMoy2MZQ+4J+9+pv51YE1j2XekbDaW7UR596cIAAAAADKRL0dmAAAAgMtlBVqmu4AsYmQGAAAAgCNRzAAAAABwJKaZAQAAABn4BTDNzCkYmQEAAADgSIzMAAAAABlYgXze7xS8UgAAAAAciZEZAAAAIAOumXEORmYAAAAAOBLFDAAAAABHYpoZAAAAkIEVyDQzpzBazBw5ckQffvihVqxYoQMHDkiSypUrp+bNm6tr164qXbq0ye4BAAAAyMOMFTO//PKLoqOjVbhwYbVq1UpXXXWVJOngwYN666239Nprr2n+/Plq2rRppu243W653W6vfcnJKXIFBeZY3wEAAJB/sQCAcxgrZp566il16NBB48ePl2V5/8DYtq0ePXroqaee0ooVKzJtJy4uTsOGDfPa99zjD+r57p2yvc8AAAAA8g5jxczatWs1efLk8woZSbIsS3379lWjRo0u2U5sbKz69evntS953ffZ1k8AAAAAeZOxYqZcuXJatWqVatWqdcHjq1atUtmyZS/Zjsvlksvl8tp3kilmAAAAuEwsAOAcxoqZAQMG6F//+pfi4+PVsmXL9MLl4MGDWrRokd5//3298cYbproHAAAAII8zVsz06tVLpUqV0ptvvql3331XaWlpkiR/f381adJEkydP1v3332+qewAAACigWADAOYwuzdyxY0d17NhRKSkpOnLkiCSpVKlSCgxkmhgAAACAzOWJm2YGBgYqPDzcdDcAAAAAWf6MzDiFn+kOAAAAAMDloJgBAAAA4Eh5YpoZAAAAkFf4Mc3MMRiZAQAAAOBIjMwAAAAAGVh+jMw4BSMzAAAAAByJYgYAAACAIzHNDAAAAMjA8ufzfqfglQIAAADgSJZt27bpTmS3hevcxrKvS1poLNvffcZYtt+BXcay7ZJlzGW7ChvLtgy+3tbZ08ayU0tVMJbtCQw2lh20f5ux7OTw6sayA04dN5adFhxiLPtUqLmf8yKJe41lpwWZ+54fKlrNWHbRd58zl129srHsIk+8Ziz7UlZed62x7Ot/XmUs24kYmQEAAADgSFwzAwAAAGTA0szOwcgMAAAAAEeimAEAAADgSEwzAwAAADLw82eamVMwMgMAAADAkRiZAQAAADKwGJlxDEZmAAAAADgSxQwAAAAAR2KaGQAAAJCB5cfn/U7BKwUAAADAkRiZAQAAADKw/FgAwCkYmQEAAADgSHm6mNm9e7e6deuW6Tlut1uJiYleW3KyO5d6CAAAgPzGz98ytsE3ebqYOXbsmKZMmZLpOXFxcQoLC/Papn8wIpd6CAAAAMAUo9fMzJ49O9Pj27dvv2QbsbGx6tevn9e+5X9eUbcAAAAAOIDRYiYmJkaWZcm27YueY1mZD7e5XC65XC6vfUFBTDMDAADA5WEBAOcwOs0sPDxcs2bNksfjueC2evVqk90DAAAAkIcZLWaaNGmi+Pj4ix6/1KgNAAAAkN0sPz9jG3xjdJrZwIEDdfr06Ysej4yM1OLFi3OxRwAAAACcwmgx06JFi0yPh4SEKCoqKpd6AwAAAMBJjBYzAAAAQF7DAgDOwcQ8AAAAAI7EyAwAAACQgZ8/IzNOwcgMAAAAAEdiZAYAAADIgGtmnIORGQAAAACORDEDAAAAwJGYZgYAAABkYPnxeb9T5MtiprH1i7Fsf/cZY9lprsLGsu2K1Y1le/wDjWWb5GeZm89rGfxZM8kv+ayx7LTiZYxlW2kpxrI9QcHGsi3bYyw75MxhY9kmn7d/8mlj2WUT/jSWHdSunbFsOzDIWDaQHfJlMQMAAABcLhYAcA7G0AAAAAA4EsUMAAAAAEdimhkAAACQAdPMnIORGQAAAACOxMgMAAAAkAEjM87ByAwAAAAAR6KYAQAAAOBITDMDAAAAMrD8+LzfKXilAAAAADgSIzMAAABABn7+LADgFIzMAAAAAHAkRmYAAACADFia2TkYmQEAAADgSMaLmbNnz2r58uX6448/zjuWlJSkqVOnZvp4t9utxMREr82dnJxT3QUAAACQRxgtZv7880/Vrl1bN910k+rVq6eoqCjt378//XhCQoIeffTRTNuIi4tTWFiY1zZ64sc53XUAAADkU5afn7ENvjH6HXv22WdVt25dHTp0SJs3b1bRokV1ww03aNeuXVluIzY2VgkJCV5bn8cfzsFeAwAAAMgLjC4A8NNPP2nhwoUqVaqUSpUqpTlz5uiJJ55QixYttHjxYoWEhFyyDZfLJZfL5bUvNSgop7oMAACAfI4FAJzD6MjM2bNnFRDw//WUZVkaN26c2rVrp6ioKP35558GewcAAAAgLzM6MlOrVi39+uuvql27ttf+t99+W5J01113megWAAAAAAcwOjJz991369NPP73gsbffflsPPvigbNvO5V4BAACgILP8LGMbfGO0mImNjdU333xz0ePvvvuuPB5PLvYIAAAAgFMYnWYGAAAA5DUskewcvFIAAAAAHImRGQAAACADrl1xDkZmAAAAADgSxQwAAAAAR2KaGQAAAJABCwA4B68UAAAAAEdiZAYAAADIyGIBAKfIl8XMcnczY9mtD4w2lm1XrG4s2+MfaCw74NQxY9me4BBj2f6nEoxl6+xpY9Fp5asZy7b9/I1lB/21zVh2WtW6xrL9zyQay/a4ChvLPhVS1lh20cQ9xrLTAoKNZR8KjTSWHfTuWGPZpRrVNJatG+41l418g2lmAAAAABwpX47MAAAAAJeL+8w4ByMzAAAAAByJkRkAAAAgA5Zmdg5eKQAAAACOxMgMAAAAkAHXzDgHIzMAAAAAHIliBgAAAIAjMc0MAAAAyIAFAJyDVwoAAACAIzEyAwAAAGTAAgDOwcgMAAAAAEeimAEAAADgSEwzAwAAADJgmplzGC9mNm7cqJUrV6pZs2aqVauWNm3apDFjxsjtduvhhx/Wrbfemunj3W633G63176U5AAFBrlystsAAAAADDM6zWzevHlq2LChBgwYoEaNGmnevHm66aabtHXrVu3cuVOtW7fW999/n2kbcXFxCgsL89o+n/xaLj0DAAAA5Dt+fuY2+MTod+yll17SwIEDdfToUU2aNEkPPfSQunfvrgULFmjRokUaOHCgXnst88IkNjZWCQkJXtt9XZ/LpWcAAAAAwBSjxcyGDRvUtWtXSdL999+vkydP6r777ks/3qlTJ61bty7TNlwul0JDQ702ppgBAADgclmWZWyDb4yPZZ170fz8/BQcHKywsLD0Y0WLFlVCQoKprgEAAAB5VlxcnK655hoVLVpUZcqUUUxMjDZv3ux1TlJSknr16qWSJUuqSJEiuvfee3Xw4EGvc3bt2qU77rhDhQsXVpkyZTRw4EClpqbm5lO5bEaLmSpVqmjLli3pX69YsUKVK1dO/3rXrl0KDw830TUAAAAgT1u6dKl69eqllStXasGCBUpJSVHr1q11+vTp9HP69u2rOXPmaObMmVq6dKn27dune+65J/14Wlqa7rjjDiUnJ+unn37SlClTNHnyZA0ZMsTEU/KZ0dXMevbsqbS0tPSv69at63X822+/veRqZgAAAEB2shxyIf68efO8vp48ebLKlCmj+Ph43XTTTUpISNAHH3ygadOmpf9NPWnSJNWuXVsrV67U9ddfr++++05//PGHFi5cqLJly6phw4Z6+eWX9eyzz+rFF19UUFCQiaeWZUaLmR49emR6/NVXX82lngAAAADmXei2Iy6XSy7Xpa8JP3d5RokSJSRJ8fHxSklJUatWrdLPqVWrlipXrqwVK1bo+uuv14oVK1SvXj2VLVs2/Zzo6Gj17NlTGzZsUKNGjbLjaeUYZ5SdAAAAQC6x/Cxj24VuOxIXF3fJPns8HvXp00c33HBD+mynAwcOKCgoSMWKFfM6t2zZsjpw4ED6ORkLmXPHzx3L64zfNBMAAADA32JjY9WvXz+vfVkZlenVq5fWr1+v5cuX51TX8iSKGQAAACCPyOqUsoyefPJJzZ07Vz/88IMqVqyYvr9cuXJKTk7WiRMnvEZnDh48qHLlyqWfs2rVKq/2zq12du6cvIxpZgAAAEBGfn7mNh/Ytq0nn3xSX3zxhb7//ntVrVrV63iTJk0UGBioRYsWpe/bvHmzdu3apWbNmkmSmjVrpt9//12HDh1KP2fBggUKDQ1VnTp1ruCbmDsYmQEAAAAcqFevXpo2bZq++uorFS1aNP0al7CwMBUqVEhhYWF67LHH1K9fP5UoUUKhoaF66qmn1KxZM11//fWSpNatW6tOnTp65JFHNGLECB04cEAvvPCCevXq5fMIkQkUMwAAAEAGlp9lugtZMm7cOEnSzTff7LV/0qRJ6tq1qyTpzTfflJ+fn+6991653W5FR0fr3XffTT/X399fc+fOVc+ePdWsWTOFhISoS5cueumll3LraVwRihkAAADAgWzbvuQ5wcHBeuedd/TOO+9c9JyIiAh988032dm1XEMxAwAAAGRgWVxW7hT5sphpHLrRWLZdsoyxbI9/oLHsgFPHjGWnFilhLNsvLcVYdmpoSWPZ/kHBxrJNsrLwCVhO8ZSuYCzbMvhz7gkOMZZt+/kbyy5y+qCxbJP8U84ayy5xdq+x7CKtWxjLVlhxc9lANqDsBAAAAOBI+XJkBgAAALhsDlkAAIzMAAAAAHAoRmYAAACADCwfb14Jc3ilAAAAADgSxQwAAAAAR2KaGQAAAJCBxQIAjsHIDAAAAABHYmQGAAAAyMji836n4JUCAAAA4EiMzAAAAAAZcM2Mc+S5kRnbtk13AQAAAIAD5LlixuVyaePGjaa7AQAAACCPMzbNrF+/fhfcn5aWptdee00lS5aUJI0aNSrTdtxut9xut/e+5GS5goKyp6MAAAAoWPzy3Of9uAhjxczo0aPVoEEDFStWzGu/bdvauHGjQkJCZFmXnq8YFxenYcOGee3r+2RP9XvqiezsLgAAAIA8xlgx8+qrr2rChAkaOXKkbr311vT9gYGBmjx5surUqZOldmJjY88b5Tm8a2u29hUAAAAFR1Y+UEfeYGwM7bnnntNnn32mnj17asCAAUpJSbmsdlwul0JDQ702ppgBAAAA+Z/RCYHXXHON4uPjdfjwYTVt2lTr16+nEgYAAACQJcbvM1OkSBFNmTJF06dPV6tWrZSWlma6SwAAACjIWADAMYwXM+c88MADuvHGGxUfH6+IiAjT3QEAAACQx+WZYkaSKlasqIoVK5ruBgAAAAowy4/LHpyCMTQAAAAAjkQxAwAAAMCR8tQ0MwAAAMA4i8/7nYJXCgAAAIAjMTIDAAAAZMQCAI7ByAwAAAAAR2JkBgAAAMjA4poZx+CVAgAAAOBIFDMAAAAAHClfTjNzpZ4xlm27ChvLNskTHGIs2y8txVi2xz/QWLZ/SpKxbE+Qy1i25Ukzlm2b/PjHMngxqm0bi7YM/n7bBXWaicHnbXnMvd6ByaeNZatYSWPRdkC+/FPwyrEAgGMU0HdqAAAAAE5HOQ4AAABkYPnxeb9T8EoBAAAAcCSKGQAAAACOxDQzAAAAICOTi67AJ4zMAAAAAHAkRmYAAACAjFgAwDF4pQAAAAA4EiMzAAAAQEZcM+MYjMwAAAAAcCSKGQAAAACOxDQzAAAAIAOLBQAcg1cKAAAAgCPlqZGZ06dPa8aMGdq6davCw8P14IMPqmTJkpk+xu12y+12e+9LTpYrKCgnuwoAAID8yuLzfqcw+krVqVNHx44dkyTt3r1bdevWVd++fbVgwQINHTpUderU0Y4dOzJtIy4uTmFhYV7bmPen5Eb3AQAAABhktJjZtGmTUlNTJUmxsbEqX768du7cqVWrVmnnzp2qX7++Bg0alGkbsbGxSkhI8Np6d++SG90HAAAAYFCemWa2YsUKjR8/XmFhYZKkIkWKaNiwYXrggQcyfZzL5ZLL5fLa52aKGQAAAC6XH/eZcQrjEwKt/92UKCkpSeHh4V7HKlSooMOHD5voFgAAAIA8zvjITMuWLRUQEKDExERt3rxZdevWTT+2c+fOSy4AAAAAAGQniwUAHMNoMTN06FCvr4sUKeL19Zw5c9SiRYvc7BIAAAAAh8hTxcw/vf7667nUEwAAAOB/uGbGMRhDAwAAAOBIFDMAAAAAHMn4AgAAAABAnsICAI7BKwUAAADAkRiZAQAAADKyWADAKRiZAQAAAOBIFDMAAAAAHIlpZgAAAEBGfnze7xS8UgAAAAAcKV+OzJx0lTSWXcS92Vi2n8GL1fxPJRjLTg0193r7pyQZy04LDDaW7X/2pLHsNFeIsWzLk2Ys2y/5rLHstGBz33Pbz+A/UwbfU1MDzP1+B509YSzbExBkLDvZVdRYtuv4YWPZVlgJY9l5GkszOwavFAAAAABHypcjMwAAAMBl82NpZqdgZAYAAACAI1HMAAAAAHAkppkBAAAAGbEAgGPwSgEAAABwJEZmAAAAgIwMLs0O3zAyAwAAAMCRKGYAAAAAOBLTzAAAAICM/Pi83yl4pQAAAAA4EiMzAAAAQEYsAOAYjMwAAAAAcCSjxczq1au1Y8eO9K8/+ugj3XDDDapUqZJuvPFGTZ8+/ZJtuN1uJSYmem1utzsnuw0AAID8zPIzt8EnRr9jjz76qLZt2yZJmjhxov7973+radOmGjRokK655hp1795dH374YaZtxMXFKSwszGsbP358bnQfAAAAgEFGr5nZsmWLatSoIUl69913NWbMGHXv3j39+DXXXKNXXnlF3bp1u2gbsbGx6tevn9e+vXv25EyHAQAAAOQZRouZwoUL68iRI4qIiNDevXt17bXXeh2/7rrrvKahXYjL5ZLL5fLad+QfXwMAAABZxtLMjmH0lWrbtq3GjRsnSYqKitLnn3/udXzGjBmKjIw00TUAAAAAeZzRkZnhw4frhhtuUFRUlJo2baqRI0dqyZIlql27tjZv3qyVK1fqiy++MNlFAAAAFDQszewYRkdmypcvr99++03NmjXTvHnzZNu2Vq1ape+++04VK1bUjz/+qNtvv91kFwEAAADkUcZvmlmsWDG99tpreu2110x3BQAAAICDGC9mAAAAgDyF+704Bq8UAAAAAEdiZAYAAADIiAUAHIORGQAAAACOxMgMAAAAkBE3zXQMXikAAAAAjkQxAwAAAMCRmGYGAAAAZGCzAIBj5Mti5khKKWPZ5c+eNpZtuQoby5bB5+0fFGws2xPkMpbtf/aksey0QkWNZVtpqcayPQHmXm+5k8xlG+SXlmIs22Pwj5k0gz9rtsH7a9iWv7Hs067ixrIDtu4wlh1cNc1YNpAd8mUxAwAAAFw2bprpGLxSAAAAAByJYgYAAACAIzHNDAAAAMiIaWaOwSsFAAAAwJEYmQEAAAAyYGlm52BkBgAAAIAjUcwAAAAAcCSmmQEAAAAZsQCAY/BKAQAAAHAkRmYAAACAjFgAwDEYmQEAAADgSIzMAAAAABn58Xm/U/BKAQAAAHAko8XMU089pWXLll1RG263W4mJiV5bcrI7m3oIAAAAIK8yWsy88847uvnmm3XVVVdp+PDhOnDggM9txMXFKSwszGubMmFUDvQWAAAABYFtWcY2+Mb4NLPvvvtOt99+u9544w1VrlxZ7du319y5c+XxeLL0+NjYWCUkJHhtXf7VL4d7DQAAAMA048VMvXr1NHr0aO3bt08ff/yx3G63YmJiVKlSJQ0aNEhbt27N9PEul0uhoaFeW1CQK5d6DwAAgHzH8jO3wSd55jsWGBio+++/X/PmzdP27dvVvXt3ffLJJ6pZs6bprgEAAADIg/JMMZNR5cqV9eKLL2rHjh2aN2+e6e4AAAAAyIOM3mcmIiJC/v7+Fz1uWZZuu+22XOwRAAAACjqb6V6OYbSY2bFjh8l4AAAAAA5G2QkAAABkZFnmNh+98847qlKlioKDg3Xddddp1apVOfANybsoZgAAAAAH+uyzz9SvXz8NHTpUq1evVoMGDRQdHa1Dhw6Z7lquoZgBAAAAMrAtP2ObL0aNGqXu3bvr0UcfVZ06dTR+/HgVLlxYH374YQ59Z/IeihkAAAAgj3C73UpMTPTa3G73eeclJycrPj5erVq1St/n5+enVq1aacWKFbnZZaMoZgAAAIA8Ii4uTmFhYV5bXFzceecdOXJEaWlpKlu2rNf+smXL6sCBA7nVXeOMrmYGAAAA5DmXcSF+domNjVW/fv289rlcLkO9yfsoZgAAAIA8wuVyZal4KVWqlPz9/XXw4EGv/QcPHlS5cuVyqnt5DtPMAAAAgIwsP3NbFgUFBalJkyZatGhR+j6Px6NFixapWbNmOfFdyZPy5chMrdMrjWWnlqpgLNuktPLVTHfBCMuTZiw7zRViLNtKSzWWbZJ/yllj2aklzH3KZvLn3BMQZCzb9vM3lh18+oixbMv2GMs2+TtW6uhmY9nWrdHGstO4072j9evXT126dFHTpk117bXXavTo0Tp9+rQeffRR013LNfmymAEAAADyu44dO+rw4cMaMmSIDhw4oIYNG2revHnnLQqQn1m2bdumO5HdEuPnG8v2S002lm2SyU8wTSqon1gX1JEZk59Yy+RbtcELYU3ifQ25xeT33Nf7mmSnItffZSz7Ukz+LRnaxNxInRMxtggAAADAkZhmBgAAAGTEtUSOwSsFAAAAwJEYmQEAAAAysFUwrxV0IkZmAAAAADgSxQwAAAAAR2KaGQAAAJCBySWr4RteKQAAAACOxMgMAAAAkBEjM47BKwUAAADAkShmAAAAADgS08wAAACADGyL+8w4hfGRmbfffludO3fW9OnTJUkfffSR6tSpo1q1aun5559Xampqpo93u91KTEz02tzJybnRdQAAAAAG+VzMbNq06aLH5s+f71Nb//nPf/T888/rzJkz6tu3r4YPH66+ffuqU6dO6tKliyZOnKiXX3450zbi4uIUFhbmtY2a9JlP/QAAAADOsS0/Yxt8Y9m2bfvygMKFC+v1119Xr1690ve53W71799fEydOVFJSUpbbioyM1IgRI3TPPfdo7dq1atKkiaZMmaJOnTpJkr744gs988wz2rJly0XbcLvdcrvd3vs2LJUrKMiXp5Vt/FIL5qiQ7edvugtGWJ40Y9meADM/45JkpWU+YppfWbbHXLhvb9XZq4BOt+B9DbnF5Pfc5B/PRa6/y1j2pRxbt8xYdon6LYxlO5HP18xMnjxZPXv21Ndff61JkyZp//79euihh+TxeLRsmW8v/L59+9S0aVNJUoMGDeTn56eGDRumH2/cuLH27duXaRsul0sul8trX6KhQgYAAAD5QAH9EMeJfC7H77//fq1du1YpKSm6+uqr1axZM0VFRWn16tW65pprfGqrXLly+uOPPyRJW7ZsUVpaWvrXkrRhwwaVKVPG1y4CAAAAKAAuezWz5ORkpaWlKS0tTeHh4QoODva5jU6dOqlz585q3769Fi1apGeeeUYDBgzQ0aNHZVmWXnnlFd13332X20UAAAAA+ZjPxcz06dPVs2dPtWjRQn/++afWrFmjRx99VPPnz9dHH32katWqZbmtYcOGqVChQlqxYoW6d++u5557Tg0aNNAzzzyjM2fOqF27dpdcAAAAAADITlyI7xw+LwAQEhKiN954Qz179kzfd/z4cf373//WvHnzlJiYmO2d9FVivG+rqmUnFgAoWFgAoGBhAYCChfc15BYWAMh7jq7/yVh2ybrNjWU7kc8jM6tXr1bNmjW99hUvXlwzZszQRx99lG0dAwAAAEywVTA/xHEin8vxmjVrKjU1VQsXLtR7772nkydPSvp7ZbK777472zsIAAAAABfi88jMzp071aZNG+3atUtut1u33XabihYtquHDh8vtdmv8+PE50U8AAAAA8OLzyEzv3r3VtGlTHT9+XIUKFUrff/fdd2vRokXZ2jkAAAAgt9mWn7ENvvF5ZGbZsmX66aefFPSPG1NWqVJFe/fuzbaOAQAAAMhf0tLSNHnyZC1atEiHDh2Sx+O9uM7333/vU3s+FzMej0dpaeevurFnzx4VLVrU1+YAAACAvKWAruKYG3r37q3JkyfrjjvuUN26dWVd4ffa52KmdevWGj16tCZMmCBJsixLp06d0tChQ3X77bdfUWcAAAAA5F/Tp0/XjBkzsq1u8LmYGTlypKKjo1WnTh0lJSXpoYce0pYtW1SqVCl9+umn2dIpAAAAwBTb98vKkUVBQUGKjIzMtvZ8vmmmJKWmpmr69Olat26dTp06pcaNG6tTp05eCwKYdPZ7c/e7SSlWzli2X/JZY9meQJexbMvkzQRN3kTRII+/uRt2+qeY+zk3emGmwRs4+qUkGcvm5rAFi8kblfqlpRjLdoeUNJYdkGru9zu08W3Gsi/l0B+/GssuU6epsezcMHLkSG3fvl1vv/32FU8xky5jZEaSAgIC9PDDD19xOAAAAICCY/ny5Vq8eLG+/fZbXX311QoMDPQ6PmvWLJ/ay1IxM3v27Cw3eNddd/nUAQAAACAvsVkAIMcUK1ZMd999d7a1l6ViJiYmxutry7L0z9lp54aJLrTSGQAAAABMmjQpW9vL0gRwj8eTvn333Xdq2LChvv32W504cUInTpzQt99+q8aNG2vevHnZ2jkAAAAgt3HTzJx3+PBhLV++XMuXL9fhw4cvux2fr5np06ePxo8frxtvvDF9X3R0tAoXLqx//etf2rhx42V3BgAAAED+dfr0aT311FOaOnVq+g0z/f391blzZ40dO1aFCxf2qT2fy79t27apWLFi5+0PCwvTX3/95WtzAAAAAAqIfv36aenSpZozZ076LK+vvvpKS5cuVf/+/X1uz+di5pprrlG/fv108ODB9H0HDx7UwIEDde211/rcAQAAACAvsWUZ2/K7//73v/rggw/Utm1bhYaGKjQ0VLfffrvef/99ff755z6353Mx8+GHH2r//v2qXLmyIiMjFRkZqcqVK2vv3r364IMPfO4AAAAAgILhzJkzKlu27Hn7y5QpozNnzvjcns/XzERGRmrdunVasGCBNm3aJEmqXbu2WrVqlS03vgEAAABMKkgX4ue2Zs2aaejQoZo6daqCg4MlSWfPntWwYcPUrFkzn9u7rJtmWpal1q1bq3Xr1pfzcAAAAAAF0JgxYxQdHa2KFSuqQYMGkqS1a9cqODhY8+fP97m9yypmFi1apEWLFunQoUPpqxCc8+GHH15OkwAAAECewE0zc07dunW1ZcsWffLJJ+mzvB588EF16tRJhQoV8rk9n4uZYcOG6aWXXlLTpk0VHh7O1DIAAAAAWVa4cGF17949W9ryuZgZP368Jk+erEceeeSKw/fv369x48Zp+fLl2r9/v/z8/FStWjXFxMSoa9eu8vf3v+IMAAAAAObMnj1bbdu2VWBgoGbPnp3puXfddZdPbftczCQnJ6t58+a+Puw8v/76q1q1aqXIyEgVKlRIW7Zs0UMPPaTk5GQNGDBAH374oebNm6eiRYtm2o7b7Zbb7fba50lOkSso8Ir7CAAAgIKnICyRnJtiYmJ04MABlSlTRjExMRc9z7IspaWl+dS2z0s1PP7445o2bZqvDztPnz591LdvX/36669atmyZJk+erD///FPTp0/X9u3bdebMGb3wwguXbCcuLk5hYWFe2+ufzrni/gEAAAC4ch6PR2XKlEn//4ttvhYykmTZtm378oDevXtr6tSpql+/vurXr6/AQO8RkFGjRmWpncKFC2v9+vWqVq2apL+fWHBwsHbv3q2yZctqwYIF6tq1q/bu3ZtpOxccmfnpc2MjMynFyhnJlSS/5LPGsj2BLmPZlm8/wtnL9lz6nHzI4x9kLNs/xdzPudGlOv3MTbv1S0kylu0JMPezZqWlGssuqGyTP+dpKcay3SEljWUHpJr7/Q5tfJux7EvZveUPY9mVatQxlp0bpk6dqo4dO8rl8v7bMTk5WdOnT1fnzp19as/naWbr1q1Tw4YNJUnr16/3OubLYgBlypTR/v3704uZgwcPKjU1VaGhoZKkGjVq6NixY5dsx+VynffNOMsUMwAAACDPefTRR9WmTZv0kZpzTp48qUcffTTni5nFixf7+pALiomJUY8ePfT666/L5XLp5ZdfVlRUVPqSbJs3b1aFChWyJQsAAACAebZtX3AAZM+ePQoLC/O5vcu6z0x2+M9//qP9+/erXbt2SktLU7NmzfTxxx+nH7csS3Fxcaa6BwAAgAKKBQCyX6NGjWRZlizLUsuWLRUQ8P9lSFpamnbs2KE2bdr43G6Wi5l77rknS+fNmjUrS+cVKVJEn332mZKSkpSamqoiRYp4HW/dunVWuwYAAAAgDzu3itmaNWsUHR3t9bd/UFCQqlSponvvvdfndrNczFzOsE9WBAcH50i7AAAAwOUwuuBLPjV06FBJUpUqVdSxY8dsqwGyXMxMmjQpWwIBAAAAFExdunTJ1vaMXTMDAAAAoGBJS0vTm2++qRkzZmjXrl1KTk72Op6V1YwzYgwNAAAAyMCWZWzL74YNG6ZRo0apY8eOSkhIUL9+/XTPPffIz89PL774os/tUcwAAAAAyBWffPKJ3n//ffXv318BAQF68MEHNXHiRA0ZMkQrV670uT2KGQAAACAD2/IztuV3Bw4cUL169ST9vbpxQkKCJOnOO+/U119/7XN7+f87BgAAACBPqFixovbv3y9Jql69ur777jtJ0i+//CKXy+Vze5dVzHz00Ue64YYbVL58ee3cuVOSNHr0aH311VeX0xwAAACQZ3DNTM65++67tWjRIknSU089pcGDB6tGjRrq3LmzunXr5nN7Pq9mNm7cOA0ZMkR9+vTRK6+8orS0NElSsWLFNHr0aLVv397nTgAAAADI/1577bX0/+/YsaMqV66sFStWqEaNGmrXrp3P7fk8MjN27Fi9//77GjRokPz9/dP3N23aVL///rvPHQAAAABQMDVr1kz9+vW7rEJGuoyRmR07dqhRo0bn7Xe5XDp9+vRldSK7rSx9n7HsZrumGMtOK17GWHbQX9uMZXtKVzCWLcvccLBf8llj2XInGYtOLVHOWLYyfICT6zxp5qIDgoxl+yeZ+3clzVXYWPapUHPva0VOHTCWbRt8Tz1Yspax7CIThxnLDo6sYixbjW8zl30JJn8W86PZs2dn+dy77rrLp7Z9LmaqVq2qNWvWKCIiwmv/vHnzVLt2bV+bAwAAAJCPxcTEZOk8y7LSL2HJKp+LmX79+qlXr15KSkqSbdtatWqVPv30U8XFxWnixIm+NgcAAADkKbbNyEx28ng8Oda2z8XM448/rkKFCumFF17QmTNn9NBDD6l8+fIaM2aMHnjggZzoIwAAAIB8JikpScHBwVfUhk8LAKSmpmrq1Klq1aqVtmzZolOnTunAgQPas2ePHnvssSvqCAAAAID8LS0tTS+//LIqVKigIkWKaPv27ZKkwYMH64MPPvC5PZ+KmYCAAPXo0UNJSX9f/Fu4cGGVKWPuonMAAAAgu9nyM7bld6+88oomT56sESNGKCjo/xeYqVu37mVdsuLzd+zaa6/Vb7/95nMQAAAAgIJt6tSpmjBhgjp16uR1m5cGDRpo06ZNPrfn8zUzTzzxhPr37689e/aoSZMmCgkJ8Tpev359nzsBAAAA5BW2WAAgp+zdu1eRkZHn7fd4PEpJSfG5PZ+LmXMX+T/99NPp+yzLkm3bl7WcGgAAAICCoU6dOlq2bNl5t3n5/PPPL3gvy0u5rJtmAgAAAPkVIzM5Z8iQIerSpYv27t0rj8ejWbNmafPmzZo6darmzp3rc3s+FzP/rKIAAAAAICvat2+vOXPm6KWXXlJISIiGDBmixo0ba86cObrtttt8bs/nYmbq1KmZHu/cubNP7SUnJ+vLL7/UihUrdODAAUlSuXLl1Lx5c7Vv395rlQMAAAAAzpSamqpXX31V3bp104IFC7KlTZ+Lmd69e3t9nZKSojNnzigoKEiFCxf2qZjZunWroqOjtW/fPl133XUqW7asJOm3337T+PHjVbFiRX377bcXvEgIAAAAyAlMM8sZAQEBGjFihM+DH5m26esDjh8/ft6+LVu2qGfPnho4cKBPbfXs2VP16tXTb7/9ptDQUK9jiYmJ6ty5s3r16qX58+f72k0AAAAAeUzLli21dOlSValSJVva87mYuZAaNWrotdde08MPP+zT+tA//vijVq1adV4hI0mhoaF6+eWXdd1112Xahtvtltvt9tqXnOxRUJAry/0AAAAAzmFkJue0bdtWzz33nH7//fcL3ublrrvu8qm9bClmpL+Hjfbt2+fTY4oVK6a//vpLdevWveDxv/76S8WKFcu0jbi4OA0bNsxrX+cez6vrEy/41BcAAAAAOeuJJ56QJI0aNeq8Y5dzmxefi5nZs2d7fW3btvbv36+3335bN9xwg09tPf744+rcubMGDx6sli1bpl8zc/DgQS1atEj/+c9/9NRTT2XaRmxsrPr16+e1b8UWj0/9AAAAAJDzPJ7s/Tvd52ImJibG62vLslS6dGndeuutGjlypE9tnVuS7fXXX1f//v1lWX8P6dm2rXLlyunZZ5/VM888k2kbLpdLLpf3lLKgoLM+9QMAAAA4x7aZZpYTUlJSVKhQIa1Zs+aiM7N85XMxk93V1LPPPqtnn31WO3bs8FqauWrVqtmaAwAAAMCcwMBAVa5c2eepZJnx8/UBL730ks6cOXPe/rNnz+qll1667I5UrVpVzZo1U7NmzdILmd27d6tbt26X3SYAAADgK1uWsS2/GzRokJ5//nkdO3YsW9rzuZgZNmyYTp06dd7+M2fOnHch/pU6duyYpkyZkq1tAgAAADDj7bff1g8//KDy5curZs2aaty4sdfmK5+nmdm2nX5tS0Zr165ViRIlfGrrn4sJ/NP27dt9ag8AAAC4UgVhhMSUf15/f6WyXMwUL15clmXJsixdddVVXgVNWlqaTp06pR49evgUHhMTI8uyZNv2Rc+5UOEEAAAAwHmGDh2are1luZgZPXq0bNtWt27dNGzYMIWFhaUfCwoKUpUqVdSsWTOfwsPDw/Xuu++qffv2Fzy+Zs0aNWnSxKc2AQAAAORt8fHx2rhxoyTp6quvVqNGjS6rnSwXM126dJH094X6zZs3V2Bg4GUFZtSkSRPFx8dftJi51KgNAAAAkN2YZpZzDh06pAceeEBLlixRsWLFJEknTpzQLbfcounTp6t06dI+tefzAgBRUVHphUxSUpISExO9Nl8MHDhQzZs3v+jxyMhILV682NcuAgAAAMiDnnrqKZ08eVIbNmzQsWPHdOzYMa1fv16JiYl6+umnfW7P5wUAzpw5o2eeeUYzZszQ0aNHzzvuy7rRLVq0yPR4SEiIoqKifO0iAAAAcNm4aWbOmTdvnhYuXKjatWun76tTp47eeecdtW7d2uf2fB6ZGThwoL7//nuNGzdOLpdLEydO1LBhw1S+fHlNnTrV5w4AAAAAKBg8Hs8FL1cJDAyUx+PxuT2fi5k5c+bo3Xff1b333quAgAC1aNFCL7zwgl599VV98sknPncAAAAAQMFw6623qnfv3tq3b1/6vr1796pv375q2bKlz+35XMwcO3ZM1apVkySFhoam373zxhtv1A8//OBzBwAAAIC8xCPL2Jbfvf3220pMTFSVKlVUvXp1Va9eXVWrVlViYqLGjh3rc3s+XzNTrVo17dixQ5UrV1atWrU0Y8YMXXvttZozZ076igQAAAAA8E+VKlXS6tWrtXDhQm3atEmSVLt2bbVq1eqy2rNsH9c+fvPNN+Xv76+nn35aCxcuVLt27WTbtlJSUjRq1Cj17t37sjqSnc4u+dRYdkpoKWPZVlqKsWzb/8qX6r5cJp+3TC4dXkBvKGt5sr7ISH7K9gQEGcsusD/nBp+3n8H3NZM/ayZ/x/ySk4xlp4QUN5ZtWz5P0sk2xRrdaiz7Un7bcsRYdqMa5v6WzEnff/+9nnzySa1cuVKhoaFexxISEtS8eXONHz/+kguE/ZPPIzN9+/ZN//9WrVpp06ZNio+PV2RkpOrXr+9rcwAAAADyudGjR6t79+7nFTKSFBYWpn//+98aNWpUzhczGSUlJSkiIkIRERFX0gwAAACQZ7A0c/Zbu3athg8fftHjrVu31htvvOFzuz6PLaalpenll19WhQoVVKRIEW3fvl2SNHjwYH3wwQc+dwAAAABA/nbw4MELLsl8TkBAgA4fPuxzuz4XM6+88oomT56sESNGKCjo/+fV1q1bVxMnTvS5AwAAAADytwoVKmj9+vUXPb5u3TqFh4f73K7PxczUqVM1YcIEderUSf7+/un7GzRokL4iAQAAAOBUtixjW351++23a/DgwUpKOn+xjbNnz2ro0KG68847fW7X52tm9u7dq8jIyPP2ezwepaQYXFUKAAAAQJ70wgsvaNasWbrqqqv05JNPqmbNmpKkTZs26Z133lFaWpoGDRrkc7s+FzN16tTRsmXLzrvo//PPP1ejRo187gAAAACQl7AAQPYrW7asfvrpJ/Xs2VOxsbE6d3cYy7IUHR2td955R2XLlvW5XZ+LmSFDhqhLly7au3evPB6PZs2apc2bN2vq1KmaO3euzx0AAAAAkP9FRETom2++0fHjx7V161bZtq0aNWqoePHLv9eSz8VM+/btNWfOHL300ksKCQnRkCFD1LhxY82ZM0e33XbbZXcEAAAAQP5XvHhxXXPNNdnSVpaLme3bt6tq1aqyLEstWrTQggULsqUDmTl48KDee+89DRkyJMezAAAAAEn5+kL8/CbLq5nVqFHDa+3njh076uDBgznSqXMOHDigYcOG5WgGAAAAAGfK8sjMuYt0zvnmm28UFxd3ReHr1q3L9PjmzZuvqH0AAADAVywA4Bw+XzOTnRo2bCjLss4rlCSl77csfpgAAAAAnC/LxYxlWecVFldaaJQoUUIjRoxQy5YtL3h8w4YNateuXaZtuN1uud1ur32e5BS5ggKvqG8AAAAomDymO4As82maWdeuXeVyuSRJSUlJ6tGjh0JCQrzOmzVrVpbDmzRpon379p13z5pzTpw4ccFRm4zi4uLOu67m+S736IWu92W5HwAAAACcJ8vFTJcuXby+fvjhh684vEePHjp9+vRFj1euXFmTJk3KtI3Y2Fj169fPa59n5ZdX3DcAAAAAeZtlX2row4HOLvnUWHZKaClj2VZairFs29/ctD6Tz1smf30K6PVklietQGZ7AoKMZRfYn3ODz9vP4PuayZ81k79jfslJxrJTQi7/hoFXyrayvLBttivW6FZj2ZeyYmOisexmtUONZTuRuZ/gLNi9e7e6detmuhsAAAAA8qA8XcwcO3ZMU6ZMMd0NAAAAFCC2LGMbfGN0aebZs2dnenz79u251BMAAAAATmO0mImJibnofWbO4T4zAAAAAC7E6DSz8PBwzZo1Sx6P54Lb6tWrTXYPAAAABZBtW8Y2+MZoMdOkSRPFx8df9PilRm0AAAAAFFxGp5kNHDgw0/vMREZGavHixbnYIwAAABR0XIjvHEaLmRYtWmR6PCQkRFFRUbnUGwAAAABOkqeXZgYAAABym8c2t+WEv/76S4899piqVq2qQoUKqXr16ho6dKiSk5O9zlu3bp1atGih4OBgVapUSSNGjDivrZkzZ6pWrVoKDg5WvXr19M033+RMp7OIYgYAAADIxzZt2iSPx6P33ntPGzZs0Jtvvqnx48fr+eefTz8nMTFRrVu3VkREhOLj4/X666/rxRdf1IQJE9LP+emnn/Tggw/qscce02+//aaYmBjFxMRo/fr1Jp6WJMmy8+EV9meXfGosOyW0lLFsKy3FWLbtH2gs2+TzlslfnwK6bLnlSSuQ2Z6AIGPZBfbn3ODz9jP4vmbyZ83k75hfcpKx7JSQ4saybcvc59rFGt1qLPtSfthw8Wu6c9pNV4fkSs7rr7+ucePGpd/Xcdy4cRo0aJAOHDigoKC/3weee+45ffnll9q0aZMkqWPHjjp9+rTmzp2b3s7111+vhg0bavz48bnS739iZAYAAADIwJZlbHO73UpMTPTa3G53tj/HhIQElShRIv3rFStW6KabbkovZCQpOjpamzdv1vHjx9PPadWqlVc70dHRWrFiRbb3L6uMLgCQUxLLXmUsu+jhbcayPUHBxrL9zyQay/YE584nGBdidDTMz9yvb0H91Nhktn+SuU8J0wz+jhkdFTIo1VXEWLZfavb/0ZRVJn/H0gILGct2Hd5pLNsONDjqiwuKi4vTsGHDvPYNHTpUL774YrZlbN26VWPHjtUbb7yRvu/AgQOqWrWq13lly5ZNP1a8eHEdOHAgfV/Gcw4cOJBtffMVIzMAAABABiZvmhkbG6uEhASvLTY29oL9fO6552RZVqbbuSli5+zdu1dt2rRRhw4d1L1799z4duaofDkyAwAAADiRy+WSy+XK0rn9+/dX165dMz2nWrVq6f+/b98+3XLLLWrevLnXhf2SVK5cOR08eNBr37mvy5Url+k5546bQDEDAAAAOFDp0qVVunTpLJ27d+9e3XLLLWrSpIkmTZokPz/vCVrNmjXToEGDlJKSosDAvxd2WrBggWrWrKnixYunn7No0SL16dMn/XELFixQs2bNsucJXQammQEAAAAZ2La5LSfs3btXN998sypXrqw33nhDhw8f1oEDB7yudXnooYcUFBSkxx57TBs2bNBnn32mMWPGqF+/funn9O7dW/PmzdPIkSO1adMmvfjii/r111/15JNP5kzHs4CRGQAAACAfW7BggbZu3aqtW7eqYsWKXsfO3aUlLCxM3333nXr16qUmTZqoVKlSGjJkiP71r3+ln9u8eXNNmzZNL7zwgp5//nnVqFFDX375perWrZurzyejfHmfmYMb441lF9TVzPzcZ4xls5pZ7iuoq5nZfv7GslnNrGDxBGRtvnxOMLmame1v7n3Ntsz9fhfU1cwK3fqIsexLWfS7ufsOtaxn7u85J2KaGQAAAABHopgBAAAA4EhcMwMAAABkYNuW6S4gixiZAQAAAOBIeaKY2bNnj06dOnXe/pSUFP3www8GegQAAICCKr8tzZyfGS1m9u/fr2uvvVYREREqVqyYOnfu7FXUHDt2TLfccovBHgIAAADIq4wWM88995z8/Pz0888/a968efrjjz90yy236Pjx4+nn5MOVowEAAJCH2bKMbfCN0WJm4cKFeuutt9S0aVO1atVKP/74o8LDw3Xrrbfq2LFjkiTL4kUFAAAAcD6jxUxCQoKKFy+e/rXL5dKsWbNUpUoV3XLLLTp06NAl23C73UpMTPTa3MnJOdltAAAAAHmA0WKmWrVqWrdunde+gIAAzZw5U9WqVdOdd955yTbi4uIUFhbmtb01YVJOdRkAAAD5nMc2t8E3RouZtm3basKECeftP1fQNGzY8JLXzMTGxiohIcFre/pfj+ZUlwEAAADkEUZvmvnKK6/ozJkzFzwWEBCg//73v9q7d2+mbbhcLrlcLq99Z4OCsq2PAAAAKFi4aaZzGB2ZCQgIUGho6EWP79+/X8OGDcvFHgEAAABwijxx08yLOXbsmKZMmWK6GwAAAADyIKPTzGbPnp3p8e3bt+dSTwAAAIC/cZtD5zBazMTExMiyrEwv8uc+MwAAAAAuxOg0s/DwcM2aNUsej+eC2+rVq012DwAAAAWQR5axDb4xWsw0adJE8fHxFz1+qVEbAAAAAAWX0WlmAwcO1OnTpy96PDIyUosXL87FHgEAAKCg47N05zBazLRo0SLT4yEhIYqKisql3gAAAABwkjy9NDMAAAAAXIzRkRkAAAAgr7FtLsR3CkZmAAAAADgSIzMAAABABh4WAHCMfFnMhCTuM5adFhxiLNuyPcayPa7CxrJtP39z2ZbBwU2DN5T1GMw2+XpbaanGstMM/o6xrE/u808+Yyzb6L8lBt9b/DzJxrJTipczlg04HdPMAAAAADhSvhyZAQAAAC4XA9LOwcgMAAAAAEdiZAYAAADIwBZLMzsFIzMAAAAAHImRGQAAACADlmZ2DkZmAAAAADgSxQwAAAAAR2KaGQAAAJABSzM7ByMzAAAAABzJ+MjM0aNHtW7dOjVo0EAlSpTQkSNH9MEHH8jtdqtDhw6qXbu26S4CAACgAGFkxjmMFjOrVq1S69atlZiYqGLFimnBggXq0KGDAgIC5PF49Nprr2n58uVq3LixyW4CAAAAyIOMTjMbNGiQOnTooISEBD3//POKiYlRy5Yt9eeff2rr1q164IEH9PLLL5vsIgAAAIA8ymgxEx8fr379+qlo0aLq3bu39u3bp+7du6cff/LJJ/XLL78Y7CEAAAAKGo9tGdvgG6PTzJKTk1WoUCFJUmBgoAoXLqxSpUqlHy9VqpSOHj2aaRtut1tut9trX0pyilxBgdnfYQAAAAB5htGRmUqVKmn79u3pX0+fPl3h4eHpX+/fv9+ruLmQuLg4hYWFeW0jp8zMsT4DAAAgf7Ntcxt8Y3Rk5oEHHtChQ4fSv77jjju8js+ePVvXXnttpm3ExsaqX79+XvtS1i7Mvk4CAAAAyJMs2867NeCZM2fk7+8vl8vl0+NO/Twnh3p0abZlbrDLsj3Gsk1+lGD7G6zJTf76WObm1VqeNGPZtp+/sWwrLdVYtsnnjdxn8nfM5L8lnoAgY9lG39cM/u1gUtFr77j0SYZMW27u3/eHbuS6GV/k6d+eo0ePqmfPnqa7AQAAACAPytPFzLFjxzRlyhTT3QAAAACQBxm9Zmb27NmZHs+4OAAAAACQGzx59iIM/JPRYiYmJkaWZSmzy3Ysg9cFAAAAAMi7jE4zCw8P16xZs+TxeC64rV692mT3AAAAUADZtmVsg2+MFjNNmjRRfHz8RY9fatQGAAAAQMFldJrZwIEDdfr06Ysej4yM1OLFi3OxRwAAAACcwmgx06JFi0yPh4SEKCoqKpd6AwAAAJi9jRx8k6eXZgYAAACAizE6MgMAAADkNSzN7ByMzAAAAABwJEZmAAAAgAy4ZsY58mUxk7xgrrFs9309jWWHnDlsLPtUSFlj2UVOHzSWbVJqQLCx7LQAl7Hs4NNHjGWbdCq0grHs0ON/GctOdRUxlu2ffMZYtu3nbyxbHnPRfilJxrL9d24xlv1rgyeNZZcISjCWXcdYMvITppkBAAAAcKR8OTIDAAAAXC6mmTkHIzMAAAAAHImRGQAAACADlmZ2DkZmAAAAADgSxQwAAAAAR2KaGQAAAJABCwA4ByMzAAAAAByJkRkAAAAgA4/Bm8fCN4zMAAAAAHCkPFnMVKtWTVu2bDHdDQAAABRAtm1ug2+MTjN76623Lrh/165dmjRpksqVKydJevrpp3OzWwAAAAAcwGgx06dPH1WoUEEBAd7d8Hg8mjp1qgIDA2VZFsUMAAAAgPMYLWb+9a9/6eeff9a0adNUu3bt9P2BgYH67rvvVKdOnUu24Xa75Xa7vfelpskV4J/t/QUAAED+x3Qv5zB6zcz48eM1ZMgQRUdH6+23376sNuLi4hQWFua1jf7ht2zuKQAAAIC8xvgCAHfffbdWrFihL774Qm3bttWBAwd8enxsbKwSEhK8tj43Ncqh3gIAACC/89jmNvjGeDEjSRUqVNDChQt10003qVGjRrJ9GNtzuVwKDQ312phiBgAAAOR/eeammZZlKTY2Vq1bt9by5csVHh5uuksAAAAA8rA8MTKTUZMmTdS7d28VL15cu3fvVrdu3Ux3CQAAAAWIbdvGNvgmzxUzGR07dkxTpkwx3Q0AAAAAeZDRaWazZ8/O9Pj27dtzqScAAADA3xggcQ6jxUxMTIwsy8p0SM2yrFzsEQAAAACnMDrNLDw8XLNmzZLH47ngtnr1apPdAwAAAJCHGS1mmjRpovj4+Isev9SoDQAAAJDdPB5zG3xjdJrZwIEDdfr06Ysej4yM1OLFi3OxRwAAAACcwmgx06JFi0yPh4SEKCoqKpd6AwAAALAAgJPk6aWZAQAAAOBijI7MAAAAAHmNh5EZx2BkBgAAAIAjUcwAAAAAcKR8Oc0sqNUdxrIDE/cay7Zsc+v5FU3cYyzbKMvc5wFBZ08Yy7YNPm+TP+cmFTl1wFi2JyDIWLZfqttYttGfNYPRtp+/uXCDN8pOrVbHWHaDU0uMZcvk663yBrMzxwIAzsHIDAAAAABHypcjMwAAAMDlso2uAGBuhNKJGJkBAAAA4EgUMwAAAAAciWlmAAAAQAbcZ8Y5GJkBAAAA4EiMzAAAAAAZsDSzczAyAwAAAMCRGJkBAAAAMvBw0YxjMDIDAAAAwJHy1MiMbdtasmSJtm7dqvDwcEVHRyswMNB0twAAAADkQUaLmdtvv12ffvqpwsLCdOzYMd1+++1atWqVSpUqpaNHj+qqq67SDz/8oNKlS5vsJgAAAAoQFgBwDqPTzObNmye32y1JeuGFF3Ty5Elt27ZNhw4d0s6dOxUSEqIhQ4aY7CIAAACAPCrPTDP7/vvvNWLECFWtWlWSVLFiRQ0fPlzdu3fP9HFutzu9IDonJTlFriCmpwEAAMB3jMw4h/EFACzLkiQdP35c1atX9zoWGRmpffv2Zfr4uLg4hYWFeW0jp36eY/0FAAAAkDcYH5np2rWrXC6XUlJStGPHDl199dXpxw4cOKBixYpl+vjY2Fj169fPa1/KmgU50VUAAAAAeYjRYqZLly7p/9++fXudOXPG6/h///tfNWzYMNM2XC6XXC6X175TTDEDAADAZfIwz8wxjBYzkyZNyvT40KFD5e/vn0u9AQAAAOAkxq+ZycyxY8f0xBNPmO4GAAAAChDbY26Db/J8MTNlyhTT3QAAAACQBxmdZjZ79uxMj2/fvj2XegIAAAD8zeaaGccwWszExMTIsqxMf2DOLd0MAAAAABkZnWYWHh6uWbNmyePxXHBbvXq1ye4BAAAAyMOMFjNNmjRRfHz8RY9fatQGAAAAyG4ej7kNvjE6zWzgwIE6ffr0RY9HRkZq8eLFudgjAAAAAE5htJhp0aJFpsdDQkIUFRWVS70BAAAAWADASfL00swAAAAAcDEUMwAAAAAcyeg0MwAAACCv8TDLzDEYmQEAAAAKCLfbrYYNG8qyLK1Zs8br2Lp169SiRQsFBwerUqVKGjFixHmPnzlzpmrVqqXg4GDVq1dP33zzTS71/MLy5chMwNF9xrKTKtQylu2ffPGV4XJaWkCwsWz/lLPGsi1PirFsT0CQsWzb8jeWbfL1tv3MPW/b4A2E/dLSjGWb/Dn3mPyepyQZy5bB521b5j5j9U8y929oYpkaxrKDkxKMZedldj4emnnmmWdUvnx5rV271mt/YmKiWrdurVatWmn8+PH6/fff1a1bNxUrVkz/+te/JEk//fSTHnzwQcXFxenOO+/UtGnTFBMTo9WrV6tu3bomng4jMwAAAEBB8O233+q7777TG2+8cd6xTz75RMnJyfrwww919dVX64EHHtDTTz+tUaNGpZ8zZswYtWnTRgMHDlTt2rX18ssvq3Hjxnr77bdz82l4oZgBAAAAMrBtc5vb7VZiYqLX5na7r/g5HTx4UN27d9dHH32kwoULn3d8xYoVuummmxQU9P8j4tHR0dq8ebOOHz+efk6rVq28HhcdHa0VK1Zccf8uF8UMAAAAkEfExcUpLCzMa4uLi7uiNm3bVteuXdWjRw81bdr0guccOHBAZcuW9dp37usDBw5kes654yZQzAAAAAB5RGxsrBISEry22NjYC5773HPPybKsTLdNmzZp7NixOnny5EXbcbJ8uQAAAAAAcLk8BhcAcLmC5XK5snRu//791bVr10zPqVatmr7//nutWLHivHabNm2qTp06acqUKSpXrpwOHjzodfzc1+XKlUv/74XOOXfcBIoZAAAAwIFKly6t0qVLX/K8t956S//5z3/Sv963b5+io6P12Wef6brrrpMkNWvWTIMGDVJKSooCAwMlSQsWLFDNmjVVvHjx9HMWLVqkPn36pLe1YMECNWvWLBuflW8oZgAAAIAMbDt/Lc1cuXJlr6+LFCkiSapevboqVqwoSXrooYc0bNgwPfbYY3r22We1fv16jRkzRm+++Wb643r37q2oqCiNHDlSd9xxh6ZPn65ff/1VEyZMyL0n8w9cMwMAAAAUcGFhYfruu++0Y8cONWnSRP3799eQIUPS7zEjSc2bN9e0adM0YcIENWjQQJ9//rm+/PJLY/eYkRiZAQAAAAqUKlWqXHD0qX79+lq2bFmmj+3QoYM6dOiQU13zmdGRmT179ujIkSPpXy9btkydOnVSixYt9PDDDxtdsxoAAAAFk+0xt8E3RouZe++9VytXrpQkffXVV7r55pt16tQp3XDDDTpz5oyioqI0d+5ck10EAAAAkEcZnWa2YcMGXX311ZL+vkHQq6++qmeffTb9+Ntvv60hQ4bozjvvNNVFAAAAFDCefLYAQH5mdGQmICBAJ0+elCTt2LFDbdu29Tretm1bbd682UTXAAAAAORxRouZqKgoffrpp5KkRo0aacmSJV7HFy9erAoVKmTahtvtVmJiotfmTknJqS4DAAAgn7Nt29gG3xidZvbaa6+pRYsW2rdvn2688UYNGjRIv/zyi2rXrq3Nmzfrs88+0/jx4zNtIy4uTsOGDfPaN+jBO/RCJ6amAQAAAPmZZRsuAbdt26YXXnhBX3/9tU6dOiXp7+ln11xzjQYOHKiYmJhMH+92u+V2u7322d9Pket/dy7NbUkVahnJlST/5NPGstMCCxnL9k85ayzbMrjsiCcgyFi2bfkbyzb5etuWucFs29/cZ09+qcnGsk3+nFueNGPZfilJxrJtfzP/fkpmf8f8k8z9G5pYpoax7OCkBGPZxRrebCz7Uvq/a+7nYeQTIcayncj4fWaqV6+uTz/9VLZt69ChQ/J4PCpVqpQCs1iMuFwuuVwur31JhgoZAAAAOJ/Hw3QvpzB6zUxGlmWpbNmyCg8PTy9kdu/erW7duhnuGQAAAIC8KM8UMxdy7NgxTZkyxXQ3AAAAUIDYtrkNvjE6zWz27NmZHt++fXsu9QQAAACA0xgtZmJiYmRZVqbL0FmWlYs9AgAAAOAURqeZhYeHa9asWfJ4PBfcVq9ebbJ7AAAAKIBsj21sg2+MFjNNmjRRfHz8RY9fatQGAAAAQMFldJrZwIEDdfr0xdfxjoyM1OLFi3OxRwAAACjoPHyY7hhGi5kWLVpkejwkJERRUVG51BsAAAAATmL8ppkAAABAXsK1K86Rp+8zAwAAAAAXQzEDAAAAwJGYZgYAAABkwDQz58iXxcyWcZ8Zy3aNmWIsu2zCn8ayD4VGGssucXavsezA5IuvxpfTkl1FjWWfdhU3ll3q6GZj2X5pKcayD5asZSy73D5z9/xKCyxkLNvPk2ws23/nFmPZqdXqGMv2TzL3npoWHGIse86eRsaya4ab+57fZCwZ+Um+LGYAAACAy8XAjHNwzQwAAAAAR6KYAQAAAOBITDMDAAAAMmABAOdgZAYAAACAIzEyAwAAAGRg24zMOAUjMwAAAAAciWIGAAAAgCMxzQwAAADIwMMCAI7ByAwAAAAARzJazIwcOVI7d+402QUAAADAi23bxjb4xmgxM3DgQFWvXl233XabPvvsMyUnJ5vsDgAAAAAHMT7NbOLEiQoJCdEjjzyi8uXLq0+fPlq/fr3pbgEAAKCAsj22sQ2+MV7M3H777fryyy+1Z88ePfPMM5o/f74aNGiga6+9Vu+//75OnjxpuosAAAAA8iDjxcw5ZcqU0TPPPKONGzdqyZIlqlOnjvr27avw8PBMH+d2u5WYmOi1Jad5cqnXAAAAAEwxWsxYlnXB/S1atNDkyZO1b98+vfnmm5m2ERcXp7CwMK9t4jYWFQAAAMDlYZqZcxgtZi61YkNoaKi6d++e6TmxsbFKSEjw2h6vHpGd3QQAAACQBxm9aabHc+XTwVwul1wul9e+IP88M3sOAAAADuNhiWTHyNN/9e/evVvdunUz3Q0AAAAAeVCeLmaOHTumKVOmmO4GAAAAgDzI6DSz2bNnZ3p8+/btudQTAAAA4G9ciO8cRouZmJgYWZaV6UIAF1vxDAAAAEDBZnSaWXh4uGbNmiWPx3PBbfXq1Sa7BwAAgALItm1jG3xjtJhp0qSJ4uPjL3r8UqM2AAAAAAouo9PMBg4cqNOnT1/0eGRkpBYvXpyLPQIAAEBB5+GaGccwWsy0aNEi0+MhISGKiorKpd4AAAAAcJI8vTQzAAAAAFyM0ZEZAAAAIK9haWbnYGQGAAAAgCMxMgMAAABkwGq6zpEvi5me1ovGsme++5yx7KB27cxlvzvWWHaR1pkvJJGjipU0Fu06fthYdsDWHcayrVujjWW7Q8y93kUmDjOWndKhs7Fs1+GdxrJTipczlv1rgyeNZTc4tcRYdmKZGsay5+xpZCz7gdCvjWX/0LCfsWylbDaXjXyDaWYAAAAAHClfjswAAAAAl8v2eEx3AVnEyAwAAAAAR2JkBgAAAMjAw9LMjsHIDAAAAABHYmQGAAAAyIClmZ2DkRkAAAAAjkQxAwAAAMCRmGYGAAAAZGCzAIBjMDIDAAAAwJGMFzNz587VkCFD9OOPP0qSvv/+e91+++1q06aNJkyYYLh3AAAAKGhsj21sg2+MFjPvvfee7r77bn3zzTe6/fbb9fHHHysmJkYVKlRQlSpV1KdPH40ZM8ZkFwEAAADkUUavmXnrrbf07rvvqnv37lq8eLFuv/12jRw5Uk888YQk6frrr9eIESPUu3dvk90EAAAAkAcZHZnZsWOHoqOjJUm33HKL0tLSdNNNN6Ufv/nmm7Vz505T3QMAAEAB5LE9xjb4xmgxU7JkyfRiZd++fUpNTdWuXbvSj+/cuVMlSpTItA23263ExESvzZOWnKP9BgAAAGCe0Wlm7du312OPPaYuXbpo9uzZ6ty5s/r37y8/Pz9ZlqWBAweqdevWmbYRFxenYcOGee2rVKOLKtd8NCe7DgAAgHyKC/Gdw+jIzPDhw3XzzTdr+vTpatiwoSZMmKDHHntM7du3V9u2bVWyZEnFxcVl2kZsbKwSEhK8toqRnXLpGQAAAAAwxejITEhIyHnLLw8YMEBPPvmkUlJSVLRo0Uu24XK55HK5vPb5+Qdlaz8BAABQcDAy4xzG7zNzIcHBwSpatKh2796tbt26me4OAAAAgDwoTxYz5xw7dkxTpkwx3Q0AAAAAeZDRaWazZ8/O9Pj27dtzqScAAADA32ybaWZOYbSYiYmJkWVZmf7AWJaViz0CAAAA4BRGp5mFh4dr1qxZ8ng8F9xWr15tsnsAAAAogC72t2lubPCN0WKmSZMmio+Pv+jxS43aAAAAACi4jE4zGzhwoE6fPn3R45GRkVq8eHEu9ggAAACAUxgtZlq0aJHp8ZCQEEVFReVSbwAAAADuM+MkeXppZgAAAAC4GKMjMwAAAEBeY9tciO8UjMwAAAAAcCRGZgAAAIAMuGbGOfJlMRPzaOYLC+Skoge+NZZtBwYZyy7VqKaxbIUVNxZtB5j7FbLCShjLDq6aZiw7zTI3oByQmmQsOziyirHsZIPfc5PvayaVCEowF+7nbyw6OMnc864ZfvHVVXPaDw37Gcu+6cdRxrKB7MA0MwAAAACOlC9HZgAAAIDLxTQz52BkBgAAAIAjMTIDAAAAZOBhaWbHYGQGAAAAgCNRzAAAAABwJKaZAQAAABmwAIBzMDIDAAAAwJEYmQEAAAAysD0sAOAUjMwAAAAAcCTjIzNnz57Vp59+quXLl2v//v3y8/NTtWrVFBMTo5YtW5ruHgAAAAoYrplxDqMjM1u3blXt2rUVGxurhQsXav78+bIsS7/88ouio6N1//33KzU11WQXAQAAAORRRouZp59+Wm3atNGBAwe0a9cuxcXFyePxaOXKldq4caN++eUX/ec//zHZRQAAAAB5lNFiZunSperfv78sy5Ik9e3bVwsXLtTRo0dVo0YNjR49WlOmTMm0DbfbrcTERK8tNcWdG90HAABAPmTbHmMbfGO0mClWrJhOnjyZ/vWZM2eUmpqqoKAgSVL9+vW1f//+TNuIi4tTWFiY1/b9rNdytN8AAAAAzDNazNx2223q16+fNm3apB07dqhHjx5q2LChihYtKknatWuXypQpk2kbsbGxSkhI8Npuvee53Og+AAAA8iGPxza2wTdGVzMbMWKE2rdvrzp16siyLFWqVElffPFF+vHDhw9r4MCBmbbhcrnkcrm89gUEMkQHAAAA5HdGi5kyZcpoxYoV2rJli9xut2rVqqWAgP/v0n333WewdwAAAADysjxx08waNWqobt26XoWMJO3evVvdunUz1CsAAAAURLbHY2yDb/JEMXMxx44du+RqZgAAAAAKJqPTzGbPnp3p8e3bt+dSTwAAAIC/2VyI7xhGi5mYmBhZliXbvvgPzLl70AAAAABARkanmYWHh2vWrFnyeDwX3FavXm2yewAAAADyMKPFTJMmTRQfH3/R45catQEAAACym217jG3wjdFpZgMHDtTp06cvejwyMlKLFy/OxR4BAAAAcAqjxUyLFi0yPR4SEqKoqKhc6g0AAADAAgBOkqeXZgYAAACAizE6MgMAAADkNdy80jkYmQEAAADgSBQzAAAAAJzJhpekpCR76NChdlJSEtlkk0022WSTTTbZeTAbOMeybW7kklFiYqLCwsKUkJCg0NBQsskmm2yyySabbLLzWDZwDtPMAAAAADgSxQwAAAAAR6KYAQAAAOBIFDP/4HK5NHToULlcLrLJJptssskmm2yy82A2cA4LAAAAAABwJEZmAAAAADgSxQwAAAAAR6KYAQAAAOBIFDMAAAAAHIliJoN33nlHVapUUXBwsK677jqtWrUqV3J/+OEHtWvXTuXLl5dlWfryyy9zJTcuLk7XXHONihYtqjJlyigmJkabN2/Olexx48apfv36Cg0NVWhoqJo1a6Zvv/02V7L/6bXXXpNlWerTp0+OZ7344ouyLMtrq1WrVo7nnrN37149/PDDKlmypAoVKqR69erp119/zZXsKlWqnPfcLctSr169cjQ3LS1NgwcPVtWqVVWoUCFVr15dL7/8snJr7ZOTJ0+qT58+ioiIUKFChdS8eXP98ssvOZJ1qfcS27Y1ZMgQhYeHq1ChQmrVqpW2bNmS47mzZs1S69atVbJkSVmWpTVr1lxxZlbzU1JS9Oyzz6pevXoKCQlR+fLl1blzZ+3bty/Hs6W/f+dr1aqlkJAQFS9eXK1atdLPP/+cK9kZ9ejRQ5ZlafTo0bmS3bVr1/N+19u0aZPjuZK0ceNG3XXXXQoLC1NISIiuueYa7dq164qzs5J/ofc4y7L0+uuv52juqVOn9OSTT6pixYoqVKiQ6tSpo/Hjx19RZlazDx48qK5du6p8+fIqXLiw2rRpky3vK0BWUcz8z2effaZ+/fpp6NChWr16tRo0aKDo6GgdOnQox7NPnz6tBg0a6J133snxrIyWLl2qXr16aeXKlVqwYIFSUlLUunVrnT59OsezK1asqNdee03x8fH69ddfdeutt6p9+/basGFDjmdn9Msvv+i9995T/fr1cy3z6quv1v79+9O35cuX50ru8ePHdcMNNygwMFDffvut/vjjD40cOVLFixfPlfxffvnF63kvWLBAktShQ4cczR0+fLjGjRunt99+Wxs3btTw4cM1YsQIjR07Nkdzz3n88ce1YMECffTRR/r999/VunVrtWrVSnv37s32rEu9l4wYMUJvvfWWxo8fr59//lkhISGKjo5WUlJSjuaePn1aN954o4YPH35FOZeTf+bMGa1evVqDBw/W6tWrNWvWLG3evFl33XVXjmdL0lVXXaW3335bv//+u5YvX64qVaqodevWOnz4cI5nn/PFF19o5cqVKl++/BVn+pLdpk0br9/5Tz/9NMdzt23bphtvvFG1atXSkiVLtG7dOg0ePFjBwcFXnJ2V/IzPd//+/frwww9lWZbuvffeHM3t16+f5s2bp48//lgbN25Unz599OSTT2r27NlXlHupbNu2FRMTo+3bt+urr77Sb7/9poiICLVq1SpX/pYAJEk2bNu27Wuvvdbu1atX+tdpaWl2+fLl7bi4uFzthyT7iy++yNXMcw4dOmRLspcuXWokv3jx4vbEiRNzLe/kyZN2jRo17AULFthRUVF27969czxz6NChdoMGDXI850KeffZZ+8YbbzSSfSG9e/e2q1evbns8nhzNueOOO+xu3bp57bvnnnvsTp065Wiubdv2mTNnbH9/f3vu3Lle+xs3bmwPGjQoR7P/+V7i8XjscuXK2a+//nr6vhMnTtgul8v+9NNPcyw3ox07dtiS7N9++y3b8nzJP2fVqlW2JHvnzp25np2QkGBLshcuXJgr2Xv27LErVKhgr1+/3o6IiLDffPPNbM29WHaXLl3s9u3bZ3vWpXI7duxoP/zwwzmam1n+P7Vv396+9dZbczz36quvtl966SWvfTnxPvPP7M2bN9uS7PXr16fvS0tLs0uXLm2///772ZoNXAwjM5KSk5MVHx+vVq1ape/z8/NTq1attGLFCoM9y10JCQmSpBIlSuRqblpamqZPn67Tp0+rWbNmuZbbq1cv3XHHHV6ve27YsmWLypcvr2rVqqlTp07ZNv3hUmbPnq2mTZuqQ4cOKlOmjBo1aqT3338/V7L/KTk5WR9//LG6desmy7JyNKt58+ZatGiR/vzzT0nS2rVrtXz5crVt2zZHcyUpNTVVaWlp530qXKhQoVwbkTtnx44dOnDggNfPe1hYmK677roC9T4n/f1eZ1mWihUrlqu5ycnJmjBhgsLCwtSgQYMcz/N4PHrkkUc0cOBAXX311Tme909LlixRmTJlVLNmTfXs2VNHjx7N0TyPx6Ovv/5aV111laKjo1WmTBldd911uTZ1+58OHjyor7/+Wo899liOZzVv3lyzZ8/W3r17Zdu2Fi9erD///FOtW7fO0Vy32y1JXu9xfn5+crlcuf4eh4KLYkbSkSNHlJaWprJly3rtL1u2rA4cOGCoV7nL4/GoT58+uuGGG1S3bt1cyfz9999VpEgRuVwu9ejRQ1988YXq1KmTK9nTp0/X6tWrFRcXlyt551x33XWaPHmy5s2bp3HjxmnHjh1q0aKFTp48mePZ27dv17hx41SjRg3Nnz9fPXv21NNPP60pU6bkePY/ffnllzpx4oS6du2a41nPPfecHnjgAdWqVUuBgYFq1KiR+vTpo06dOuV4dtGiRdWsWTO9/PLL2rdvn9LS0vTxxx9rxYoV2r9/f47nZ3Tuvawgv89JUlJSkp599lk9+OCDCg0NzZXMuXPnqkiRIgoODtabb76pBQsWqFSpUjmeO3z4cAUEBOjpp5/O8ax/atOmjaZOnapFixZp+PDhWrp0qdq2bau0tLQcyzx06JBOnTql1157TW3atNF3332nu+++W/fcc4+WLl2aY7kXM2XKFBUtWlT33HNPjmeNHTtWderUUcWKFRUUFKQ2bdronXfe0U033ZSjubVq1VLlypUVGxur48ePKzk5WcOHD9eePXty/T0OBVeA6Q4gb+jVq5fWr1+fq5+k1KxZU2vWrFFCQoI+//xzdenSRUuXLs3xgmb37t3q3bu3FixYkG3zqLMq42hA/fr1dd111ykiIkIzZszI8U/vPB6PmjZtqldffVWS1KhRI61fv17jx49Xly5dcjT7nz744AO1bds2W+fwX8yMGTP0ySefaNq0abr66qu1Zs0a9enTR+XLl8+V5/3RRx+pW7duqlChgvz9/dW4cWM9+OCDio+Pz/FseEtJSdH9998v27Y1bty4XMu95ZZbtGbNGh05ckTvv/++7r//fv38888qU6ZMjmXGx8drzJgxWr16dY6Pfl7IAw88kP7/9erVU/369VW9enUtWbJELVu2zJFMj8cjSWrfvr369u0rSWrYsKF++uknjR8/XlFRUTmSezEffvihOnXqlCv/zowdO1YrV67U7NmzFRERoR9++EG9evVS+fLlc3T2QWBgoGbNmqXHHntMJUqUkL+/v1q1aqW2bdvm2iIrACMzkkqVKiV/f38dPHjQa//BgwdVrlw5Q73KPU8++aTmzp2rxYsXq2LFirmWGxQUpMjISDVp0kRxcXFq0KCBxowZk+O58fHxOnTokBo3bqyAgAAFBARo6dKleuuttxQQEJCjnxz+U7FixXTVVVdp69atOZ4VHh5+XqFYu3btXJvmds7OnTu1cOFCPf7447mSN3DgwPTRmXr16umRRx5R3759c21Urnr16lq6dKlOnTql3bt3a9WqVUpJSVG1atVyJf+cc+9lBfV97lwhs3PnTi1YsCDXRmUkKSQkRJGRkbr++uv1wQcfKCAgQB988EGOZi5btkyHDh1S5cqV09/ndu7cqf79+6tKlSo5mn0h1apVU6lSpXL0va5UqVIKCAjIE+9zy5Yt0+bNm3Plfe7s2bN6/vnnNWrUKLVr107169fXk08+qY4dO+qNN97I8fwmTZpozZo1OnHihPbv36958+bp6NGjuf4eh4KLYkZ//1HdpEkTLVq0KH2fx+PRokWLcvUajtxm27aefPJJffHFF/r+++9VtWpVo/3xeDzp829zUsuWLfX7779rzZo16VvTpk3VqVMnrVmzRv7+/jneh3NOnTqlbdu2KTw8PMezbrjhhvOW3v7zzz8VERGR49kZTZo0SWXKlNEdd9yRK3lnzpyRn5/3W52/v3/6p7i5JSQkROHh4Tp+/Ljmz5+v9u3b52p+1apVVa5cOa/3ucTERP3888/5+n1O+v9CZsuWLVq4cKFKlixptD+58V73yCOPaN26dV7vc+XLl9fAgQM1f/78HM2+kD179ujo0aM5+l4XFBSka665Jk+8z33wwQdq0qRJrlwblZKSopSUFOPvc2FhYSpdurS2bNmiX3/9Ndff41BwMc3sf/r166cuXbqoadOmuvbaazV69GidPn1ajz76aI5nnzp1yuvTqh07dmjNmjUqUaKEKleunGO5vXr10rRp0/TVV1+paNGi6fPmw8LCVKhQoRzLlaTY2Fi1bdtWlStX1smTJzVt2jQtWbIkV/6RLVq06HnXBYWEhKhkyZI5fr3QgAED1K5dO0VERGjfvn0aOnSo/P399eCDD+ZoriT17dtXzZs316uvvqr7779fq1at0oQJEzRhwoQczz7H4/Fo0qRJ6tKliwICcuftp127dnrllVdUuXJlXX311frtt980atQodevWLVfy58+fL9u2VbNmTW3dulUDBw5UrVq1cuS95VLvJX369NF//vMf1ahRQ1WrVtXgwYNVvnx5xcTE5GjusWPHtGvXrvR7u5z7Y7NcuXLZMiqUWX54eLjuu+8+rV69WnPnzlVaWlr6e12JEiUUFBSUY9klS5bUK6+8orvuukvh4eE6cuSI3nnnHe3duzdbliS/1Pf9n0VbYGCgypUrp5o1a+ZodokSJTRs2DDde++9KleunLZt26ZnnnlGkZGRio6OzrHcypUra+DAgerYsaNuuukm3XLLLZo3b57mzJmjJUuWXFFuVvOlvz8kmDlzpkaOHJktmVnJjYqK0sCBA1WoUCFFRERo6dKlmjp1qkaNGpXj2TNnzlTp0qVVuXJl/f777+rdu7diYmJyfPEBIJ3RtdTymLFjx9qVK1e2g4KC7GuvvdZeuXJlruQuXrzYlnTe1qVLlxzNvVCmJHvSpEk5mmvbtt2tWzc7IiLCDgoKskuXLm23bNnS/u6773I892Jya2nmjh072uHh4XZQUJBdoUIFu2PHjvbWrVtzPPecOXPm2HXr1rVdLpddq1Yte8KECbmWbdu2PX/+fFuSvXnz5lzLTExMtHv37m1XrlzZDg4OtqtVq2YPGjTIdrvduZL/2Wef2dWqVbODgoLscuXK2b169bJPnDiRI1mXei/xeDz24MGD7bJly9oul8tu2bJltrwWl8qdNGnSBY8PHTr0irMvlX9uOegLbYsXL87R7LNnz9p33323Xb58eTsoKMgODw+377rrLnvVqlVX/qQvkX0h2bk0c2bZZ86csVu3bm2XLl3aDgwMtCMiIuzu3bvbBw4cyNHccz744AM7MjLSDg4Oths0aGB/+eWXV5zrS/57771nFypUKFt/zy+Vu3//frtr1652+fLl7eDgYLtmzZr2yJEjs2Xp+0tljxkzxq5YsaIdGBhoV65c2X7hhRdy7f0VsG3btmybK7QAAAAAOA/XzAAAAABwJIoZAAAAAI5EMQMAAADAkShmAAAAADgSxQwAAAAAR6KYAQAAAOBIFDMAAAAAHIliBgAAAIAjUcwAAC7Isix9+eWXprsBAMBFUcwAKFC6du0qy7LO27Zu3Zot7U+ePFnFihXLlrYuV9euXRUTE2O0DwAA5IYA0x0AgNzWpk0bTZo0yWtf6dKlDfXm4lJSUhQYGGi6GwAA5FmMzAAocFwul8qVK+e1+fv7S5K++uorNW7cWMHBwapWrZqGDRum1NTU9MeOGjVK9erVU0hIiCpVqqQnnnhCp06dkiQtWbJEjz76qBISEtJHfF588UVJF56yVaxYMU2ePFmS9Ndff8myLH322WeKiopScHCwPvnkE0nSxIkTVbt2bQUHB6tWrVp69913fXq+N998s55++mk988wzKlGihMqVK5fer3O2bNmim266ScHBwapTp44WLFhwXju7d+/W/fffr2LFiqlEiRJq3769/vrrL0nSpk2bVLhwYU2bNi39/BkzZqhQoUL6448/fOovAABZRTEDAP+zbNkyde7cWb1799Yff/yh9957T5MnT9Yrr7ySfo6fn5/eeustbdiwQVOmTNH333+vZ555RpLUvHlzjR49WqGhodq/f7/279+vAQMG+NSH5557Tr1799bGjRsVHR2tTz75REOGDNErr7yijRs36tVXX9XgwYM1ZcoUn9qdMmWKQkJC9PPPP2vEiBF66aWX0gsWj8eje+65R0FBQfr55581fvx4Pfvss16PT0lJUXR0tIoWLaply5bpxx9/VJEiRdSmTRslJyerVq1aeuONN/TEE09o165d2rNnj3r06KHhw4erTp06PvUVAIAsswGgAOnSpYvt7+9vh4SEpG/33Xefbdu23bJlS/vVV1/1Ov+jjz6yw8PDL9rezJkz7ZIlS6Z/PWnSJDssLOy88yTZX3zxhde+sLAwe9KkSbZt2/aOHTtsSfbo0aO9zqlevbo9bdo0r30vv/yy3axZs0yfY/v27dO/joqKsm+88Uavc6655hr72WeftW3btufPn28HBATYe/fuTT/+7bffevX5o48+smvWrGl7PJ70c9xut12oUCF7/vz56fvuuOMOu0WLFnbLli3t1q1be50PAEB245oZAAXOLbfconHjxqV/HRISIklau3atfvzxR6+RmLS0NCUlJenMmTMqXLiwFi5cqLi4OG3atEmJiYlKTU31On6lmjZtmv7/p0+f1rZt2/TYY4+pe/fu6ftTU1MVFhbmU7v169f3+jo8PFyHDh2SJG3cuFGVKlVS+fLl0483a9bM6/y1a9dq69atKlq0qNf+pKQkbdu2Lf3rDz/8UFdddZX8/Py0YcMGWZblUz8BAPAFxQyAAickJESRkZHn7T916pSGDRume+6557xjwcHB+uuvv3TnnXeqZ8+eeuWVV1SiRAktX75cjz32mJKTkzMtZizLkm3bXvtSUlIu2LeM/ZGk999/X9ddd53Xeeeu8cmqfy4kYFmWPB5Plh9/6tQpNWnSJP06nowyLp6wdu1anT59Wn5+ftq/f7/Cw8N96icAAL6gmAGA/2ncuLE2b958wUJHkuLj4+XxeDRy5Ej5+f19yeGMGTO8zgkKClJaWtp5jy1durT279+f/vWWLVt05syZTPtTtmxZlS9fXtu3b1enTp18fTpZVrt2be3evdur+Fi5cqXXOY0bN9Znn32mMmXKKDQ09ILtHDt2TF27dtWgQYO0f/9+derUSatXr1ahQoVyrO8AgIKNBQAA4H+GDBmiqVOnatiwYdqwYYM2btyo6dOn64UXXpAkRUZGKiUlRWPHjtX27dv10Ucfafz48V5tVKlSRadOndKiRYt05MiR9ILl1ltv1dtvv63ffvtNv/76q3r06JGlZZeHDRumuLg4vfXWW/rzzz/1+++/a9KkSRo1alS2Pe9WrVrpqquuUpcuXbR27VotW7ZMgwYN8jqnU6dOKlWqlNq3b69ly5Zpx44dWrJkiZ5++mnt2bNHktSjRw9VqlRJL7zwgkaNGqW0tDSfF0AAAMAXFDMA8D/R0dGaO3euvvvuO11zzTW6/vrr9eabbyoiIkKS1KBBA40aNUrDhw9X3bp19cknnyguLs6rjebNm6tHjx7q2LGjSpcurREjRkiSRo4cqUqVKqlFixZ66KGHNGDAgCxdY/P4449r4sSJmjRpkurVq6eoqChNnjxZVatWzbbn7efnpy+++EJnz57Vtddeq8cff9zruiFJKly4sH744QdVrlxZ99xzj2rXrq3HHntMSUlJCg0N1dSpU/XNN9/oo48+UkBAgEJCQvTxxx/r/fff17fffpttfQUAICPL/uckbgAAAABwAEZmAAAAADgSxQwAAAAAR6KYAQAAAOBIFDMAAAAAHIliBgAAAIAjUcwAAAAAcCSKGQAAAACORDEDAAAAwJEoZgAAAAA4EsUMAAAAAEeimAEAAADgSP8HsjdSqOIZP+4AAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "from sklearn.metrics import mean_squared_error\n",
        "from sklearn.datasets import make_spd_matrix\n",
        "from scipy.optimize import nnls\n",
        "# Import the precise package or specific functions if necessary\n",
        "# from precise import <list_of_port_functions>\n",
        "\n",
        "def nnls_regression(X, y):\n",
        "    beta_nnls, rnorm = nnls(X, y)\n",
        "    return beta_nnls\n",
        "\n",
        "\n",
        "def test_scaled_regression():\n",
        "    np.random.seed(0)\n",
        "    # Setup and data generation as specified in your script\n",
        "\n",
        "    # Generate correlated X\n",
        "    n_samples =  50\n",
        "    train_proportion = 0.5\n",
        "    n_features = 20\n",
        "    n_informative = 5\n",
        "    beta_true = np.concatenate([\n",
        "    np.random.normal(loc=0.0, scale=0.1, size=n_features-n_informative),\n",
        "    np.random.normal(loc=1.0, scale=0.5, size=n_informative)\n",
        "    ])\n",
        "\n",
        "\n",
        "    # Initialize the covariance matrix\n",
        "    cov_matrix = np.zeros((n_features, n_features))\n",
        "\n",
        "    # Initialize the covariance matrix\n",
        "    cov_matrix = make_spd_matrix(n_features)\n",
        "    cov_matrix = (cov_matrix + cov_matrix.T) / 2\n",
        "\n",
        "\n",
        "    # Generate X data with hierarchical covariance structure\n",
        "    X = np.random.multivariate_normal(mean=np.zeros(n_features), cov=cov_matrix, size=n_samples)\n",
        "\n",
        "\n",
        "    # Generate y with some noise\n",
        "    y = X @ beta_true + np.random.randn(n_samples) * 5.\n",
        "\n",
        "    # Split the data into training and testing sets\n",
        "\n",
        "    split_idx = int(n_samples * train_proportion)\n",
        "    X_train, X_test = X[:split_idx], X[split_idx:]\n",
        "    y_train, y_test = y[:split_idx], y[split_idx:]\n",
        "\n",
        "\n",
        "\n",
        "    # Placeholder for storing MSE results\n",
        "    mse_results = {}\n",
        "\n",
        "    # Fit models on the training set\n",
        "    beta_custom = scaled_regression(X_train, y_train)\n",
        "    beta_port = unit_regression(X_train, y_train)\n",
        "    beta_hrp = hrp_regression(X_train, y_train)\n",
        "    beta_schur = schur_regression(X_train, y_train)\n",
        "    beta_nnls = nnls_regression(X_train, y_train)\n",
        "\n",
        "    # Regularized regressions\n",
        "    ridge_model = Ridge(alpha=1.0, fit_intercept=False)\n",
        "    ridge_model.fit(X_train, y_train)\n",
        "    beta_ridge = ridge_model.coef_\n",
        "\n",
        "    lasso_model = Lasso(alpha=0.1, fit_intercept=False)\n",
        "    lasso_model.fit(X_train, y_train)\n",
        "    beta_lasso = lasso_model.coef_\n",
        "\n",
        "    # Scale X for comparison with scikit-learn\n",
        "    scaler = StandardScaler()\n",
        "    X_train_scaled = scaler.fit_transform(X_train)\n",
        "    model = LinearRegression(fit_intercept=False)\n",
        "    model.fit(X_train_scaled, y_train)\n",
        "    beta_sklearn = model.coef_ / scaler.scale_\n",
        "\n",
        "    # Evaluate model performance on the test set\n",
        "    y_pred_sklearn = X_test @ beta_sklearn\n",
        "    y_pred_custom = X_test @ beta_custom\n",
        "    y_pred_port = X_test @ beta_port\n",
        "    y_pred_hrp = X_test @ beta_hrp\n",
        "    y_pred_schur = X_test @ beta_schur\n",
        "    y_pred_nnls = X_test @ beta_nnls\n",
        "\n",
        "    # Calculate and print out-of-sample mean squared error\n",
        "    mse_zero = mean_squared_error(y_test,np.zeros_like(y_test))\n",
        "    mse_custom = mean_squared_error(y_test, y_pred_custom)\n",
        "    mse_sklearn = mean_squared_error(y_test, y_pred_sklearn)\n",
        "    mse_port = mean_squared_error(y_test, y_pred_port)\n",
        "    mse_hrp = mean_squared_error(y_test, y_pred_hrp)\n",
        "    mse_schur = mean_squared_error(y_test, y_pred_schur)\n",
        "    mse_nnls = mean_squared_error(y_test, y_pred_nnls)\n",
        "    mse_ridge = mean_squared_error(y_test, X_test @ beta_ridge)\n",
        "    mse_lasso = mean_squared_error(y_test, X_test @ beta_lasso)\n",
        "\n",
        "\n",
        "    mse_results = {\n",
        "    'Sklearn lin reg      ': mse_sklearn,\n",
        "    'Zero prediction      ': mse_zero,\n",
        "    'Ridge                ': mse_ridge,\n",
        "    'Lasso                ': mse_lasso,\n",
        "    'Custom Implementation': mse_custom,\n",
        "    'Port Regression.     ': mse_port,\n",
        "    'HRP Regression.      ': mse_hrp,\n",
        "    'Schur Regression.    ': mse_schur,\n",
        "    'NNLS Regression.     ': mse_nnls\n",
        "    }\n",
        "\n",
        "\n",
        "    # Get the list of all port methods from precise\n",
        "    from precise.skaters.portfoliostatic.allstaticport import PORT as port_functions\n",
        "\n",
        "    for k, port_function in enumerate(port_functions):\n",
        "        print(port_function.__name__)\n",
        "        # Run regression with the current port function\n",
        "        beta_port = port_regression(X_train, y_train, port_function)\n",
        "\n",
        "        # Evaluate model performance on the test set\n",
        "        y_pred_port = X_test @ beta_port\n",
        "        mse_port = mean_squared_error(y_test, y_pred_port)\n",
        "\n",
        "        # Store the MSE result\n",
        "        mse_results[port_function.__name__] = mse_port\n",
        "\n",
        "        if (k % 20 ==1) or k==len(port_functions)-1:\n",
        "          # Sort and display MSE results from all regressions\n",
        "          sorted_mse = sorted(mse_results.items(), key=lambda item: item[1])\n",
        "          print(' ')\n",
        "          print(' ')\n",
        "          print(\"MSE values sorted from lowest to highest:\")\n",
        "          for name, mse in sorted_mse:\n",
        "              print(f\"{name} MSE: {mse}\")\n",
        "\n",
        "# Call the test function\n",
        "test_scaled_regression()\n",
        "\n",
        "print('Done')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QTx6InoA4jWH",
        "outputId": "b0ee22ee-1be2-49ad-9d99-042426b809ca"
      },
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "ppo_sharpe_long_port\n",
            "ppo_vol_long_port\n",
            " \n",
            " \n",
            "MSE values sorted from lowest to highest:\n",
            "HRP Regression.       MSE: 20.529032445089566\n",
            "ppo_sharpe_long_port MSE: 21.027663141648112\n",
            "ppo_vol_long_port MSE: 21.027663141648112\n",
            "Schur Regression.     MSE: 21.468808596555878\n",
            "Zero prediction       MSE: 22.084882881358794\n",
            "Ridge                 MSE: 23.08802200551705\n",
            "Lasso                 MSE: 26.324642233612632\n",
            "NNLS Regression.      MSE: 28.993185344984393\n",
            "Custom Implementation MSE: 84.5315967717654\n",
            "Sklearn lin reg       MSE: 84.53159677176784\n",
            "Port Regression.      MSE: 2910.8321083339983\n",
            "ppo_quad_long_port\n",
            "weak_long_port\n",
            "weak_h125_long_port\n",
            "weak_h150_long_port\n",
            "weak_h200_long_port\n",
            "weak_h400_long_port\n",
            "diag_long_port\n",
            "hrp_weak_weak_s5_long_port\n",
            "hrp_diag_diag_s5_long_port\n",
            "hrp_diag_weak_s5_long_port\n",
            "hrp_weak_diag_s5_long_port\n",
            "schur_diag_diag_s5_g100_long_port\n",
            "schur_weak_weak_s5_g100_long_port\n",
            "schur_diag_diag_s5_g050_long_port\n",
            "schur_weak_weak_s5_g050_long_port\n",
            "schur_weak_diag_s5_g100_long_port\n",
            "schur_weak_diag_s5_g050_long_port\n",
            "schur_weak_weak_s2_g100_long_port\n",
            "schur_weak_weak_s2_g050_long_port\n",
            "schur_weak_weak_s2_gcombo_long_port\n",
            " \n",
            " \n",
            "MSE values sorted from lowest to highest:\n",
            "schur_weak_weak_s2_gcombo_long_port MSE: 19.134836718960184\n",
            "schur_weak_weak_s5_g050_long_port MSE: 19.311820986394476\n",
            "schur_weak_weak_s5_g100_long_port MSE: 19.7675871916133\n",
            "schur_weak_weak_s2_g050_long_port MSE: 19.76821580863897\n",
            "schur_weak_diag_s5_g100_long_port MSE: 19.768538063518232\n",
            "schur_weak_weak_s2_g100_long_port MSE: 19.789400307404755\n",
            "hrp_weak_diag_s5_long_port MSE: 19.858717030060003\n",
            "schur_weak_diag_s5_g050_long_port MSE: 19.864045764936314\n",
            "diag_long_port MSE: 20.13437928466365\n",
            "hrp_weak_weak_s5_long_port MSE: 20.434951749560117\n",
            "HRP Regression.       MSE: 20.529032445089566\n",
            "hrp_diag_diag_s5_long_port MSE: 20.634545367699005\n",
            "weak_long_port MSE: 20.87391134278893\n",
            "weak_h125_long_port MSE: 20.87391134278893\n",
            "weak_h150_long_port MSE: 20.87391134278893\n",
            "weak_h200_long_port MSE: 20.87391134278893\n",
            "weak_h400_long_port MSE: 20.87391134278893\n",
            "schur_diag_diag_s5_g050_long_port MSE: 20.911259752347476\n",
            "hrp_diag_weak_s5_long_port MSE: 20.988040934668785\n",
            "ppo_sharpe_long_port MSE: 21.027663141648112\n",
            "ppo_vol_long_port MSE: 21.027663141648112\n",
            "ppo_quad_long_port MSE: 21.027663141648112\n",
            "schur_diag_diag_s5_g100_long_port MSE: 21.3061222194207\n",
            "Schur Regression.     MSE: 21.468808596555878\n",
            "Zero prediction       MSE: 22.084882881358794\n",
            "Ridge                 MSE: 23.08802200551705\n",
            "Lasso                 MSE: 26.324642233612632\n",
            "NNLS Regression.      MSE: 28.993185344984393\n",
            "Custom Implementation MSE: 84.5315967717654\n",
            "Sklearn lin reg       MSE: 84.53159677176784\n",
            "Port Regression.      MSE: 2910.8321083339983\n",
            "schur_weak_weak_s5_gcombo_long_port\n",
            "schur_diag_diag_s2_gcombo_long_port\n",
            "schur_diag_diag_s5_gcombo_long_port\n",
            "schur_diag_weak_s5_gcombo_long_port\n",
            "schur_weak_diag_s5_gcombo_long_port\n",
            "equal_long_port\n",
            "rp_port_p100\n",
            "rp_port_p70\n",
            "rp_port_p80\n",
            "rp_port_p90\n",
            "rp_port_p40\n",
            "rp_port_p50\n",
            "rp_port_p59\n",
            "rp_port_p60\n",
            "rp_port_p61\n",
            "rp_port_p10\n",
            "rp_port_p20\n",
            "rp_port_p30\n",
            "rp_port_p0\n",
            "ppo_sharpe_port\n",
            " \n",
            " \n",
            "MSE values sorted from lowest to highest:\n",
            "schur_weak_weak_s5_gcombo_long_port MSE: 18.759812868684357\n",
            "schur_weak_weak_s2_gcombo_long_port MSE: 19.134836718960184\n",
            "rp_port_p100 MSE: 19.226872311624806\n",
            "schur_weak_weak_s5_g050_long_port MSE: 19.311820986394476\n",
            "rp_port_p90 MSE: 19.314999300031683\n",
            "rp_port_p80 MSE: 19.395374542692107\n",
            "rp_port_p70 MSE: 19.469665198245412\n",
            "rp_port_p61 MSE: 19.53411136975739\n",
            "rp_port_p60 MSE: 19.541264829187686\n",
            "rp_port_p59 MSE: 19.548432189547736\n",
            "rp_port_p50 MSE: 19.61426049430657\n",
            "rp_port_p40 MSE: 19.69286045420432\n",
            "schur_weak_weak_s5_g100_long_port MSE: 19.7675871916133\n",
            "schur_weak_weak_s2_g050_long_port MSE: 19.76821580863897\n",
            "schur_weak_diag_s5_g100_long_port MSE: 19.768538063518232\n",
            "schur_weak_diag_s5_gcombo_long_port MSE: 19.776074595001926\n",
            "rp_port_p30 MSE: 19.781259289768762\n",
            "schur_weak_weak_s2_g100_long_port MSE: 19.789400307404755\n",
            "hrp_weak_diag_s5_long_port MSE: 19.858717030060003\n",
            "schur_weak_diag_s5_g050_long_port MSE: 19.864045764936314\n",
            "rp_port_p20 MSE: 19.88351999811004\n",
            "rp_port_p10 MSE: 20.002598843541012\n",
            "rp_port_p0 MSE: 20.134379284333324\n",
            "diag_long_port MSE: 20.13437928466365\n",
            "equal_long_port MSE: 20.134379284959596\n",
            "schur_diag_weak_s5_gcombo_long_port MSE: 20.26218140424147\n",
            "schur_diag_diag_s2_gcombo_long_port MSE: 20.284525637758012\n",
            "hrp_weak_weak_s5_long_port MSE: 20.434951749560117\n",
            "HRP Regression.       MSE: 20.529032445089566\n",
            "schur_diag_diag_s5_gcombo_long_port MSE: 20.607364656218515\n",
            "hrp_diag_diag_s5_long_port MSE: 20.634545367699005\n",
            "weak_long_port MSE: 20.87391134278893\n",
            "weak_h125_long_port MSE: 20.87391134278893\n",
            "weak_h150_long_port MSE: 20.87391134278893\n",
            "weak_h200_long_port MSE: 20.87391134278893\n",
            "weak_h400_long_port MSE: 20.87391134278893\n",
            "schur_diag_diag_s5_g050_long_port MSE: 20.911259752347476\n",
            "hrp_diag_weak_s5_long_port MSE: 20.988040934668785\n",
            "ppo_sharpe_long_port MSE: 21.027663141648112\n",
            "ppo_vol_long_port MSE: 21.027663141648112\n",
            "ppo_quad_long_port MSE: 21.027663141648112\n",
            "schur_diag_diag_s5_g100_long_port MSE: 21.3061222194207\n",
            "Schur Regression.     MSE: 21.468808596555878\n",
            "ppo_sharpe_port MSE: 21.97015533644057\n",
            "Zero prediction       MSE: 22.084882881358794\n",
            "Ridge                 MSE: 23.08802200551705\n",
            "Lasso                 MSE: 26.324642233612632\n",
            "NNLS Regression.      MSE: 28.993185344984393\n",
            "Custom Implementation MSE: 84.5315967717654\n",
            "Sklearn lin reg       MSE: 84.53159677176784\n",
            "Port Regression.      MSE: 2910.8321083339983\n",
            "ppo_vol_port\n",
            "ppo_quad_port\n",
            "hrp_unit_unit_s5_port\n",
            "hrp_unit_weak_s5_port\n",
            "hrp_unit_diag_s5_port\n",
            "schur_unit_unit_s5_g100_port\n",
            "schur_unit_unit_s5_g050_port\n",
            "unit_port\n",
            "unit_port_p000\n",
            "unit_port_p010\n",
            "unit_port_p020\n",
            "unit_port_p030\n",
            "unit_port_p040\n",
            "unit_port_p050\n",
            "unit_port_p060\n",
            "unit_port_p070\n",
            "unit_port_p080\n",
            "unit_port_p090\n",
            "unit_port_p100\n",
            " \n",
            " \n",
            "MSE values sorted from lowest to highest:\n",
            "schur_weak_weak_s5_gcombo_long_port MSE: 18.759812868684357\n",
            "schur_weak_weak_s2_gcombo_long_port MSE: 19.134836718960184\n",
            "rp_port_p100 MSE: 19.226872311624806\n",
            "schur_weak_weak_s5_g050_long_port MSE: 19.311820986394476\n",
            "rp_port_p90 MSE: 19.314999300031683\n",
            "rp_port_p80 MSE: 19.395374542692107\n",
            "rp_port_p70 MSE: 19.469665198245412\n",
            "rp_port_p61 MSE: 19.53411136975739\n",
            "rp_port_p60 MSE: 19.541264829187686\n",
            "rp_port_p59 MSE: 19.548432189547736\n",
            "rp_port_p50 MSE: 19.61426049430657\n",
            "rp_port_p40 MSE: 19.69286045420432\n",
            "schur_weak_weak_s5_g100_long_port MSE: 19.7675871916133\n",
            "schur_weak_weak_s2_g050_long_port MSE: 19.76821580863897\n",
            "schur_weak_diag_s5_g100_long_port MSE: 19.768538063518232\n",
            "schur_weak_diag_s5_gcombo_long_port MSE: 19.776074595001926\n",
            "rp_port_p30 MSE: 19.781259289768762\n",
            "schur_weak_weak_s2_g100_long_port MSE: 19.789400307404755\n",
            "unit_port_p040 MSE: 19.823230835739455\n",
            "unit_port_p030 MSE: 19.828162699950475\n",
            "hrp_weak_diag_s5_long_port MSE: 19.858717030060003\n",
            "schur_weak_diag_s5_g050_long_port MSE: 19.864045764936314\n",
            "unit_port_p050 MSE: 19.875159170302013\n",
            "unit_port_p020 MSE: 19.883196093946296\n",
            "rp_port_p20 MSE: 19.88351999811004\n",
            "hrp_unit_diag_s5_port MSE: 19.95865054505531\n",
            "unit_port_p010 MSE: 19.985773527421195\n",
            "unit_port_p060 MSE: 19.997947237908228\n",
            "rp_port_p10 MSE: 20.002598843541012\n",
            "rp_port_p0 MSE: 20.134379284333324\n",
            "unit_port_p000 MSE: 20.134379284367704\n",
            "diag_long_port MSE: 20.13437928466365\n",
            "equal_long_port MSE: 20.134379284959596\n",
            "unit_port_p070 MSE: 20.21686340327325\n",
            "schur_diag_weak_s5_gcombo_long_port MSE: 20.26218140424147\n",
            "schur_diag_diag_s2_gcombo_long_port MSE: 20.284525637758012\n",
            "hrp_weak_weak_s5_long_port MSE: 20.434951749560117\n",
            "HRP Regression.       MSE: 20.529032445089566\n",
            "unit_port_p080 MSE: 20.572723179792543\n",
            "schur_diag_diag_s5_gcombo_long_port MSE: 20.607364656218515\n",
            "hrp_diag_diag_s5_long_port MSE: 20.634545367699005\n",
            "hrp_unit_unit_s5_port MSE: 20.65816147726662\n",
            "weak_long_port MSE: 20.87391134278893\n",
            "weak_h125_long_port MSE: 20.87391134278893\n",
            "weak_h150_long_port MSE: 20.87391134278893\n",
            "weak_h200_long_port MSE: 20.87391134278893\n",
            "weak_h400_long_port MSE: 20.87391134278893\n",
            "schur_diag_diag_s5_g050_long_port MSE: 20.911259752347476\n",
            "hrp_diag_weak_s5_long_port MSE: 20.988040934668785\n",
            "ppo_sharpe_long_port MSE: 21.027663141648112\n",
            "ppo_vol_long_port MSE: 21.027663141648112\n",
            "ppo_quad_long_port MSE: 21.027663141648112\n",
            "unit_port_p090 MSE: 21.117010554894932\n",
            "hrp_unit_weak_s5_port MSE: 21.249247040657075\n",
            "schur_diag_diag_s5_g100_long_port MSE: 21.3061222194207\n",
            "Schur Regression.     MSE: 21.468808596555878\n",
            "schur_unit_unit_s5_g050_port MSE: 21.571861290767217\n",
            "ppo_sharpe_port MSE: 21.97015533644057\n",
            "ppo_vol_port MSE: 21.97015533644057\n",
            "ppo_quad_port MSE: 21.97015533644057\n",
            "unit_port MSE: 21.970161424505562\n",
            "unit_port_p100 MSE: 21.970161424505562\n",
            "Zero prediction       MSE: 22.084882881358794\n",
            "schur_unit_unit_s5_g100_port MSE: 22.617404426826905\n",
            "Ridge                 MSE: 23.08802200551705\n",
            "Lasso                 MSE: 26.324642233612632\n",
            "NNLS Regression.      MSE: 28.993185344984393\n",
            "Custom Implementation MSE: 84.5315967717654\n",
            "Sklearn lin reg       MSE: 84.53159677176784\n",
            "Port Regression.      MSE: 2910.8321083339983\n",
            "Done\n"
          ]
        }
      ]
    }
  ]
}