Merge pull request #98 from rzyu45/dev-on-paper

To prepare the release of SolMuseum and standardize namespaces

.github/workflows/run-tests.yml

@@ -24,4 +24,4 @@ jobs:
         pip install -r requirements.txt
         cd ..
     - run: | # run both independent pytest and doctest
-        pytest -vv
+        pytest

Solverz/__init__.py

@@ -1,10 +1,9 @@
-from Solverz.equation.eqn import Eqn, Ode, HyperbolicPde
+from Solverz.equation.eqn import Eqn, Ode
 from Solverz.equation.equations import AE, FDAE, DAE
 from Solverz.equation.param import Param, IdxParam, TimeSeriesParam
 from Solverz.sym_algebra.symbols import idx, Para, iVar, iAliasVar
 from Solverz.sym_algebra.functions import (Sign, Abs, transpose, exp, Diag, Mat_Mul, sin, cos, Min, AntiWindUp,
                                            Saturation, heaviside, ln)
-from Solverz.num_api.custom_function import minmod_flag, minmod
 from Solverz.variable.variables import Vars, TimeVars, as_Vars
 from Solverz.solvers import *
 from Solverz.code_printer import made_numerical, module_printer

Solverz/code_printer/python/inline/inline_printer.py

@@ -4,7 +4,7 @@
 
 from Solverz.variable.variables import Vars
 from Solverz.code_printer.python.utilities import *
-
+from Solverz.num_api.module_parser import modules
 
 # %%
 
@@ -162,12 +162,11 @@ def made_numerical(eqs: SymEquations,
                              sparse)
         code['HVP'] = code_HVP
     custom_func = dict()
-    custom_func.update(numerical_interface)
     custom_func.update(parse_trigger_func(eqs.PARAM))
-    F = Solverzlambdify(code_F, 'F_', modules=[custom_func, 'numpy'])
-    J = Solverzlambdify(code_J, 'J_', modules=[custom_func, 'numpy'])
+    F = Solverzlambdify(code_F, 'F_', modules=[custom_func]+modules)
+    J = Solverzlambdify(code_J, 'J_', modules=[custom_func]+modules)
     if make_hvp:
-        HVP = Solverzlambdify(code_HVP, 'Hvp_', modules=[custom_func, 'numpy'])
+        HVP = Solverzlambdify(code_HVP, 'Hvp_', modules=[custom_func]+modules)
     p = parse_p(eqs.PARAM)
     print('Complete!')
     if isinstance(eqs, SymAE) and not isinstance(eqs, SymFDAE):

Solverz/code_printer/python/inline/tests/test_inline_printer.py

@@ -5,22 +5,22 @@
 from numpy.testing import assert_allclose
 import re
 import pytest
-from sympy import sin
 from sympy.codegen.ast import Assignment, AddAugmentedAssignment
 
-from Solverz import Model, Var, AliasVar
+from Solverz import Model, Var, AliasVar, sin
 from Solverz.equation.eqn import Eqn, Ode
 from Solverz.equation.equations import DAE, AE
 from Solverz.equation.param import Param
 from Solverz.sym_algebra.symbols import iVar, idx, Para
 from Solverz.variable.variables import combine_Vars, as_Vars
-from Solverz.equation.jac import JacBlock, Ones, Jac
+from Solverz.equation.jac import JacBlock, Jac
 from Solverz.equation.hvp import Hvp
-from Solverz.sym_algebra.functions import Diag, sin, cos, exp
+from Solverz.sym_algebra.functions import *
 from Solverz.code_printer.python.inline.inline_printer import print_J_block, extend, SolList, print_J_blocks, print_J, \
     print_F, print_Hvp, made_numerical, Solverzlambdify
 from Solverz.utilities.address import Address
-from Solverz.num_api.custom_function import numerical_interface
+from Solverz.num_api.module_parser import modules
+
 
 # %%
 row = iVar('row', internal_use=True)
@@ -104,8 +104,8 @@ def test_jbs_printer():
     h = y_[0:1]
     v = y_[1:2]
     g = p_["g"]
-    J_ = zeros((2, 2))
-    J_[0:1,1] += ones(1)
+    J_ = np.zeros((2, 2))
+    J_[0:1,1] += np.ones(1)
     return J_
 """.strip()
 
@@ -118,15 +118,15 @@ def test_jbs_printer():
     data = []
     row.extend([0])
     col.extend([1])
-    data.extend(ones(1))
-    return coo_array((data, (row, col)), (2, 2)).tocsc()
+    data.extend(np.ones(1))
+    return sps.coo_array((data, (row, col)), (2, 2)).tocsc()
 """.strip()
 
 expected_F = """def F_(t, y_, p_):
     h = y_[0:1]
     v = y_[1:2]
     g = p_["g"]
-    _F_ = zeros((2, ))
+    _F_ = np.zeros((2, ))
     _F_[0:1] = v
     _F_[1:2] = g
     return _F_
@@ -141,6 +141,7 @@ def test_print_F_J():
     m.f1 = Ode('f1', f=m.v, diff_var=m.h)
     m.f2 = Ode('f2', f=m.g, diff_var=m.v)
     bball, y0 = m.create_instance()
+    bball.FormJac(y0)
     assert print_J(bball.__class__.__name__,
                    bball.jac,
                    bball.a,
@@ -175,11 +176,11 @@ def test_print_F_J():
     q = y_[3:6]
     p_tag_0 = y_0[0:3]
     q_tag_0 = y_0[3:6]
-    J_ = zeros((6, 6))
-    J_[0:2,1:3] += diagflat(ones(2))
-    J_[2:4,0:2] += diagflat(-ones(2))
-    J_[4:5,2] += -ones(1)
-    J_[5:6,2] += ones(1)
+    J_ = np.zeros((6, 6))
+    J_[0:2,1:3] += np.diagflat(np.ones(2))
+    J_[2:4,0:2] += np.diagflat(-np.ones(2))
+    J_[4:5,2] += -np.ones(1)
+    J_[5:6,2] += np.ones(1)
     return J_
 """.strip()
 
@@ -193,25 +194,25 @@ def test_print_F_J():
     data = []
     row.extend([0, 1])
     col.extend([1, 2])
-    data.extend(ones(2))
+    data.extend(np.ones(2))
     row.extend([2, 3])
     col.extend([0, 1])
-    data.extend(-ones(2))
+    data.extend(-np.ones(2))
     row.extend([4])
     col.extend([2])
-    data.extend(-ones(1))
+    data.extend(-np.ones(1))
     row.extend([5])
     col.extend([2])
-    data.extend(ones(1))
-    return coo_array((data, (row, col)), (6, 6)).tocsc()
+    data.extend(np.ones(1))
+    return sps.coo_array((data, (row, col)), (6, 6)).tocsc()
 """.strip()
 
 expected_F_fdae = """def F_(t, y_, p_, y_0):
     p = y_[0:3]
     q = y_[3:6]
     p_tag_0 = y_0[0:3]
     q_tag_0 = y_0[3:6]
-    _F_ = zeros((6, ))
+    _F_ = np.zeros((6, ))
     _F_[0:2] = -p_tag_0[0:2] + p[1:3]
     _F_[2:4] = p_tag_0[1:3] - p[0:2]
     _F_[4:5] = p_tag_0[0] - p[2]
@@ -237,6 +238,7 @@ def test_print_F_J_FDAE():
     m.ae4 = Eqn('cha4',
                 m.p[2] - m.p0[0])
     fdae, y0 = m.create_instance()
+    fdae.FormJac(y0)
     assert print_J(fdae.__class__.__name__,
                    fdae.jac,
                    fdae.a,
@@ -334,7 +336,7 @@ def test_hvp_printer():
                      var_addr,
                      dict())
 
-    HVP = Solverzlambdify(code, 'Hvp_', [numerical_interface, 'numpy'])
+    HVP = Solverzlambdify(code, 'Hvp_', modules)
     np.testing.assert_allclose(HVP(np.array([1, 2]), dict(), np.array([1, 1])).toarray(),
                                np.array([[2.71828183, -0.90929743],
                                          [0.,  2.]]))

Solverz/code_printer/python/module/module_generator.py

@@ -11,7 +11,7 @@
 from Solverz.code_printer.python.module.module_printer import print_F, print_inner_F, print_sub_inner_F, \
     print_J, print_inner_J, print_Hvp, print_inner_Hvp
 from Solverz.equation.equations import Equations as SymEquations
-from Solverz.num_api.custom_function import numerical_interface
+from Solverz.num_api.module_parser import modules
 from Solverz.variable.variables import Vars, combine_Vars
 from Solverz.equation.hvp import Hvp
 
@@ -67,8 +67,7 @@ def render_modules(eqs: SymEquations,
         code_dict["inner_Hvp"] = inner_Hvp['code_inner_Hvp']
         code_dict["sub_inner_Hvp"] = inner_Hvp['code_sub_inner_Hvp']
 
-    custom_func = dict()
-    custom_func.update(numerical_interface)
+
 
     def print_trigger_func_code():
         code_tfuc = dict()
@@ -117,7 +116,7 @@ def print_trigger_func_code():
                       p,
                       eqn_parameter,
                       y,
-                      [custom_func, 'numpy'],
+                      modules,
                       numba,
                       directory)
     print('Complete!')
@@ -265,16 +264,24 @@ def print_module_code(code_dict: Dict[str, str], numba=False):
 
     return code
 
+#
+code_from_SolMuseum="""
+try:
+    import SolMuseum.num_api as SolMF
+except ImportError as e:
+    pass
 
+"""
 def print_dependency_code(modules):
     code = "import os\n"
     code += "current_module_dir = os.path.dirname(os.path.abspath(__file__))\n"
     code += 'from Solverz import load\n'
     code += 'auxiliary = load(f"{current_module_dir}\\\\param_and_setting.pkl")\n'
     code += 'from numpy import *\n'
-    code += 'from numpy import abs\n'
-    code += 'from Solverz.num_api.custom_function import *\n'  # import Solverz built-in func
-    code += 'from scipy.sparse import *\n'
+    code += 'import numpy as np\n'
+    code += 'import Solverz.num_api.custom_function as SolCF\n'  # import Solverz built-in func
+    code += code_from_SolMuseum
+    code += 'import scipy.sparse as sps\n'
     code += 'from numba import njit\n'
     code += 'setting = auxiliary["eqn_param"]\n'
     code += 'row = setting["row"]\n'

Solverz/code_printer/python/module/test/test_module_generator.py

@@ -18,9 +18,15 @@
 from Solverz import load
 auxiliary = load(f"{current_module_dir}\\param_and_setting.pkl")
 from numpy import *
-from numpy import abs
-from Solverz.num_api.custom_function import *
-from scipy.sparse import *
+import numpy as np
+import Solverz.num_api.custom_function as SolCF
+
+try:
+    import SolMuseum.num_api as SolMF
+except ImportError as e:
+    pass
+
+import scipy.sparse as sps
 from numba import njit
 setting = auxiliary["eqn_param"]
 row = setting["row"]
@@ -99,7 +105,7 @@ def test_AE_module_printer():
         inner_F) == '@njit(cache=True)\ndef inner_F(_F_, x):\n    _F_[0:1] = inner_F0(x)\n    _F_[1:2] = inner_F1(x)\n    return _F_\n'
     assert inspect.getsource(inner_F0.func_code) == '@njit(cache=True)\ndef inner_F0(x):\n    return 2*x[0] + x[1]\n'
     assert inspect.getsource(
-        inner_F1.func_code) == '@njit(cache=True)\ndef inner_F1(x):\n    return x[0]**2 + sin(x[1])\n'
+        inner_F1.func_code) == '@njit(cache=True)\ndef inner_F1(x):\n    return x[0]**2 + np.sin(x[1])\n'
     assert inspect.getsource(
         inner_J) == '@njit(cache=True)\ndef inner_J(_data_, x):\n    _data_[2:3] = inner_J0(x)\n    _data_[3:4] = inner_J1(x)\n    return _data_\n'
 
@@ -125,7 +131,7 @@ def test_AE_module_printer():
     assert inspect.getsource(
         inner_F) == 'def inner_F(_F_, x):\n    _F_[0:1] = inner_F0(x)\n    _F_[1:2] = inner_F1(x)\n    return _F_\n'
     assert inspect.getsource(inner_F0) == 'def inner_F0(x):\n    return 2*x[0] + x[1]\n'
-    assert inspect.getsource(inner_F1) == 'def inner_F1(x):\n    return x[0]**2 + sin(x[1])\n'
+    assert inspect.getsource(inner_F1) == 'def inner_F1(x):\n    return x[0]**2 + np.sin(x[1])\n'
     assert inspect.getsource(
         inner_J) == 'def inner_J(_data_, x):\n    _data_[2:3] = inner_J0(x)\n    _data_[3:4] = inner_J1(x)\n    return _data_\n'
 
@@ -136,7 +142,7 @@ def test_AE_module_printer():
     x = y_[0:2]
     c = p_["c"]
     data_hvp = inner_Hvp(_data_hvp, v_, x, c)
-    return coo_array((data_hvp, (row_hvp, col_hvp)), (2, 2)).tocsc()
+    return sps.coo_array((data_hvp, (row_hvp, col_hvp)), (2, 2)).tocsc()
 """
 
 expected_inner_Hvp = """@njit(cache=True)
@@ -149,17 +155,17 @@ def inner_Hvp(_data_hvp, v_, x, c):
 
 expected_inner_Hvp0 = """@njit(cache=True)
 def inner_Hvp0(v_, x):
-    return v_[0]*exp(x[0])*ones(1)
+    return v_[0]*np.exp(x[0])*np.ones(1)
 """
 
 expected_inner_Hvp1 = """@njit(cache=True)
 def inner_Hvp1(v_, x):
-    return -v_[1]*sin(x[1])*ones(1)
+    return -v_[1]*sin(x[1])*np.ones(1)
 """
 
 expected_inner_Hvp2 = """@njit(cache=True)
 def inner_Hvp2(c, v_):
-    return v_[1]*c*ones(1)
+    return v_[1]*c*np.ones(1)
 """
 
 
@@ -236,7 +242,7 @@ def inner_F(_F_, p, q, p_tag_0, q_tag_0, pb, qb):
     pb = p_["pb"].get_v_t(t)
     qb = p_["qb"]
     data = inner_J(_data_, p, q, p_tag_0, q_tag_0, pb, qb)
-    return coo_array((data, (row, col)), (164, 164)).tocsc()
+    return sps.coo_array((data, (row, col)), (164, 164)).tocsc()
 """
 
 expected_inner_J = """@njit(cache=True)
@@ -333,7 +339,7 @@ def inner_F(_F_, h, v):
     h = y_[0:1]
     v = y_[1:2]
     data = inner_J(_data_, h, v)
-    return coo_array((data, (row, col)), (2, 2)).tocsc()
+    return sps.coo_array((data, (row, col)), (2, 2)).tocsc()
 """
 
 expected_inner_J1 = """@njit(cache=True)

Solverz/code_printer/python/module/test/test_module_printer.py

@@ -26,7 +26,7 @@
     lam = p_["lam"]
     ax = ax_trigger_func(x)
     data_hvp = inner_Hvp(_data_hvp, v_, omega, delta, x, y, ax, lam)
-    return coo_array((data_hvp, (row_hvp, col_hvp)), (25, 25)).tocsc()
+    return sps.coo_array((data_hvp, (row_hvp, col_hvp)), (25, 25)).tocsc()
 """.strip()
 expected_FDAE_Hvp = """def Hvp_(t, y_, p_, v_, y_0):
     omega = y_[0:10]
@@ -41,7 +41,7 @@
     lam = p_["lam"]
     ax = ax_trigger_func(x)
     data_hvp = inner_Hvp(_data_hvp, v_, omega, delta, x, y, omega_tag_0, delta_tag_0, x_tag_0, y_tag_0, ax, lam)
-    return coo_array((data_hvp, (row_hvp, col_hvp)), (25, 25)).tocsc()
+    return sps.coo_array((data_hvp, (row_hvp, col_hvp)), (25, 25)).tocsc()
 """.strip()
 
 
@@ -87,13 +87,13 @@ def test_print_Hvp():
     return _data_hvp
 """.strip()
 expected_inner_Hvp0 = """def inner_Hvp0(v_, x):
-    return v_[0]*exp(x[0])*ones(1)
+    return v_[0]*np.exp(x[0])*np.ones(1)
 """.strip()
 expected_inner_Hvp1 = """def inner_Hvp1(v_, x):
-    return -v_[1]*sin(x[1])*ones(1)
+    return -v_[1]*sin(x[1])*np.ones(1)
 """.strip()
 expected_inner_Hvp2 = """def inner_Hvp2(c, v_):
-    return v_[1]*c*ones(1)
+    return v_[1]*c*np.ones(1)
 """.strip()
 
 
@@ -169,7 +169,7 @@ def test_print_inner_Hvp():
     G6 = p_["G6"].get_v_t(t)
     ax = ax_trigger_func(x)
     data = inner_J(_data_, omega, delta, x, y, ax, lam, G6)
-    return coo_array((data, (row, col)), (25, 25)).tocsc()
+    return sps.coo_array((data, (row, col)), (25, 25)).tocsc()
 """.strip()
 expected1 = """def J_(t, y_, p_, y_0):
     omega = y_[0:10]
@@ -185,7 +185,7 @@ def test_print_inner_Hvp():
     G6 = p_["G6"].get_v_t(t)
     ax = ax_trigger_func(x)
     data = inner_J(_data_, omega, delta, x, y, omega_tag_0, delta_tag_0, x_tag_0, y_tag_0, ax, lam, G6)
-    return coo_array((data, (row, col)), (25, 25)).tocsc()
+    return sps.coo_array((data, (row, col)), (25, 25)).tocsc()
 """.strip()
 
 
@@ -234,7 +234,7 @@ def test_print_J():
     return _data_
 """.strip()
 expected3 = """def inner_J0(y):
-    return y*ones(3)
+    return y*np.ones(3)
 """.strip()
 expected4 = """def inner_J1(y):
     return y**2