Feature: Symbolic outputs of tendencies

documentation/source/files/technical_description.rst

@@ -94,6 +94,20 @@ The computational flow is as follows:
 
     Sketch of the computational flow.
 
+
+Generating Symbolic Equations
+-----------------------------
+Using `Sympy`_ qgs offers the functionality to return the ordinary differential equations of the projected model as a string, with any parameters the user chooses to be returned as varibales in the equations. In addition, the resulting equations can be returned already formatted in the programming language of the users' choice. This allows the qgs framework to feed directly into pipelines in other programming languages. Currently the framework can return the model equations in the following languages:
+
+* Python
+* Julia
+* Fortran-90
+* Mathematica
+* AUTO-p07 continuation software
+
+This also allows the user to specify their own integration method for solving the model equations in python.
+
+
 Additional technical information
 --------------------------------
 
@@ -103,6 +117,8 @@ Additional technical information
 * qgs has a `tangent linear model`_ optimized to run ensembles of initial conditions as well, with a broadcast
   integration of the tangent model thanks to `Numpy`_.
 
+* The symbolic output functionality of the qgs model relies on `Sympy`_ to perform the tensor calculations. This library is significantly slower than the numerical equivilent and as a result it is currently only feasible to generate the symbolic model equations for model resolutions up to :math:`4x4`. 
+
 References
 ----------
 

documentation/source/files/user_guide.rst

@@ -462,7 +462,33 @@ This concludes the initialization of qgs, the function :meth:`f` hence produced
 An example of the construction exposed here, along with plots of the trajectories generated, can be found in the section :ref:`files/examples/manual:Manual setting of the basis`.
 See also the following section for the possible usages.
 
-4. Using qgs (once configured)
+
+4. Symbolic outputs 
+-------------------
+
+If the user wants to generate the model tendencies with non-fixed parameters, use the tendencies in another programming language, or use their own integrator in python, the qgs framework can use `Sympy`_ to proform the above calculations symbolically rather than numerically.  
+
+To return the tendencies including specified parameter values, or to format the tendencies in another programming language, the qgs model is configured as shown in section :ref:`files/user_guide:Configuration of qgs`, however the **Symbolic** inner product is always used in this pipeline. To create the symbolic tendencies, the instance of :class:`.QgParams` is passed to the function :func:`.create_symbolic_tendencies`:
+
+.. code:: ipython3
+
+    from qgs.functions.symbolic_tendencies import create_symbolic_equations
+
+    parameters = [model_parameters.par1, model_parameters.par2, model_parameters.par3]
+
+    f = create_symbolic_equations(model_parameters, continuation_variables=parameters, language='python')
+
+The varibale :meth:`f` is a string of the model tendencies, formatted in the programming language specified by the keyword :meth:`language`. The model tendencies will contain the specified :meth:`continuation_variables` as free parameters. 
+
+Currently the framework can format the equations in the following programming languages:
+* :meth:`python`
+* :meth:`julia`
+* :meth:`fortran`
+* :meth:`mathematica`
+* :meth:`auto`
+
+
+1. Using qgs (once configured)
 ---------------------------------
 
 Once the function :math:`\boldsymbol{f}` giving the model's tendencies has been obtained, it is possible to use it with
@@ -482,7 +508,7 @@ the qgs built-in integrator to obtain the model's trajectories:
 
 Note that it is also possible to use other ordinary differential equations integrators available on the market, see for instance the :ref:`files/examples/diffeq_example:Example of DiffEqPy usage`.
 
-4.1 Analyzing the model's output with diagnostics
+5.1 Analyzing the model's output with diagnostics
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 .. include:: diagnostic.rst
@@ -524,7 +550,7 @@ Note that it is also possible to use other ordinary differential equations integ
 
 More diagnostics will be implemented soon.
 
-4.2 Toolbox
+5.2 Toolbox
 ^^^^^^^^^^^
 
 The toolbox regroups submodules to make specific analysis with the model and are available in the :mod:`qgs.toolbox` module.
@@ -536,7 +562,7 @@ Presently, the list of submodules available is the following:
 
 More submodules will be implemented soon.
 
-4.2.1 Lyapunov toolbox
+5.2.1 Lyapunov toolbox
 """"""""""""""""""""""
 
 This module allows you to integrate the model and simultaneously obtain the *local* `Lyapunov vectors`_ and `exponents`_.
@@ -555,10 +581,10 @@ See also :cite:`user-KP2012` for a description of these methods.
 
 Some example notebooks on how to use this module are available in the `notebooks/lyapunov <../../../../notebooks/lyapunov>`_ folder.
 
-5. Developers information
+6. Developers information
 -------------------------
 
-5.1 Running the test
+6.1 Running the test
 ^^^^^^^^^^^^^^^^^^^^^
 
 The model core tensors can be tested by running `pytest`_ in the main folder: ::

model_test/test_aotensor_sym.py

@@ -0,0 +1,103 @@
+
+import sys
+import os
+
+path = os.path.abspath('./')
+base = os.path.basename(path)
+if base == 'model_test':
+    sys.path.extend([os.path.abspath('../')])
+else:
+    sys.path.extend([path])
+
+
+import unittest
+import numpy as np
+
+from qgs.params.params import QgParams
+from qgs.inner_products import symbolic
+from qgs.tensors.qgtensor import QgsTensor
+from qgs.tensors.symbolic_qgtensor import SymbolicQgsTensor
+
+from model_test.test_base_symbolic import TestBaseSymbolic
+
+real_eps = np.finfo(np.float64).eps
+
+
+class TestSymbolicAOTensor(TestBaseSymbolic):
+    '''
+        Test class for the Linear Symbolic Tensor
+        The tensor is tested against the reference file, and then the numerical tensor calculated in qgs.
+    '''
+
+    filename = 'test_aotensor.ref'
+
+    def test_sym_against_ref(self):
+        self.check_lists_flt()
+
+    def test_sym_against_num(self):
+        self.check_numerical_lists_flt()
+
+    def symbolic_outputs(self, output_func=None):
+
+        if output_func is None:
+            self.symbolic_values.clear()
+            tfunc = self.save_ip_symbolic
+        else:
+            tfunc = output_func
+
+        params = QgParams({'rr': 287.e0, 'sb': 5.6e-8})
+        params.set_atmospheric_channel_fourier_modes(2, 2, mode="symbolic")
+        params.set_oceanic_basin_fourier_modes(2, 4, mode="symbolic")
+
+        # Setting MAOOAM default parameters
+        params.set_params({'kd': 0.04, 'kdp': 0.04, 'n': 1.5})
+
+        aip = symbolic.AtmosphericSymbolicInnerProducts(params, return_symbolic=True, make_substitution=True)
+        oip = symbolic.OceanicSymbolicInnerProducts(params, return_symbolic=True, make_substitution=True)
+
+        aip.connect_to_ocean(oip)
+
+        sym_aotensor = SymbolicQgsTensor(params=params, atmospheric_inner_products=aip, oceanic_inner_products=oip)
+
+        subbed_tensor = sym_aotensor.sub_tensor()
+
+        for coo, val in zip(subbed_tensor.keys(), subbed_tensor.values()):
+            _ip_string_format(tfunc, 'sym_aotensor', coo, val)
+
+    def numerical_outputs(self, output_func=None):
+
+        if output_func is None:
+            self.numerical_values.clear()
+            tfunc = self.save_ip_numeric
+        else:
+            tfunc = output_func
+
+        params = QgParams({'rr': 287.e0, 'sb': 5.6e-8})
+        params.set_atmospheric_channel_fourier_modes(2, 2, mode="symbolic")
+        params.set_oceanic_basin_fourier_modes(2, 4, mode="symbolic")
+
+        # Setting MAOOAM default parameters
+        params.set_params({'kd': 0.04, 'kdp': 0.04, 'n': 1.5})
+
+        aip = symbolic.AtmosphericSymbolicInnerProducts(params, return_symbolic=False)
+        oip = symbolic.OceanicSymbolicInnerProducts(params, return_symbolic=False)
+
+        aip.connect_to_ocean(oip)
+
+        num_aotensor = QgsTensor(params=params, atmospheric_inner_products=aip, oceanic_inner_products=oip)
+
+        for coo, val in zip(num_aotensor.tensor.coords.T, num_aotensor.tensor.data):
+            _ip_string_format(tfunc, 'num_aotensor', coo, val)
+
+
+def _ip_string_format(func, symbol, indices, value):
+    if abs(value) >= real_eps:
+        s = symbol
+        for i in indices:
+            s += "["+str(i)+"]"
+        s += " = % .5E" % value
+        func(s)
+
+
+if __name__ == "__main__":
+    unittest.main()

model_test/test_aotensor_sym_dynT.py

@@ -0,0 +1,97 @@
+
+import sys
+import os
+
+path = os.path.abspath('./')
+base = os.path.basename(path)
+if base == 'model_test':
+    sys.path.extend([os.path.abspath('../')])
+else:
+    sys.path.extend([path])
+
+
+import unittest
+import numpy as np
+
+from qgs.params.params import QgParams
+from qgs.inner_products import symbolic
+from qgs.tensors.qgtensor import QgsTensorDynamicT
+from qgs.tensors.symbolic_qgtensor import SymbolicQgsTensorDynamicT
+
+from model_test.test_base_symbolic import TestBaseSymbolic
+
+real_eps = np.finfo(np.float64).eps
+
+class TestSymbolicAOTensorDynT(TestBaseSymbolic):
+    '''
+        Test class for the Dynamic T Symbolic Tensor
+        The tensor is tested against the reference file, and then the numerical tensor calculated in qgs.
+    '''
+
+    def test_sym_against_num(self):
+        self.check_numerical_lists_flt()
+
+    def symbolic_outputs(self, output_func=None):
+
+        if output_func is None:
+            self.symbolic_values.clear()
+            tfunc = self.save_ip_symbolic
+        else:
+            tfunc = output_func
+
+        params = QgParams({'rr': 287.e0, 'sb': 5.6e-8}, dynamic_T=True)
+        params.set_atmospheric_channel_fourier_modes(2, 2, mode="symbolic")
+        params.set_oceanic_basin_fourier_modes(2, 4, mode="symbolic")
+
+        # Setting MAOOAM default parameters
+        params.set_params({'kd': 0.04, 'kdp': 0.04, 'n': 1.5})
+
+        aip = symbolic.AtmosphericSymbolicInnerProducts(params, return_symbolic=True, make_substitution=True)
+        oip = symbolic.OceanicSymbolicInnerProducts(params, return_symbolic=True, make_substitution=True)
+
+        if not aip.connected_to_ocean:
+            aip.connect_to_ocean(oip)
+
+        sym_aotensor = SymbolicQgsTensorDynamicT(params=params, atmospheric_inner_products=aip, oceanic_inner_products=oip)
+
+        subbed_tensor = sym_aotensor.sub_tensor()
+
+        for coo, val in zip(subbed_tensor.keys(), subbed_tensor.values()):
+            _ip_string_format(tfunc, 'sym_aotensor', coo, val)
+
+    def numerical_outputs(self, output_func=None):
+
+        if output_func is None:
+            self.numerical_values.clear()
+            tfunc = self.save_ip_numeric
+        else:
+            tfunc = output_func
+
+        params = QgParams({'rr': 287.e0, 'sb': 5.6e-8}, dynamic_T=True)
+        params.set_atmospheric_channel_fourier_modes(2, 2, mode="symbolic")
+        params.set_oceanic_basin_fourier_modes(2, 4, mode="symbolic")
+
+        # Setting MAOOAM default parameters
+        params.set_params({'kd': 0.04, 'kdp': 0.04, 'n': 1.5})
+
+        aip = symbolic.AtmosphericSymbolicInnerProducts(params, return_symbolic=False)
+        oip = symbolic.OceanicSymbolicInnerProducts(params, return_symbolic=False)
+
+        if not aip.connected_to_ocean:
+            aip.connect_to_ocean(oip)
+
+        num_aotensor = QgsTensorDynamicT(params=params, atmospheric_inner_products=aip, oceanic_inner_products=oip)
+
+        for coo, val in zip(num_aotensor.tensor.coords.T, num_aotensor.tensor.data):
+            _ip_string_format(tfunc, 'num_aotensor', coo, val)
+
+def _ip_string_format(func, symbol, indices, value):
+    if abs(value) >= real_eps:
+        s = symbol
+        for i in indices:
+            s += "["+str(i)+"]"
+        s += " = % .5E" % value
+        func(s)
+
+if __name__ == "__main__":
+    unittest.main()