Merge pull request #30 from modelon-community/dev-PM-4738

Added C version of Radau5 solver
Improved tests to use nose.tools_assert_x

CHANGELOG

@@ -2,6 +2,8 @@
 
 --- Assimulo-3.2.8 ---
     * Sundials 5.x port
+    * Improved asserts in tests
+    * Added C version of the Radau5 solver
 
 --- Assimulo-3.2.7 ---
     * Resolved deprecation warnings visible when using numpy 1.20.

doc/sphinx/source/changelog.rst

@@ -3,6 +3,11 @@
 Changelog
 ==========
 
+--- Assimulo-3.2.8 ---
+    * Sundials 5.x port
+    * Improved asserts in tests
+    * Added C version of the Radau5 solver
+
 --- Assimulo-3.2.7 ---
     * Resolved deprecation warnings visible when using numpy 1.20, related to deprecation of the alias numpy.float.
     * Resolved deprecation warnings visible when creating an ndarray from ragged nested sequences.

doc/sphinx/source/conf.py

@@ -50,9 +50,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = '3.2.7'
+version = '3.2.8'
 # The full version, including alpha/beta/rc tags.
-release = '3.2.7'
+release = '3.2.8'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -109,7 +109,7 @@
 
 # The name for this set of Sphinx documents.  If None, it defaults to
 # "<project> v<release> documentation".
-html_title = 'Assimulo 3.2.7 documentation'
+html_title = 'Assimulo 3.2.8 documentation'
 
 # A shorter title for the navigation bar.  Default is the same as html_title.
 #html_short_title = None

examples/cvode_stability.py

@@ -89,7 +89,7 @@ def f(t,y):
 
     #Basic test
     x1 = y[:,0]
-    assert N.abs(x1[-1]-1.8601438) < 1e-1 #For test purpose
+    nose.tools.assert_less(N.abs(float(x1[-1]) - 1.8601438), 1e-1)
 
     return exp_mod, exp_sim
 

examples/euler_vanderpol.py

@@ -86,7 +86,7 @@ def jac(t,y):
 
     #Basic test
     x1 = y[:,0]
-    assert N.abs(x1[-1]-1.8601438) < 1e-1 #For test purpose
+    nose.tools.assert_less(N.abs(float(x1[-1]) - 1.8601438), 1e-1)
 
     return exp_mod, exp_sim
 

examples/glimda_vanderpol.py

@@ -86,7 +86,7 @@ def f(t,y,yd):
 
     #Basic test
     x1 = y[:,0]
-    assert N.abs(x1[-1]-1.706168035) < 1e-3 #For test purpose
+    nose.tools.assert_almost_equal(float(x1[-1]), 1.706168035, 3)
 
     return imp_mod, imp_sim
 

examples/lsodar_vanderpol.py

@@ -75,7 +75,7 @@ def f(t,y):
 
     #Basic test
     x1 = y[:,0]
-    assert N.abs(x1[-1]-1.706168035) < 1e-3 #For test purpose
+    nose.tools.assert_less(N.abs(x1[-1] - 1.706168035), 1e-3)
 
     return exp_mod, exp_sim
 

examples/lsodar_with_disc.py

@@ -67,7 +67,6 @@ def state_events(self,t,y,sw):
 
         return N.array([event_0,event_1,event_2])
 
-
     #Responsible for handling the events.
     def handle_event(self, solver, event_info):
         """
@@ -121,8 +120,6 @@ def init_mode(self, solver):
         solver.y[1] = (-1.0 if solver.sw[1] else 3.0)
         solver.y[2] = (0.0 if solver.sw[2] else 2.0)
 
-
-
 def run_example(with_plots=True):
     r"""
     Example of the use of Euler's method for a differential equation
@@ -154,15 +151,12 @@ def run_example(with_plots=True):
         P.xlabel('Time')
         P.show()
 
-    return exp_mod, exp_sim
-
     #Basic test
     nose.tools.assert_almost_equal(y[-1][0],8.0)
     nose.tools.assert_almost_equal(y[-1][1],3.0)
     nose.tools.assert_almost_equal(y[-1][2],2.0)
+
+    return exp_mod, exp_sim
 
 if __name__=="__main__":
     mod,sim = run_example()
-
-
-

examples/mech_system_pendulum.py

@@ -19,6 +19,7 @@
 import assimulo.problem as ap
 import assimulo.special_systems as ass
 import numpy as N
+import nose
 import scipy.linalg as sl
 from assimulo.solvers import IDA, ODASSL
 
@@ -63,7 +64,7 @@ def run_example(index="ind1", with_plots=True, with_test=False):
     print(final_residual, 'Norm:  ', sl.norm(final_residual))
 
     if with_test:
-        assert(sl.norm(final_residual) < 1.5e-1)
+        nose.tools.assert_less(sl.norm(final_residual), 1.5e-1)
     if with_plots:
         dae_pend.plot(mask=[1,1]+(len(my_pend.y0)-2)*[0]) 
     return my_pend, dae_pend
@@ -75,9 +76,3 @@ def run_example(index="ind1", with_plots=True, with_test=False):
     mod={}
     for ind in index_values:
         mod[ind], sim[ind]=run_example(index=ind)
-
-
-
-
-
-

examples/radau5dae_time_events.py

@@ -49,7 +49,7 @@ def time_events(self, t,y,yd,sw):
     def handle_event(self, solver, event_info):
         self.my *= 1e-1
 
-def run_example(with_plots=True):
+def run_example(with_plots=True, solver='c'):
 
     y0 = [2.0,-0.6] #Initial conditions
     yd0 = [-.6,-200000.]
@@ -59,6 +59,7 @@ def run_example(with_plots=True):
 
     #Define an explicit solver
     imp_sim = Radau5DAE(imp_mod) #Create a Radau5 solver
+    imp_sim.solver = solver
 
     #Simulate
     t, y, yd = imp_sim.simulate(8.) #Simulate 8 seconds
@@ -74,8 +75,8 @@ def run_example(with_plots=True):
 
     #Basic test
     x1 = y[:,0]
-    assert N.abs(x1[-1]-1.14330840983) < 1e-3 #For test purpose
+    nose.tools.assert_less(N.abs(float(x1[-1]) - 1.14330840983), 1e-3)
+
     return imp_mod, imp_sim
 if __name__=='__main__':
     mod,sim = run_example()
-

examples/radau5dae_vanderpol.py

@@ -20,7 +20,7 @@
 from assimulo.solvers import Radau5DAE
 from assimulo.problem import Implicit_Problem
 
-def run_example(with_plots=True):
+def run_example(with_plots=True,solver='c'):
     r"""
     Example for the use of Radau5DAE to solve
     Van der Pol's equation
@@ -61,6 +61,7 @@ def f(t,y,yd):
 
     #Define an explicit solver
     imp_sim = Radau5DAE(imp_mod) #Create a Radau5 solver
+    imp_sim.solver = solver
 
     #Sets the parameters
     imp_sim.atol = 1e-4 #Default 1e-6
@@ -86,9 +87,9 @@ def f(t,y,yd):
 
     #Basic test
     x1 = y[:,0]
-    assert N.abs(x1[-1]-1.706168035) < 1e-3 #For test purpose
+    nose.tools.assert_less(N.abs(float(x1[-1]) - 1.706168035), 1e-3)
+
     return imp_mod, imp_sim
 
 if __name__=='__main__':
     mod,sim = run_example()
-

examples/radau5ode_vanderpol.py

@@ -20,7 +20,7 @@
 from assimulo.solvers import Radau5ODE
 from assimulo.problem import Explicit_Problem
 
-def run_example(with_plots=True):
+def run_example(with_plots=True,solver='c'):
     r"""
     Example for the use of the implicit Euler method to solve
     Van der Pol's equation
@@ -56,6 +56,7 @@ def f(t,y):
 
     #Define an explicit solver
     exp_sim = Radau5ODE(exp_mod) #Create a Radau5 solver
+    exp_sim.solver = solver
 
     #Sets the parameters
     exp_sim.atol = 1e-4 #Default 1e-6
@@ -76,10 +77,9 @@ def f(t,y):
 
     #Basic test
     x1 = y[:,0]
-    assert N.abs(x1[-1]-1.706168035) < 1e-3 #For test purpose
+    nose.tools.assert_less(N.abs(float(x1[-1]) - 1.706168035), 1e-3)
 
     return exp_mod, exp_sim
 
 if __name__=='__main__':
     mod,sim = run_example()
-

examples/radau5ode_with_disc.py

@@ -123,11 +123,12 @@ def init_mode(self, solver):
 
 
 
-def run_example(with_plots=True):
+def run_example(with_plots=True,solver='c'):
     #Create an instance of the problem
     exp_mod = Extended_Problem() #Create the problem
 
     exp_sim = Radau5ODE(exp_mod) #Create the solver
+    exp_sim.solver = solver
 
     exp_sim.verbosity = 0
     exp_sim.report_continuously = True
@@ -153,6 +154,4 @@ def run_example(with_plots=True):
 
 if __name__=="__main__":
     mod,sim = run_example()
-
 
-

examples/rodasode_vanderpol.py

@@ -86,7 +86,8 @@ def jac(t,y):
 
     #Basic test
     x1 = y[:,0]
-    assert N.abs(x1[-1]-1.706168035) < 1e-3 #For test purpose
+    nose.tools.assert_less(N.abs(float(x1[-1]) - 1.706168035), 1e-3)
+
     return exp_mod, exp_sim
 
 

examples/rungekutta34_with_disc.py

@@ -153,6 +153,3 @@ def run_example(with_plots=True):
 
 if __name__=="__main__":
     mod,sim = run_example()
-
-
-

setup.py

@@ -209,7 +209,6 @@ def _set_directories(self):
         self.desThirdParty=dict([(thp,os.path.join(self.curdir,self.build_assimulo_thirdparty,thp)) 
                                           for thp in self.thirdparty_methods])
         # filelists
-
         self.fileSrc     = os.listdir("src")
         self.fileLib     = os.listdir(os.path.join("src","lib"))
         self.fileSolvers = os.listdir(os.path.join("src","solvers"))
@@ -459,12 +458,13 @@ def check_LAPACK(self):
 
     def cython_extensionlists(self):
         extra_link_flags = self.static_link_gcc + self.flag_32bit
-
+        
         #Cythonize main modules
         ext_list = cythonize(["assimulo"+os.path.sep+"*.pyx"], 
                              include_path=[".","assimulo"])
-         #Cythonize Solvers
-         # Euler
+
+        #Cythonize Solvers
+        # Euler
         ext_list += cythonize(["assimulo"+os.path.sep+"solvers"+os.path.sep+"euler.pyx"], 
                              include_path=[".","assimulo"])
         for el in ext_list:
@@ -493,7 +493,6 @@ def cython_extensionlists(self):
                 ext_list[-1].include_dirs.append(self.SLUincdir)
                 ext_list[-1].library_dirs.append(self.SLUlibdir)
                 ext_list[-1].libraries.extend(self.superLUFiles)
-
 
             #Kinsol
             ext_list += cythonize(["assimulo"+os.path.sep+"solvers"+os.path.sep+"kinsol.pyx"], 
@@ -507,6 +506,16 @@ def cython_extensionlists(self):
                 ext_list[-1].include_dirs.append(self.SLUincdir)
                 ext_list[-1].library_dirs.append(self.SLUlibdir)
                 ext_list[-1].libraries.extend(self.superLUFiles)
+
+        ## Radau
+        ext_list += cythonize([os.path.join("assimulo","thirdparty","hairer","radau5_c_py.pyx")],
+                              include_path=[".", "assimulo", os.path.join("assimulo", "lib")],
+                              force = True)
+        ext_list[-1].include_dirs = [np.get_include(), "assimulo", os.path.join("assimulo", "lib"),
+                                     os.path.join("assimulo","thirdparty","hairer"),
+                                     self.incdirs]
+        ext_list[-1].sources = ext_list[-1].sources + [os.path.join("assimulo","thirdparty","hairer","radau_decsol_c.c")]
+        ext_list[-1].name = "assimulo.lib.radau5_c_py"
 
         for el in ext_list:
             #Debug
@@ -597,10 +606,8 @@ def fortran_extensionlists(self):
         else:
             L.warning("Could not find Blas or Lapack, disabling support for the solver GLIMDA.")
 
-
         return config.todict()["ext_modules"]
-
-
+
 prepare=Assimulo_prepare(args, thirdparty_methods)
 curr_dir=os.getcwd()
 if not os.path.isdir("assimulo"):

src/lib/radau_core.py

@@ -434,3 +434,41 @@ def _set_maxsteps(self, max_steps):
         self.options["maxsteps"] = max_steps
 
     maxsteps = property(_get_maxsteps, _set_maxsteps)
+
+    def _get_solver(self):
+        """
+        Solver implementation used, "f" for Fortran, "c" for C
+        
+            Parameters::
+            
+                solver
+                            - Default "f"
+                            
+                            - needs to be either "f" (Fotran) or "c" (C)
+        """
+        return self.options["solver"]
+
+    def _set_solver(self, solver):
+        try:
+            solver.lower()
+        except:
+            raise Radau_Exception("'solver' parameters needs to be the STRING 'c' or 'f'. Set value: {}, type: {}".format(solver, type(solver)))
+        if solver.lower() == "f": ## Fortran
+            try:
+                from assimulo.lib import radau5 as radau5_f
+                self.radau5 = radau5_f
+                self.solver_module_imported = True
+            except:
+                raise Radau_Exception("Failed to import the Fotran based Radau5 solver. Try using solver = 'c' for the C based solver instead.")
+        elif solver.lower() == "c":
+            try:
+                from assimulo.lib import radau5_c_py as radau5_c
+                self.radau5 = radau5_c
+                self.solver_module_imported = True
+            except:
+                raise Radau_Exception("Failed to import the C based Radau5 solver. Try using solver = 'f' for the Fortran based solver instead.")
+        else:
+            raise Radau_Exception("Solver parameters needs to be either 'f' or 'c'. Set value: {}".format(solver))
+        self.options["solver"] = solver.lower()
+
+    solver = property(_get_solver, _set_solver)