Updated to work with Python3

README.md

@@ -11,20 +11,37 @@ and then adding "spy" as in the English word *spy*.
 
 ### How do I install Odespy?
 
-#### Pip
 
-Odespy requires Python version 2.7.  The simplest procedure for
-installing Odespy is to use `pip`:
+Odespy works with Python 2.7, and Python 3 (tested with v3.6). 
+
+<!--
+#### Pip
+The simplest procedure for installing Odespy is to use `pip`:
 
 
 ```
 Terminal> sudo pip install -e git+https://github.com/hplgit/odespy.git#egg=odespy
 ```
+commented out this section as it does not work currently
+-->
+
+
 
 #### Clone/fork repo
 
-Alternatively, you can check out this repo and run `setup.py`:
+You can check out this repo and run `setup.py`:
+
+
+##### For Python 3:
+
+```
+Terminal> git clone --single-branch -b python3 https://github.com/briandrawert/odespy.git
+Terminal> cd odespy
+Terminal> sudo python setup.py install
+```
+
 
+##### For Python 2:
 
 ```
 Terminal> git clone [email protected]:hplgit/odespy.git

odespy/PyDSTool.py

@@ -1,7 +1,7 @@
 # Author: Liwei Wang
 """
 """
-from solvers import *
+from .solvers import *
 import numpy as np
 
 class Pyds(Solver):
@@ -24,10 +24,10 @@ class Pyds(Solver):
 
     def initialize(self):
         try:
-            import PyDSTool
+            from . import PyDSTool
         except ImportError:
-            raise ImportError,'''
-        PyDSTool is not installed - required for solvers from PyDSTool'''
+            raise ImportError('''
+        PyDSTool is not installed - required for solvers from PyDSTool''')
 
     def solve(self, time_points, terminate=None):
         # Common parts as superclass
@@ -45,7 +45,7 @@ def solve(self, time_points, terminate=None):
         # through Python dictionaries with string keys.
 
         # Start setting for PyDSTool
-        import PyDSTool
+        from . import PyDSTool
         neq, f, u0 = self.neq, self.f, self.U0
 
         # Initialize variables as trajectories in PyDSTOOL
@@ -139,7 +139,7 @@ def initialize_for_solve(self):
         if not hasattr(self,'init_step'):
             self.init_step = self.t[1] - self.t[0]
         self.params_pydstool = dict(\
-            (key,value) for key,value in self.__dict__.items() \
+            (key,value) for key,value in list(self.__dict__.items()) \
             if key in self._params_pydstool)    
 
 if __name__ == '__main__':
@@ -148,7 +148,7 @@ def initialize_for_solve(self):
     method = Vode_pyds(f)
     method.set_initial_condition([0.,1.])
     u,t = method.solve(np.linspace(0.,10.,50))
-    print u
+    print(u)
     import scitools.std as st
     st.plot(t,u[:,0])
-    print max(u[:,0]-np.sin(t))
+    print(max(u[:,0]-np.sin(t)))

odespy/RungeKutta.py

@@ -1,4 +1,5 @@
-from solvers import Solver, Adaptive
+from __future__ import print_function
+from .solvers import Solver, Adaptive
 import numpy as np
 
 def _calculate_order_1_level(coefficients):
@@ -134,7 +135,7 @@ def middle(x,y,z):    # Auxilary function
         k = np.zeros((k_len, self.neq), self.dtype)  # intern stages
 
         if self.verbose > 0:
-            print 'advance solution in [%s, %s], h=%g' % (t_n, t_next, h)
+            print('advance solution in [%s, %s], h=%g' % (t_n, t_next, h))
 
         # Loop until next time point is reached
         while (abs(t - t_n) < abs(t_next - t_n)):
@@ -150,7 +151,7 @@ def middle(x,y,z):    # Auxilary function
 
             self.info['rejected'] += 1  # reduced below if accepted
             if self.verbose > 0:
-                print '  u(t=%g)=%g: ' % (t+h, u_new),
+                print('  u(t=%g)=%g: ' % (t+h, u_new), end=' ')
 
             # local error between 2 levels
             error = h*np.abs(np.dot(factors_error, k))
@@ -171,18 +172,18 @@ def middle(x,y,z):    # Auxilary function
                 self.info['rejected'] -= 1
 
                 if self.verbose > 0:
-                    print 'accepted, ',
+                    print('accepted, ', end=' ')
             else:
                 if self.verbose > 0:
-                    print 'rejected, ',
+                    print('rejected, ', end=' ')
 
             if self.verbose > 0:
-                print 'err=%s, ' % str(error),
+                print('err=%s, ' % str(error), end=' ')
                 if hasattr(self, 'u_exact') and callable(self.u_exact):
-                    print 'exact-err=%s, ' % \
-                          (np.asarray(self.u_exact(t+h))-u_new),
+                    print('exact-err=%s, ' % \
+                          (np.asarray(self.u_exact(t+h))-u_new), end=' ')
                 if h <= self.min_step:
-                    print 'h=min_step!! ',
+                    print('h=min_step!! ', end=' ')
 
 
            # Replace 0 values by 1e-16 since we will divide by error
@@ -209,7 +210,7 @@ def middle(x,y,z):    # Auxilary function
             h = min(h, t_next - t_intermediate[-1])
 
             if self.verbose > 0:
-                print 'new h=%g' % h
+                print('new h=%g' % h)
 
             if h == 0:
                 break
@@ -367,16 +368,16 @@ def validate_data(self):
         # Check for dimension of user-defined butcher table.
         array_shape = self.butcher_tableau.shape
         if len(array_shape) is not 2:
-            raise ValueError,'''
-        Illegal input! Your input butcher_tableau should be a 2d-array!'''
+            raise ValueError('''
+        Illegal input! Your input butcher_tableau should be a 2d-array!''')
         else:
             m,n = array_shape
             if m not in (n, n + 1):
-                raise ValueError, '''\
+                raise ValueError('''\
         The dimension of 2d-array <method_yours_array> should be:
         1. Either (n, n), --> For 1-level RungeKutta methods
         2. Or (n+1, n),   --> For 2-levels RungeKutta methods
-        The shape of your input array is (%d, %d).''' % (m,n)
+        The shape of your input array is (%d, %d).''' % (m,n))
         self._butcher_tableau = self.butcher_tableau
 
         # Check for user-defined order,
@@ -393,19 +394,19 @@ def validate_data(self):
             if array_shape[0] == array_shape[1] + 1:
                 # 2-level RungeKutta methods
                 if type(self.method_order) is int:
-                    raise ValueError, error_2level
+                    raise ValueError(error_2level)
                 try:
                     order1, order2 = self.method_order
                     if abs(order1-order2) != 1 or \
                             order1 < 1 or order2 < 1:
-                        raise ValueError, error_2level
+                        raise ValueError(error_2level)
                 except:
-                    raise ValueError,error_2level
+                    raise ValueError(error_2level)
             else:
                 # 1-level RungeKutta methods
                 if type(self.method_order) is not int or \
                         self.method_order < 1:
-                    raise ValueError,error_1level
+                    raise ValueError(error_1level)
             self._method_order = self.method_order
 
         else:   # method_order is not specified
@@ -418,14 +419,14 @@ def validate_data(self):
         for i in range(1,array_shape[1] - 1):
             if not np.allclose(self.butcher_tableau[i][0],\
                                sum(self.butcher_tableau[i][1:])):
-                raise ValueError, '''
+                raise ValueError('''
         Inconsistent data in Butcher_Tableau!
         In each lines of stage-coefficients, first number should be
         equal to the sum of other numbers.
         That is, for a butcher_table with K columns,
             a[i][0] == a[i][1] + a[i][2] + ... + a[i][K - 1]
             where 1 <= i <= K - 1
         Your input for line %d is :%s
-        ''' % (i,str(self.butcher_tableau[i]))
+        ''' % (i,str(self.butcher_tableau[i])))
 
         return True

odespy/__init__.py

@@ -1255,17 +1255,17 @@ def f(u, t):
 with :math:`f(u,t)` implemented in Fortran.
 '''
 
-from solvers import *
-from RungeKutta import *
-from rkc import *
-from rkf45 import *
-from odepack import *
-from radau5 import *
-import problems
+from .solvers import *
+from .RungeKutta import *
+from .rkc import *
+from .rkf45 import *
+from .odepack import *
+from .radau5 import *
+from . import problems
 
 # Update doc strings with common info
 class_, doc_str, classname = None, None, None
-classnames = [name for name, obj in locals().items() \
+classnames = [name for name, obj in list(locals().items()) \
                if inspect.isclass(obj)]
 
 toc = []
@@ -1283,7 +1283,7 @@ def f(u, t):
 
 # Do not pollute namespace
 del class_, doc_str, classname, classnames, toc, typeset_toc, \
-    table_of_parameters, name, obj, inspect
+    table_of_parameters, inspect
 
 if __name__ == '__main__':
     from os.path import join

odespy/demos/demo_Lsodi_1.py

@@ -52,7 +52,7 @@ def jac(u, t, s):
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'b-', title="Lsodi with Python functions",
         legend="with res, adda, ydoti & jac", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 # Test case 2: Lsodi, with res, ydoti & adda
 m = method(res=res, rtol=rtol, atol=atol, ydoti=ydoti,
@@ -61,6 +61,6 @@ def jac(u, t, s):
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r*', title="Lsodi with Python functions",
         legend="with res, adda & ydoti", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 

odespy/demos/demo_Lsodi_1_fortran.py

@@ -83,7 +83,7 @@
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r-', title="Lsodi with Fortran subroutines",
         legend="with res, adda, ydoti & jac", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 # Test case 2: Lsodi, with res, ydoti & adda
 m = method(res_f77=res_f77, rtol=rtol, atol=atol, ydoti=ydoti,
@@ -92,6 +92,6 @@
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'g*', title="Lsodi with Fortran subroutines",
         legend="with res, adda & ydoti", hold="on")
-print 'Max error for test case 2 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 2 is %g' % max(u[-1] - exact_final))
 
 os.remove('callback.so')

odespy/demos/demo_Lsodi_2.py

@@ -113,7 +113,7 @@ def jac_banded(u, t, s, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r-', title="Lsodi with Python functions",
         legend="with res, adda_full & jac_full", hold="on")
-print 'Max error with test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error with test case 1 is %g' % max(u[-1] - exact_final))
 
 
 # Test case 2: Lsodi, with res & adda_full
@@ -122,7 +122,7 @@ def jac_banded(u, t, s, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'b*', title="Lsodi with Python functions",
         legend="with res & adda_full", hold="on")
-print 'Max error with test case 2 is %g' % max(u[-1] - exact_final)
+print('Max error with test case 2 is %g' % max(u[-1] - exact_final))
 
 # Test case 3: Lsodi, with res, adda_banded, ml, mu, jac_banded
 m = method(res=res, rtol=rtol, atol=atol, 
@@ -133,7 +133,7 @@ def jac_banded(u, t, s, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'go', title="Lsodi with Python functions",
         legend="with res, adda_banded, jac_banded, ml, mu", hold="on")
-print 'Max error with test case 3 is %g' % max(u[-1] - exact_final)
+print('Max error with test case 3 is %g' % max(u[-1] - exact_final))
 
 # Test case 4: Lsodi, with res, adda_banded, ml, mu
 m = method(res=res, rtol=rtol, atol=atol, 
@@ -143,7 +143,7 @@ def jac_banded(u, t, s, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'y-', title="Lsodi with Python functions",
         legend="with res, adda_banded, ml, mu", hold="on")
-print 'Max error with test case 4 is %g' % max(u[-1] - exact_final)
+print('Max error with test case 4 is %g' % max(u[-1] - exact_final))
 
 
 

odespy/demos/demo_Lsodis_1.py

@@ -73,11 +73,11 @@ def jac(y, t, s, j, ia, ja):
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol, jac_lsodis=jac)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 1 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(y[-1] - exact_final))
 
 # Test case 2:  With res & adda
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol,lrw=4000,liw=100)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 2 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 2 is %g' % max(y[-1] - exact_final))
 

odespy/demos/demo_Lsodis_2.py

@@ -69,25 +69,25 @@ def jac(y, t, s, j, ia, ja):
            ia=ia, ja=ja, ic=ic, jc=jc)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 1 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(y[-1] - exact_final))
 
 # Test case 2: with res, adda, ia, ja & jac
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol, jac_lsodis=jac,
            ia=ia, ja=ja)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 2 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 2 is %g' % max(y[-1] - exact_final))
 
 # Test case 3: with res, adda & jac
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol, jac_lsodis=jac)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 3 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 3 is %g' % max(y[-1] - exact_final))
 
 # Test case 4:  With res & adda
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol,lrw=4000,liw=100)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 4 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 4 is %g' % max(y[-1] - exact_final))
 
 

odespy/demos/demo_Lsoibt_2.py

@@ -209,7 +209,7 @@ def jac(y, t, s):
 u_final = u[-1].reshape(99,3)
 u1, u2, u3 = u_final[:, 0], u_final[:, 1], u_final[:, 2]
 max_error = max(max(u1 - u1_exact), max(u2 - u2_exact), max(u3 - u3_exact))
-print 'Max error with Test case 1 is %g' % max_error
+print('Max error with Test case 1 is %g' % max_error)
 
 # Test case 2: Lsoibt, with res, adda, mb, nb
 m = method(rtol=rtol, atol=atol, res=res, adda_lsoibt=adda, 
@@ -222,4 +222,4 @@ def jac(y, t, s):
 u_final = u[-1].reshape(99,3)
 u1, u2, u3 = u_final[:, 0], u_final[:, 1], u_final[:, 2]
 max_error = max(max(u1 - u1_exact), max(u2 - u2_exact), max(u3 - u3_exact))
-print 'Max error with Test case 2 is %g' % max_error
+print('Max error with Test case 2 is %g' % max_error)

odespy/demos/demo_MyRungeKutta.py

@@ -38,14 +38,14 @@ def f(u,t):
 time_points = np.linspace(t0, tn, n_points)
 
 orders = []  # list for calculated orders
-for m in bt.keys():
+for m in list(bt.keys()):
     # user-defined method, without order suplied 
     method = MyRungeKutta(f, butcher_tableau=bt[m]['table'])
     orders += [method.get_order()]
     method.set_initial_condition(u0)
     u,t = method.solve(time_points)
     error = abs((u[-1] - np.exp(-1.))/np.exp(-1.))
-    print 'Error is %g with solver %s' % (error, m)
+    print('Error is %g with solver %s' % (error, m))