Replace pk1Comp.py

pk_optimizer/pk1Comp.py

 #!/usr/bin/env python
 # coding: utf-8
 
-# In[1]:
-
-
 # Import commands
 from scipy.stats import gamma
 import numpy as np
@@ -16,7 +13,7 @@ class pk1Comp:
     
     """The pk1Comp object is a one compartment PK model that outputs graphs of mass of tracer over time."""
 
-    def __init__ (self, numParam = 4, Flow = 1, Vp = 0.1, Visf = 0.5, PS = 0.15):
+    def __init__ (self, numParam = 4, Flow = 1, Vol = 0.5, PS = 0.15):
         
         """Initializes the model with default parameter values for flow, Vp, Visf, and PS.
         Parameters
@@ -36,24 +33,26 @@ class pk1Comp:
         PS : double
             PS is the permeability-surface area constant in mL/(g*min). Defaults to 0.15.      
         """
-        
+        if numParam <= 0 or Flow <= 0 or Vol < 0 or PS<0 or PS>10 or Vol>10:
+            raise ValueError("Input values are incorrect.")
+            
         # Declare Variables for initial conditions
         self.numParam = numParam
         self.Flow = Flow
-        self.Vp = Vp
-        self.Visf = Visf
+        self.Vol = Vol
         self.PS = PS
-        C0 = 0 # Initial concentration of tracer in plasma
-        tmax = 10 #Time in seconds
-        dt = 0.1 #Time step
-        a = 2. # Alpha for gamma distribution
-        rv = gamma(a, loc = 2, scale = 0.65) #input function
+        self.C0 = 0 # Initial concentration of tracer in plasma
+        self.tmax = 10 #Time in seconds
+        self.dt = 1 #Time step
+        self.a = 2. # Alpha for gamma distribution
+        self.rv = gamma(self.a, loc = 2, scale = 0.65) #input function
 
+        self.sol = []
         # Define the time array
-        time = np.arange(0, tmax + dt, dt)
+        self.time = np.arange(0, self.tmax + self.dt, self.dt)
         
     # Derivative function
-    def derivs(curr_vals, time):
+    def derivs(self, curr_vals, time):
         """Finds derivatives of ODEs.
         
         Parameters
@@ -71,13 +70,13 @@ class pk1Comp:
         """
 
         # Define value of input function Cin
-        Cin = rv.pdf(time)    
+        Cin = self.rv.pdf(time)    
 
         # Unpack the current values of the variables we wish to "update" from the curr_vals list
         C = curr_vals
 
         # Right-hand side of odes, which are used to computer the derivative
-        dC_dt = flow*(Cin - C)/Vol
+        dC_dt = self.Flow*(Cin - C)/self.Vol
         #Cout = C
         return dC_dt
     
@@ -93,8 +92,8 @@ class pk1Comp:
         # Plot the results using the values stored in the solution variable, "sol"
         # Plot Cp using the "0" element from the solution
         plt.figure(1)
-        plt.plot(time, rv.pdf(time), color = 'blue', label = 'Input Function')
-        plt.plot(time, sol[:,0],color="green", label = 'Cout')
+        plt.plot(self.time, self.rv.pdf(self.time), color = 'blue', label = 'Input Function')
+        plt.plot(self.time, self.sol[:,0],color="green", label = 'Cout')
 
         # Plot Cisf using the "1" element from the solution
         #plt.plot(time, sol[:,1],color="purple", label = 'Cisf')
@@ -107,34 +106,7 @@ class pk1Comp:
         """Main function to run and solve ODEs"""
         
         # Store the initial values in a list
-        init = [C0]
+        init = [self.C0]
 
         # Solve the odes with odeint
-        sol = odeint(derivs, init, time)
-
-        #Mass_plasma = Vp * sol[:,0] #mass of tracer in plasma
-        #Mass_isf = Visf * sol[:,1] #mass of tracer in isf
-        #Tp = Vp/(flow + PS) # mean transit time
-        #E = 1 - np.exp(-PS/flow) #extraction fraction
-        #Q = Mass_plasma + Mass_isf
-
-        #print('The mean transit time is ' + str(Tp))
-        #print('The extraction fraction is ' + str(E))
-
-        # Plot mass of tracer using the "2" element from the solution
-        #plt.figure(2)
-        #plt.plot(time, Mass_plasma,color="red", label = 'Plasma')
-        # Plot mass of tracer in tissue using the "3" element from the solution
-        #plt.plot(time, Mass_isf,color="black", label = 'Interstitial Space')
-        #plt.plot(time, Q, color="blue", label = 'Total mass')
-        #plt.xlabel('Time [s]')
-        #plt.ylabel('Mass [mg]')
-        #plt.legend(loc = 'best')
-        #plt.grid()
-
-
-# In[ ]:
-
-
-
-
+        self.sol = odeint(self.derivs, init, self.time).round(4)