Removing placenta from humanHG.json and creating (uncomplete) humanHG_pregnancy.json (to avoid zero mass divergent transfer compartments); allowing for interpolation of parameter specific derivatives to qt.txt, sOut.txt; adding write3D and read3D to runSolver to store 3D arrays; some cleanup of setup parsing in runSolver

models/humanHG.json

@@ -3,7 +3,7 @@
    "compartments":[
       "redBloodCells","plasma","venous","kidney","urine",
       "richlyPerfused","fat","slowlyPerfused","hair",
-      "brainBlood","brain","placenta",
+      "brainBlood","brain",
       "liver","gut","intestine","feces","inorganicMercury"],
    "commentParameters":"every parameter must have a value, name is optional and used in SE computation",
    "volumeComment":"in a compartment w/o volume, value is mass",
@@ -15,7 +15,6 @@
       "fat":"fatVolume",
       "slowlyPerfused":"slowlyPerfusedVolume",
       "brainBlood":"brainBloodVolume",
-      "placenta":"placentaVolume",
       "gut":"gutVolume",
       "intestine":"intestineVolume",
       "redBloodCells":"redBloodCellsVolume",
@@ -29,7 +28,6 @@
       "(plasma,venous):fat":"fatFlow",
       "(plasma,venous):slowlyPerfused":"slowlyPerfusedFlow",
       "(plasma,venous):brainBlood":"brainBloodFlow",
-      "(plasma,venous):placenta":"placentaFlow",
       "plasma:liver":"liverInFlow",
       "liver:venous":"liverOutFlow",
       "(plasma,liver):gut":"gutFlow",
@@ -42,9 +40,9 @@
    "bindings": {
       "flowComments":"group common sources or targets with brackets",
       "flow":[
-         "(kidney,richlyPerfused,fat,slowlyPerfused,brainBlood,placenta,liver)->venous",
+         "(kidney,richlyPerfused,fat,slowlyPerfused,brainBlood,liver)->venous",
          "gut->liver",
-         "plasma->(kidney,richlyPerfused,fat,slowlyPerfused,brainBlood,placenta,liver,gut)",
+         "plasma->(kidney,richlyPerfused,fat,slowlyPerfused,brainBlood,liver,gut)",
          "venous->plasma"],
       "diffusion":{
          "plasma->redBloodCells":"kRBC",
@@ -66,11 +64,9 @@
        "fat":"fatPC",
        "slowlyPerfused":"slowlyPerfusedPC",
        "brainBlood":"brainBloodPC",
-       "placenta":"placentaPC",
        "liver":"liverPC",
        "gut":"gutPC",
        "brainBlood->brain":"brainPC",
        "hair->venous":"hairPC",
-       "redBloodCells->plasma":"rbcPC"},
-     "parameters":"humanHG_parameters1.json"
+       "redBloodCells->plasma":"rbcPC"}
 }

models/humanHG_parameters1.json

@@ -40,9 +40,6 @@
         "liverVolume":{"derived":"liverVolume"},
         "gutVolume":{"derived":"gutVolume"},
         "intestineVolume":{"derived":"intestineVolume"},
-        "placentaVolume":{"derived":"placentaVolume"},
-        "fetusWeight":{"function":"fetusWeight"},
-        "fetusWeightMax":{"value":4.0,"unit":"kg","dist":"normal","cv":0.1},
         "commentFractions":"X",
         "plasmaVolumeFraction":{"value":0.024,"dist":"normal","cv":0.14},
         "redBloodCellsVolumeFraction":{"value":0.024,"dist":"normal","cv":0.25},
@@ -54,8 +51,6 @@
         "brainBloodVolumeFraction":{"value":0.007,"dist":"normal","cv":0.3},
         "brainVolumeFraction":{"value":0.02,"dist":"normal","cv":0.3},
         "liverVolumeFraction":{"value":0.026,"dist":"normal","cv":0.25},
-        "placentaVolumeFraction":{"value":0.18,"dist":"normal","cv":0.13},
-        "placentaVolumeFractionSurce":"https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3530253",
         "gutVolumeFraction":{"value":0.017,"dist":"normal","cv":0.15},
         "intestineVolumeFraction":{"value":0.14,"dist":"normal","cv":0.3},
         "commentDiff":"X",
@@ -75,10 +70,7 @@
         "heavisideTransition":{"value":100},
         "hairGrowthStart":{"value":0e6,"unit":"min"},
         "hairGrowthStop":{"derived":"hairGrowthStop"},
-        "hairGrowth":{"derived":"hairGrowth"},
-        "fetusGrowthStart":{"value":8.6e5,"unit":"min"},
-        "fetusGrowthDuration":{"value":4e5,"unit":"min"},
-        "fetusGrowthStop":{"derived":"fetusGrowthStop"}},
+        "hairGrowth":{"derived":"hairGrowth"}},
     "functions":{
        "commentHairGrowth":"2yr hair growth to a total mass of 140g",
        "hairGrowth":{
@@ -87,14 +79,7 @@
            "W1":"heavisideTransition",
            "t2":"hairGrowthStop",
            "W2":"heavisideTransition",
-           "C":"hairVolumeMax"},
-       "fetusWeight":{
-           "type":"linearGrowth",
-           "t1":"fetusGrowthStart",
-           "t2":"fetusGrowthStop",
-           "W1":"heavisideTransition",
-           "W2":"heavisideTransition",
-           "C":"fetusWeightMax"}},
+           "C":"hairVolumeMax"}},
     "derivedParameters":{
         "plasmaVolume":{
            "type":"product",
@@ -140,10 +125,6 @@
            "type":"product",
            "a":"intestineVolumeFraction",
            "b":"bodyWeight"},
-        "placentaVolume":{
-           "type":"product",
-           "a":"placentaVolumeFraction",
-           "b":"fetusWeight"},
        "hairVolumeMax":{
            "type":"product",
            "a":"hairVolumeFraction",
@@ -160,10 +141,6 @@
             "type":"ratio",
             "a":"kH",
             "b":"hairPC"},
-       "fetusGrowthStop":{
-            "type":"sum",
-            "a":"fetusGrowthStart",
-            "b":"fetusGrowthDuration"},
        "bwPower0p75":{
           "type":"power",
           "a":"bodyWeight",

models/humanHG_parameters3.json

@@ -0,0 +1,179 @@
+{"parameters":{
+        "bodyWeight":{"value":70,"units":"kg"},
+        "bodyWeight0p75":{"derived":"bwPower0p75"},
+        "kidneyPC":{"value":4.0},
+        "richlyPerfusedPC":{"value":1},
+        "fatPC":{"value":0.15},
+        "slowlyPerfusedPC":{"value":2},
+        "brainBloodPC":{"value":1},
+        "brainPC":{"value":3},
+        "hairPCBlood":{"value":248.7,"dist":"lognormal","cv":0.7},
+        "hairPC":{"derived":"hairPC"},
+        "placentaPC":{"value":2},
+	     "liverPC":{"value":5},
+        "gutPC":{"value":1},
+        "rbcPC":{"value":12},
+        "flowNotes":"assuming 70kg individual, fetus at 0.1kg",
+        "plasmaFlow":{"value":8.067,"unit":"l/min"},
+        "plasmaFlowScaled":{"value":0.333,"unit":"l/min"},
+        "kidneyFlow":{"value":1.412,"unit":"l/min"},
+        "richlyPerfusedFlow":{"value":1.484,"unit":"l/min"},
+        "fatFlow":{"value":0.419,"unit":"l/min"},
+        "slowlyPerfusedFlow":{"value":2.01,"unit":"l/min"},
+        "brainBloodFlow":{"value":0.920,"unit":"l/min"},
+        "placentaFlow":{"value":0.173,"unit":"l/min"},
+        "liverInFlow":{"value":0.391,"unit":"l/min"},
+        "liverOutFlow":{"value":1.928,"unit":"l/min"},
+        "gutFlow":{"value":1.538,"unit":"l/min"},
+        "commentVolumes":"X",
+        "plasmaVolume":{"derived":"plasmaVolume"},
+        "redBloodCellsVolume":{"derived":"redBloodCellsVolume"},
+        "kidneyVolume":{"derived":"kidneyVolume"},
+        "richlyPerfusedVolume":{"derived":"richlyPerfusedVolume"},
+        "fatVolume":{"derived":"fatVolume"},
+        "hairVolume":{"function":"hairGrowth"},
+        "hairGrowthDuration":{"derived":"hairGrowthDuration"},
+        "hairVolumeMax":{"derived":"hairVolumeMax"},
+        "slowlyPerfusedVolume":{"derived":"slowlyPerfusedVolume"},
+        "brainBloodVolume":{"derived":"brainBloodVolume"},
+        "brainVolume":{"derived":"brainVolume"},
+        "liverVolume":{"derived":"liverVolume"},
+        "gutVolume":{"derived":"gutVolume"},
+        "intestineVolume":{"derived":"intestineVolume"},
+        "placentaVolume":{"derived":"placentaVolume"},
+        "fetusWeight":{"function":"fetusWeight"},
+        "fetusWeightMax":{"value":4.0,"unit":"kg","dist":"normal","cv":0.1},
+        "commentFractions":"X",
+        "plasmaVolumeFraction":{"value":0.024,"dist":"normal","cv":0.14},
+        "redBloodCellsVolumeFraction":{"value":0.024,"dist":"normal","cv":0.25},
+        "kidneyVolumeFraction":{"value":0.004,"dist":"normal","cv":0.3},
+        "richlyPerfusedVolumeFraction":{"value":0.1,"dist":"normal","cv":0.3},
+        "fatVolumeFraction":{"value":0.273,"dist":"normal","cv":0.24},
+        "hairVolumeFraction":{"value":0.002,"dist":"normal","cv":0.5},
+        "slowlyPerfusedVolumeFraction":{"value":0.35,"dist":"normal","cv":0.16},
+        "brainBloodVolumeFraction":{"value":0.007,"dist":"normal","cv":0.3},
+        "brainVolumeFraction":{"value":0.02,"dist":"normal","cv":0.3},
+        "liverVolumeFraction":{"value":0.026,"dist":"normal","cv":0.25},
+        "placentaVolumeFraction":{"value":0.18,"dist":"normal","cv":0.13},
+        "placentaVolumeFractionSurce":"https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3530253",
+        "gutVolumeFraction":{"value":0.017,"dist":"normal","cv":0.15},
+        "intestineVolumeFraction":{"value":0.14,"dist":"normal","cv":0.3},
+        "commentDiff":"X",
+        "kRBC":{"value":0.60,"unit":"l/min"},
+        "kU":{"value":0,"unit":"l/min"},
+        "commentkH":"assume pc-times larger input flow, 2.823e-6*248.7",
+        "kH":{"value":2.823e-6,"unit":"l/min"},
+        "kBR":{"value":4.033e-3,"unit":"l/min"},
+        "kI":{"value":4.033e-6,"unit":"l/min"},
+        "kB":{"value":4.033e-5,"unit":"l/min"},
+        "kD":{"value":4.033e-5,"unit":"l/min"},
+        "kF":{"value":8.066e-5,"unit":"l/min"},
+        "kF1":{"value":0,"unit":"l/min"},
+        "kR":{"value":2.016e-3,"unit":"l/min"},
+        "zero":{"value":0},
+        "threeQuarters":{"value":0.75},
+        "heavisideTransition":{"value":100},
+        "hairGrowthStart":{"value":0e6,"unit":"min"},
+        "hairGrowthStop":{"derived":"hairGrowthStop"},
+        "hairGrowth":{"derived":"hairGrowth"},
+        "fetusGrowthStart":{"value":8.6e5,"unit":"min"},
+        "fetusGrowthDuration":{"value":4e5,"unit":"min"},
+        "fetusGrowthStop":{"derived":"fetusGrowthStop"}},
+    "functions":{
+       "commentHairGrowth":"2yr hair growth to a total mass of 140g",
+       "hairGrowth":{
+           "type":"linearGrowth",
+           "t1":"hairGrowthStart",
+           "W1":"heavisideTransition",
+           "t2":"hairGrowthStop",
+           "W2":"heavisideTransition",
+           "C":"hairVolumeMax"},
+       "fetusWeight":{
+           "type":"linearGrowth",
+           "t1":"fetusGrowthStart",
+           "t2":"fetusGrowthStop",
+           "W1":"heavisideTransition",
+           "W2":"heavisideTransition",
+           "C":"fetusWeightMax"}},
+    "derivedParameters":{
+        "plasmaVolume":{
+           "type":"product",
+           "a":"plasmaVolumeFraction",
+           "b":"bodyWeight"},
+        "redBloodCellsVolume":{
+           "type":"product",
+           "a":"redBloodCellsVolumeFraction",
+           "b":"bodyWeight"},
+        "kidneyVolume":{
+           "type":"product",
+           "a":"kidneyVolumeFraction",
+           "b":"bodyWeight"},
+        "richlyPerfusedVolume":{
+           "type":"product",
+           "a":"richlyPerfusedVolumeFraction",
+           "b":"bodyWeight"},
+        "fatVolume":{
+           "type":"product",
+           "a":"fatVolumeFraction",
+           "b":"bodyWeight"},
+        "slowlyPerfusedVolume":{
+           "type":"product",
+           "a":"slowlyPerfusedVolumeFraction",
+           "b":"bodyWeight"},
+        "brainBloodVolume":{
+           "type":"product",
+           "a":"brainBloodVolumeFraction",
+           "b":"bodyWeight"},
+        "brainVolume":{
+           "type":"product",
+           "a":"brainVolumeFraction",
+           "b":"bodyWeight"},
+        "liverVolume":{
+           "type":"product",
+           "a":"liverVolumeFraction",
+           "b":"bodyWeight"},
+        "gutVolume":{
+           "type":"product",
+           "a":"gutVolumeFraction",
+           "b":"bodyWeight"},
+        "intestineVolume":{
+           "type":"product",
+           "a":"intestineVolumeFraction",
+           "b":"bodyWeight"},
+        "placentaVolume":{
+           "type":"product",
+           "a":"placentaVolumeFraction",
+           "b":"fetusWeight"},
+       "hairVolumeMax":{
+           "type":"product",
+           "a":"hairVolumeFraction",
+           "b":"bodyWeight"},
+        "hairGrowthStop":{
+           "type":"sum",
+           "a":"hairGrowthStart",
+           "b":"hairGrowthDuration"},
+        "hairGrowthDuration":{
+            "type":"ratio",
+            "a":"hairVolumeMax",
+            "b":"hairGrowth"},
+        "hairGrowth":{
+            "type":"ratio",
+            "a":"kH",
+            "b":"hairPC"},
+       "fetusGrowthStop":{
+            "type":"sum",
+            "a":"fetusGrowthStart",
+            "b":"fetusGrowthDuration"},
+       "bwPower0p75":{
+          "type":"power",
+          "a":"bodyWeight",
+          "n":"threeQuarters"},
+       "plasmaFlow":{
+          "type":"product",
+          "a":"plasmaFlowScaled",
+          "b":"bodyWeight0p75"},
+       "hairPC":{
+          "type":"ratio",
+          "a":"hairPCBlood",
+          "b":"slowlyPerfusedPC"}}
+}

models/humanHG_pregnancy.json

@@ -0,0 +1,76 @@
+{
+   "timeUnit":"min",
+   "compartments":[
+      "redBloodCells","plasma","venous","kidney","urine",
+      "richlyPerfused","fat","slowlyPerfused","hair",
+      "brainBlood","brain","placenta",
+      "liver","gut","intestine","feces","inorganicMercury"],
+   "commentParameters":"every parameter must have a value, name is optional and used in SE computation",
+   "volumeComment":"in a compartment w/o volume, value is mass",
+   "volumes":{
+      "plasma":"plasmaVolume",
+      "venous":"plasmaVolume",
+      "kidney":"kidneyVolume",
+      "richlyPerfused":"richlyPerfusedVolume",
+      "fat":"fatVolume",
+      "slowlyPerfused":"slowlyPerfusedVolume",
+      "brainBlood":"brainBloodVolume",
+      "placenta":"placentaVolume",
+      "gut":"gutVolume",
+      "intestine":"intestineVolume",
+      "redBloodCells":"redBloodCellsVolume",
+      "hair":"hairVolume",
+      "brain":"brainVolume",
+      "liver":"liverVolume"
+      },
+   "flows":{
+      "(plasma,venous):kidney":"kidneyFlow",
+      "(plasma,venous):richlyPerfused":"richlyPerfusedFlow",
+      "(plasma,venous):fat":"fatFlow",
+      "(plasma,venous):slowlyPerfused":"slowlyPerfusedFlow",
+      "(plasma,venous):brainBlood":"brainBloodFlow",
+      "(plasma,venous):placenta":"placentaFlow",
+      "plasma:liver":"liverInFlow",
+      "liver:venous":"liverOutFlow",
+      "(plasma,liver):gut":"gutFlow",
+      "plasma:venous":"plasmaFlow"},
+   "commentSources":"constants in concentration units per s",
+   "commentHeavyside":"value is mass of exogene/volume of container/duration",
+   "sources":{
+      "intestine":{"name":"constant","value":0.0486,"unit":"ug/min"}},
+   "commentBindings":"split to flow and diffusion dominated",
+   "bindings": {
+      "flowComments":"group common sources or targets with brackets",
+      "flow":[
+         "(kidney,richlyPerfused,fat,slowlyPerfused,brainBlood,placenta,liver)->venous",
+         "gut->liver",
+         "plasma->(kidney,richlyPerfused,fat,slowlyPerfused,brainBlood,placenta,liver,gut)",
+         "venous->plasma"],
+      "diffusion":{
+         "plasma->redBloodCells":"kRBC",
+         "redBloodCells->plasma":"kRBC",
+         "kidney->urine":"kU",
+         "slowlyPerfused->hair":"kH",
+         "hair->venous":"kH",
+         "brainBlood->brain":"kBR",
+         "liver->intestine":"kB",
+         "liver->inorganicMercury":"kI",
+         "intestine->gut":"kR",
+         "intestine->inorganicMercury":"kD",
+         "intestine->feces":"kF",
+         "inorganicMercury->feces":"kF1"}},
+    "pcComments":"always related to outflow of organs", 
+    "partitionCoefficients":{
+       "kidney":"kidneyPC",
+       "richlyPerfused":"richlyPerfusedPC",
+       "fat":"fatPC",
+       "slowlyPerfused":"slowlyPerfusedPC",
+       "brainBlood":"brainBloodPC",
+       "placenta":"placentaPC",
+       "liver":"liverPC",
+       "gut":"gutPC",
+       "brainBlood->brain":"brainPC",
+       "hair->venous":"hairPC",
+       "redBloodCells->plasma":"rbcPC"},
+     "parameters":"humanHG_parameters1.json"
+}

pythonScripts/ivp.py

@@ -73,7 +73,7 @@ def solveSimultaneous(sys,tmax,atol,rtol,method='LSODA'):
     ysol=sFull[:,:,0]
     se=sys.calculateUncertainty(ysol,s1)
     print('Done simultaneous LSODA SE')
-    return t,ysol,se
+    return t,ysol,se,s1
 
 def solveSequential(sys,tmax,atol,rtol,method='LSODA'):
    y0=numpy.zeros(sys.n)
@@ -98,7 +98,7 @@ def solveSequential(sys,tmax,atol,rtol,method='LSODA'):
          s1[:,i,j]=scipy.interpolate.splev(t, tck, der=0)
 
    se=sys.calculateUncertainty(sol,s1)
-   return t,sol,se
+   return t,sol,se,s1
 
 def solveSimultaneousOdeint(sys,tmax,nt=201):
    t = numpy.linspace(0,tmax, nt)
@@ -113,7 +113,7 @@ def solveSimultaneousOdeint(sys,tmax,nt=201):
    sol=sFull[:,:,0]
    se=sys.calculateUncertainty(sol,s1)
    print('Done simultaneous SE')
-   return t,sol,se
+   return t,sol,se,s1
 
 def solveSequentialOdeint(sys,tmax,nt=201):
    t = numpy.linspace(0,tmax, nt)
@@ -128,6 +128,6 @@ def solveSequentialOdeint(sys,tmax,nt=201):
    s1=numpy.reshape(solSE,(len(t),sys.n,sys.m))
    se=sys.calculateUncertainty(sol,s1)
    print('Done sequential SE')
-   return t,sol,se
+   return t,sol,se,s1
 
 

pythonScripts/runSolver.py

@@ -5,9 +5,9 @@ import sys
 import json
 import numpy
 import time
+import scipy.interpolate
 
-def get(setup,par):
-   defaultValues={\
+defaultValues={\
       'method':'LSODA',\
       'atol':1e-4,\
       'rtol':1e-4,\
@@ -15,52 +15,111 @@ def get(setup,par):
       'mode':'IVPSimultaneous',\
       'nt':201,
       'tUnit':'day'}
+
+def get(setup,par):
    try:
       return setup[par]
    except KeyError:
       pass
    return defaultValues[par]
 
-def main(parFiles,jobDir):
-   sys=cModel.model()
-   setupFile=parFiles[1]
-   parameterFile=parFiles[2]
-   sys.parse(setupFile,parameterFile)
-   with open(parFiles[0],'r') as f:
+def parseSetup(setupFile):
+   with open(setupFile,'r') as f:
       setup=json.load(f)
-   tmax=get(setup,'tmax')
-   atol=get(setup,'atol')
-   rtol=get(setup,'rtol')
-   mode=get(setup,'mode')
-   method=get(setup,'method')
-   nt=get(setup,'nt')
-   tUnit=get(setup,'tUnit')
+   out={}
+   for p in defaultValues:
+      out[p]=get(setup,p)
+   return out
+
+def getScale(setup):
+   tUnit=setup['tUnit']
    if tUnit=='min':
-      pass
+      return 1
    if tUnit=='hour':
-      tmax*=60
+      return 60
    if tUnit=='day':
-      tmax*=24*60
+      return 24*60
    if tUnit=='month':
-      tmax*=24*60*30
+      return 24*60*30
    if tUnit=='year':
-      tmax*=24*60*30*365
+      return 24*60*30*365
+   return 1
+
+
+
+def main(parFiles,jobDir):
+   model=cModel.model()
+   setupFile=parFiles[0]
+   modelFile=parFiles[1]
+   parameterFile=parFiles[2]
+   model.parse(modelFile,parameterFile)
+   setup=parseSetup(setupFile)
+   scale=getScale(setup)
+   tmax=setup['tmax']*scale
 
    start_time=time.time()
-   if mode=='SequentialOdeint':
-      t,sol,se=ivp.solveSequentialOdeint(sys,tmax,nt)
-   if mode=='SimultaneousOdeint':
-      t,sol,se=ivp.solveSimultaneousOdeint(sys,tmax,nt)
+   if setup['mode']=='SequentialOdeint':
+      t,sol,se,s1=ivp.solveSequentialOdeint(model,tmax,setup['nt'])
+   if setup['mode']=='SimultaneousOdeint':
+      t,sol,se,s1=ivp.solveSimultaneousOdeint(model,tmax,setup['nt'])
 
-   if mode=='IVP':
-      t,sol,se=ivp.solveSequential(sys,tmax,atol=atol,rtol=rtol,method=method)
+   if setup['mode']=='IVP':
+      t,sol,se,s1=ivp.solveSequential(model,tmax,atol=setup['atol'],
+         rtol=setup['rtol'],method=setup['method'])
         
-   if mode=='IVPSimultaneous':
-      t,sol,se=ivp.solveSimultaneous(sys,tmax,atol=atol,rtol=rtol,method=method)
+   if setup['mode']=='IVPSimultaneous':
+      t,sol,se,s1=ivp.solveSimultaneous(model,tmax,atol=setup['atol'],
+         rtol=setup['rtol'],method=setup['method'])
 
    end_time=time.time()
    print('Time: {:.3f} s'.format(end_time-start_time))
    #store
+   #interpolate on s1
+   #s1 has shape ft x n x m
+   #where ft is LSODE driven number of points
+   qt,sOut=interpolate(setup,model,t,s1,tmax)
+
    numpy.savetxt(os.path.join(jobDir,'t.txt'),t)
    numpy.savetxt(os.path.join(jobDir,'sol.txt'),sol)
    numpy.savetxt(os.path.join(jobDir,'se.txt'),se)
+   numpy.savetxt(os.path.join(jobDir,'qt.txt'),qt)
+   #this is 3D, so new routines
+   write3D(os.path.join(jobDir,'sOut.txt'),sOut)
+
+def interpolate(setup,model,t,s1,tmax):
+   #interpolate on s1
+   #s1 has shape ft x n x m
+   #where ft is LSODE driven number of points
+   sOut=numpy.zeros((setup['nt'],model.n,model.m))
+   qt=numpy.linspace(0,tmax,setup['nt'])
+   for i in range(model.n):
+      for j in range(model.m):
+         y=s1[:,i,j]
+         f = scipy.interpolate.interp1d(t, y, kind='cubic')
+         sOut[:,i,j]=f(qt)
+   return qt,sOut
+
+
+
+def write3D(fName,a):
+   with open(fName,'w') as f:
+      f.write('#Shape {}\n'.format(a.shape))
+      i=0
+      for data_slice in a:
+         f.write('#\t Slice {}\n'.format(i))
+         numpy.savetxt(f,data_slice)
+         i+=1
+      
+def read3D(fName):
+   new_data=numpy.loadtxt(fName)
+   #read and parse first row to get shape
+   with open(fName,'r') as f:
+      firstLine=f.readline()
+
+   shape=firstLine.replace('#Shape ','').\
+      replace('\n','').\
+      replace('(','').\
+      replace(')','').\
+      split(',')
+   shape=[int(x) for x in shape]
+   return new_data.reshape(shape)