GenerateData_ParallelVersion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#############################################################################################\n",
    "### This file is used to generate Data needed for \"PDE example.py\" or \"PDE example.ipynb\" ###\n",
    "#############################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-19T20:02:49.208911Z",
     "start_time": "2019-11-19T20:02:49.122284Z"
    }
   },
   "outputs": [],
   "source": [
    "import ipyparallel as ipp\n",
    "rc = ipp.Client()\n",
    "view = rc[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-19T20:02:49.962775Z",
     "start_time": "2019-11-19T20:02:49.953412Z"
    }
   },
   "outputs": [],
   "source": [
    "#Setup local environment libraries and also the list of processors on the cluster that we can access \n",
    "import os\n",
    "import numpy as np\n",
    "import dolfin as fn\n",
    "import scipy.io as sio #for the i/o\n",
    "import time #to wait between checking if jobs are done\n",
    "import numpy.polynomial.hermite_e as H \n",
    "from math import factorial\n",
    "from scipy.stats import gaussian_kde as kde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### To make it cleaner, create Directory \"Data\" to store all the data ####\n",
    "datapath = os.path.join(os.getcwd(),\"Data\")\n",
    "os.makedirs(datapath,exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-19T20:02:51.853114Z",
     "start_time": "2019-11-19T20:02:51.845109Z"
    }
   },
   "outputs": [],
   "source": [
    "procs = np.array(rc.ids) #get list of processors on cluster that local environment can access\n",
    "print(procs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-19T20:02:55.205385Z",
     "start_time": "2019-11-19T20:02:55.062102Z"
    }
   },
   "outputs": [],
   "source": [
    "%%px\n",
    "#Setup environment on each processor in cluster\n",
    "\n",
    "import numpy as np\n",
    "import dolfin as fn\n",
    "import scipy.io as sio\n",
    "import numpy.polynomial.hermite_e as H \n",
    "from scipy.stats import gaussian_kde as kde\n",
    "from scipy.stats import norm\n",
    "from math import factorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-19T20:02:56.654220Z",
     "start_time": "2019-11-19T20:02:56.466157Z"
    }
   },
   "outputs": [],
   "source": [
    "%%px\n",
    "\n",
    "def Compute_Q(proc_num, proc_max, mu1=0, mu2=0, sigma1=0.1, sigma2=0.1, gridx=50, gridy=50, p=1):\n",
    "    num_quad = 20\n",
    "    lambda1 = H.hermegauss(num_quad)[0]\n",
    "    lambda2 = H.hermegauss(num_quad)[0]\n",
    "    \n",
    "    # Create the characteristic function class used to define the QoI\n",
    "    class AvgCharFunc(fn.UserExpression):\n",
    "        def __init__(self, region, **kwargs):\n",
    "            self.a = region[0]\n",
    "            self.b = region[1]\n",
    "            self.c = region[2]\n",
    "            self.d = region[3]\n",
    "            super().__init__(**kwargs)\n",
    "        def eval(self, v, x):\n",
    "            v[0] = 0\n",
    "            if (x[0] >= self.a) & (x[0] <= self.b) & (x[1] >= self.c) & (x[1] <= self.d):\n",
    "                v[0] = 1./( (self.b-self.a) * (self.d-self.c) )\n",
    "            return v\n",
    "        def value_shape(self):\n",
    "            return ()\n",
    "        \n",
    "    def QoI_FEM(lam1,lam2,pointa,pointb,gridx,gridy,p):\n",
    "        aa = pointa[0]\n",
    "        bb = pointb[0]\n",
    "        cc = pointa[1]\n",
    "        dd = pointb[1]\n",
    "\n",
    "        mesh = fn.UnitSquareMesh(gridx, gridy)\n",
    "        V = fn.FunctionSpace(mesh, \"Lagrange\", p)\n",
    "\n",
    "        # Define diffusion tensor (here, just a scalar function) and parameters\n",
    "        A = fn.Expression((('exp(lam1)','a'),\n",
    "                    ('a','exp(lam2)')), a = fn.Constant(0.0), lam1 = lam1, lam2 = lam2, degree=3) \n",
    "\n",
    "        u_exact = fn.Expression(\"sin(lam1*pi*x[0])*cos(lam2*pi*x[1])\", lam1 = lam1, lam2 = lam2, degree=2+p)\n",
    "\n",
    "        # Define the mix of Neumann and Dirichlet BCs\n",
    "        class LeftBoundary(fn.SubDomain):\n",
    "            def inside(self, x, on_boundary):\n",
    "                return (x[0] < fn.DOLFIN_EPS)\n",
    "        class RightBoundary(fn.SubDomain):\n",
    "            def inside(self, x, on_boundary):\n",
    "                return (x[0] > 1.0 - fn.DOLFIN_EPS)\n",
    "        class TopBoundary(fn.SubDomain):\n",
    "            def inside(self, x, on_boundary):\n",
    "                return (x[1] > 1.0 - fn.DOLFIN_EPS)\n",
    "        class BottomBoundary(fn.SubDomain):\n",
    "            def inside(self, x, on_boundary):\n",
    "                return (x[1] < fn.DOLFIN_EPS)\n",
    "\n",
    "        # Create a mesh function (mf) assigning an unsigned integer ('uint')\n",
    "        # to each edge (which is a \"Facet\" in 2D)\n",
    "        mf = fn.MeshFunction('size_t', mesh, 1)\n",
    "        mf.set_all(0) # initialize the function to be zero\n",
    "        # Setup the boundary classes that use Neumann boundary conditions\n",
    "        NTB = TopBoundary() # instatiate\n",
    "        NTB.mark(mf, 1) # set all values of the mf to be 1 on this boundary\n",
    "        NBB = BottomBoundary()\n",
    "        NBB.mark(mf, 2) # set all values of the mf to be 2 on this boundary\n",
    "        NRB = RightBoundary()\n",
    "        NRB.mark(mf, 3)\n",
    "\n",
    "        # Define Dirichlet boundary conditions\n",
    "        Gamma_0 = fn.DirichletBC(V, u_exact, LeftBoundary())\n",
    "        bcs = [Gamma_0]\n",
    "\n",
    "        # Define data necessary to approximate exact solution\n",
    "        f = ( fn.exp(lam1)*(lam1*fn.pi)**2 + fn.exp(lam2)*(lam2*fn.pi)**2 ) * u_exact\n",
    "        #g1:#pointing outward unit normal vector, pointing upaward (0,1)\n",
    "        g1 = fn.Expression(\"-exp(lam2)*lam2*pi*sin(lam1*pi*x[0])*sin(lam2*pi*x[1])\", lam1=lam1, lam2=lam2, degree=2+p)    #pointing outward unit normal vector, pointing upaward (0,1)\n",
    "        #g2:pointing downward (0,1)\n",
    "        g2 = fn.Expression(\"exp(lam2)*lam2*pi*sin(lam1*pi*x[0])*sin(lam2*pi*x[1])\", lam1=lam1, lam2=lam2, degree=2+p)     #pointing downward (0,1)\n",
    "        g3 = fn.Expression(\"exp(lam1)*lam1*pi*cos(lam1*pi*x[0])*cos(lam2*pi*x[1])\", lam1=lam1, lam2=lam2, degree=2+p)\n",
    "\n",
    "        fn.ds = fn.ds(subdomain_data=mf)\n",
    "        # Define variational problem\n",
    "        u = fn.TrialFunction(V)\n",
    "        v = fn.TestFunction(V)\n",
    "        a = fn.inner(A*fn.grad(u), fn.grad(v))*fn.dx\n",
    "        L = f*v*fn.dx + g1*v*fn.ds(1) + g2*v*fn.ds(2) + g3*v*fn.ds(3)  #note the 1, 2 and 3 correspond to the mf\n",
    "\n",
    "        # Compute solution\n",
    "        u = fn.Function(V)\n",
    "        fn.solve(a == L, u, bcs)\n",
    "        psi = AvgCharFunc([aa, bb, cc, dd], degree=0)\n",
    "        Q = fn.assemble(fn.project(psi * u, V) * fn.dx)\n",
    "        \n",
    "        return Q\n",
    "\n",
    "    Q_FEM = np.zeros(400)\n",
    "    num_Q_per_proc = 400//proc_max\n",
    "    if proc_num != procs.size -1:       \n",
    "        for i in range(proc_num*num_Q_per_proc, (proc_num+1)*num_Q_per_proc):\n",
    "            Q_FEM[i] = QoI_FEM(mu1+sigma1*lambda1[i%num_quad],mu2+sigma2*lambda2[i//num_quad],[0.4,0.4],[0.6,0.6],gridx,gridy,p)\n",
    "    else:\n",
    "        for i in range(proc_num*num_Q_per_proc,400):\n",
    "            Q_FEM[i] = QoI_FEM(mu1+sigma1*lambda1[i%num_quad],mu2+sigma2*lambda2[i//num_quad],[0.4,0.4],[0.6,0.6],gridx,gridy,p)\n",
    "    filename = os.path.join(os.getcwd(), \"Data\", \"Q_FEM_quad_\") + str(proc_num) + \".mat\"\n",
    "    data_dict = {'Q_FEM': Q_FEM}\n",
    "    sio.savemat(filename, data_dict)\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-19T20:03:58.552528Z",
     "start_time": "2019-11-19T20:02:57.445014Z"
    }
   },
   "outputs": [],
   "source": [
    "#########################################################\n",
    "##### Generate datafiles Data/Q_FEM_quad_[0-24].mat #####\n",
    "#########################################################\n",
    "\n",
    "view.push(dict(procs=procs)) #This pushes procs to all processors on the cluster\n",
    "\n",
    "async_process = view.map_async(lambda proc_num: Compute_Q(proc_num, proc_max=procs.size), range(procs.size) ) \n",
    "\n",
    "time.sleep(1) #give the process time to start and see if any errors occur\n",
    "if async_process.error[0] is None:\n",
    "    done = False    \n",
    "    while done == False:\n",
    "        if async_process.done():\n",
    "            done = True\n",
    "            #print('Now we can load in the data')\n",
    "        else:\n",
    "            time.sleep(1)\n",
    "else:\n",
    "    raise RuntimeError(async_process.error[0])    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-19T20:04:49.871851Z",
     "start_time": "2019-11-19T20:04:49.047837Z"
    }
   },
   "outputs": [],
   "source": [
    "%%px\n",
    "def r(nn, proc_num, proc_max):\n",
    "    mu1 = 0\n",
    "    mu2 = 0\n",
    "    sigma1 = 0.1\n",
    "    sigma2 = 0.1\n",
    "    N_size = int(1E4)\n",
    "    np.random.seed(123456)\n",
    "    lam1 = np.random.normal(mu1,sigma1,N_size)\n",
    "    lam2 = np.random.normal(mu2,sigma2,N_size)\n",
    "    \n",
    "    def Hermite_2d(i,j,x,y):\n",
    "        c = np.zeros((20,20))\n",
    "        c[i,j] = 1\n",
    "        return H.hermeval2d(x, y, c)\n",
    "\n",
    "    Q_FEM_quad = np.zeros(int(400))   #already include information of mu1, mu2, sigma1, sigma2\n",
    "    for i in range(procs.size):         \n",
    "        filename = os.path.join(os.getcwd(), \"Data\", \"Q_FEM_quad_\") + str(i) + '.mat' \n",
    "        partial_data = sio.loadmat(filename)\n",
    "        Q_FEM_quad += partial_data['Q_FEM'].reshape(int(400))\n",
    "\n",
    "    def Phi(n):\n",
    "        #define H_n\n",
    "        coeffs = [0]*(n+1)\n",
    "        coeffs[n] = 1\n",
    "        return coeffs\n",
    "\n",
    "    def q(i,j):\n",
    "        x, w=H.hermegauss(20)     \n",
    "        Q=sum([w[ldx]*sum([w[kdx] * Q_FEM_quad[ldx*20+kdx] * H.hermeval(x[kdx],Phi(i)) for kdx in range(20)])*H.hermeval(x[ldx],Phi(j)) for ldx in range(20)])     \n",
    "        q= Q/(2*np.pi*factorial(i)*factorial(j))\n",
    "\n",
    "        return q\n",
    "\n",
    "    qij = np.zeros((10,10))\n",
    "    for i in range(10):\n",
    "        for j in range(10):\n",
    "            qij[i,j] = q(i,j)\n",
    "\n",
    "    def Q(n,x,y):\n",
    "        result = 0 \n",
    "        for i in range(n+1):\n",
    "            for j in range(n+1):\n",
    "                if i+j <=n:\n",
    "                    result += qij[i,j]*Hermite_2d(i,j,(x-mu1)/sigma1,(y-mu2)/sigma2)\n",
    "        return result\n",
    "\n",
    "    def Qexact(x,y,a=0.4,b=0.6,c=0.4,d=0.6):\n",
    "        sol = (np.cos(x*np.pi*a)-np.cos(x*np.pi*b))*(np.sin(y*np.pi*d)-np.sin(y*np.pi*c))/((b-a)*(d-c)*x*y*np.pi**2)\n",
    "        return sol\n",
    "    qexact = Qexact(lam1,lam2)\n",
    "    pfprior_dens = kde(qexact)\n",
    "\n",
    "    def pfprior_dens_n(n,x):\n",
    "        pfprior_sample_n =  Q(n,lam1,lam2)\n",
    "        pdf = kde(pfprior_sample_n)\n",
    "        return pdf(x)\n",
    "     \n",
    "    def pdf_obs(x):\n",
    "        return norm.pdf(x, loc=0.3, scale=0.1)\n",
    "\n",
    "    r = np.zeros(int(1E4))    \n",
    "    num_r_per_proc = int(1E4)//proc_max\n",
    "    if proc_num != procs.size -1: \n",
    "        for i in range(proc_num*num_r_per_proc, (proc_num+1)*num_r_per_proc):\n",
    "            if nn == 0:\n",
    "                r[i] = pdf_obs(qexact[i])/pfprior_dens(qexact[i])\n",
    "            else:\n",
    "                q = Q(nn,lam1[i],lam2[i])\n",
    "                r[i] = pdf_obs(q)/pfprior_dens_n(nn,q)\n",
    "    elif proc_num == procs.size-1:\n",
    "        for i in range(proc_num*num_r_per_proc,int(1E4)):\n",
    "            if nn == 0:\n",
    "                r[i] = pdf_obs(qexact[i])/pfprior_dens(qexact[i])\n",
    "            else:\n",
    "                q = Q(nn,lam1[i],lam2[i])\n",
    "                r[i] = pdf_obs(q)/pfprior_dens_n(nn,q)\n",
    "    if nn == 0:\n",
    "        filename = os.path.join(os.getcwd(),\"Data\",\"r_\") + str(proc_num) + '.mat'\n",
    "    else:\n",
    "        filename = os.path.join(os.getcwd(),\"Data\",\"r\") + str(nn) + '_' + str(proc_num) + '.mat' \n",
    "    data_dict = {'r': r}\n",
    "    sio.savemat(filename, data_dict)\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-19T20:29:17.267346Z",
     "start_time": "2019-11-19T20:04:58.001776Z"
    }
   },
   "outputs": [],
   "source": [
    "#########################################################\n",
    "######### Generate datafiles Data/r_[0-24].mat ##########\n",
    "######### Generate datafiles Data/r1_[0-24].mat #########\n",
    "######### Generate datafiles Data/r2_[0-24].mat #########\n",
    "######### Generate datafiles Data/r3_[0-24].mat #########\n",
    "######### Generate datafiles Data/r4_[0-24].mat #########\n",
    "######### Generate datafiles Data/r5_[0-24].mat #########\n",
    "######### Generate datafiles Data/r6_[0-24].mat #########\n",
    "#########################################################\n",
    "nn = [0,1,2,3,4,5,6]\n",
    "\n",
    "view.push(dict(procs=procs, nn=nn)) #This pushes procs to all processors on the cluster\n",
    "for j in range(len(nn)):\n",
    "    view.push(dict(j=j))\n",
    "    async_process = view.map_async(lambda proc_num: r(nn[j], proc_num, proc_max=procs.size), range(procs.size) ) \n",
    "\n",
    "    time.sleep(1) #give the process time to start and see if any errors occur\n",
    "    if async_process.error[0] is None:\n",
    "        done = False    \n",
    "        while done == False:\n",
    "            if async_process.done():\n",
    "                done = True\n",
    "                #print('Now we can load in the data')\n",
    "            else:\n",
    "                time.sleep(1)\n",
    "    else:\n",
    "        raise RuntimeError(async_process.error[0])"
   ]
  }
 ],
 "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.8.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}