Changing output to nt

pythonScripts/compartmentModel.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 220,
+   "execution_count": 223,
    "metadata": {},
    "outputs": [
     {
@@ -12,7 +12,7 @@
       "Using /home/studen/temp/humanHGplusI\n",
       "At t=58.68\n",
       "Done simultaneous LSODA SE\n",
-      "Time: 2.736 s\n"
+      "Time: 2.793 s\n"
      ]
     }
    ],
@@ -33,6 +33,7 @@
     "import runSolver\n",
     "importlib.reload(runSolver)\n",
     "\n",
+    "  \n",
     "#run solver\n",
     "fh=os.path.expanduser('~')\n",
     "jobDir=os.path.join(fh,'temp','humanHGplusI')\n",
@@ -55,7 +56,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "jobDir=os.path.join(fh,'temp','current')\n",
+    "jobDir=os.path.join(fh,'temp','humanHGplusI')\n",
     "model=cModel.model()\n",
     "#sys.parse(os.path.join(fh,'software','src','Integra','models','cDiazepam.json'))\n",
     "setupFile=os.path.join(fh,'software','src','PBPK','setup','setupFastLoad.json')\n",

pythonScripts/ivp.py

@@ -66,31 +66,43 @@ def jacobiSEFull(t,S,system):
         fJ[i*system.n:(i+1)*system.n,i*system.n:(i+1)*system.n]=system.M(t)
     return fJ
 
-def solveSimultaneous(model,tmax,atol,rtol,method='LSODA',t0=0,y0=numpy.array([]), sIn=numpy.array([])):
-    s0=numpy.zeros((model.n,model.m+1))
-    if sIn.size!=0:
-       s0[:,1:]=sIn
-    if y0.size==0:
-       y0=numpy.zeros(model.n)
-    #set initial condition
-    s0[:,0]=y0
-    s0=s0.ravel()
-    sol=scipy.integrate.solve_ivp(dfdySFull,[t0, tmax],s0, args=(model,), jac=jacobiSEFull,
-                                  method=method, atol=atol, rtol=rtol)
-    t=sol.t
-    sFull=numpy.reshape(numpy.transpose(sol.y),(len(t),model.n,model.m+1))
-    s1=sFull[:,:,1:]
-    ysol=sFull[:,:,0]
-    se=model.calculateUncertainty(s1)
-    print('Done simultaneous LSODA SE')
-    return t,ysol,se,s1
-
-def solveSequential(model,tmax,atol,rtol,method='LSODA',t0=0,y0=numpy.array([]), sIn=numpy.array([])):
+def solveSimultaneous(model,tmax,atol,rtol,method='LSODA',nt=201,\
+      t0=0,y0=numpy.array([]), sIn=numpy.array([])):
+
+   t_eval=numpy.linspace(t0,tmax,nt)
+
+   s0=numpy.zeros((model.n,model.m+1))
+
+   if sIn.size!=0:
+      s0[:,1:]=sIn
+   if y0.size==0:
+      y0=numpy.zeros(model.n)
+
+   #set initial condition
+   s0[:,0]=y0
+   s0=s0.ravel()
+   sol=scipy.integrate.solve_ivp(dfdySFull,[t0, tmax],s0, args=(model,),\
+                                 jac=jacobiSEFull, method=method, \
+                                 atol=atol, rtol=rtol, t_eval=t_eval)
+   t=sol.t
+   sFull=numpy.reshape(numpy.transpose(sol.y),(len(t),model.n,model.m+1))
+   s1=sFull[:,:,1:]
+   ysol=sFull[:,:,0]
+   se=model.calculateUncertainty(s1)
+   print('Done simultaneous LSODA SE')
+   return t,ysol,se,s1
+
+def solveSequential(model,tmax,atol,rtol,method='LSODA',nt=201,\
+      t0=0,y0=numpy.array([]), sIn=numpy.array([])):
+
+   t_eval=numpy.linspace(t0,tmax,nt)
+
    if y0.size==0:
       y0=numpy.zeros(model.n)
-   model.iPrint=0
+   
    solIVP=scipy.integrate.solve_ivp(dfdy,[t0, tmax],y0, args=(model,), jac=jacobi,
-                                  method=method, atol=atol, rtol=rtol)
+                                  method=method, atol=atol, rtol=rtol, t_eval=t_eval)
+
    #y is n x nt (odeint nt x n)
    sol=numpy.transpose(solIVP.y)
    t=solIVP.t
@@ -102,7 +114,8 @@ def solveSequential(model,tmax,atol,rtol,method='LSODA',t0=0,y0=numpy.array([]),
    s0=s0.ravel()
     
    solIVPSE=scipy.integrate.solve_ivp(dfdyS,[0, tmax],s0, args=(model,), jac=jacobiSE,
-                               method=method, atol=atol, rtol=rtol)
+                               method=method, atol=atol, rtol=rtol,t_eval=t_eval)
+
    sraw=numpy.reshape(numpy.transpose(solIVPSE.y),(len(solIVPSE.t),model.n,model.m))
    #interpolate on t
    s1=numpy.zeros((len(t),model.n,model.m))

pythonScripts/runSolver.py

@@ -81,12 +81,16 @@ def main(parFiles,jobDir):
       t,sol,se,s1=ivp.solveSimultaneousOdeint(model,tmax,setup['nt'],t0,y0,sIn)
 
    if setup['mode']=='IVP':
-      t,sol,se,s1=ivp.solveSequential(model,tmax,atol=setup['atol'],\
+      t,sol,se,s1=ivp.solveSequential(model,tmax,\
+         nt=setup['nt'],\
+         atol=setup['atol'],\
          rtol=setup['rtol'],method=setup['method'],\
          t0=t0,y0=y0,sIn=sIn)
         
    if setup['mode']=='IVPSimultaneous':
-      t,sol,se,s1=ivp.solveSimultaneous(model,tmax,atol=setup['atol'],\
+      t,sol,se,s1=ivp.solveSimultaneous(model,tmax,\
+         nt=setup['nt'],\
+         atol=setup['atol'],\
          rtol=setup['rtol'],method=setup['method'],\
          t0=t0,y0=y0,sIn=sIn)