Surrogates/PolynomialChaosExpansion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "source": [
    "import numpy as np\n",
    "\n",
    "factorial = np.math.factorial\n",
    "\n",
    "from itertools import permutations, combinations\n",
    "import matplotlib.pyplot as plt"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Orthonormal Polynomials\n",
    "\n",
    "## _Hermite Polynomials_\n",
    "\n",
    "Recursion relation : \n",
    "${\\mathit {He}}_{n+1}(x)=x{\\mathit {He}}_{n}(x)-n{\\mathit {He}}_{n-1}(x)$\n",
    "\n",
    "alternatively,${\\mathit {He}}_{n}(x)=x{\\mathit {He}}_{n-1}(x)-(n-1){\\mathit {He}}_{n-2}(x)$\n",
    "\n",
    "Normalisation constant: $\\sqrt{n!}$\n",
    "\n",
    "Normalised Polynomial: $ \\tilde{P_{n}}(x) = \\frac{P_{n}(x)} {\\sqrt{n!}} $\n",
    "\n",
    "## _Legendre Polynomials_\n",
    " \n",
    "Recursion relation : \n",
    "$ (n+1)P_{n+1}(x)=(2n+1)xP_{n}(x)-nP_{n-1}(x) $\n",
    "\n",
    "alternatively, $ P_{n}(x)=\\frac{(2n-1)xP_{n-1}(x)-(n-1)P_{n-2}(x)}{n} $\n",
    "\n",
    "Normalisation constant: $\\frac{1}{\\sqrt{(2n+1)}}$\n",
    "\n",
    "Normalised Polynomial: $ \\tilde{P_{n}}(x) = P_{n}(x) \\cdot \\sqrt{(2n+1)} $\n",
    "\n"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "source": [
    "'''\n",
    "Input:\n",
    "    n: degree of polynomial\n",
    "    x: random variable\n",
    "\n",
    "Output:\n",
    "    out: normalised polynomial of degree 'n'\n",
    "         evaluated at 'x'\n",
    "\n",
    "'''\n",
    "\n",
    "def Hermite(n, x):\n",
    "    \n",
    "    def recurence_relation(n, x):\n",
    "    \n",
    "        if n<=0:\n",
    "            return 1\n",
    "\n",
    "        else:\n",
    "            return x*recurence_relation(n-1,x) - (n-1)*recurence_relation(n-2,x)\n",
    "    \n",
    "    # normalise\n",
    "    out = recurence_relation(n,x)/np.sqrt(factorial(n))\n",
    "    \n",
    "    return out\n",
    "\n",
    "\n",
    "def Legendre(n,x):\n",
    "                 \n",
    "    def recurence_relation(n,x):    \n",
    "        \n",
    "        if n<=0:\n",
    "            return 1\n",
    "\n",
    "        else:\n",
    "            return ((2*n-1)*x*recurence_relation(n-1,x))/n - (n-1) * (recurence_relation(n-2,x)/n)\n",
    "\n",
    "    # normalise        \n",
    "    out = recurence_relation(n,x)*np.sqrt(2*n+1)\n",
    "            \n",
    "    return out"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Combinatorics\n",
    "\n",
    "This section is divided into three simple sub-sections:\n",
    "\n",
    "1. method: __combination_for_a_given_sum__ </p>\n",
    "This function computes the sets of natural numbers (<= n) whose sum is 'n'. </p>\n",
    "Example: For n = 3, we get {{1,1,1}, {1,2}, {3}}. The function also discards set with cardinality greater than N_p. </p>\n",
    "\n",
    "2. method:  __find_permutations__ </p>\n",
    "This function finds the permutations of the sets from the above function in a set of size N_p. </p>\n",
    "Example: For set q = {1,2}, and N_p = 3, we get Q = {(0, 2, 1), (1, 2, 0), (2, 1, 0), (2, 0, 1), (0, 1, 2), (1, 0, 2)}. </p>\n",
    "    \n",
    "3. __A for loop__ </p>\n",
    "   To iterate over each q from __combination_for_a_given_sum__ to get corresponding Q and store them in the variable 'all_permutations'\n",
    "    \n",
    "recursion in first function adapted from [this source](https://www.techiedelight.com/print-all-combination-numbers-from-1-to-n/)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "source": [
    "def compute_all_permutations(n, N_p):\n",
    "    \n",
    "    '''\n",
    "    Input :\n",
    "        n : total degree of polynomial\n",
    "        N_p : Number of random variables\n",
    "        \n",
    "    Output: \n",
    "        all_permutations: [{}, {}, ..]\n",
    "                          All sets sets of size N_p containing all combinations \n",
    "                          of whole numbers whose sum is equal \n",
    "                          to 'n'.\n",
    "        comb_dict: A dictionary which lists above combinations systematically\n",
    "    \n",
    "    '''\n",
    "\n",
    "    def combination_for_a_given_sum(n):\n",
    "\n",
    "        '''\n",
    "            Input: \n",
    "               n : [natural number]\n",
    "                   sum of natural numbers\n",
    "\n",
    "            Output: \n",
    "            store_combinations: [{}, {}, ..]\n",
    "                                sets of natural numbers\n",
    "                                Each set has size <= N_p and \n",
    "                                sum of elements = n\n",
    "        '''\n",
    "\n",
    "        def compute_Combinations(i, n, out, index):\n",
    "\n",
    "            '''\n",
    "            Recursive function to compute all combinations of numbers from 'i' to 'n' having sum 'n'. \n",
    "            The 'index' denotes the next free slot in the output list 'out'.\n",
    "\n",
    "            '''\n",
    "\n",
    "            # if the sum becomes 'n', store the combination\n",
    "            if n == 0:\n",
    "                combination = out[:index]\n",
    "                store_combinations.append(combination)\n",
    "\n",
    "            # start from the previous element in the combination till 'n'\n",
    "            for j in range(i, n + 1):\n",
    "\n",
    "                # place current element at the current index\n",
    "                out[index] = j\n",
    "\n",
    "                # recur with a reduced sum\n",
    "                compute_Combinations(j, n - j, out, index + 1)\n",
    "\n",
    "        out = [None] * n\n",
    "\n",
    "        # store the combinations\n",
    "        store_combinations = []  \n",
    "\n",
    "        # store trimmed combinations\n",
    "        trimmed_combinations = []\n",
    "        \n",
    "        # compute all combinations of numbers from 1 to 'n' having sum 'n'\n",
    "        compute_Combinations(1, n, out, 0)\n",
    "\n",
    "        # Remove combinations with number of elements greater than N_p\n",
    "        for combination in store_combinations:\n",
    "            if len(combination) <= N_p:\n",
    "                trimmed_combinations.append(combination)\n",
    "        \n",
    "        return trimmed_combinations\n",
    "\n",
    "\n",
    "    def find_permutations(store_combinations):\n",
    "\n",
    "        store = set() # an empty set to store all permutations for a given sum 'p'\n",
    "\n",
    "        for i, combination in enumerate(store_combinations):\n",
    "            \n",
    "            # create an array 'C' of size 'N_p' filled with 'combination' elements \n",
    "            # and padded with zeros if necessary\n",
    "            C = np.zeros(N_p, dtype = int)\n",
    "            C[0:len(combination)] = combination\n",
    "            \n",
    "            # find all permutations of 'combination' elements\n",
    "            # 'set' method removes duplicates, as 'combination' sometimes contains repeating elements        \n",
    "            perm = set(permutations(C))\n",
    "\n",
    "            # add to previous set and store the permutation \n",
    "            store = store.union(store, perm)\n",
    "\n",
    "        return store\n",
    "\n",
    "\n",
    "    comb_dict = {}\n",
    "    \n",
    "    all_permutations = set()\n",
    "    \n",
    "    for i in range(1, n+1):\n",
    "\n",
    "        store_combinations = combination_for_a_given_sum(i)\n",
    "\n",
    "        all_permuatations_for_given_sum = find_permutations(store_combinations)\n",
    "\n",
    "        # store in dictionary \n",
    "        comb_dict[i] = all_permuatations_for_given_sum\n",
    "\n",
    "        # store in set\n",
    "        all_permutations = all_permutations.union(all_permutations,all_permuatations_for_given_sum)\n",
    "\n",
    "    return all_permutations, comb_dict"
   ],
   "outputs": [],
   "metadata": {
    "tags": []
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# PCE Surrogate Class"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "source": [
    "class PCE_surrogate():\n",
    "    \"\"\"\n",
    "    Create PCE surrogate\n",
    "\n",
    "    Refer UQLab Manual for mathematical details\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    total_polynomial_degree : int\n",
    "                              Total polynomial degree for PCE surrogate\n",
    "    polynomial_classes_of_random_variables : list of strings\n",
    "                                              List specifying polynomial class of each random variable\n",
    "                                              Uniform: 'Legendre'\n",
    "                                              Normal: 'Hermite'\n",
    "                                              (Only above RV types supported)\n",
    "    isoprob_transform : method\n",
    "                        Isoprobabilistic transform required for PCE surrogate\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, total_polynomial_degree, polynomial_classes_of_random_variables, isoprob_transform):\n",
    "\n",
    "        self.n = total_polynomial_degree\n",
    "        self.pcrv = polynomial_classes_of_random_variables\n",
    "        self.N_p = len(self.pcrv)\n",
    "        self.isoprob_transform = isoprob_transform \n",
    "\n",
    "        #total number of polynomial terms in the PCE\n",
    "        self.number_of_PCE_terms = int (factorial(self.N_p+self.n) / ( factorial(self.N_p) * factorial(self.n) ))\n",
    "        \n",
    "        self.all_permutations, self.comb_dict = compute_all_permutations(self.n, self.N_p)\n",
    "\n",
    "    def print_combinations(self):\n",
    "        \n",
    "        print(\"Number of polynomial terms:\", self.number_of_PCE_terms)\n",
    "        \n",
    "        for key, value in self.comb_dict.items():\n",
    "\n",
    "            print(\"Sum =\",key,\":\",value, '\\n')\n",
    "\n",
    "\n",
    "    def phi(self, alpha, X):\n",
    "        \"\"\"\n",
    "        Orthogonal Multivariate Polynomial Basis\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        alpha : list\n",
    "                powers of random variables\n",
    "        X : ndarray\n",
    "            [N, N_p]\n",
    "            N : Number of MC samples\n",
    "            N_p : Number of random variables\n",
    "            Matrix of random variables\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        ndarray\n",
    "            Evaluations of X in the PCE basis\n",
    "        \"\"\"\n",
    "\n",
    "        phi_value = 1\n",
    "\n",
    "        for i in range(self.N_p):\n",
    "\n",
    "            phi_value *= self.pcrv[i](alpha[i],X[:,i])\n",
    "\n",
    "        return phi_value\n",
    "\n",
    "\n",
    "    def compute_vanderMonde_matrix(self, any_X):\n",
    "\n",
    "        # first column must contain ones\n",
    "        # other columns will be replaced by polynomial evaluations\n",
    "        vanderMonde_matrix =  np.ones((any_X.shape[0],self.number_of_PCE_terms))\n",
    "\n",
    "        for i, alpha in enumerate (self.all_permutations):\n",
    "\n",
    "            vanderMonde_matrix[:,i+1] = self.phi(alpha, any_X)\n",
    "\n",
    "        return vanderMonde_matrix\n",
    "\n",
    "\n",
    "    def Moore_Penrose_inverse(self, X, Y):\n",
    "\n",
    "        N = np.matmul(X.T,X)\n",
    "        M = np.matmul(X.T,Y)\n",
    "        beta = np.linalg.solve(N,M)\n",
    "\n",
    "        return beta \n",
    "\n",
    "\n",
    "    def find_coefficients(self, X, Y):\n",
    "        \"\"\"\n",
    "        Find coefficients of the model in the polynomial basis surrogate\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            [N, N_p]\n",
    "            N : Number of MC samples\n",
    "            N_p : Number of random variables\n",
    "            Matrix of random variables\n",
    "        Y : ndarray\n",
    "            [N, N_p]\n",
    "            N : Number of MC samples\n",
    "            Model evaluations for 'N' MC samples\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        beta: ndarray\n",
    "              Vector of coefficients of model in the polynomial basis\n",
    "        \"\"\"\n",
    "\n",
    "        # iso_probabilistic transform\n",
    "        # apply transform to the samples \n",
    "        # to convert to the range [-1,1]\n",
    "        X = self.isoprob_transform(X)\n",
    "\n",
    "        vanderMonde_train =  self.compute_vanderMonde_matrix(X) \n",
    "\n",
    "        # compute coefficients\n",
    "        self.beta = self.Moore_Penrose_inverse(vanderMonde_train, Y)\n",
    "\n",
    "        return self.beta\n",
    "\n",
    "\n",
    "    def evaluate(self, X):\n",
    "        \"\"\"\n",
    "        Evaluate random variables in the polynomial basis surrogate created\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            [N, N_p]\n",
    "            N : Number of MC samples\n",
    "            N_p : Number of random variables\n",
    "            Matrix of random variables\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        Y_hat : ndarray\n",
    "                [N, 1]\n",
    "                N : Number of MC samples\n",
    "                Surrogate evaluations for 'N' MC samples\n",
    "        \"\"\"\n",
    "        X_test = self.isoprob_transform(X)\n",
    "\n",
    "        vanderMonde_test = self.compute_vanderMonde_matrix(X_test)\n",
    "\n",
    "        Y_hat = vanderMonde_test @ self.beta\n",
    "\n",
    "        return Y_hat\n",
    "\n",
    "\n",
    "    def coefficient_pickers(self):\n",
    "        \"\"\"\n",
    "        Coefficients necessary for Sobol Indices\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        first_order_index_coefficient_picker : ndarray\n",
    "                                               N_p, number_of_PCE_terms\n",
    "        total_index_coefficient_pickerndarray : ndarray\n",
    "                                                N_p, number_of_PCE_terms\n",
    "        \"\"\"\n",
    "\n",
    "        coefficient_store = np.zeros((self.N_p, self.number_of_PCE_terms), dtype = int)\n",
    "\n",
    "        for i, alpha in enumerate (self.all_permutations):\n",
    "            coefficient_store[:,i+1] = list(alpha)\n",
    "\n",
    "        # create coefficient pickers : Matrices which multiply coefficients by 0 or 1 based on\n",
    "        # first order or total Sobol index\n",
    "\n",
    "        first_order_index_coefficient_picker = copy.deepcopy(coefficient_store)\n",
    "\n",
    "        # Coefficient picker for First order index \n",
    "        # set to zero if other coefficients are present\n",
    "        for i in range(self.number_of_PCE_terms):\n",
    "            column_i = first_order_index_coefficient_picker[:,i]\n",
    "\n",
    "            if np.count_nonzero(column_i) > 1:\n",
    "                first_order_index_coefficient_picker[:,i] = np.zeros((self.N_p))\n",
    "\n",
    "        # set all non-zero elements to 1\n",
    "        first_order_index_coefficient_picker = np.where(first_order_index_coefficient_picker == 0, 0, 1)\n",
    "        total_index_coefficient_picker = np.where(coefficient_store == 0, 0, 1)\n",
    "\n",
    "        return first_order_index_coefficient_picker, total_index_coefficient_picker\n",
    "\n",
    "\n",
    "    def LeaveOneOut(self, X, Y):\n",
    "        \"\"\"\n",
    "        Leave One Out error for PCE surrogate\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            [N, N_p]\n",
    "            N : Number of MC samples\n",
    "            N_p : Number of random variables\n",
    "            Matrix of random variables\n",
    "        Y : ndarray\n",
    "            [description]\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        float\n",
    "        Error in PCE approximation\n",
    "\n",
    "        Notes\n",
    "        -----\n",
    "        [1] UQLab Manual: PolynomialChaos Section 1.4.3.2\n",
    "\n",
    "        \"\"\"\n",
    "\n",
    "        # total variance \n",
    "        denominator = np.var(Y)\n",
    "\n",
    "        # H matrix\n",
    "        A = self.compute_vanderMonde_matrix(X)\n",
    "        N = A.T @ A\n",
    "        H = A @ np.linalg.solve(N,A.T)\n",
    "        h = H.diagonal()\n",
    "\n",
    "        # evalute model\n",
    "        Y_hat = self.evaluate(X)\n",
    "\n",
    "        # numerator\n",
    "        numerator = ((Y - Y_hat) / (1-h))**2\n",
    "        numerator = numerator.sum()/Y.size\n",
    "\n",
    "        return numerator/denominator"
   ],
   "outputs": [],
   "metadata": {}
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3613jvsc74a57bd081f210cd5dcd35cb0a15aa45415a139281dabdc726128a552d7e7a882a257bd3",
   "display_name": "Python 3.6.13 64-bit ('uq_matlab': conda)"
  },
  "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.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}