GenData.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "36b6cc32",
   "metadata": {},
   "source": [
    "See the final report appendix or the final project data assignment for more information on the data generation for this project"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6ce70021",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Surrogates,Flux,Statistics,DifferentialEquations,Plots,Suppressor,MLDataUtils,DelimitedFiles\n",
    "using Flux.Data: DataLoader\n",
    "using DelimitedFiles\n",
    "include(\"soluteperm.jl\");\n",
    "include(\"utilityfunc.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3a27e836",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float32}:\n",
       "  5.0f-6\n",
       "  5.0f-6\n",
       " 30.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#set lower and upper bounds on Lpt and Kt (Lpt is first element and Kt is second element)\n",
    "lb=Float32[5e-7,5e-7,5]\n",
    "ub=Float32[5e-6,5e-6,30]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3f3390b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "inv_oscaler (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conc_ub=.6\n",
    "i_lb=Float32[5e-7,7e-7,5.0,0.0,0.0]\n",
    "i_ub=Float32[3.5e-6,4e-6,30,(300/3600),conc_ub]\n",
    "iscaler(x)=(x-i_lb)./(i_ub-i_lb)\n",
    "inv_iscaler(x)=x.*(i_ub.-i_lb).+i_lb\n",
    "oscaler(x)= x ./ conc_ub\n",
    "inv_oscaler(x)=x.*conc_ub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dc038a7a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "iscalerbatch (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function iscalerbatch(x)\n",
    "    nfeat,sl = size(x)\n",
    "    out=Array{Float32,2}(undef,nfeat,sl)\n",
    "    for i=1:sl\n",
    "        out[:,i]=iscaler(x[:,i])\n",
    "    end\n",
    "    return out\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "62018489",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "inv_iscalerbatch (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function inv_iscalerbatch(x)\n",
    "    nfeat,sl = size(x)\n",
    "    out=Array{Float32,2}(undef,nfeat,sl)\n",
    "    for i=1:sl\n",
    "        out[:,i]=inv_iscaler(x[:,i])\n",
    "    end\n",
    "    return out\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "206ac29a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "prepdata (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function Diffusion(rs)\n",
    "    Diam_parc = 2*rs\n",
    "    return D = 1.981e-06*Diam_parc^(-1.157) + 2.221e-08\n",
    "end\n",
    "\n",
    "function Blood_half_life(rs)\n",
    "    Diam_parc = 2*rs\n",
    "    a1 = 1081\n",
    "    b1 = -16.63\n",
    "    c1 = 84.82\n",
    "    a2 = 517.4\n",
    "    b2 = 65.61\n",
    "    c2 = 996.6\n",
    "    return kd = (a1*exp(-((Diam_parc-b1)/c1)^2) + a2*exp(-((Diam_parc-b2)/c2)^2))*60\n",
    "end\n",
    "function a_pressure(r,α)\n",
    "    return r == 0 ? 1 - csch(α) : 1 - sinh(α*r)/(r*sinh(α))\n",
    "    #return 1 - sinh(α*r)/(r*sinh(α))\n",
    "end\n",
    "function a_pressure_prime(r,α)\n",
    "    return r == 0 ? 0 : (sinh(α*r)-α*r*cosh(r*α))/(sinh(α)*r^2)\n",
    "    #return (sinh(α*r)-α*r*cosh(r*α))/(sinh(α)*r^2)\n",
    "end\n",
    "function MSTHighPeclet!(dc,c,p,t)\n",
    "    #dc = preallocated vector corresponding to f(c,t)= dC/dt\n",
    "    #c = concentration vector at all spatial points\n",
    "    #p = vector of parameters as define in modelanalysis\n",
    "    #t = current time\n",
    "    #Unpackage parameters\n",
    "    P,N,sigma,Perm,Lpt,Svt,Kt,Pvv,Pv,D,r,dr,kd=p\n",
    "    co=1 #initial concentration\n",
    "    tspan=3600 #timespan\n",
    "    cv=co*exp(-t*tspan/kd) #time-dependent exponential decay term\n",
    "    Pe=Lpt*Pv*(Pvv-P[1])*(1-sigma)/Perm #Initial Peclet value\n",
    "    #interior boundary condition\n",
    "    dc[1]=tspan*(2*D*(c[2]-c[1])/dr^2 +Lpt*Svt*Pv*(Pvv-P[1])*(1-sigma)*cv)\n",
    "    for j=2:N-1\n",
    "        dc[j]=tspan*(((2*D/r[j])*((c[j+1]-c[j])/dr))+(D*(c[j+1]-2*c[j]+c[j-1])/dr^2)+ Kt*((P[j+1]-P[j-1])/(2*dr))*((c[j+1]-c[j])/dr)+ Lpt*Svt*Pv*(Pvv-P[j])*(1-sigma)*cv)\n",
    "    end\n",
    "    dc[N]=0 # concentration at the outtermost boundary\n",
    "end\n",
    "function MST!(dc,c,p,t)\n",
    "    #dc = preallocated vector corresponding to f(c,t)= dC/dt\n",
    "    #c = concentration vector at all spatial points\n",
    "    #p = vector of parameters as define in modelanalysis\n",
    "    #t = current time\n",
    "    #Unpackage parameters\n",
    "    P,N,sigma,Perm,Lpt,Svt,Kt,Pvv,Pv,D,r,dr,kd=p\n",
    "    co=1 #initial concentration\n",
    "    tspan=3600 #timespan\n",
    "    cv=co*exp(-t*tspan/kd) #time-dependent exponential decay term\n",
    "    Pe=Lpt*Pv*(Pvv-P[1])*(1-sigma)/Perm #Initial Peclet value\n",
    "    dc[1]= tspan*(2*D*(c[2]-c[1])*dr^-2 + Lpt*Svt*Pv*(Pvv-P[1])*(1-sigma)*(cv*exp(Pe)-c[1])/(exp(Pe)-1)) # R=0 boundary condition\n",
    "    for j=2:N-1\n",
    "        #Define interior nodes\n",
    "        Pe=Lpt*Pv*(Pvv-P[j])*(1-sigma)/Perm #Peclet number\n",
    "        dc[j]=tspan*(((2*D/r[j])*((c[j+1]-c[j])/dr))+(D*(c[j+1]-2*c[j]+c[j-1])/dr^2)+ Kt*((P[j+1]-P[j-1])/(2*dr))*((c[j+1]-c[j])/dr)+ Lpt*Svt*Pv*(Pvv-P[j])*(1-sigma)*(cv*exp(Pe)-c[j])/(exp(Pe)-1))\n",
    "    end\n",
    "    dc[N]=0 # concentration at the outtermost boundary\n",
    "end\n",
    "function Accumulation_Model(sol,n_spatial,n_time,n_nodes)\n",
    "    c_model=zeros(n_nodes)\n",
    "    spacingfactor = n_time ÷ n_nodes\n",
    "    for j=1:n_nodes\n",
    "        c_model[j]=mean(sol[:,spacingfactor*j])\n",
    "    end\n",
    "    return c_model\n",
    "end\n",
    "function SepBatch(Lpt,Kt,rs,t,accum_dat,n)\n",
    "    idat = zeros(5,n-1)\n",
    "    odat=zeros(n-1)\n",
    "    for i=1:n-1\n",
    "        idat[:,i]=[Lpt,Kt,rs,t[i],accum_dat[i]]\n",
    "        odat[i]=accum_dat[i+1]\n",
    "    end\n",
    "    return idat,odat\n",
    "end\n",
    "function accumlation_intermediate(x) #x is a vector containing Lpt and Kt\n",
    "    n_spatial=100\n",
    "    n_time=21\n",
    "    n_nodes = 21 # number of nodes desired to be outputted for accumulation data\n",
    "    Lpt,Kt,rs = x\n",
    "    Svt = 200;      #tumor vascular density\n",
    "    D = Diffusion(rs)\n",
    "    kd = Blood_half_life(rs)\n",
    "    Perm,sigma=soluteperm(Lpt,rs) #Vascular permeability and solute reflection coefficient respectively\n",
    "    R=1. # Tumor Radius\n",
    "    Pv=25. #\n",
    "    Pvv=1.\n",
    "    att=R*sqrt(Lpt*Svt/Kt) #parameter alpha for tumor, ignoring lymphatics\n",
    "    Press(r)= a_pressure(r,att)\n",
    "    co=1.\n",
    "    tdomain=300\n",
    "\n",
    "    r= (1/R)*(range(0,stop=R,length=n_spatial))\n",
    "    dr=1/(n_spatial-1)\n",
    "\n",
    "    P= broadcast(Press,r)\n",
    "\n",
    "    c0=zeros(n_spatial,1)\n",
    "    c0[end]=0.\n",
    "    #define timespan\n",
    "    time_end =(300/3600)\n",
    "    tspan=(0.,time_end)\n",
    "    #package parameters\n",
    "    p=P,n_spatial,sigma,Perm,Lpt,Svt,Kt,Pvv,Pv,D,r,dr,kd\n",
    "    #Define and solve system of ODEs using time constant of 1/n_nodes and the QNDF solver for stiff systems\n",
    "    dt=time_end/(n_time-1)\n",
    "    t = dt*[i for i=0:n_nodes-1]\n",
    "    @suppress_err begin\n",
    "        prob=ODEProblem(MST!,c0,tspan,p)\n",
    "        sol= solve(prob,QNDF(),saveat=dt)\n",
    "        if length(sol.t)==1\n",
    "            prob=ODEProblem(MSTHighPeclet!,c0,tspan,p)\n",
    "            sol= solve(prob,QNDF(),saveat=dt)\n",
    "        end\n",
    "        accum_dat = Accumulation_Model(sol,n_spatial,n_time,n_nodes)\n",
    "        return SepBatch(Lpt,Kt,rs,t,accum_dat,n_nodes)\n",
    "    end\n",
    "end\n",
    "function prepdata(s)\n",
    "   x,y= accumlation_intermediate(s)\n",
    "    return Float32.(iscalerbatch(x)),reshape(Float32.(oscaler.(y)),(1,20))\n",
    "end\n",
    "#accum(x) =  output_scaler(accumlation_intermediate(invscaler(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ae89d5e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samp=Int64(1e4); # number of samples\n",
    "nbatch = n_samp; # number of batches\n",
    "nfeat=5; #5 features\n",
    "ntarg=1; #1 target values\n",
    "sl=20; #sequence length of 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "616599af",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create Data Sets\n",
    "#First Randomly Sample the biological parameters using a Sobol Sequence\n",
    "S=Float32.(reshape([x[j] for x in (sample(n_samp, lb,ub, SobolSample())) for j in eachindex(x)],(3,n_samp)));\n",
    "#Define empty arrays to store inputs and outputs into\n",
    "X = Array{Float32, 2}(undef,5, n_samp*20);\n",
    "Y= Array{Float32, 2}(undef,1, n_samp*20);\n",
    "#Fill in the arrays using the numerical model\n",
    "for i=1:n_samp\n",
    "    X[:,(20*(i-1)+1):(20*i)],Y[1,(20*(i-1)+1):(20*i)]= accumlation_intermediate(S[:,i][:])\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "27bf092d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×200000 Matrix{Float32}:\n",
       "  5.00824f-7   5.00824f-7   5.00824f-7  …  8.01437f-7  8.01437f-7  8.01437f-7\n",
       "  1.69998f-6   1.69998f-6   1.69998f-6     4.0468f-6   4.0468f-6   4.0468f-6\n",
       " 11.3431      11.3431      11.3431         7.64359     7.64359     7.64359\n",
       "  0.0          0.00416667   0.00833333     0.0708333   0.075       0.0791667\n",
       "  0.0          0.00512395   0.0102461      0.191997    0.203264    0.214528\n",
       "  0.00512395   0.0102461    0.0153663   …  0.203264    0.214528    0.225788"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = [X ; Y]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d661f552",
   "metadata": {},
   "outputs": [],
   "source": [
    "#writedlm( \"MLFinalData_v.csv\",  data, ',')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.1",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}