update

equadratures/__init__.py

@@ -9,6 +9,7 @@
 from equadratures.weight import Weight
 from equadratures.poly import evaluate_model, evaluate_model_gradients, vector_to_2D_grid
 from equadratures.polytree import PolyTree
+from equadratures.plot import Plot
 from equadratures.scalers import *
 import numpy as np
 import os, sys

equadratures/parameter.py

@@ -17,10 +17,11 @@
 from equadratures.distributions.gumbel import Gumbel
 from equadratures.distributions.chi import Chi
 from equadratures.distributions.analytical import Analytical
+from equadratures.plot import Plot
 import numpy as np
 import scipy as sc
 
-class Parameter(object):
+class Parameter(Plot):
     """
     This class defines a univariate parameter. Below are details of its constructor.
 

equadratures/plot.py

@@ -0,0 +1,134 @@
+import seaborn as sns
+import matplotlib.pyplot as plt
+from equadratures.datasets import score
+import numpy as np
+sns.set(font_scale=1.5)
+sns.set_style("white")
+sns.set_style("ticks")
+
+class Plot:
+    """
+    Plotting utilities.
+    """
+    def plot_pdf(self, data=None, save=False, xlim=None, ylim=None):
+        """
+        Plots the probability density function for a Parameter.
+        """
+        s_values, pdf = self.get_pdf()
+        fig = plt.figure(figsize=(8,6))
+        ax = fig.add_subplot(1,1,1)
+        plt.fill_between(s_values,  pdf*0.0, pdf, color="gold" , label='Density', interpolate=True, hatch="\\\\\\\\", edgecolor="grey",  linewidth=0.5,alpha=0.5)
+        if data is not None:
+            plt.hist(data, 50, density=True, facecolor='dodgerblue', alpha=0.7, label='Data', edgecolor='white')
+        plt.xlabel(self.variable.capitalize())
+        plt.ylabel('PDF')
+        if xlim is not None:
+            plt.xlim([xlim[0], xlim[1]])
+        if ylim is not None:
+            plt.ylim([ylim[0], ylim[1]])
+        plt.legend()
+        sns.despine(offset=10, trim=True)
+        if save:
+            plt.savefig('pdf_plot.png', dpi=140, bbox_inches='tight')
+        else:
+            plt.show()
+    def plot_orthogonal_polynomials(self, order_limit=None, number_of_points=200, save=False, xlim=None, ylim=None):
+        """
+        Plots the first K orthogonal polynomials.
+
+        :param Parameter self: An instance of the Parameter class.
+        """
+        Xi = np.linspace(self.distribution.x_range_for_pdf[0], \
+                    self.distribution.x_range_for_pdf[-1], number_of_points).reshape(number_of_points, 1)
+        P, _, _ = self._get_orthogonal_polynomial(Xi)
+        fig = plt.figure(figsize=(8,6))
+        ax = fig.add_subplot(1,1,1)
+        if order_limit is None:
+            max_order = P.shape[0]
+        else:
+            max_order = order_limit
+        for i in range(0, max_order):
+            plt.plot(Xi, P[i,:], '-', lw=2, label='order %d'%(i))
+        if xlim is not None:
+            ax.set_xlim([xlim[0], xlim[1]])
+        if ylim is not None:
+            ax.set_ylim([ylim[0], ylim[1]])
+        plt.legend()
+        sns.despine(offset=10, trim=True)
+        plt.xlabel(self.name.capitalize()+' parameter')
+        plt.ylabel('Orthogonal polynomials')
+        if save:
+            plt.savefig('polyfit_1D_plot.png', dpi=140, bbox_inches='tight')
+        else:
+            plt.show()
+    def plot_polyfit_1D(self, uncertainty=True, output_variances=None, number_of_points=200, save=False, xlim=None, ylim=None):
+        """
+        Plots a univariate polynomial.
+        """
+        if self.dimensions != 1:
+            raise(ValueError, 'plot_polyfit_1D is only meant for univariate polynomials.')
+        Xi = np.linspace(self.parameters[0].distribution.x_range_for_pdf[0], \
+                    self.parameters[0].distribution.x_range_for_pdf[-1], number_of_points).reshape(number_of_points, 1)
+        if uncertainty:
+            if output_variances is None:
+                y, ystd = self.get_polyfit(Xi,uq=True)
+            else:
+                self.output_variances = output_variances
+                y, ystd = self.get_polyfit(Xi,uq=True)
+            ystd = ystd.squeeze()
+        else:
+            y = self.get_polyfit(Xi)
+        y = y.squeeze()
+        X = self.get_points()
+        y_truth = self._model_evaluations
+        fig = plt.figure(figsize=(8,6))
+        ax = fig.add_subplot(1,1,1)
+        plt.plot(Xi, y, '-',label='Polynomial fit', color='navy')
+        plt.plot(X.flatten(), y_truth.flatten(), 'o', color='dodgerblue', ms=10, markeredgecolor='k',lw=1, alpha=0.6, label='Data')
+        if uncertainty:
+            ax.fill_between(Xi.flatten(), y+ystd, y-ystd, alpha=.10, color='deepskyblue',label='Polynomial $\sigma$')
+        if xlim is not None:
+            ax.set_xlim([xlim[0], xlim[1]])
+        if ylim is not None:
+            ax.set_ylim([ylim[0], ylim[1]])
+        plt.legend()
+        sns.despine(offset=10, trim=True)
+        plt.xlabel(self.parameters[0].variable.capitalize())
+        plt.ylabel('Polynomial fit')
+        if save:
+            plt.savefig('polyfit_1D_plot.png', dpi=140, bbox_inches='tight')
+        else:
+            plt.show()
+
+    def plot_model_vs_data(self, sample_data=None, metric='adjusted_r2', save=False, xlim=None, ylim=None):
+        """
+        Plots the polynomial approximation against the true data.
+
+        :param Poly self: An instance of the Poly class.
+
+        """
+        if sample_data is None:
+            X = self.get_points()
+            y_truth = self._model_evaluations
+            y_model = self.get_polyfit(X)
+        else:
+            X, y_truth = sample_data[0], sample_data[1]
+            y_model = self.get_polyfit(X)
+        R2score = score(y_truth, y_model, metric, X)
+        fig = plt.figure(figsize=(8,6))
+        ax = fig.add_subplot(1,1,1)
+        plt.plot(y_model, y_truth, 'o', color='dodgerblue', ms=10, markeredgecolor='k',lw=1, alpha=0.6)
+        plt.xlabel('Polynomial model')
+        plt.ylabel('True data')
+        displaytext = '$R^2$ = '+str(np.round(float(R2score), 2))
+        if xlim is not None:
+            plt.xlim([xlim[0], xlim[1]])
+        if ylim is not None:
+            plt.ylim([ylim[0], ylim[1]])
+        plt.text(0.3, 0.9, displaytext, transform=ax.transAxes, \
+            horizontalalignment='center', verticalalignment='center', fontsize=14, color='grey')
+        sns.despine(offset=10, trim=True)
+        if save:
+            plt.savefig('model_vs_data_plot.png', dpi=140, bbox_inches='tight')
+        else:
+            plt.show()

equadratures/poly.py

@@ -6,11 +6,13 @@
 from equadratures.subsampling import Subsampling
 from equadratures.quadrature import Quadrature
 from equadratures.datasets import score
+from equadratures.plot import Plot
 import scipy.stats as st
 import numpy as np
 from copy import deepcopy
 MAXIMUM_ORDER_FOR_STATS = 8
-class Poly(object):
+
+class Poly(Plot):
     """
     Definition of a polynomial object.
 
@@ -815,7 +817,7 @@ def get_poly_hess(self, stack_of_points):
         return H
     def get_polyscore(self,X_test=None,y_test=None,metric='adjusted_r2'):
         """
-        Evaluates the accuracy of the polynomial approximation using the selected accuracy metric. Training accuracy is evaluated on the data used for fitting the polynomial. Testing accuracy is evaluated on new data if it is provided by the ``X_test`` and ``y_test`` arguments (both must be provided together). 
+        Evaluates the accuracy of the polynomial approximation using the selected accuracy metric. Training accuracy is evaluated on the data used for fitting the polynomial. Testing accuracy is evaluated on new data if it is provided by the ``X_test`` and ``y_test`` arguments (both must be provided together).
 
         :param Poly self:
             An instance of the Poly class.
@@ -824,7 +826,7 @@ def get_polyscore(self,X_test=None,y_test=None,metric='adjusted_r2'):
         :param numpy.ndarray y_test:
             An ndarray with shape (number_of_observations, 1) containing new ``y_test`` data (optional).
         :param string metric:
-            An optional string containing the scoring metric to use. Avaliable options are: ``adjusted_r2``, ``r2``, ``mae``, ``rmse``, or ``normalised_mae`` (default: ``adjusted_r2``). 
+            An optional string containing the scoring metric to use. Avaliable options are: ``adjusted_r2``, ``r2``, ``mae``, ``rmse``, or ``normalised_mae`` (default: ``adjusted_r2``).
 
         :return:
             **score_train**: The training score of the model, output as a float.
@@ -840,7 +842,6 @@ def get_polyscore(self,X_test=None,y_test=None,metric='adjusted_r2'):
             return train_score, test_score
         else:
             return train_score
-
     def _get_polystd(self, stack_of_points):
         """
         Private function to evaluate the uncertainty of the polynomial approximation at prescribed points, following the approach from [7].
@@ -853,8 +854,8 @@ def _get_polystd(self, stack_of_points):
             **y_std**: A numpy.ndarray of shape (number_of_observations,1) corresponding to the uncertainty (one standard deviation) of the polynomial approximation at each point.
         """
         # Training data
-        X_train = self.inputs
-        y_train = self.outputs
+        X_train = self._quadrature_points
+        y_train = self._model_evaluations
 
         # Define covariance matrix - TODO: allow non-diagonal matrix?
         # Empirical variance
@@ -882,7 +883,7 @@ def _get_polystd(self, stack_of_points):
 
         # Propagate the uncertainties
         Sigma_X = np.dot( np.dot(Q, Sigma), Q.T)
-        Sigma_F = np.dot( np.dot(Ao, Sigma_X), Ao.T) 
+        Sigma_F = np.dot( np.dot(Ao, Sigma_X), Ao.T)
         std_F = np.sqrt( np.diag(Sigma_F) )
         return std_F.reshape(-1,1)
 

source/_documentation/tutorial_1.txt

@@ -1,6 +1,6 @@
 Foundations I: Parameter
 ===========================
-A **parameter** is one of the main building blocks in equadratures. Let :math:`s` be a parameter defined on a domain :math:`\mathcal{D} \in \mathbb{R}`. The support of the domain :math:`\mathcal{D}` may be:
+A **parameter** is one of the main building blocks in :code:`equadratures`. Let :math:`s` be a parameter defined on a domain :math:`\mathcal{D} \in \mathbb{R}`. The support of the domain :math:`\mathcal{D}` may be:
 
 * closed :math:`[a,b]`;
 * semi-infinite :math:`(-\infty, b)` or :math:`[a, \infty)`;
@@ -12,6 +12,7 @@ Further, let us assume that this parameter is characterized by a positive weight
 
 	\int_{\mathcal{D}}\rho\left(s\right)ds=1.
 
+
 We now demonstrate some basic functionality of this parameter. First consider the case where :math:`\rho(s) = \mathcal{N} (0, 1)` is a standard Gaussian distribution with a mean of 0.0 and a variance of 1.0. We then plot its PDF and cumulative density function (CDF) and demonstrate how we can generate random samples from this distribution.
 
 .. code::

source/index.txt

@@ -2,7 +2,11 @@
 equadratures
 =============
 
-:code:`equadratures` is an open-source python code tailored for tackling problems in **uncertainty quantification**, **surrogate-based optimisation**, **numerical integration**, and **data-driven dimension reduction**.
+.. meta::
+   :description: equadratures is an open-source python code for uncertainty quantification, surrogate-based optimisation, numerical integration and data-driven dimension reduction. It is built on orthogonal polynomial approximations.
+   :keywords: polynomials, polynomial chaos, optimisation, orthogonal polynomials, dimension reduction
+
+:code:`equadratures` is an open-source python code tailored for **uncertainty quantification**, **surrogate-based optimisation**, **numerical integration**, and **data-driven dimension reduction**.
 
 It does this by constructing orthogonal polynomials. It is particularly useful for problems where output quantities of interest are smooth and continuous and to this extent it has found widespread applications in computational engineering models.