step0-cython-step0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTES\n",
    "=====\n",
    "\n",
    "Uses the cython based on the unmodified source \n",
    "(The optimization done by numpy cannot be exploited by cython)\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",
    "import math\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",
    "\n",
    "k = 0.8  # 1/s\n",
    "v = 0.325  # v= v0 + v1*beta\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",
    "kf = 1.0\n",
    "\n",
    "sampleFreq = 1.0 # s\n",
    "duration = 5\n",
    "\n",
    "h = 0.02\n",
    "N = int(duration/h) + 1\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",
    "\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 f( X,  Y):\n",
    "    return VM2*X/(K2+X) - VM3*Y*X*X/((KR+Y)*(KA*KA+X*X)) \\\n",
    "            - kf*Y\n",
    "\n",
    "cdef dXdt( X,  Y,  flux):\n",
    "    return v - f(X,Y) - k*X + flux\n",
    "\n",
    "cdef dYdt( X, Y):\n",
    "    return f(X,Y)\n",
    "\n",
    "cpdef timestep(X, Y, cylinderConcPerCell, index, \n",
    "                ring,  t,  tau,  omega,  Am):\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 + math.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
}