LU+Decomposition.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Imports \n",
    "import pprint\n",
    "import scipy\n",
    "import scipy.linalg   # SciPy Linear Algebra Library\n",
    "import numpy as np\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 25. ,  50. ,  -1.3],\n",
       "       [ 75. ,  80. ,   2.5],\n",
       "       [ 10. ,  25. ,   3. ],\n",
       "       [ 95. ,  15. ,  -2.7]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Assign Data Values - assume 4 data locations\n",
    "ndata = 4\n",
    "data = np.zeros((ndata,3))  # x, y, value\n",
    "data[0,0] = 25.0;  data[0,1] = 50.0;  data[0,2] = -1.3\n",
    "data[1,0] = 75.0;  data[1,1] = 80.0;  data[1,2] = 2.5\n",
    "data[2,0] = 10.0;  data[2,1] = 25.0;  data[2,2] = 3.0\n",
    "data[3,0] = 95.0;  data[3,1] = 15.0;  data[3,2] = -2.7\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.2244834 ,  0.57506938,  0.06574285],\n",
       "       [ 0.2244834 ,  1.        ,  0.03145365,  0.13715671],\n",
       "       [ 0.57506938,  0.03145365,  1.        ,  0.0296663 ],\n",
       "       [ 0.06574285,  0.13715671,  0.0296663 ,  1.        ]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate Symmetric Covariance Array - assuming variogram with spherical structure with range specified\n",
    "cov = np.zeros((ndata,ndata))\n",
    "var_range = 100.0\n",
    "for i in range(0, ndata):\n",
    "    for j in range(0, ndata):\n",
    "        distance = math.sqrt(math.pow((data[i,0]-data[j,0]),2) + math.pow((data[i,1]-data[j,1]),2))\n",
    "        cova = 0.0\n",
    "        if(distance < var_range):\n",
    "            hr = distance / var_range\n",
    "            cova = 1.0 - hr*(1.5 - 0.5*hr*hr)  # spherical structure, no nugget\n",
    "            cov[i,j] = cova\n",
    "cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.          0.          0.          0.        ]\n",
      " [ 0.2244834   1.          0.          0.        ]\n",
      " [ 0.57506938 -0.10282133  1.          0.        ]\n",
      " [ 0.06574285  0.12889386  0.00674212  1.        ]]\n",
      "[[ 1.          0.2244834   0.57506938  0.06574285]\n",
      " [ 0.          0.9496072  -0.09763988  0.12239854]\n",
      " [ 0.          0.          0.65925575  0.00444478]\n",
      " [ 0.          0.          0.          0.97987149]]\n"
     ]
    }
   ],
   "source": [
    "# LU Decomposition using scipy (used tutorial at www.quantstart.com) \n",
    "P, L, U = scipy.linalg.lu(cov)\n",
    "print(L);  print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "          0.00000000e+00],\n",
       "       [  0.00000000e+00,   0.00000000e+00,   1.38777878e-17,\n",
       "          0.00000000e+00],\n",
       "       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "          0.00000000e+00],\n",
       "       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,\n",
       "          0.00000000e+00]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Test the LU Decomposition\n",
    "Test = cov - np.matmul(L,U)\n",
    "Test  # should be zeros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.84674634  0.90973758  2.40580429  1.60080233]\n"
     ]
    }
   ],
   "source": [
    "# Unconditional Realization at the Specified Locations\n",
    "rand = np.zeros((ndata))\n",
    "for i in range(0, ndata):\n",
    "    rand[i] = np.random.normal()\n",
    "realization = np.matmul(L,rand)\n",
    "print(realization)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}