step0-cython.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTES\n",
    "=====\n",
    "\n",
    "Uses the cython heavily, based on the unmodified source \n",
    "(The optimization done by numpy cannot be exploited by cython)\n",
    "\n",
    "The complete simulation has been moved to cython and optimized.\n",
    "\n",
    "typed memory views\n",
    "------------\n",
    "The X, Y, and ConcPerCell Arrays have been used as a typed memory view.\n",
    "\n",
    "type annotation\n",
    "--------\n",
    "Every single type annotation has been introduced, otherwise no changes\n",
    "\n",
    "Boundschecking off\n",
    "-----\n",
    "Bounds checking off gives a speedup\n",
    "\n",
    "Installing Cython\n",
    "-----\n",
    "pip3 install Cython\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from  matplotlib import pyplot \n",
    "%reload_ext Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "cimport cython\n",
    "\n",
    "from libc.math cimport cos\n",
    "\n",
    "cdef double VM2 = 20.0 # 10^-6 M/s\n",
    "cdef double VM3 = 23.0 # 10^-6 M/s\n",
    "cdef double K2 = 1.0 # 10^-6 M/s\n",
    "cdef double KR = 0.8 # 10^-6 M/s\n",
    "cdef double KA = 0.9 # 10^-6 M/s\n",
    "cdef double kf = 1.0  # 1/s\n",
    "\n",
    "cdef double k = 0.8  # 1/s\n",
    "cdef double v = 0.325  # v= v0 + v1*beta\n",
    "\n",
    "cdef double Dxl = 2.0  # 1/s\n",
    "cdef double Dxr = 6.0  # 1/s\n",
    "cdef double Dxd = 1.0  # 1/s\n",
    "cdef double kt = 0.1 # 1/s\n",
    "\n",
    "cdef double sampleFreq = 1.0 # s\n",
    "cdef int duration = 5\n",
    "\n",
    "cdef double h = 0.02\n",
    "cdef int N = int(duration/h) + 1\n",
    "#######################################\n",
    "# Initial conditions\n",
    "#######################################\n",
    "\n",
    "cdef int ringsize=50\n",
    "cdef int rings=200\n",
    "#fileX = fromfile(\"X0.dat\", sep=';')\n",
    "#fileY = fromfile(\"Y0.dat\", sep=';')\n",
    "\n",
    "\n",
    "\n",
    "#X0 = np.zeros((ringsize,rings), double)\n",
    "#Y0 = np.zeros((ringsize,rings), double)\n",
    "#Amp = 0.55\n",
    "#index=0\n",
    "#for ring in range(0,3):\n",
    "#    for cell in range(24,27):\n",
    "#        X0[cell][ring] = Amp\n",
    "#        Y0[cell][ring] = Amp \n",
    "#        index += 1          \n",
    "  \n",
    "#for ring in range(rings):\n",
    "#    for cell in range(ringsize):\n",
    "#        X0[cell][ring] = fileX[index]\n",
    "#        Y0[cell][ring] = fileY[index] \n",
    "#        index += 1  \n",
    "\n",
    "\n",
    "cdef double f(double X, double Y):\n",
    "    return VM2*X/(K2+X) - VM3*Y*X*X/((KR+Y)*(KA*KA+X*X)) \\\n",
    "            - kf*Y\n",
    "\n",
    "cdef double dXdt(double X, double Y, double flux):\n",
    "    return v - f(X,Y) - k*X + flux\n",
    "\n",
    "cdef double dYdt(double X,double Y):\n",
    "    return f(X,Y)\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "cpdef timestep(double[:,:,:] X, double[:,:,:] Y, double[:,:] cylinderConcPerCell, int index, \n",
    "               int ring, int t, int tau, double omega, double Am):\n",
    "    cdef double flux\n",
    "    cdef double x, y\n",
    "    cdef double k11, k12, k21, k22, k31, k32, k41, k42\n",
    "    cdef int cell\n",
    "    for cell in range(0,ringsize): \n",
    "            flux = 0.0\n",
    "# Sinus-Exitation\n",
    "            if 0 <= ring <= 3:\n",
    "                if 24 <= cell <= 27:\n",
    "                    X[index][ring][cell] = Am*0.5*(1 + cos(omega*t*h))\n",
    "#                   Y[index][ring][cell] = Amp*(math.cos(omega*t*h))**2\n",
    "# horizontal diffusion\n",
    "            if cell == 0:\n",
    "                flux = Dxr*(X[index][ring][ringsize-1] - X[index][ring][cell]) \\\n",
    "                   + Dxl*(X[index][ring][1] - X[index][ring][cell])\n",
    "            elif cell==ringsize-1:\n",
    "                flux = Dxr*(X[index][ring][cell-1] - X[index][ring][cell]) \\\n",
    "                       + Dxl*(X[index][ring][0] - X[index][ring][cell])\n",
    "            else :\n",
    "                flux = Dxr*(X[index][ring][cell-1] - X[index][ring][cell]) \\\n",
    "                       + Dxl*(X[index][ring][cell+1] - X[index][ring][cell])\n",
    "# vertical diffusion\n",
    "            if ring == 0:\n",
    "                if X[index][ring][cell] > X[index][ring+1][cell]:\n",
    "                    flux -= Dxd*(X[index][ring][cell] - X[index][ring+1][cell])\n",
    "            elif ring==rings-1:\n",
    "                if X[index][ring-1][cell] > X[index][ring][cell]:\n",
    "                    flux += Dxd*(X[index][ring-1][cell] - X[index][ring][cell])\n",
    "            else:\n",
    "                if X[index][ring-1][cell] > X[index][ring][cell]:\n",
    "                    flux += Dxd*(X[index][ring-1][cell] - X[index][ring][cell])\n",
    "                if X[index][ring][cell] > X[index][ring+1][cell]:\n",
    "                    flux -= Dxd*(X[index][ring][cell] - X[index][ring+1][cell])\n",
    "\n",
    "# Feedback across ring, tau elements back in time\n",
    "            if t>=tau:\n",
    "                flux += kt*(X[(t-tau)%tau][ring][(cell + ringsize//2)%ringsize]\\\n",
    "                              - X[index][ring][cell])\n",
    "\n",
    "\n",
    "# Perform Runge-Kutta integration on current cell\n",
    "            x = X[index][ring][cell]\n",
    "            y = Y[index][ring][cell]\n",
    "            k11 = h*dXdt(x,y,flux)\n",
    "            k12 = h*dYdt(x,y)\n",
    "\n",
    "            k21 = h*dXdt(x+0.5*k11, y+0.5*k12,flux)\n",
    "            k22 = h*dYdt(x+0.5*k11, y+0.5*k12)\n",
    "\n",
    "            k31 = h*dXdt(x+0.5*k21, y+0.5*k22,flux)\n",
    "            k32 = h*dYdt(x+0.5*k21, y+0.5*k22)\n",
    "\n",
    "            k41 = h*dXdt(x+k31, y+k32,flux)\n",
    "            k42 = h*dYdt(x+k31, y+k32)\n",
    "\n",
    "            Xret = x + (k11 + 2*k21 + 2*k31 + k41)/6\n",
    "            Yret = y + (k12 + 2*k22 + 2*k32 + k42)/6\n",
    "            X[(index+1)%tau][ring][cell] = Xret\n",
    "            Y[(index+1)%tau][ring][cell] = Yret\n",
    "            cylinderConcPerCell[cell][t] += X[index][ring][cell]\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original Code\n",
    "#######################################\n",
    "#######################################\n",
    "## 50x50Nodes.py\n",
    "##\n",
    "## Runge-Kutta simulation of 50 cells/ring x 50 rings with feedback\n",
    "## Create image contour plots of full X matrix at given intervals\n",
    "## Calculate sum through column and difference in sum across ring\n",
    "##\n",
    "##Store minimal amount of data in RAM\n",
    "#######################################\n",
    "#######################################\n",
    "\n",
    "\n",
    "#######################################\n",
    "# Library imports\n",
    "#######################################\n",
    "\n",
    "\n",
    "import time as mytime\n",
    "import math\n",
    "import numpy as np\n",
    "\n",
    "np.random.seed(1234)\n",
    "#######################################\n",
    "# Constants\n",
    "#######################################\n",
    "\n",
    "VM2 = 20.0 # 10^-6 M/s\n",
    "VM3 = 23.0 # 10^-6 M/s\n",
    "K2 = 1.0 # 10^-6 M/s\n",
    "KR = 0.8 # 10^-6 M/s\n",
    "KA = 0.9 # 10^-6 M/s\n",
    "kf = 1.0  # 1/s\n",
    "\n",
    "k = 0.8  # 1/s\n",
    "v = 0.325  # v= v0 + v1*beta\n",
    "\n",
    "tau = 16 # s\n",
    "\n",
    "Dxl = 2.0  # 1/s\n",
    "Dxr = 6.0  # 1/s\n",
    "Dxd = 1.0  # 1/s\n",
    "kt = 0.1 # 1/s\n",
    "\n",
    "sampleFreq = 1.0 # s\n",
    "duration = 5\n",
    "\n",
    "\n",
    "#######################################\n",
    "# Initial conditions\n",
    "#######################################\n",
    "\n",
    "ringsize=50\n",
    "rings=200\n",
    "#fileX = fromfile(\"X0.dat\", sep=';')\n",
    "#fileY = fromfile(\"Y0.dat\", sep=';')\n",
    "omega = 1.047\n",
    "Am =  1.2\n",
    "\n",
    "#X0 = np.zeros((ringsize,rings), double)\n",
    "#Y0 = np.zeros((ringsize,rings), double)\n",
    "#Amp = 0.55\n",
    "#index=0\n",
    "#for ring in range(0,3):\n",
    "#    for cell in range(24,27):\n",
    "#        X0[cell][ring] = Amp\n",
    "#        Y0[cell][ring] = Amp \n",
    "#        index += 1          \n",
    "  \n",
    "#for ring in range(rings):\n",
    "#    for cell in range(ringsize):\n",
    "#        X0[cell][ring] = fileX[index]\n",
    "#        Y0[cell][ring] = fileY[index] \n",
    "#        index += 1  \n",
    "        \n",
    "X0 = (np.random.rand(ringsize,rings)-0.5)*0.1 + 0.406\n",
    "Y0 = (np.random.rand(ringsize,rings)-0.5)*0.2 + 2.76\n",
    "\n",
    "#######################################\n",
    "# Functions\n",
    "#######################################\n",
    "\n",
    "# Log data to file\n",
    "def writeRing(filename,data,ring):\n",
    "    f= open (filename + str(ring) + '.dat', 'w')\n",
    "    for t in range(0,N):\n",
    "        string = \"\"\n",
    "        for cell in range(0,ringsize-1):\n",
    "            string += str(data[t][ring][cell]) + ';'\n",
    "        string += str(data[t][ring][ringsize-1]) + '\\n'\n",
    "        f.write(string) \n",
    "\n",
    "def writeCylinder(filename,data,ring):\n",
    "    f= open (filename + str(ring) + '.dat', 'w')\n",
    "    for ring in range(0,rings):\n",
    "        string = \"\"\n",
    "        for cell in range(0,ringsize-1):\n",
    "            string += str(data[time%tau][ring][cell]) + ';'\n",
    "        string += str(data[time%tau][ring][ringsize-1]) + '\\n'\n",
    "        f.write(string) \n",
    "        \n",
    "        \n",
    "#######################################\n",
    "# Perform simulation\n",
    "#######################################\n",
    "\n",
    "# Set step size \n",
    "h = 0.02\n",
    "N = int(duration/h) + 1\n",
    "\n",
    "sampleFreq = sampleFreq/h\n",
    "tau = int(tau/h)\n",
    "\n",
    "cylinderConcPerCell = np.zeros((ringsize,N),float)\n",
    "totalConcVector = np.zeros((2,N),float)\n",
    "angle = 2*math.pi/ringsize\n",
    "\n",
    "# initialize X and Y matrices\n",
    "# X[time][ring][cell]\n",
    "X = np.zeros((tau,rings, ringsize),float)\n",
    "Y = np.zeros((tau,rings, ringsize),float)\n",
    "for ring in range(rings):\n",
    "    for cell in range(ringsize):\n",
    "        X[0][ring][cell] = X0[cell][ring]\n",
    "        Y[0][ring][cell] = Y0[cell][ring]\n",
    "\n",
    "time = np.arange(0,h*N+h,h)\n",
    "\n",
    "# Loop over all cells and time points\n",
    "startTime = mytime.time()\n",
    "for t in range(0,N):\n",
    "    index = t%tau\n",
    "    for ring in range(0,rings):\n",
    "        timestep(X, Y, cylinderConcPerCell, index, ring, t, tau, omega, Am)\n",
    "            \n",
    "# Generate progress output to standard out        \n",
    "    if t%5==0 and t!=0:\n",
    "        elapsed = mytime.time()-startTime\n",
    "        ratio = 1.0*t/N\n",
    "        print ('%d/%d, \\t %.1f%%, \\t elapsed: %.1f \\t ETA: %.1f'\\\n",
    "               % (t, N, ratio*100, elapsed, elapsed/ratio-elapsed))\n",
    "   \n",
    "# Write data output to files\n",
    "    if t%sampleFreq==0:\n",
    "# Write raw data to files\n",
    "#        writeCylinder('outputX_F_Dxr_' + str(Dxr) + '_k_' + ':', X, t)\n",
    "#        writeCylinder('outputY_F_Dxr_' + str(Dxr) + '_k_' + ':', Y, t)\n",
    "\n",
    "# Write X matrix to image file\n",
    "        fig = pyplot.figure(figsize=(5,16))\n",
    "        pyplot.clf()\n",
    "        fig = pyplot.contourf(X[(index+1)%tau,:,:])\n",
    "        pyplot.xlabel('Cell')\n",
    "        pyplot.ylabel('Ring')\n",
    "        pyplot.colorbar()\n",
    "        pyplot.title('Dxr = ' + str(Dxr) + ', Dxl = ' + str(Dxl) + ', time = ' + str(t*h) + ' s')\n",
    "        \n",
    "        \n",
    "#######################################\n",
    "# Post processing\n",
    "#######################################\n",
    "\n",
    "# calculate total concentration totalConcVector \n",
    "    x = 0\n",
    "    y = 0\n",
    "    for cell in range(0,ringsize):\n",
    "        x += cylinderConcPerCell[cell][t]*math.cos(angle*cell)\n",
    "        y += cylinderConcPerCell[cell][t]*math.sin(angle*cell)\n",
    "        \n",
    "# Angle\n",
    "    totalConcVector[0][t] = math.atan(y/x)\n",
    "    if totalConcVector[0][t]<0:\n",
    "        totalConcVector[0][t] += 2*math.pi\n",
    "    totalConcVector[0][t] = totalConcVector[0][t]/2/math.pi*ringsize\n",
    "# Length\n",
    "    totalConcVector[1][t] = math.sqrt(x*x+y*y)\n",
    "    \n",
    "    \n",
    "    \n",
    "# Calculate total diff concentration vector across cylinderConcPerCell\n",
    "totalDiffConcVector = np.zeros((ringsize,N),float)\n",
    "for t in range(N):\n",
    "    for cell in range(ringsize):\n",
    "        totalDiffConcVector[cell][t] = cylinderConcPerCell[cell][t] - cylinderConcPerCell[(cell + ringsize//2)%ringsize][t]\n",
    "        \n",
    "        \n",
    "        \n",
    "# Plot Cylinder Sum\n",
    "fig = pyplot.figure(figsize=(5,16))\n",
    "pyplot.clf()\n",
    "fig = pyplot.contourf(cylinderConcPerCell)\n",
    "pyplot.xlabel('Iterations (' + str(1/h) + '/s)')\n",
    "pyplot.ylabel('Cell')\n",
    "pyplot.colorbar()\n",
    "pyplot.title('Vertical sum through cylinder\\n Dxr = ' + str(Dxr) + '_Dxl = ' + str(Dxl) +\\\n",
    "      ', k='+ str(k))\n",
    "\n",
    "\n",
    "fig = pyplot.figure(figsize=(12,4))\n",
    "pyplot.clf()\n",
    "pyplot.plot(time[1:], totalConcVector[0][:], label='Angle [Cell Number]')\n",
    "pyplot.plot(time[1:], totalConcVector[1][:], label='Length [uM]')\n",
    "pyplot.legend()\n",
    "pyplot.xlabel('Time [s]')\n",
    "pyplot.title('Total Ca++ Concentration vector for cylinder')\n",
    "        \n",
    "        \n",
    "# Plot Diff Cylinder Sum\n",
    "fig = pyplot.figure(figsize=(5,16))\n",
    "pyplot.clf()\n",
    "fig = pyplot.contourf(totalDiffConcVector)\n",
    "pyplot.xlabel('Iterations (' + str(1/h) + '/s)')\n",
    "pyplot.ylabel('Cell')\n",
    "pyplot.colorbar()\n",
    "pyplot.title('Vertical diff sum through cylinder\\n Dxr = ' + str(Dxr) + ', Dxl = ' +  str(Dxl) +\\\n",
    "      ', k='+ str(k))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}