diff --git a/README.md b/README.md index 60b9746..f7907cd 100644 --- a/README.md +++ b/README.md @@ -115,7 +115,7 @@ Note that the Jupyter notebooks are not concerned by this recommendation and wor #### Activating DifferentialEquations.jl optional support In addition to the qgs builtin Runge-Kutta integrator, the qgs model can alternatively be integrated with a package called [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) written in [Julia](https://julialang.org/), and available through the -[diffeqpy](https://github.com/SciML/diffeqpy) python package. +[diffeqpy](https://github.com/SciML/diffeqpy) Python package. The diffeqpy package first installation step is done by Anaconda in the qgs environment but then you must [install Julia](https://julialang.org/downloads/) and follow the final manual installation instruction found in the [diffeqpy README](https://github.com/SciML/diffeqpy). These can be summed up as opening a terminal and doing: @@ -208,7 +208,7 @@ Non-exhaustive list: * [Q-GCM](http://q-gcm.org/): A mid-latitude grid based ocean-atmosphere model like MAOOAM. Code in Fortran, interface is in Python. -* [pyqg](https://github.com/pyqg/pyqg): A pseudo-spectral python solver for quasi-geostrophic systems. +* [pyqg](https://github.com/pyqg/pyqg): A pseudo-spectral Python solver for quasi-geostrophic systems. * [Isca](https://execlim.github.io/IscaWebsite/index.html): Research GCM written in Fortran and largely configured with Python scripts, with internal coding changes required for non-standard cases. diff --git a/documentation/source/files/general_information.rst b/documentation/source/files/general_information.rst index 2626422..be6a53d 100644 --- a/documentation/source/files/general_information.rst +++ b/documentation/source/files/general_information.rst @@ -161,7 +161,7 @@ A review of your pull request will follow with possibly suggestions of changes b Please consider the following guidelines before submitting: * Before submitting a pull request, double check that the branch to be merged contains only changes you wish to add to the master branch. This will save time in reviewing the code. -* For any changes to the core model files, please check your submission by running the tests found in the folder `model_test <../../../../model_test>`_ to ensure that the model tensors are still valid (see the section :ref:`files/user_guide:5. Developers information` of the :ref:`files/user_guide:User guide`). Please do not make changes to existing test cases. +* For any changes to the core model files, please check your submission by running the tests found in the folder `model_test <../../../../model_test>`_ to ensure that the model tensors are still valid (see the section :ref:`files/user_guide:6. Developers information` of the :ref:`files/user_guide:User guide`). Please do not make changes to existing test cases. * For substantial additions of code, including a test case (using `unittest`_) in the folder `model_test <../../../../model_test>`_ is recommended. * Please document the new functionalities in the documentation. Code addition without documentation addition will not be accepted. The documentation is done with `sphinx`_ and follows the Numpy conventions. Please take a look to the actual code to get an idea about how to document the code. * If your addition can be considered as a tool not directly related to the core of the model, please develop it in the toolbox folder. @@ -182,7 +182,7 @@ Non-exhaustive list: * `Q-GCM `_: A mid-latitude grid based ocean-atmosphere model like MAOOAM. Code in Fortran, interface is in Python. -* `pyqg `_: A pseudo-spectral python solver for quasi-geostrophic systems. +* `pyqg `_: A pseudo-spectral Python solver for quasi-geostrophic systems. * `Isca `_: Research GCM written in Fortran and largely configured with Python scripts, with internal coding changes required for non-standard cases. diff --git a/documentation/source/files/technical/functions.rst b/documentation/source/files/technical/functions.rst index fed6003..065b867 100644 --- a/documentation/source/files/technical/functions.rst +++ b/documentation/source/files/technical/functions.rst @@ -10,5 +10,11 @@ The functions module contains utility functions used by the model, as well as so .. automodule:: qgs.functions.tendencies :members: +.. automodule:: qgs.functions.symbolic_mul + :members: + +.. automodule:: qgs.functions.symbolic_tendencies + :members: + .. automodule:: qgs.functions.util :members: diff --git a/documentation/source/files/technical/tensors.rst b/documentation/source/files/technical/tensors.rst index 3a20dee..486e49c 100644 --- a/documentation/source/files/technical/tensors.rst +++ b/documentation/source/files/technical/tensors.rst @@ -11,3 +11,7 @@ Module holding the model's tendencies tensor encoding each of their additive ter .. automodule:: qgs.tensors.atmo_thermo_tensor :show-inheritance: :members: + +.. automodule:: qgs.tensors.symbolic_qgtensor + :show-inheritance: + :members: diff --git a/documentation/source/files/technical_description.rst b/documentation/source/files/technical_description.rst index 7270c02..d403ea6 100644 --- a/documentation/source/files/technical_description.rst +++ b/documentation/source/files/technical_description.rst @@ -101,11 +101,11 @@ Using `Sympy`_ qgs offers the functionality to return the ordinary differential * Python * Julia -* Fortran-90 -* Mathematica -* AUTO-p07 continuation software +* Fortran 90 +* `AUTO-p07 `_ continuation software +* `Mathematica `_ (still being tested) -This also allows the user to specify their own integration method for solving the model equations in python. +This also allows the user to specify their own integration method for solving the model equations in Python. Additional technical information @@ -117,7 +117,9 @@ 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`. +* The symbolic output functionality of the qgs model relies on `Sympy`_ to perform the tensor calculations. + This library is significantly slower than the numerical equivalent and as a result it is currently only feasible to + generate the symbolic model equations for model resolutions up to :math:`4x4`. References ---------- @@ -131,3 +133,4 @@ References .. _multiprocessing: https://docs.python.org/3.7/library/multiprocessing.html#module-multiprocessing .. _tangent linear model: http://glossary.ametsoc.org/wiki/Tangent_linear_model .. _ordinary differential equations: https://en.wikipedia.org/wiki/Ordinary_differential_equation +.. _Sympy: https://www.sympy.org/ diff --git a/documentation/source/files/user_guide.rst b/documentation/source/files/user_guide.rst index f73e18a..64c3e7d 100644 --- a/documentation/source/files/user_guide.rst +++ b/documentation/source/files/user_guide.rst @@ -18,7 +18,7 @@ This guide provides: 1. The different ways to configure qgs in order to obtain the function :math:`\boldsymbol{f}` are explained in the section :ref:`files/user_guide:2. Configuration of qgs`. 2. Some specific ways to configure qgs detailed in the section :ref:`files/user_guide:3. Using user-defined symbolic basis`. -3. Examples of usages of the model's tendencies function :math:`\boldsymbol{f}` are given in the section :ref:`files/user_guide:4. Using qgs (once configured)`. +3. Examples of usages of the model's tendencies function :math:`\boldsymbol{f}` are given in the section :ref:`files/user_guide:5. Using qgs (once configured)`. 2. Configuration of qgs ----------------------- @@ -227,7 +227,7 @@ Jacobian matrix :math:`\boldsymbol{\mathrm{Df}}`. Just pass it to the function : f, Df = create_tendencies(model_parameters) The function :meth:`f` hence produced can be used to generate the model's trajectories. -See the section :ref:`files/user_guide:4. Using qgs (once configured)` for the possible usages. +See the section :ref:`files/user_guide:5. Using qgs (once configured)` for the possible usages. 2.5 Saving your model ^^^^^^^^^^^^^^^^^^^^^^^ @@ -466,29 +466,37 @@ See also the following section for the possible usages. 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. +If the user wants to generate the model tendencies with non-fixed parameters, or use the tendencies in another programming language, +the qgs framework can use `Sympy`_ to perform 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`: +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:2. 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 + from qgs.functions.symbolic_tendencies import create_symbolic_tendencies parameters = [model_parameters.par1, model_parameters.par2, model_parameters.par3] - f = create_symbolic_equations(model_parameters, continuation_variables=parameters, language='python') + f = create_symbolic_tendencies(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. +The variable :code:`f` is a string of the model tendencies, formatted in the programming language specified by the +keyword :code:`language`. The model tendencies will contain the specified :code:`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` +* :code:`python` +* :code:`julia` +* :code:`fortran` +* :code:`auto` +* :code:`mathematica` (still under test) -1. Using qgs (once configured) +Some example on how to use this functionality are provided in the +`notebooks/symbolic_outputs <../../../../notebooks/symbolic_outputs>`_ folder. + +5. 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 diff --git a/model_test/test_aotensor_sym_dynT.py b/model_test/test_aotensor_sym_dynT.py index 852b2a6..2f7904c 100644 --- a/model_test/test_aotensor_sym_dynT.py +++ b/model_test/test_aotensor_sym_dynT.py @@ -22,6 +22,7 @@ real_eps = np.finfo(np.float64).eps + class TestSymbolicAOTensorDynT(TestBaseSymbolic): ''' Test class for the Dynamic T Symbolic Tensor @@ -85,6 +86,7 @@ def numerical_outputs(self, output_func=None): 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 @@ -93,5 +95,6 @@ def _ip_string_format(func, symbol, indices, value): s += " = % .5E" % value func(s) + if __name__ == "__main__": unittest.main() diff --git a/model_test/test_base_symbolic.py b/model_test/test_base_symbolic.py index b540410..05b56c3 100644 --- a/model_test/test_base_symbolic.py +++ b/model_test/test_base_symbolic.py @@ -13,6 +13,7 @@ real_eps = np.finfo(np.float64).eps + class TestBaseSymbolic(unittest.TestCase): reference = list() diff --git a/notebooks/symbolic_outputs/symbolic_output_land_atmosphere-AUTO.ipynb b/notebooks/symbolic_outputs/symbolic_output_land_atmosphere-AUTO.ipynb index 06495d6..a7866e6 100644 --- a/notebooks/symbolic_outputs/symbolic_output_land_atmosphere-AUTO.ipynb +++ b/notebooks/symbolic_outputs/symbolic_output_land_atmosphere-AUTO.ipynb @@ -65,7 +65,7 @@ "outputs": [], "source": [ "import sys, os\n", - "sys.path.extend([os.path.abspath('../')])" + "sys.path.extend([os.path.abspath('../../')])" ] }, { @@ -114,7 +114,7 @@ "metadata": {}, "outputs": [], "source": [ - "from qgs.functions.symbolic_tendencies import create_symbolic_equations" + "from qgs.functions.symbolic_tendencies import create_symbolic_tendencies" ] }, { @@ -238,7 +238,7 @@ "metadata": {}, "outputs": [], "source": [ - "funcs, = create_symbolic_equations(model_parameters, continuation_variables=[model_parameters.gotemperature_params.C[0], model_parameters.atemperature_params.C[0], model_parameters.atmospheric_params.kd, model_parameters.atmospheric_params.kdp], language='auto')" + "funcs, = create_symbolic_tendencies(model_parameters, continuation_variables=[model_parameters.gotemperature_params.C[0], model_parameters.atemperature_params.C[0], model_parameters.atmospheric_params.kd, model_parameters.atmospheric_params.kdp], language='auto')" ] }, { diff --git a/notebooks/symbolic_outputs/symbolic_output_land_atmosphere.ipynb b/notebooks/symbolic_outputs/symbolic_output_land_atmosphere.ipynb index 822fe14..5cf7ed6 100644 --- a/notebooks/symbolic_outputs/symbolic_output_land_atmosphere.ipynb +++ b/notebooks/symbolic_outputs/symbolic_output_land_atmosphere.ipynb @@ -54,7 +54,7 @@ "outputs": [], "source": [ "import sys, os\n", - "sys.path.extend([os.path.abspath('../')])" + "sys.path.extend([os.path.abspath('../../')])" ] }, { @@ -113,7 +113,7 @@ "metadata": {}, "outputs": [], "source": [ - "from qgs.functions.symbolic_tendencies import create_symbolic_equations\n", + "from qgs.functions.symbolic_tendencies import create_symbolic_tendencies\n", "from qgs.functions.tendencies import create_tendencies" ] }, @@ -141,7 +141,7 @@ "outputs": [], "source": [ "# Time parameters\n", - "dt = 0.1\n" + "dt = 0.1" ] }, { @@ -260,7 +260,15 @@ "outputs": [], "source": [ "%%time\n", - "funcs, = create_symbolic_equations(model_parameters, continuation_variables=[model_parameters.gotemperature_params.C[0], model_parameters.atemperature_params.C[0], model_parameters.atmospheric_params.kd, model_parameters.atmospheric_params.kdp], language='python')" + "funcs, = create_symbolic_tendencies(\n", + " model_parameters, \n", + " continuation_variables=[\n", + " model_parameters.gotemperature_params.C[0],\n", + " model_parameters.atemperature_params.C[0],\n", + " model_parameters.atmospheric_params.kd,\n", + " model_parameters.atmospheric_params.kdp\n", + " ], \n", + " language='python')" ] }, { @@ -476,10 +484,7 @@ " plt.ylabel('$'+model_parameters.latex_var_string[vary]+'$');\n", "\n", " plt.legend()\n", - " k+=1\n", - "\n", - "\n", - "\n" + " k+=1" ] } ], @@ -499,7 +504,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.16" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/qgs/functions/symbolic_mul.py b/qgs/functions/symbolic_mul.py index 3707ea4..3f32d92 100644 --- a/qgs/functions/symbolic_mul.py +++ b/qgs/functions/symbolic_mul.py @@ -8,20 +8,7 @@ """ from sympy import tensorproduct, tensorcontraction - - -def add_to_dict(dic, loc, value): - """Adds `value` to dictionary `dic`, with the dictionary key of `loc`. - If the dictionary did not have a key of `loc` before, a new key is made. - - # Jonathan: Add parameters descriptions - """ - - try: - dic[loc] += value - except: - dic[loc] = value - return dic +from qgs.functions.util import add_to_dict def symbolic_tensordot(a, b, axes=2): @@ -29,6 +16,8 @@ def symbolic_tensordot(a, b, axes=2): This is based on `Numpy`_ :meth:`~numpy.tensordot` . + .. _Numpy: https://numpy.org/ + Parameters ---------- a, b: ~sympy.tensor.array.DenseNDimArray or ~sympy.tensor.array.SparseNDimArray @@ -40,11 +29,9 @@ def symbolic_tensordot(a, b, axes=2): Returns ------- - output: sympy tensor + output: Sympy tensor The tensor dot product of the input. - .. _Numpy: https://numpy.org/ - """ as_ = a.shape nda = len(as_) diff --git a/qgs/functions/symbolic_tendencies.py b/qgs/functions/symbolic_tendencies.py index 4fd9160..0e37fb8 100644 --- a/qgs/functions/symbolic_tendencies.py +++ b/qgs/functions/symbolic_tendencies.py @@ -7,8 +7,7 @@ """ -import numpy as np -from sympy import Symbol, lambdify +from sympy import Symbol import warnings from qgs.functions.symbolic_mul import symbolic_sparse_mult2, symbolic_sparse_mult3, symbolic_sparse_mult4, \ @@ -19,19 +18,16 @@ python_lang_translation = { 'sqrt': 'math.sqrt', - 'pi': 'math.pi', - 'lambda': 'lmda' # Remove conflict for lambda function in python + 'lambda': 'lmda' # Remove conflict for lambda function in python } fortran_lang_translation = { 'conjugate': 'CONJG', 'epsilon': 'eps' # Remove conflict for EPSILON function in fortran - # TODO: may need to add variable for pi } julia_lang_translation = { '**': '^', - 'pi': 'pi()', 'conjugate': 'conj' } @@ -40,33 +36,35 @@ } -def create_symbolic_equations(params, atm_ip=None, ocn_ip=None, gnd_ip=None, continuation_variables=None, - language='python', return_inner_products=False, return_jacobian=False, - return_symbolic_eqs=False, return_symbolic_qgtensor=False): - """Function to output the raw symbolic functions of the qgs model. +def create_symbolic_tendencies(params, continuation_variables, atm_ip=None, ocn_ip=None, gnd_ip=None, + language='python', return_inner_products=False, return_jacobian=False, + return_symbolic_eqs=False, return_symbolic_qgtensor=False): + """Function to output the raw symbolic tendencies of the qgs model. Parameters ---------- params: QgParams The parameters fully specifying the model configuration. - atm_ip: SymbolicAtmosphericInnerProducts, optional + continuation_variables: list(Parameter, ScalingParameter or ParametersArray) or None + The variables to not substitute by their numerical value and to leave in the equations. + If `None`, no variables are substituted. + If an empty list is provided, then all variables are substituted, providing fully numerical tendencies. + atm_ip: AtmosphericSymbolicInnerProducts, optional Allows for stored inner products to be input. - ocn_ip: SymbolicOceanInnerProducts, optional + ocn_ip: OceanSymbolicInnerProducts, optional Allows for stored inner products to be input. - gnd_ip: SymbolicGroundInnerProducts, optional + gnd_ip: GroundSymbolicInnerProducts, optional Allows for stored inner products to be input. - continuation_variables: Iterable(Parameter, ScalingParameter, ParametersArray) - The variables to not substitute and to leave in the equations, if `None` all variables are substituted. - language: str + language: str, optional Options for the output language syntax: 'python', 'julia', 'fortran', 'auto', 'mathematica'. Default to 'python'. - return_inner_products: bool + return_inner_products: bool, optional If `True`, return the inner products of the model. Default to `False`. - return_jacobian: bool + return_jacobian: bool, optional If `True`, return the Jacobian of the model. Default to `False`. - return_symbolic_eqs: bool + return_symbolic_eqs: bool, optional If `True`, return the substituted symbolic equations. - return_symbolic_qgtensor: bool + return_symbolic_qgtensor: bool, optional If `True`, return the symbolic tendencies tensor of the model. Default to `False`. Returns @@ -75,7 +73,7 @@ def create_symbolic_equations(params, atm_ip=None, ocn_ip=None, gnd_ip=None, con The substituted functions in the language syntax specified, as a string. dict_eq_simplified: dict(~sympy.core.expr.Expr) Dictionary of the substituted Jacobian matrix. - inner_products: tuple(SymbolicAtmosphericInnerProducts, SymbolicOceanicInnerProducts, SymbolicGroundInnerProducts) + inner_products: tuple(AtmosphericSymbolicInnerProducts, OceanicSymbolicInnerProducts, GroundSymbolicInnerProducts) If `return_inner_products` is `True`, the inner products of the system. eq_simplified: dict(~sympy.core.expr.Expr) If `return_symbolic_eqs` is `True`, dictionary of the model tendencies symbolic functions. @@ -86,7 +84,11 @@ def create_symbolic_equations(params, atm_ip=None, ocn_ip=None, gnd_ip=None, con make_ip_subs = True if continuation_variables is None: - make_ip_subs = False # TODO: check this !!! For me it should be True because if no cv then all subs are done. + make_ip_subs = False + + # Generates list of all available parameters + continuation_variables = params._all_items + else: for cv in continuation_variables: try: @@ -198,15 +200,17 @@ def translate_equations(equations, language='python'): Parameters ---------- - equations: dict(string), list, string + equations: dict(str), list, str Dictionary, list, or string of the symbolic model equations. - language: string + language: str, optional Language syntax that the equations are returned in. Options are: + - `python` - `fortran` - `julia` - `auto` - `mathematica` + Default to `python`. Returns @@ -243,7 +247,7 @@ def translate_equations(equations, language='python'): for k in translator.keys(): str_eq = str_eq.replace(k, translator[k]) else: - raise warnings.warn("Expected a dict, list, or string input") + raise ValueError("Expected a dict, list, or string input") return str_eq @@ -258,17 +262,19 @@ def format_equations(equations, params, save_loc=None, language='python', print_ Dictionary of symbolic model equations. params: QgParams The parameters fully specifying the model configuration. - save_loc: str + save_loc: str, optional Location to save the outputs as a .txt file. - language: str + language: str, optional Language syntax that the equations are returned in. Options are: + - `python` - `fortran` - `julia` - `auto` - `mathematica` + Default to `python`. - print_equations: bool + print_equations: bool, optional If `True`, equations are printed by the function, if `False`, equation strings are returned by the function. Defaults to `False` @@ -335,7 +341,7 @@ def equations_to_string(equations): Returns ------- - dict(~string) + dict(str) Dictionary of the substituted symbolic model equations. """ @@ -345,7 +351,7 @@ def equations_to_string(equations): return str_eq -def equation_as_function(equations, params, language='python', continuation_variables=None): +def equation_as_function(equations, params, continuation_variables, language='python'): """Converts the symbolic equations to a function in string format in the language syntax specified, or a lambdified python function. @@ -355,23 +361,28 @@ def equation_as_function(equations, params, language='python', continuation_vari Dictionary of the substituted symbolic model equations. params: QgParams The parameters fully specifying the model configuration. - language: str + continuation_variables: list(Parameter, ScalingParameter, ParametersArray) or None + The variables to not substitute by their numerical value and to leave in the equations. + If `None`, no variables are substituted. + If an empty list is provided, then all variables are substituted, providing fully numerical tendencies. + language: str, optional Language syntax that the equations are returned in. Options are: + - `python` - `fortran` - `julia` - `auto` - `mathematica` + Default to `python`. - continuation_variables: Iterable(Parameter, ScalingParameter, ParametersArray) - Variables that are not substituted with numerical values. If `None`, all symbols are substituted. Returns ------- - f_output: callable or str + f_output: str Output is a function as a string in the specified language syntax. """ + if continuation_variables is None: continuation_variables = list() @@ -446,6 +457,7 @@ def equation_as_function(equations, params, language='python', continuation_vari f_output = ['\n'.join(auto_file), '\n'.join(auto_config)] if language == 'mathematica': + raise NotImplemented("Mathematica code output is not yet available.") # TODO: This function needs testing before release f_output.append('F = Array[' + str(len(eq_dict)) + ']') @@ -460,7 +472,7 @@ def equation_as_function(equations, params, language='python', continuation_vari return f_output -def jacobian_as_function(equations, params, language='python', continuation_variables=None): +def jacobian_as_function(equations, params, continuation_variables, language='python'): """Converts the symbolic equations of the jacobain to a function in string format in the language syntax specified, or a lambdified python function. @@ -470,23 +482,28 @@ def jacobian_as_function(equations, params, language='python', continuation_vari Dictionary of the substituted symbolic model equations. params: QgParams The parameters fully specifying the model configuration. - language: str + continuation_variables: list(Parameter, ScalingParameter, ParametersArray) or None + The variables to not substitute by their numerical value and to leave in the equations. + If `None`, no variables are substituted. + If an empty list is provided, then all variables are substituted, providing fully numerical tendencies. + language: str, optional Language syntax that the equations are returned in. Options are: + - `python` - `fortran` - `julia` - `auto` - `mathematica` + Default to `python`. - continuation_variables: Iterable(Parameter, ScalingParameter, ParametersArray) - Variables that are not substituted with numerical values. If `None`, all symbols are substituted. Returns ------- - f_output: callable or str + f_output: str Output is a function as a string in the specified language syntax """ + if continuation_variables is None: continuation_variables = list() @@ -562,6 +579,7 @@ def jacobian_as_function(equations, params, language='python', continuation_vari if language == 'mathematica': # TODO: This function needs testing before release + raise NotImplemented("Mathematica code output is not yet available.") f_output.append('jac = Array[' + str(len(eq_dict)) + ']') @@ -578,7 +596,7 @@ def jacobian_as_function(equations, params, language='python', continuation_vari def create_auto_file(equations, params, continuation_variables, auto_main_template=None, auto_c_template=None, initialize_params=False, initialize_solution=False): - """Creates the auto configuration file and the model file. + """Creates the AUTO configuration file and the model file. Saves files to specified folder. Parameters @@ -587,8 +605,9 @@ def create_auto_file(equations, params, continuation_variables, auto_main_templa Dictionary of the substituted symbolic model equations params: QgParams The parameters fully specifying the model configuration. - continuation_variables: Iterable(Parameter, ScalingParameter, ParametersArray) - Variables that are not substituted with numerical values. If `None`, all symbols are substituted + continuation_variables: list(Parameter, ScalingParameter, ParametersArray) + The variables to not substitute by their numerical value and to leave in the equations. + There must be at least one variable in this list and, due to AUTO constraints, its maximum length is 10. auto_main_template: str, optional The template to be used to generate the main AUTO file. If not provided, use the default template. @@ -602,20 +621,13 @@ def create_auto_file(equations, params, continuation_variables, auto_main_templa Returns ------- - auto_file: Str + auto_file: str The auto model file as a string - auto_config: Str + auto_config: str Auto configuration file as a string """ - # TODO: There is some weird double tab spacings in the output, and I am not sure why - - # User passes the equations, with the variables to leave as variables. - # The existing model parameters are used to populate the auto file - # The variables given as `continuation_variables` remain in the equations. - # There is a limit of 1-10 remaining variables - if (len(continuation_variables) < 1) or (len(continuation_variables) > 10): raise ValueError("Too many variables for auto file") diff --git a/qgs/functions/util.py b/qgs/functions/util.py index d6ae059..c78bdae 100644 --- a/qgs/functions/util.py +++ b/qgs/functions/util.py @@ -11,6 +11,26 @@ from numba import njit +def add_to_dict(dic, loc, value): + """Adds `value` to dictionary `dic`, with the dictionary key of `loc`. + If the dictionary did not have a key of `loc` before, a new key is made. + + Parameters + ---------- + dic: dict + Dictionary to add the value to. + loc: + Item of the dictionary to add the value to. + value: + Value to add. + """ + try: + dic[loc] += value + except: + dic[loc] = value + return dic + + @njit def reverse(a): """Numba-jitted function to reverse a 1D array. diff --git a/qgs/params/parameter.py b/qgs/params/parameter.py index 1917307..c741031 100644 --- a/qgs/params/parameter.py +++ b/qgs/params/parameter.py @@ -64,7 +64,6 @@ # TODO: Automatize warnings and errors -# TODO: Implement operations for arrays class ScalingParameter(float): """Class of model's dimension parameter. @@ -998,7 +997,7 @@ class ParametersArray(np.ndarray): description: str or list(str) or array(str), optional String or an iterable of strings, describing the parameters. If an iterable, should have the same length or shape as `values`. - symbols ~sympy.core.symbol.Symbol or list(~sympy.core.symbol.Symbol) or ~numpy.ndarray(~sympy.core.symbol.Symbol), optional + symbols: ~sympy.core.symbol.Symbol or list(~sympy.core.symbol.Symbol) or ~numpy.ndarray(~sympy.core.symbol.Symbol), optional A `Sympy`_ symbol or an iterable of symbols, to represent the parameters in symbolic expressions. If an iterable, should have the same length or shape as `values`. symbolic_expressions: ~sympy.core.expr.Expr or list(~sympy.core.expr.Expr) or ~numpy.ndarray(~sympy.core.expr.Expr), optional diff --git a/qgs/params/params.py b/qgs/params/params.py index 5d34b41..13c6fae 100644 --- a/qgs/params/params.py +++ b/qgs/params/params.py @@ -1329,6 +1329,48 @@ def dimensional_time(self): c = 24 * 3600 * 365 return 1 / (self.scale_params.f0 * c) + def _parameter_values(self, obj=None): + """Function produces a list of the values in the parameters class, or any class passed""" + + subs = list() + iter_vals = obj.__dict__.values() if obj is not None else self.__dict__.values() + for val in iter_vals: + if isinstance(val, Parameter): + subs.append(val) + + if isinstance(val, ScalingParameter): + subs.append(val) + + if isinstance(val, ParametersArray): + for v in val: + if v.symbol != 0: + subs.append(v) + return subs + + @property + def _all_items(self): + """ + Function to return a list of all values in the parameter class, + along with their current numerical values. + + Returns + ------- + list + """ + + subs = self._parameter_values() + + for _, obj in self.__dict__.items(): + if issubclass(obj.__class__, Params): + for v in self._parameter_values(obj): + subs.append(v) + + # Manually add properties from scaling class + subs.append(self.scale_params.L) + subs.append(self.scale_params.beta) + + return subs + # ------------------------------------------------------------------- # Config setters to be used with symbolic inner products # ------------------------------------------------------------------- diff --git a/qgs/tensors/symbolic_qgtensor.py b/qgs/tensors/symbolic_qgtensor.py index 64b0111..4af20bb 100644 --- a/qgs/tensors/symbolic_qgtensor.py +++ b/qgs/tensors/symbolic_qgtensor.py @@ -6,7 +6,8 @@ equations. """ -from qgs.functions.symbolic_mul import add_to_dict, symbolic_tensordot +from qgs.functions.symbolic_mul import symbolic_tensordot +from qgs.functions.util import add_to_dict from qgs.params.params import Parameter, ScalingParameter, ParametersArray, Params import numpy as np @@ -54,9 +55,9 @@ class SymbolicQgsTensor(object): ground_inner_products: None or GroundInnerProducts The inner products of the ground basis functions on which the model's PDE ground equations are projected. If `None`, disable the ground tendencies. Default to `None`. - tensor: sparse.COO(float) + tensor: ~sympy.tensor.array.ImmutableSparseNDimArray The tensor :math:`\\mathcal{T}_{i,j,k}` :math:`i`-th components. - jacobian_tensor: sparse.COO(float) + jacobian_tensor: ~sympy.tensor.array.ImmutableSparseNDimArray The jacobian tensor :math:`\\mathcal{T}_{i,j,k} + \\mathcal{T}_{i,k,j}` :math:`i`-th components. """ @@ -800,16 +801,20 @@ def sub_tensor(self, tensor=None, continuation_variables=None): Parameters ---------- - tensor: dict or ~sympy.tensor.array.ImmutableSparseNDimArray + tensor: dict(~sympy.core.expr.Expr or float) or ~sympy.tensor.array.ImmutableSparseNDimArray + Tensor of model tendencies, either as - continuation_variables: Iterable(Parameter, ScalingParameter, ParametersArray) + - a dictionary with keys of non-zero coordinates, and values of Sympy expressions or floats + - or a sparse Sympy tensor + + continuation_variables: list(Parameter, ScalingParameter or ParametersArray) or `None` Variables which remain symbolic, all other variables are substituted with numerical values. If `None` all variables are substituted. Returns ------- - ten_out: dict - Dictionary of the substituted tensor of the model tendencies, with coordinates and value + ten_out: dict(float) + Dictionary of the substituted tensor of the model tendencies, with coordinates and numerical values """ if continuation_variables is None: @@ -842,14 +847,21 @@ def print_tensor(self, tensor=None, dict_opp=True, tol=1e-10): Parameters ---------- - tensor: dict or ~sympy.tensor.array.ImmutableSparseNDimArray - Tensor of model tendencies, either as a dictionary with keys of non-zero coordinates, and values of ~sympy.core.symbol.Symbol or floats, or as a - ~sympy.tensor.array.ImmutableSparseNDimArray . + tensor: dict(~sympy.core.expr.Expr or float) or ~sympy.tensor.array.ImmutableSparseNDimArray or `None` + Tensor of model tendencies, either as + + - a dictionary with keys of non-zero coordinates, and values of Sympy expressions or floats + - or a sparse Sympy tensor + + If `None`, defaults to the stored tensor. + Defaults to `None`. + dict_opp: bool - ... + If `True`, returns the unsimplified symbolic expressions, if `False` the simplified expressions are returned. tol: float - ... + The tolerance to allow for numerical errors when finding non-zero values. """ + if tensor is None: if dict_opp: temp_ten = self.tensor_dic @@ -910,9 +922,9 @@ class SymbolicQgsTensorDynamicT(SymbolicQgsTensor): ground_inner_products: None or GroundInnerProducts The inner products of the ground basis functions on which the model's PDE ground equations are projected. If `None`, disable the ground tendencies. Default to `None`. - tensor: sparse.COO(float) + tensor: ~sympy.tensor.array.ImmutableSparseNDimArray The tensor :math:`\\mathcal{T}_{i,j,k}` :math:`i`-th components. - jacobian_tensor: sparse.COO(float) + jacobian_tensor: ~sympy.tensor.array.ImmutableSparseNDimArray The jacobian tensor :math:`\\mathcal{T}_{i,j,k} + \\mathcal{T}_{i,k,j}` :math:`i`-th components. """ @@ -1220,7 +1232,6 @@ def _compute_non_stored_full_dict(self): def compute_tensor(self): """Routine to compute the tensor.""" - # gathering if self.params.T4: # TODO: Make a proper error message for here raise ValueError("Parameters are set for T4 version, set dynamic_T=True") @@ -1238,9 +1249,6 @@ def compute_tensor(self): class SymbolicQgsTensorT4(SymbolicQgsTensor): - # TODO: this takes a long time (>1hr) to run. I think we need a better way to run the non-stored z, v, Z, V IPs. Maybe do not allow `n` as a continuation parameter for this version? - # TODO: Create a warning about long run-times. - """qgs dynamical temperature first order (linear) symbolic tendencies tensor class. Parameters @@ -1270,9 +1278,9 @@ class SymbolicQgsTensorT4(SymbolicQgsTensor): ground_inner_products: None or GroundInnerProducts The inner products of the ground basis functions on which the model's PDE ground equations are projected. If `None`, disable the ground tendencies. Default to `None`. - tensor: sparse.COO(float) + tensor: ~sympy.tensor.array.ImmutableSparseNDimArray The tensor :math:`\\mathcal{T}_{i,j,k}` :math:`i`-th components. - jacobian_tensor: sparse.COO(float) + jacobian_tensor: ~sympy.tensor.array.ImmutableSparseNDimArray The jacobian tensor :math:`\\mathcal{T}_{i,j,k} + \\mathcal{T}_{i,k,j}` :math:`i`-th components. """ @@ -1420,7 +1428,6 @@ def _compute_non_stored_full_dict(self): def compute_tensor(self): """Routine to compute the tensor.""" - # gathering if not self.params.T4: raise ValueError("Parameters are not set for T4 version") @@ -1437,10 +1444,8 @@ def compute_tensor(self): def _kronecker_delta(i, j): - if i == j: return 1 - else: return 0 @@ -1451,9 +1456,11 @@ def _shift_dict_keys(dic, shift): Parameters ---------- - dic: dictionary + dic: dict + Dictionary representing a tensor. - shift: Tuple + shift: tuple + A tuple that represents the shift in the tensor indices. """ shifted_dic = dict() @@ -1465,47 +1472,32 @@ def _shift_dict_keys(dic, shift): def _parameter_substitutions(params, continuation_variables): - - subs = _parameter_values(params) - for _, obj in params.__dict__.items(): - if issubclass(obj.__class__, Params): - subs.update(_parameter_values(obj)) + """ + Returns a dict of parameters values that are to be substituted, + removing the parameters given in `continuation_variables`. + """ + + subs = params._all_items - # Manually add properties from class - subs[params.scale_params.L.symbol] = params.scale_params.L - subs[params.scale_params.beta.symbol] = params.scale_params.beta + if continuation_variables is None: + continuation_variables = list() # Remove variables in continuation variables for cv in continuation_variables: if isinstance(cv, ParametersArray): - for cv_i in cv.symbols: - subs.pop(cv_i) - elif hasattr(cv, "symbol"): - subs.pop(cv.symbol) + for cv_i in cv: + subs.remove(cv_i) + elif isinstance(cv, Parameter): + subs.remove(cv) else: # Try ... who knows... - subs.pop(cv) - - return subs - - -def _parameter_values(pars): - """Function takes a parameter class and produces a dictionary of the symbol and the corresponding numerical value""" - - subs = dict() - for val in pars.__dict__.values(): - if isinstance(val, Parameter): - if val.symbol is not None: - subs[val.symbol] = val - - if isinstance(val, ScalingParameter): - if val.symbol is not None: - subs[val.symbol] = val - - if isinstance(val, ParametersArray): - for v in val: - if v.symbol is not None or v.symbol != 0: - subs[v.symbol] = v - return subs + subs.remove(cv) + + # make the remaining items into a dict to pass to sympy subs function + sub_dic = {} + for p in subs: + if p.symbol is not None: + sub_dic[p.symbol] = float(p) + return sub_dic if __name__ == "__main__":