From 0ae244c293bf847ce5f64d95246d6fd510906df1 Mon Sep 17 00:00:00 2001 From: Alec Bills Date: Sun, 13 Nov 2022 15:56:23 -0500 Subject: [PATCH 01/10] cmd-shift-f pass --- .github/workflows/test_on_push.yml | 10 - .gitignore | 3 - docs/source/util.rst | 2 - pybamm/__init__.py | 6 - .../operations/evaluate_julia.py | 1084 ----------------- pybamm/models/base_model.py | 77 -- requirements.txt | 1 - setup.py | 3 - .../test_operations/test_evaluate_julia.py | 366 ------ tests/unit/test_models/test_base_model.py | 36 - .../test_base_model_generate_julia_diffeq.py | 131 -- tests/unit/test_solvers/test_julia_mtk.py | 240 ---- tox.ini | 12 +- 13 files changed, 1 insertion(+), 1970 deletions(-) delete mode 100644 pybamm/expression_tree/operations/evaluate_julia.py delete mode 100644 tests/unit/test_expression_tree/test_operations/test_evaluate_julia.py delete mode 100644 tests/unit/test_models/test_base_model_generate_julia_diffeq.py delete mode 100644 tests/unit/test_solvers/test_julia_mtk.py diff --git a/.github/workflows/test_on_push.yml b/.github/workflows/test_on_push.yml index 203af0e771..c892aa3121 100644 --- a/.github/workflows/test_on_push.yml +++ b/.github/workflows/test_on_push.yml @@ -55,12 +55,6 @@ jobs: with: python-version: ${{ matrix.python-version }} - - name: Set up Julia - if: matrix.os != 'windows-latest' - uses: julia-actions/setup-julia@v1 - with: - version: 1.7 - - name: Install Linux system dependencies if: matrix.os == 'ubuntu-latest' run: | @@ -92,10 +86,6 @@ jobs: if: matrix.os == 'ubuntu-latest' run: tox -e pybamm-requires - - name: Install Julia - if: matrix.os != 'windows-latest' - run: tox -e julia - - name: Run unit tests for GNU/Linux with Python 3.8 if: matrix.os == 'ubuntu-latest' && matrix.python-version != 3.9 run: python -m tox -e unit diff --git a/.gitignore b/.gitignore index 4986d7f0a5..e0b2d092f2 100644 --- a/.gitignore +++ b/.gitignore @@ -112,9 +112,6 @@ test.json # tox .tox/ -# julia -Manifest.toml - # vcpkg vcpkg_installed/ diff --git a/docs/source/util.rst b/docs/source/util.rst index f6074eacf4..937771d0ad 100644 --- a/docs/source/util.rst +++ b/docs/source/util.rst @@ -22,8 +22,6 @@ Utility functions .. autofunction:: pybamm.get_parameters_filepath -.. autofunction:: pybamm.have_julia - .. autofunction:: pybamm.have_jax .. autofunction:: pybamm.is_jax_compatible diff --git a/pybamm/__init__.py b/pybamm/__init__.py index cf2b690c2b..699984b162 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -49,7 +49,6 @@ have_jax, install_jax, is_jax_compatible, - have_julia, get_git_commit_info, ) from .logger import logger, set_logging_level @@ -96,11 +95,6 @@ from .expression_tree.operations.jacobian import Jacobian from .expression_tree.operations.convert_to_casadi import CasadiConverter from .expression_tree.operations.unpack_symbols import SymbolUnpacker -from .expression_tree.operations.evaluate_julia import ( - get_julia_function, - get_julia_mtk_model, -) - # # Model classes # diff --git a/pybamm/expression_tree/operations/evaluate_julia.py b/pybamm/expression_tree/operations/evaluate_julia.py deleted file mode 100644 index f5fa656343..0000000000 --- a/pybamm/expression_tree/operations/evaluate_julia.py +++ /dev/null @@ -1,1084 +0,0 @@ -# -# Write a symbol to Julia -# -import pybamm - -import numpy as np -import scipy.sparse -from collections import OrderedDict -from pybamm.util import is_constant_and_can_evaluate - -import numbers - - -def id_to_julia_variable(symbol_id, prefix): - """ - This function defines the format for the julia variable names used in find_symbols - and to_julia. Variable names are based on a nodes' id to make them unique - """ - var_format = prefix + "_{:05d}" - # Need to replace "-" character to make them valid julia variable names - return var_format.format(symbol_id).replace("-", "m") - - -def find_symbols( - symbol, - constant_symbols, - variable_symbols, - variable_symbol_sizes, - round_constants=True, -): - """ - This function converts an expression tree to a dictionary of node id's and strings - specifying valid julia code to calculate that nodes value, given y and t. - - The function distinguishes between nodes that represent constant nodes in the tree - (e.g. a pybamm.Matrix), and those that are variable (e.g. subtrees that contain - pybamm.StateVector). The former are put in `constant_symbols`, the latter in - `variable_symbols` - - Note that it is important that the arguments `constant_symbols` and - `variable_symbols` be and *ordered* dict, since the final ordering of the code lines - are important for the calculations. A dict is specified rather than a list so that - identical subtrees (which give identical id's) are not recalculated in the code - - Parameters - ---------- - symbol : :class:`pybamm.Symbol` - The symbol or expression tree to convert - - constant_symbol : collections.OrderedDict - The output dictionary of constant symbol ids to lines of code - - variable_symbol : collections.OrderedDict - The output dictionary of variable (with y or t) symbol ids to lines of code - - variable_symbol_sizes : collections.OrderedDict - The output dictionary of variable (with y or t) symbol ids to size of that - variable, for caching - - """ - # ignore broadcasts for now - if isinstance(symbol, pybamm.Broadcast): - symbol = symbol.child - if is_constant_and_can_evaluate(symbol): - value = symbol.evaluate() - if round_constants: - value = np.round(value, decimals=11) - if not isinstance(value, numbers.Number): - if scipy.sparse.issparse(value): - # Create Julia SparseArray - row, col, data = scipy.sparse.find(value) - if round_constants: - data = np.round(data, decimals=11) - m, n = value.shape - # Set print options large enough to avoid ellipsis - # at least as big is len(row) = len(col) = len(data) - np.set_printoptions( - threshold=max(np.get_printoptions()["threshold"], len(row) + 10) - ) - # increase precision for printing - np.set_printoptions(precision=20) - # add 1 to correct for 1-indexing in Julia - # use array2string so that commas are included - constant_symbols[symbol.id] = "sparse({}, {}, {}, {}, {})".format( - np.array2string(row + 1, separator=","), - np.array2string(col + 1, separator=","), - np.array2string(data, separator=","), - m, - n, - ) - elif value.shape == (1, 1): - # Extract value if array has only one entry - constant_symbols[symbol.id] = value[0, 0] - variable_symbol_sizes[symbol.id] = 1 - elif value.shape[1] == 1: - # Set print options large enough to avoid ellipsis - # at least as big as len(row) = len(col) = len(data) - np.set_printoptions( - threshold=max( - np.get_printoptions()["threshold"], value.shape[0] + 10 - ) - ) - # Flatten a 1D array - constant_symbols[symbol.id] = np.array2string( - value.flatten(), separator="," - ) - variable_symbol_sizes[symbol.id] = symbol.shape[0] - else: - constant_symbols[symbol.id] = value - # No need to save the size as it will not need to be used - return - - # process children recursively - for child in symbol.children: - find_symbols( - child, - constant_symbols, - variable_symbols, - variable_symbol_sizes, - round_constants=round_constants, - ) - - # calculate the variable names that will hold the result of calculating the - # children variables - children_vars = [] - for child in symbol.children: - if isinstance(child, pybamm.Broadcast): - child = child.child - if is_constant_and_can_evaluate(child): - child_eval = child.evaluate() - if isinstance(child_eval, numbers.Number): - children_vars.append(str(child_eval)) - else: - children_vars.append(id_to_julia_variable(child.id, "const")) - else: - children_vars.append(id_to_julia_variable(child.id, "cache")) - - if isinstance(symbol, pybamm.BinaryOperator): - # TODO: we can pass through a dummy y and t to get the type and then hardcode - # the right line, avoiding these checks - if isinstance(symbol, pybamm.MatrixMultiplication): - symbol_str = "{0} @ {1}".format(children_vars[0], children_vars[1]) - elif isinstance(symbol, pybamm.Inner): - symbol_str = "{0} * {1}".format(children_vars[0], children_vars[1]) - elif isinstance(symbol, pybamm.Minimum): - symbol_str = "min({},{})".format(children_vars[0], children_vars[1]) - elif isinstance(symbol, pybamm.Maximum): - symbol_str = "max({},{})".format(children_vars[0], children_vars[1]) - elif isinstance(symbol, pybamm.Power): - # julia uses ^ instead of ** for power - # include dot for elementwise operations - symbol_str = children_vars[0] + " .^ " + children_vars[1] - else: - # all other operations use the same symbol - symbol_str = children_vars[0] + " " + symbol.name + " " + children_vars[1] - - elif isinstance(symbol, pybamm.UnaryOperator): - # Index has a different syntax than other univariate operations - if isinstance(symbol, pybamm.Index): - # Because of how julia indexing works, add 1 to the start, but not to the - # stop - symbol_str = "{}[{}:{}]".format( - children_vars[0], symbol.slice.start + 1, symbol.slice.stop - ) - elif isinstance(symbol, pybamm.Gradient): - symbol_str = "grad_{}({})".format(tuple(symbol.domain), children_vars[0]) - elif isinstance(symbol, pybamm.Divergence): - symbol_str = "div_{}({})".format(tuple(symbol.domain), children_vars[0]) - elif isinstance(symbol, pybamm.Broadcast): - # ignore broadcasts for now - symbol_str = children_vars[0] - elif isinstance(symbol, pybamm.BoundaryValue): - symbol_str = "boundary_value_{}({})".format(symbol.side, children_vars[0]) - else: - symbol_str = symbol.name + children_vars[0] - - elif isinstance(symbol, pybamm.Function): - # write functions directly - symbol_str = "{}({})".format(symbol.julia_name, ", ".join(children_vars)) - - elif isinstance(symbol, (pybamm.Variable, pybamm.ConcatenationVariable)): - # No need to do anything if a Variable is found - return - - elif isinstance(symbol, pybamm.Concatenation): - if isinstance(symbol, (pybamm.NumpyConcatenation, pybamm.SparseStack)): - # return a list of the children variables, which will be converted to a - # line by line assignment - # return this as a string so that other functionality still works - # also save sizes - symbol_str = "[" - for child in children_vars: - child_id = child[6:].replace("m", "-") - size = variable_symbol_sizes[int(child_id)] - symbol_str += "{}::{}, ".format(size, child) - symbol_str = symbol_str[:-2] + "]" - - # DomainConcatenation specifies a particular ordering for the concatenation, - # which we must follow - elif isinstance(symbol, pybamm.DomainConcatenation): - if symbol.secondary_dimensions_npts == 1: - all_child_vectors = children_vars - all_child_sizes = [ - variable_symbol_sizes[int(child[6:].replace("m", "-"))] - for child in children_vars - ] - else: - slice_starts = [] - all_child_vectors = [] - all_child_sizes = [] - for i in range(symbol.secondary_dimensions_npts): - child_vectors = [] - child_sizes = [] - for child_var, slices in zip( - children_vars, symbol._children_slices - ): - for child_dom, child_slice in slices.items(): - slice_starts.append(symbol._slices[child_dom][i].start) - # add 1 to slice start to account for julia indexing - child_vectors.append( - "@view {}[{}:{}]".format( - child_var, - child_slice[i].start + 1, - child_slice[i].stop, - ) - ) - child_sizes.append( - child_slice[i].stop - child_slice[i].start - ) - all_child_vectors.extend( - [v for _, v in sorted(zip(slice_starts, child_vectors))] - ) - all_child_sizes.extend( - [v for _, v in sorted(zip(slice_starts, child_sizes))] - ) - # return a list of the children variables, which will be converted to a - # line by line assignment - # return this as a string so that other functionality still works - # also save sizes - symbol_str = "[" - for child, size in zip(all_child_vectors, all_child_sizes): - symbol_str += "{}::{}, ".format(size, child) - symbol_str = symbol_str[:-2] + "]" - - else: - # A regular Concatenation for the MTK model - # We will define the concatenation function separately - symbol_str = "concatenation(" + ", ".join(children_vars) + ")" - - # Note: we assume that y is being passed as a column vector - elif isinstance(symbol, pybamm.StateVectorBase): - if isinstance(symbol, pybamm.StateVector): - name = "@view y" - elif isinstance(symbol, pybamm.StateVectorDot): - name = "@view dy" - indices = np.argwhere(symbol.evaluation_array).reshape(-1).astype(np.int32) - # add 1 since julia uses 1-indexing - indices += 1 - if len(indices) == 1: - symbol_str = "{}[{}]".format(name, indices[0]) - else: - # julia does include the final value - symbol_str = "{}[{}:{}]".format(name, indices[0], indices[-1]) - - elif isinstance(symbol, pybamm.Time): - symbol_str = "t" - - elif isinstance(symbol, pybamm.InputParameter): - symbol_str = "inputs['{}']".format(symbol.name) - - elif isinstance(symbol, pybamm.SpatialVariable): - symbol_str = symbol.name - - elif isinstance(symbol, pybamm.FunctionParameter): - symbol_str = "{}({})".format(symbol.name, ", ".join(children_vars)) - - else: - raise NotImplementedError( - "Conversion to Julia not implemented for a symbol of type '{}'".format( - type(symbol) - ) - ) - - variable_symbols[symbol.id] = symbol_str - - # Save the size of the symbol - try: - if symbol.shape == (): - variable_symbol_sizes[symbol.id] = 1 - else: - variable_symbol_sizes[symbol.id] = symbol.shape[0] - except NotImplementedError: - pass - - -def to_julia(symbol, round_constants=True): - """ - This function converts an expression tree into a dict of constant input values, and - valid julia code that acts like the tree's :func:`pybamm.Symbol.evaluate` function - - Parameters - ---------- - symbol : :class:`pybamm.Symbol` - The symbol to convert to julia code - - Returns - ------- - constant_values : collections.OrderedDict - dict mapping node id to a constant value. Represents all the constant nodes in - the expression tree - str - valid julia code that will evaluate all the variable nodes in the tree. - - """ - - constant_values = OrderedDict() - variable_symbols = OrderedDict() - variable_symbol_sizes = OrderedDict() - find_symbols( - symbol, - constant_values, - variable_symbols, - variable_symbol_sizes, - round_constants=round_constants, - ) - - return constant_values, variable_symbols, variable_symbol_sizes - - -def get_julia_function( - symbol, - funcname="f", - input_parameter_order=None, - len_rhs=None, - preallocate=True, - round_constants=True, -): - """ - Converts a pybamm expression tree into pure julia code that will calculate the - result of calling `evaluate(t, y)` on the given expression tree. - - Parameters - ---------- - symbol : :class:`pybamm.Symbol` - The symbol to convert to julia code - funcname : str, optional - The name to give to the function (default 'f') - input_parameter_order : list, optional - List of input parameter names. Defines the order in which the input parameters - are extracted from 'p' in the julia function that is created - len_rhs : int, optional - The number of ODEs in the discretized differential equations. This also - determines whether the model has any algebraic equations: if None (default), - the model is assume to have no algebraic parts and ``julia_str`` is compatible - with an ODE solver. If not None, ``julia_str`` is compatible with a DAE solver - preallocate : bool, optional - Whether to write the function in a way that preallocates memory for the output. - Default is True, which is faster. Must be False for the function to be - modelingtoolkitized. - - Returns - ------- - julia_str : str - String of julia code, to be evaluated by ``julia.Main.eval`` - - """ - if len_rhs is None: - typ = "ode" - else: - typ = "dae" - # Take away dy from the differential states - # we will return a function of the form - # out[] = .. - dy[] for the differential states - # out[] = .. for the algebraic states - symbol_minus_dy = [] - end = 0 - for child in symbol.orphans: - start = end - end += child.size - if end <= len_rhs: - symbol_minus_dy.append(child - pybamm.StateVectorDot(slice(start, end))) - else: - symbol_minus_dy.append(child) - symbol = pybamm.numpy_concatenation(*symbol_minus_dy) - constants, var_symbols, var_symbol_sizes = to_julia( - symbol, round_constants=round_constants - ) - - # extract constants in generated function - const_and_cache_str = "cs = (\n" - shorter_const_names = {} - for i_const, (symbol_id, const_value) in enumerate(constants.items()): - const_name = id_to_julia_variable(symbol_id, "const") - const_name_short = "const_{}".format(i_const) - const_and_cache_str += " {} = {},\n".format(const_name_short, const_value) - shorter_const_names[const_name] = const_name_short - - # Pop (get and remove) items from the dictionary of symbols one by one - # If they are simple operations (@view, +, -, *, /), replace all future - # occurences instead of assigning them. This "inlining" speeds up the computation - inlineable_symbols = ["@view", "+", "-", "*", "/"] - var_str = "" - input_parameters = {} - while var_symbols: - var_symbol_id, symbol_line = var_symbols.popitem(last=False) - julia_var = id_to_julia_variable(var_symbol_id, "cache") - # Look for lists in the variable symbols. These correpsond to concatenations, so - # assign the children to the right parts of the vector - if symbol_line[0] == "[" and symbol_line[-1] == "]": - # convert to actual list - symbol_line = symbol_line[1:-1].split(", ") - start = 0 - if preallocate is True or var_symbol_id == symbol.id: - for child_size_and_name in symbol_line: - child_size, child_name = child_size_and_name.split("::") - end = start + int(child_size) - # add 1 to start to account for julia 1-indexing - var_str += "@. {}[{}:{}] = {}\n".format( - julia_var, start + 1, end, child_name - ) - start = end - else: - concat_str = "{} = vcat(".format(julia_var) - for i, child_size_and_name in enumerate(symbol_line): - child_size, child_name = child_size_and_name.split("::") - var_str += "x{} = @. {}\n".format(i + 1, child_name) - concat_str += "x{}, ".format(i + 1) - var_str += concat_str[:-2] + ")\n" - # use mul! for matrix multiplications (requires LinearAlgebra library) - elif " @ " in symbol_line: - if preallocate is False: - symbol_line = symbol_line.replace(" @ ", " * ") - var_str += "{} = {}\n".format(julia_var, symbol_line) - else: - symbol_line = symbol_line.replace(" @ ", ", ") - var_str += "mul!({}, {})\n".format(julia_var, symbol_line) - # find input parameters - elif symbol_line.startswith("inputs"): - input_parameters[julia_var] = symbol_line[8:-2] - elif "minimum" in symbol_line or "maximum" in symbol_line: - var_str += "{} .= {}\n".format(julia_var, symbol_line) - else: - # don't replace the matrix multiplication cases (which will be - # turned into a mul!), since it is faster to assign to a cache array - # first in that case - # e.g. mul!(cs.cache_1, cs.cache_2, cs.cache_3) - # unless it is a @view in which case we don't - # need to cache - # e.g. mul!(cs.cache_1, cs.cache_2, @view y[1:10]) - # also don't replace the minimum() or maximum() cases as we can't - # broadcast them - any_matmul_min_max = any( - julia_var in next_symbol_line - and ( - any( - x in next_symbol_line - for x in [" @ ", "mul!", "minimum", "maximum"] - ) - and not symbol_line.startswith("@view") - ) - for next_symbol_line in var_symbols.values() - ) - # inline operation if it can be inlined - if ( - any(x in symbol_line for x in inlineable_symbols) or symbol_line == "t" - ) and not any_matmul_min_max: - found_replacement = False - # replace all other occurrences of the variable - # in the dictionary with the symbol line - for next_var_id, next_symbol_line in var_symbols.items(): - if julia_var in next_symbol_line: - if symbol_line == "t": - # no brackets needed - var_symbols[next_var_id] = next_symbol_line.replace( - julia_var, symbol_line - ) - else: - # add brackets so that the order of operations is maintained - var_symbols[next_var_id] = next_symbol_line.replace( - julia_var, "({})".format(symbol_line) - ) - found_replacement = True - if not found_replacement: - var_str += "@. {} = {}\n".format(julia_var, symbol_line) - - # otherwise assign - else: - var_str += "@. {} = {}\n".format(julia_var, symbol_line) - # Replace all input parameter names - for input_parameter_id, input_parameter_name in input_parameters.items(): - var_str = var_str.replace(input_parameter_id, input_parameter_name) - - # indent code - var_str = " " + var_str - var_str = var_str.replace("\n", "\n ") - - # add the cache variables to the cache NamedTuple - i_cache = 0 - for var_symbol_id, var_symbol_size in var_symbol_sizes.items(): - # Skip caching the result variable since this is provided as dy - # Also skip caching the result variable if it doesn't appear in the var_str, - # since it has been inlined and does not need to be assigned to - julia_var = id_to_julia_variable(var_symbol_id, "cache") - if var_symbol_id != symbol.id and julia_var in var_str: - julia_var_short = "cache_{}".format(i_cache) - var_str = var_str.replace(julia_var, julia_var_short) - i_cache += 1 - if preallocate is True: - const_and_cache_str += " {} = zeros({}),\n".format( - julia_var_short, var_symbol_size - ) - else: - # Cache variables have not been preallocated - var_str = var_str.replace( - "@. {} = ".format(julia_var_short), - "{} = @. ".format(julia_var_short), - ) - - # Shorten the name of the constants from id to const_0, const_1, etc. - for long, short in shorter_const_names.items(): - var_str = var_str.replace(long, "cs." + short) - - # close the constants and cache string - const_and_cache_str += ")\n" - - # remove the constant and cache sring if it is empty - const_and_cache_str = const_and_cache_str.replace("cs = (\n)\n", "") - - # calculate the final variable that will output the result - if symbol.is_constant(): - result_var = id_to_julia_variable(symbol.id, "const") - if result_var in shorter_const_names: - result_var = shorter_const_names[result_var] - result_value = symbol.evaluate() - if isinstance(result_value, numbers.Number): - var_str = var_str + "\n dy .= " + str(result_value) + "\n" - else: - var_str = var_str + "\n dy .= cs." + result_var + "\n" - else: - result_var = id_to_julia_variable(symbol.id, "cache") - if typ == "ode": - out = "dy" - elif typ == "dae": - out = "out" - # replace "cache_123 = ..." with "dy .= ..." (ensure we allocate to the - # variable that was passed in) - var_str = var_str.replace(f" {result_var} =", f" {out} .=") - # catch other cases for dy - var_str = var_str.replace(result_var, out) - - # add "cs." to cache names - if preallocate is True: - var_str = var_str.replace("cache", "cs.cache") - - # line that extracts the input parameters in the right order - if input_parameter_order is None: - input_parameter_extraction = "" - elif len(input_parameter_order) == 1: - # extract the single parameter - input_parameter_extraction = " " + input_parameter_order[0] + " = p[1]\n" - else: - # extract all parameters - input_parameter_extraction = " " + ", ".join(input_parameter_order) + " = p\n" - - if preallocate is False or const_and_cache_str == "": - func_def = f"{funcname}!" - else: - func_def = f"{funcname}_with_consts!" - - # add function def - if typ == "ode": - function_def = f"\nfunction {func_def}(dy, y, p, t)\n" - elif typ == "dae": - function_def = f"\nfunction {func_def}(out, dy, y, p, t)\n" - julia_str = ( - "begin\n" - + const_and_cache_str - + function_def - + input_parameter_extraction - + var_str - ) - - # close the function, with a 'nothing' to avoid allocations - julia_str += "nothing\nend\n\n" - julia_str = julia_str.replace("\n \n", "\n") - - if not (preallocate is False or const_and_cache_str == ""): - # Use a let block for the cached variables - # open the let block - julia_str = julia_str.replace("cs = (", f"{funcname}! = let cs = (") - # close the let block - julia_str += "end\n" - - # close the "begin" - julia_str += "end" - - return julia_str - - -def convert_var_and_eqn_to_str(var, eqn, all_constants_str, all_variables_str, typ): - """ - Converts a variable and its equation to a julia string - - Parameters - ---------- - var : :class:`pybamm.Symbol` - The variable (key in the dictionary of rhs/algebraic/initial conditions) - eqn : :class:`pybamm.Symbol` - The equation (value in the dictionary of rhs/algebraic/initial conditions) - all_constants_str : str - String containing all the constants defined so far - all_variables_str : str - String containing all the variables defined so far - typ : str - The type of the variable/equation pair being converted ("equation", "initial - condition", or "boundary condition") - - Returns - ------- - all_constants_str : str - Updated string of all constants - all_variables_str : str - Updated string of all variables - eqn_str : str - The string describing the final equation result, perhaps as a function of some - variables and/or constants in all_constants_str and all_variables_str - - """ - if isinstance(eqn, pybamm.Broadcast): - # ignore broadcasts for now - eqn = eqn.child - - var_symbols = to_julia(eqn)[1] - - # var_str = "" - # for symbol_id, symbol_line in var_symbols.items(): - # var_str += f"{id_to_julia_variable(symbol_id)} = {symbol_line}\n" - # Pop (get and remove) items from the dictionary of symbols one by one - # If they are simple operations (+, -, *, /), replace all future - # occurences instead of assigning them. - inlineable_symbols = [" + ", " - ", " * ", " / "] - var_str = "" - while var_symbols: - var_symbol_id, symbol_line = var_symbols.popitem(last=False) - julia_var = id_to_julia_variable(var_symbol_id, "cache") - # inline operation if it can be inlined - if "concatenation" not in symbol_line: - found_replacement = False - # replace all other occurrences of the variable - # in the dictionary with the symbol line - for next_var_id, next_symbol_line in var_symbols.items(): - if ( - symbol_line == "t" - or " " not in symbol_line - or symbol_line.startswith("grad") - or not any(x in next_symbol_line for x in inlineable_symbols) - ): - # cases that don't need brackets - var_symbols[next_var_id] = next_symbol_line.replace( - julia_var, symbol_line - ) - elif next_symbol_line.startswith("concatenation"): - var_symbols[next_var_id] = next_symbol_line.replace( - julia_var, f"\n {symbol_line}\n" - ) - else: - # add brackets so that the order of operations is maintained - var_symbols[next_var_id] = next_symbol_line.replace( - julia_var, "({})".format(symbol_line) - ) - found_replacement = True - if not found_replacement: - var_str += "{} = {}\n".format(julia_var, symbol_line) - - # otherwise assign - else: - var_str += "{} = {}\n".format(julia_var, symbol_line) - - # If we have created a concatenation we need to define it - # Hardcoded to the negative electrode, separator, positive electrode case for now - if "concatenation" in var_str and "function concatenation" not in all_variables_str: - concatenation_def = ( - "\nfunction concatenation(n, s, p)\n" - + " # A concatenation in the electrolyte domain\n" - + " IfElse.ifelse(\n" - + " x < neg_width, n, IfElse.ifelse(\n" - + " x < neg_plus_sep_width, s, p\n" - + " )\n" - + " )\n" - + "end\n" - ) - else: - concatenation_def = "" - - # Define the FunctionParameter objects that have not yet been defined - function_defs = "" - for x in eqn.pre_order(): - if ( - isinstance(x, pybamm.FunctionParameter) - and f"function {x.name}" not in all_variables_str - and typ == "equation" - ): - function_def = ( - f"\nfunction {x.name}(" - + ", ".join(x.arg_names) - + ")\n" - + " {}\n".format(str(x.callable).replace("**", "^")) - + "end\n" - ) - function_defs += function_def - - if concatenation_def + function_defs != "": - function_defs += "\n" - - var_str = concatenation_def + function_defs + var_str - - # add a comment labeling the equation, and the equation itself - if var_str == "": - all_variables_str += "" - else: - all_variables_str += f"# '{var.name}' {typ}\n" + var_str + "\n" - - # calculate the final variable that will output the result - if eqn.is_constant(): - result_var = id_to_julia_variable(eqn.id, "const") - else: - result_var = id_to_julia_variable(eqn.id, "cache") - if is_constant_and_can_evaluate(eqn): - result_value = eqn.evaluate() - else: - result_value = None - - # define the variable that goes into the equation - if eqn.is_constant() and isinstance(result_value, numbers.Number): - eqn_str = str(result_value) - else: - eqn_str = result_var - - return all_constants_str, all_variables_str, eqn_str - - -def get_julia_mtk_model(model, geometry=None, tspan=None): - """ - Converts a pybamm model into a Julia ModelingToolkit model - - Parameters - ---------- - model : :class:`pybamm.BaseModel` - The model to be converted - geometry : dict, optional - Dictionary defining the geometry. Must be provided if the model is a PDE model - tspan : array-like, optional - Time for which to solve the model. Must be provided if the model is a PDE model - - Returns - ------- - mtk_str : str - String of julia code representing a model in MTK, - to be evaluated by ``julia.Main.eval`` - """ - # Extract variables - variables = {**model.rhs, **model.algebraic}.keys() - variable_to_print_name = {} - for i, var in enumerate(variables): - if var.print_name is not None: - print_name = var._raw_print_name - else: - print_name = f"u{i+1}" - variable_to_print_name[var] = print_name - if isinstance(var, pybamm.ConcatenationVariable): - for child in var.children: - variable_to_print_name[child] = print_name - - # Extract domain and auxiliary domains - all_domains = set( - [tuple(dom) for var in variables for dom in var.domains.values() if dom != []] - ) - is_pde = bool(all_domains) - - # Check geometry and tspan have been provided if a PDE - if is_pde: - if geometry is None: - raise ValueError("must provide geometry if the model is a PDE model") - if tspan is None: - raise ValueError("must provide tspan if the model is a PDE model") - - # Read domain names - domain_name_to_symbol = {} - long_domain_symbol_to_short = {} - for dom in all_domains: - # Read domain name from geometry - domain_symbol = list(geometry[dom[0]].keys())[0] - if len(dom) > 1: - domain_symbol = domain_symbol[0] - # For multi-domain variables keep only the first letter of the domain - domain_name_to_symbol[tuple(dom)] = domain_symbol - # Record which domain symbols we shortened - for d in dom: - long = list(geometry[d].keys())[0] - long_domain_symbol_to_short[long] = domain_symbol - else: - # Otherwise keep the whole domain - domain_name_to_symbol[tuple(dom)] = domain_symbol - - # Read domain limits - domain_name_to_limits = {(): None} - for dom in all_domains: - limits = list(geometry[dom[0]].values())[0].values() - if len(limits) > 1: - lower_limit, _ = list(geometry[dom[0]].values())[0].values() - _, upper_limit = list(geometry[dom[-1]].values())[0].values() - domain_name_to_limits[tuple(dom)] = ( - lower_limit.evaluate(), - upper_limit.evaluate(), - ) - else: - # Don't record limits for variables that have "limits" of length 1 i.e. - # a zero-dimensional domain - domain_name_to_limits[tuple(dom)] = None - - # Define independent variables for each variable - var_to_ind_vars = {} - var_to_ind_vars_left_boundary = {} - var_to_ind_vars_right_boundary = {} - for var in variables: - if var.domain in [[], ["current collector"]]: - var_to_ind_vars[var] = "(t)" - else: - # all independent variables e.g. (t, x) or (t, rn, xn) - domain_symbols = ", ".join( - domain_name_to_symbol[tuple(dom)] - for dom in var.domains.values() - if domain_name_to_limits[tuple(dom)] is not None - ) - var_to_ind_vars[var] = f"(t, {domain_symbols})" - if isinstance(var, pybamm.ConcatenationVariable): - for child in var.children: - var_to_ind_vars[child] = f"(t, {domain_symbols})" - aux_domain_symbols = ", ".join( - domain_name_to_symbol[tuple(dom)] - for level, dom in var.domains.items() - if level != "primary" and domain_name_to_limits[tuple(dom)] is not None - ) - if aux_domain_symbols != "": - aux_domain_symbols = ", " + aux_domain_symbols - - limits = domain_name_to_limits[tuple(var.domain)] - # left bc e.g. (t, 0) or (t, 0, xn) - var_to_ind_vars_left_boundary[var] = f"(t, {limits[0]}{aux_domain_symbols})" - # right bc e.g. (t, 1) or (t, 1, xn) - var_to_ind_vars_right_boundary[ - var - ] = f"(t, {limits[1]}{aux_domain_symbols})" - - mtk_str = "begin\n" - # Define parameters (including independent variables) - # Makes a line of the form '@parameters t x1 x2 x3 a b c d' - ind_vars = ["t"] + [ - sym - for dom, sym in domain_name_to_symbol.items() - if domain_name_to_limits[dom] is not None - ] - for domain_name, domain_symbol in domain_name_to_symbol.items(): - if domain_name_to_limits[domain_name] is not None: - mtk_str += f"# {domain_name} -> {domain_symbol}\n" - mtk_str += "@parameters " + " ".join(ind_vars) - if len(model.input_parameters) > 0: - mtk_str += "\n# Input parameters\n@parameters" - for param in model.input_parameters: - mtk_str += f" {param.name}" - mtk_str += "\n" - - # Add a comment with the variable names - for var in variables: - mtk_str += f"# '{var.name}' -> {variable_to_print_name[var]}\n" - # Makes a line of the form '@variables u1(t) u2(t)' - dep_vars = [] - mtk_str += "@variables" - for var in variables: - mtk_str += f" {variable_to_print_name[var]}(..)" - dep_var = variable_to_print_name[var] + var_to_ind_vars[var] - dep_vars.append(dep_var) - mtk_str += "\n" - - # Define derivatives - for domain_symbol in ind_vars: - mtk_str += f"D{domain_symbol} = Differential({domain_symbol})\n" - mtk_str += "\n" - - # Define equations - all_eqns_str = "" - all_constants_str = "" - all_julia_str = "" - for var, eqn in {**model.rhs, **model.algebraic}.items(): - all_constants_str, all_julia_str, eqn_str = convert_var_and_eqn_to_str( - var, eqn, all_constants_str, all_julia_str, "equation" - ) - - if var in model.rhs: - all_eqns_str += ( - f" Dt({variable_to_print_name[var]}{var_to_ind_vars[var]}) " - + f"~ {eqn_str},\n" - ) - elif var in model.algebraic: - all_eqns_str += f" 0 ~ {eqn_str},\n" - - # Replace any long domain symbols with the short version - # e.g. "xn" gets replaced with "x" - for long, short in long_domain_symbol_to_short.items(): - # we need to add a space to avoid accidentally replacing 'exp' with 'ex' - all_julia_str = all_julia_str.replace(" " + long, " " + short) - - # Replace variables in the julia strings that correspond to pybamm variables with - # their julia equivalent - for var, julia_id in variable_to_print_name.items(): - # e.g. boundary_value_right(cache_123456789) gets replaced with u1(t, 1) - cache_var_id = id_to_julia_variable(var.id, "cache") - if f"boundary_value_right({cache_var_id})" in all_julia_str: - all_julia_str = all_julia_str.replace( - f"boundary_value_right({cache_var_id})", - julia_id + var_to_ind_vars_right_boundary[var], - ) - # e.g. cache_123456789 gets replaced with u1(t, x) - all_julia_str = all_julia_str.replace( - cache_var_id, julia_id + var_to_ind_vars[var] - ) - - # Replace independent variables (domain names) in julia strings with the - # corresponding symbol - for domain_name, domain_symbol in domain_name_to_symbol.items(): - all_julia_str = all_julia_str.replace( - f"grad_{domain_name}", f"D{domain_symbol}" - ) - # Different divergence depending on the coordinate system - coord_sys = getattr(pybamm.standard_spatial_vars, domain_symbol).coord_sys - if coord_sys == "cartesian": - all_julia_str = all_julia_str.replace( - f"div_{domain_name}", f"D{domain_symbol}" - ) - elif coord_sys == "spherical polar": - all_julia_str = all_julia_str.replace( - f"div_{domain_name}(", - f"1 / {domain_symbol}^2 * D{domain_symbol}({domain_symbol}^2 * ", - ) - - # Replace the thicknesses in the concatenation with the actual thickness from the - # geometry - if "neg_width" in all_julia_str or "neg_plus_sep_width" in all_julia_str: - var = pybamm.standard_spatial_vars - x_n = geometry["negative electrode"]["x_n"]["max"].evaluate() - x_s = geometry["separator"]["x_s"]["max"].evaluate() - all_julia_str = all_julia_str.replace("neg_width", str(x_n)) - all_julia_str = all_julia_str.replace("neg_plus_sep_width", str(x_s)) - - # Update the MTK string - mtk_str += all_constants_str + all_julia_str + "\n" + f"eqs = [\n{all_eqns_str}]\n" - - #################################################################################### - # Initial and boundary conditions - #################################################################################### - # Initial conditions - all_ic_bc_str = " # initial conditions\n" - all_ic_bc_constants_str = "" - all_ic_bc_julia_str = "" - for var, eqn in model.initial_conditions.items(): - ( - all_ic_bc_constants_str, - all_ic_bc_julia_str, - eqn_str, - ) = convert_var_and_eqn_to_str( - var, eqn, all_ic_bc_constants_str, all_ic_bc_julia_str, "initial condition" - ) - - if not is_pde: - all_ic_bc_str += f" {variable_to_print_name[var]}(t) => {eqn_str},\n" - else: - if var.domain == []: - doms = "" - else: - doms = ", " + domain_name_to_symbol[tuple(var.domain)] - - all_ic_bc_str += f" {variable_to_print_name[var]}(0{doms}) ~ {eqn_str},\n" - # Boundary conditions - if is_pde: - all_ic_bc_str += " # boundary conditions\n" - for var, eqn_side in model.boundary_conditions.items(): - if isinstance(var, (pybamm.Variable, pybamm.ConcatenationVariable)): - for side, (eqn, typ) in eqn_side.items(): - ( - all_ic_bc_constants_str, - all_ic_bc_julia_str, - eqn_str, - ) = convert_var_and_eqn_to_str( - var, - eqn, - all_ic_bc_constants_str, - all_ic_bc_julia_str, - "boundary condition", - ) - - if side == "left": - limit = var_to_ind_vars_left_boundary[var] - elif side == "right": - limit = var_to_ind_vars_right_boundary[var] - - bc = f"{variable_to_print_name[var]}{limit}" - if typ == "Dirichlet": - bc = bc - elif typ == "Neumann": - bc = f"D{domain_name_to_symbol[tuple(var.domain)]}({bc})" - all_ic_bc_str += f" {bc} ~ {eqn_str},\n" - - # Replace variables in the julia strings that correspond to pybamm variables with - # their julia equivalent - for var, julia_id in variable_to_print_name.items(): - # e.g. boundary_value_right(cache_123456789) gets replaced with u1(t, 1) - cache_var_id = id_to_julia_variable(var.id, "cache") - if f"boundary_value_right({cache_var_id})" in all_ic_bc_julia_str: - all_ic_bc_julia_str = all_ic_bc_julia_str.replace( - f"boundary_value_right({cache_var_id})", - julia_id + var_to_ind_vars_right_boundary[var], - ) - # e.g. cache_123456789 gets replaced with u1(t, x) - all_ic_bc_julia_str = all_ic_bc_julia_str.replace( - cache_var_id, julia_id + var_to_ind_vars[var] - ) - - #################################################################################### - - # Create ODESystem or PDESystem - if not is_pde: - mtk_str += "sys = ODESystem(eqs, t)\n\n" - mtk_str += ( - all_ic_bc_constants_str - + all_ic_bc_julia_str - + "\n" - + f"u0 = [\n{all_ic_bc_str}]\n" - ) - else: - # Initial and boundary conditions - mtk_str += ( - all_ic_bc_constants_str - + all_ic_bc_julia_str - + "\n" - + f"ics_bcs = [\n{all_ic_bc_str}]\n" - ) - - # Domains - mtk_str += "\n" - tpsan_str = ",".join( - map(lambda x: f"{x / model.timescale.evaluate():.3f}", tspan) - ) - mtk_str += f"t_domain = Interval({tpsan_str})\n" - domains = "domains = [\n t in t_domain,\n" - for domain, symbol in domain_name_to_symbol.items(): - limits = domain_name_to_limits[tuple(domain)] - if limits is not None: - mtk_str += f"{symbol}_domain = Interval{limits}\n" - domains += f" {symbol} in {symbol}_domain,\n" - domains += "]\n" - - mtk_str += "\n" - mtk_str += domains - - # Independent and dependent variables - mtk_str += "ind_vars = [{}]\n".format(", ".join(ind_vars)) - mtk_str += "dep_vars = [{}]\n\n".format(", ".join(dep_vars)) - - name = model.name.replace(" ", "_").replace("-", "_") - mtk_str += ( - name - + "_pde_system = PDESystem(eqs, ics_bcs, domains, ind_vars, dep_vars)\n\n" - ) - - # Replace parameters in the julia strings in the form "inputs[name]" - # with just "name" - for param in model.input_parameters: - mtk_str = mtk_str.replace(f"inputs['{param.name}']", param.name) - - # Need to add 'nothing' to the end of the mtk string to avoid errors in MTK v4 - # See https://github.com/SciML/diffeqpy/issues/82 - mtk_str += "nothing\nend\n" - - return mtk_str diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py index a84f84dd22..76d294cf2d 100644 --- a/pybamm/models/base_model.py +++ b/pybamm/models/base_model.py @@ -1079,83 +1079,6 @@ def generate( C.add(variables_fn) C.generate() - def generate_julia_diffeq( - self, - input_parameter_order=None, - get_consistent_ics_solver=None, - dae_type="semi-explicit", - **kwargs, - ): - """ - Generate a Julia representation of the model, ready to be solved by Julia's - DifferentialEquations library. - - Parameters - ---------- - input_parameter_order : list, optional - Order in which input parameters will be provided when solving the model - get_consistent_ics_solver : pybamm solver, optional - Solver to use to get consistent initial conditions. If None, the initial - guesses for boundary conditions (non-consistent) are used. - dae_type : str, optional - How to write the DAEs. Options are "semi-explicit" (default) or "implicit". - - Returns - ------- - eqn_str : str - The Julia-compatible equations for the model in string format, - to be evaluated by eval(Meta.parse(...)) - ics_str : str - The Julia-compatible initial conditions for the model in string format, - to be evaluated by eval(Meta.parse(...)) - """ - self.check_discretised_or_discretise_inplace_if_0D() - - name = self.name.replace(" ", "_") - - if self.algebraic == {}: - # ODE model: form dy[] = ... - eqn_str = pybamm.get_julia_function( - self.concatenated_rhs, - funcname=name, - input_parameter_order=input_parameter_order, - **kwargs, - ) - else: - if dae_type == "semi-explicit": - len_rhs = None - else: - len_rhs = self.concatenated_rhs.size - # DAE model: form out[] = ... - dy[] - eqn_str = pybamm.get_julia_function( - pybamm.numpy_concatenation( - self.concatenated_rhs, self.concatenated_algebraic - ), - funcname=name, - input_parameter_order=input_parameter_order, - len_rhs=len_rhs, - **kwargs, - ) - - if get_consistent_ics_solver is None or self.algebraic == {}: - ics = self.concatenated_initial_conditions - else: - get_consistent_ics_solver.set_up(self) - get_consistent_ics_solver._set_initial_conditions(self, {}, False) - ics = pybamm.Vector(self.y0.full()) - - ics_str = pybamm.get_julia_function( - ics, - funcname=name + "_u0", - input_parameter_order=input_parameter_order, - **kwargs, - ) - # Change the string to a form for u0 - ics_str = ics_str.replace("(dy, y, p, t)", "(u0, p)") - ics_str = ics_str.replace("dy", "u0") - - return eqn_str, ics_str - def latexify(self, filename=None, newline=True): # For docstring, see pybamm.expression_tree.operations.latexify.Latexify return Latexify(self, filename, newline).latexify() diff --git a/requirements.txt b/requirements.txt index 75431982d9..0988a947dc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,7 +6,6 @@ autograd >= 1.2 scikit-fem >= 0.2.0 casadi >= 3.5.0 imageio>=2.9.0 -julia>=0.5.6 jupyter # For example notebooks pybtex>=0.24.0 sympy >= 1.8 diff --git a/setup.py b/setup.py index af904285da..e044444a7b 100644 --- a/setup.py +++ b/setup.py @@ -197,9 +197,6 @@ def compile_KLU(): "scikit-fem>=0.2.0", "casadi>=3.5.0", "imageio>=2.9.0", - # Julia pip packaged can be installed even if - # julia programming language is not installed - "julia>=0.5.6", "jupyter", # For example notebooks "pybtex>=0.24.0", "sympy>=1.8", diff --git a/tests/unit/test_expression_tree/test_operations/test_evaluate_julia.py b/tests/unit/test_expression_tree/test_operations/test_evaluate_julia.py deleted file mode 100644 index eee960e069..0000000000 --- a/tests/unit/test_expression_tree/test_operations/test_evaluate_julia.py +++ /dev/null @@ -1,366 +0,0 @@ -# -# Test for the evaluate-to-Julia functions -# -import pybamm - -from tests import get_mesh_for_testing, get_1p1d_mesh_for_testing -import unittest -import numpy as np -import scipy.sparse -from platform import system - -have_julia = pybamm.have_julia() -if have_julia and system() != "Windows": - from julia.api import Julia - - Julia(compiled_modules=False) - from julia import Main - from julia.core import JuliaError - - # load julia libraries required for evaluating the strings - Main.eval("using SparseArrays, LinearAlgebra") - - -@unittest.skipIf(not have_julia, "Julia not installed") -@unittest.skipIf(system() == "Windows", "Julia not supported on windows") -class TestEvaluate(unittest.TestCase): - def evaluate_and_test_equal( - self, expr, y_tests, t_tests=0.0, inputs=None, decimal=14, **kwargs - ): - if not isinstance(y_tests, list): - y_tests = [y_tests] - if not isinstance(t_tests, list): - t_tests = [t_tests] - if inputs is None: - input_parameter_order = None - p = 0.0 - else: - input_parameter_order = list(inputs.keys()) - p = list(inputs.values()) - - pybamm_eval = expr.evaluate(t=t_tests[0], y=y_tests[0], inputs=inputs).flatten() - for preallocate in [True, False]: - kwargs["funcname"] = ( - kwargs.get("funcname", "f") + "_" + str(int(preallocate)) - ) - evaluator_str = pybamm.get_julia_function( - expr, - input_parameter_order=input_parameter_order, - preallocate=preallocate, - **kwargs, - ) - Main.eval(evaluator_str) - funcname = kwargs.get("funcname", "f") - Main.p = p - for t_test, y_test in zip(t_tests, y_tests): - Main.dy = np.zeros_like(pybamm_eval) - Main.y = y_test - Main.t = t_test - try: - Main.eval(f"{funcname}!(dy,y,p,t)") - except JuliaError as e: - # debugging - print(Main.dy, y_test, p, t_test) - print(evaluator_str) - raise e - pybamm_eval = expr.evaluate(t=t_test, y=y_test, inputs=inputs).flatten() - try: - np.testing.assert_array_almost_equal( - Main.dy.flatten(), - pybamm_eval, - decimal=decimal, - ) - except AssertionError as e: - # debugging - print(Main.dy, y_test, p, t_test) - print(evaluator_str) - raise e - - def test_exceptions(self): - a = pybamm.Symbol("a") - with self.assertRaisesRegex(NotImplementedError, "Conversion to Julia"): - pybamm.get_julia_function(a) - - def test_evaluator_julia(self): - a = pybamm.StateVector(slice(0, 1)) - b = pybamm.StateVector(slice(1, 2)) - - y_tests = [np.array([[2], [3]]), np.array([[1], [3]])] - t_tests = [1, 2] - - # test a * b - expr = a * b - self.evaluate_and_test_equal(expr, np.array([2.0, 3.0])) - self.evaluate_and_test_equal(expr, np.array([1.0, 3.0])) - - # test function(a*b) - expr = pybamm.cos(a * b) - self.evaluate_and_test_equal(expr, np.array([1.0, 3.0]), funcname="g") - - # test a constant expression - expr = pybamm.Multiplication(pybamm.Scalar(2), pybamm.Scalar(3)) - self.evaluate_and_test_equal(expr, 0.0) - - expr = pybamm.Multiplication(pybamm.Scalar(2), pybamm.Vector([1, 2, 3])) - self.evaluate_and_test_equal(expr, None, funcname="g2") - - # test a larger expression - expr = a * b + b + a**2 / b + 2 * a + b / 2 + 4 - self.evaluate_and_test_equal(expr, y_tests) - - # test something with time - expr = a * pybamm.t - self.evaluate_and_test_equal(expr, y_tests, t_tests=t_tests) - - # test something with a matrix multiplication - A = pybamm.Matrix([[1, 2], [3, 4]]) - expr = A @ pybamm.StateVector(slice(0, 2)) - self.evaluate_and_test_equal(expr, y_tests, funcname="g3") - - # test something with a heaviside - a = pybamm.Vector([1, 2]) - expr = a <= pybamm.StateVector(slice(0, 2)) - self.evaluate_and_test_equal(expr, y_tests, funcname="g4") - - # test something with a minimum or maximum - a = pybamm.Vector([1, 2]) - for expr in [ - pybamm.minimum(a, pybamm.StateVector(slice(0, 2))), - pybamm.maximum(a, pybamm.StateVector(slice(0, 2))), - ]: - self.evaluate_and_test_equal(expr, y_tests, funcname="g5") - - # test something with an index - expr = pybamm.Index(A @ pybamm.StateVector(slice(0, 2)), 0) - self.evaluate_and_test_equal(expr, y_tests, funcname="g6") - - # test something with a sparse matrix multiplication - A = pybamm.Matrix([[1, 2], [3, 4]]) - B = pybamm.Matrix(scipy.sparse.csr_matrix(np.array([[1, 0], [0, 4]]))) - expr = A @ B @ pybamm.StateVector(slice(0, 2)) - self.evaluate_and_test_equal(expr, y_tests, funcname="g7") - - expr = B @ pybamm.StateVector(slice(0, 2)) - self.evaluate_and_test_equal(expr, y_tests, funcname="g8") - - # test Inner - expr = pybamm.Inner(pybamm.Vector([1, 2]), pybamm.StateVector(slice(0, 2))) - self.evaluate_and_test_equal(expr, y_tests, funcname="g12") - - # test numpy concatenation - a = pybamm.StateVector(slice(0, 3)) - b = pybamm.Vector([2, 3, 4]) - c = pybamm.Vector([5]) - - y_tests = [np.array([[2], [3], [4]]), np.array([[1], [3], [2]])] - - expr = pybamm.NumpyConcatenation(a, b, c) - self.evaluate_and_test_equal(expr, y_tests, funcname="g9") - - expr = pybamm.NumpyConcatenation(a, c) - self.evaluate_and_test_equal(expr, y_tests, funcname="g10") - - # test sparse stack - A = pybamm.Matrix(scipy.sparse.csr_matrix(np.array([[1, 0], [0, 4]]))) - B = pybamm.Matrix(scipy.sparse.csr_matrix(np.array([[2, 0], [5, 0]]))) - c = pybamm.StateVector(slice(0, 2)) - expr = pybamm.SparseStack(A, B) @ c - self.evaluate_and_test_equal(expr, y_tests, funcname="g11") - - def test_evaluator_julia_input_parameters(self): - a = pybamm.StateVector(slice(0, 1)) - b = pybamm.StateVector(slice(1, 2)) - c = pybamm.InputParameter("c") - d = pybamm.InputParameter("d") - - # test one input parameter - expr = a * c - self.evaluate_and_test_equal(expr, np.array([2.0, 3.0]), inputs={"c": 5}) - - # test several input parameters - expr = a * c + b * d - self.evaluate_and_test_equal( - expr, np.array([2.0, 3.0]), inputs={"c": 5, "d": 6} - ) - - def test_evaluator_julia_all_functions(self): - a = pybamm.StateVector(slice(0, 3)) - y_test = np.array([1, 2, 3]) - - for function in [ - pybamm.arcsinh, - pybamm.cos, - pybamm.cosh, - pybamm.exp, - pybamm.log, - pybamm.log10, - pybamm.sin, - pybamm.sinh, - pybamm.sqrt, - pybamm.tanh, - pybamm.arctan, - ]: - expr = function(a) - self.evaluate_and_test_equal(expr, y_test) - - for function in [ - pybamm.min, - pybamm.max, - ]: - expr = function(a) - self.evaluate_and_test_equal(expr, y_test) - - # More advanced tests for min - b = pybamm.StateVector(slice(3, 6)) - concat = pybamm.NumpyConcatenation(2 * a, 3 * b) - expr = pybamm.min(concat) - self.evaluate_and_test_equal(expr, np.array([1, 2, 3, 4, 5, 6]), funcname="h1") - - v = pybamm.Vector([1, 2, 3]) - expr = pybamm.min(v * a) - self.evaluate_and_test_equal(expr, y_test, funcname="h2") - - def test_evaluator_julia_domain_concatenation(self): - c_n = pybamm.Variable("c_n", domain="negative electrode") - c_s = pybamm.Variable("c_s", domain="separator") - c_p = pybamm.Variable("c_p", domain="positive electrode") - c = pybamm.concatenation(c_n / 2, c_s / 3, c_p / 4) - # create discretisation - mesh = get_mesh_for_testing() - spatial_methods = { - "macroscale": pybamm.FiniteVolume(), - "current collector": pybamm.FiniteVolume(), - } - disc = pybamm.Discretisation(mesh, spatial_methods) - - combined_submesh = mesh.combine_submeshes(*c.domain) - nodes = combined_submesh.nodes - y_tests = [nodes**2 + 1, np.cos(nodes)] - - # discretise and evaluate the variable - disc.set_variable_slices([c_n, c_s, c_p]) - c_disc = disc.process_symbol(c) - self.evaluate_and_test_equal(c_disc, y_tests) - - def test_evaluator_julia_domain_concatenation_2D(self): - c_n = pybamm.Variable( - "c_n", - domain="negative electrode", - auxiliary_domains={"secondary": "current collector"}, - ) - c_s = pybamm.Variable( - "c_s", - domain="separator", - auxiliary_domains={"secondary": "current collector"}, - ) - c_p = pybamm.Variable( - "c_p", - domain="positive electrode", - auxiliary_domains={"secondary": "current collector"}, - ) - c = pybamm.concatenation(c_n / 2, c_s / 3, c_p / 4) - # create discretisation - mesh = get_1p1d_mesh_for_testing() - spatial_methods = {"macroscale": pybamm.FiniteVolume()} - disc = pybamm.Discretisation(mesh, spatial_methods) - - combined_submesh = mesh.combine_submeshes(*c.domain) - nodes = np.linspace( - 0, 1, combined_submesh.npts * mesh["current collector"].npts - ) - y_tests = [nodes**2 + 1, np.cos(nodes)] - - # discretise and evaluate the variable - disc.set_variable_slices([c_n, c_s, c_p]) - c_disc = disc.process_symbol(c) - - self.evaluate_and_test_equal(c_disc, y_tests) - - def test_evaluator_julia_discretised_operators(self): - whole_cell = ["negative electrode", "separator", "positive electrode"] - # create discretisation - mesh = get_mesh_for_testing() - spatial_methods = {"macroscale": pybamm.FiniteVolume()} - disc = pybamm.Discretisation(mesh, spatial_methods) - - combined_submesh = mesh.combine_submeshes(*whole_cell) - - var = pybamm.Variable("var", domain=whole_cell) - boundary_conditions = { - var: { - "left": (pybamm.Scalar(1), "Dirichlet"), - "right": (pybamm.Scalar(2), "Neumann"), - } - } - disc.bcs = boundary_conditions - disc.set_variable_slices([var]) - - # grad - grad_eqn = pybamm.grad(var) - grad_eqn_disc = disc.process_symbol(grad_eqn) - - # div - div_eqn = pybamm.div(var * grad_eqn) - div_eqn_disc = disc.process_symbol(div_eqn) - - # test - nodes = combined_submesh.nodes - y_tests = [nodes**2 + 1, np.cos(nodes)] - - for i, expr in enumerate([grad_eqn_disc, div_eqn_disc]): - self.evaluate_and_test_equal(expr, y_tests, funcname=f"f{i}", decimal=10) - - def test_evaluator_julia_discretised_microscale(self): - # create discretisation - mesh = get_1p1d_mesh_for_testing(xpts=5, rpts=5, zpts=2) - spatial_methods = {"negative particle": pybamm.FiniteVolume()} - disc = pybamm.Discretisation(mesh, spatial_methods) - - submesh = mesh["negative particle"] - - # grad - # grad(r) == 1 - var = pybamm.Variable( - "var", - domain=["negative particle"], - auxiliary_domains={ - "secondary": "negative electrode", - "tertiary": "current collector", - }, - ) - grad_eqn = pybamm.grad(var) - div_eqn = pybamm.div(var * grad_eqn) - - boundary_conditions = { - var: { - "left": (pybamm.Scalar(1), "Dirichlet"), - "right": (pybamm.Scalar(2), "Neumann"), - } - } - - disc.bcs = boundary_conditions - - disc.set_variable_slices([var]) - grad_eqn_disc = disc.process_symbol(grad_eqn) - div_eqn_disc = disc.process_symbol(div_eqn) - - # test - total_npts = ( - submesh.npts - * mesh["negative electrode"].npts - * mesh["current collector"].npts - ) - y_tests = [np.linspace(0, 1, total_npts) ** 2] - - for i, expr in enumerate([grad_eqn_disc, div_eqn_disc]): - self.evaluate_and_test_equal(expr, y_tests, funcname=f"f{i}", decimal=11) - - -if __name__ == "__main__": - print("Add -v for more debug output") - import sys - - if "-v" in sys.argv: - debug = True - pybamm.settings.debug_mode = True - unittest.main() diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py index 132e9826fe..7a5f5777c7 100644 --- a/tests/unit/test_models/test_base_model.py +++ b/tests/unit/test_models/test_base_model.py @@ -586,42 +586,6 @@ def test_generate_casadi(self): os.remove("test.c") os.remove("test.so") - @unittest.skipIf(platform.system() == "Windows", "Skipped for Windows") - def test_generate_julia_diffeq(self): - # ODE model with no input parameters - model = pybamm.BaseModel(name="ode test model") - a = pybamm.Variable("a") - b = pybamm.Variable("b") - model.rhs = {a: -a, b: a - b} - model.initial_conditions = {a: 1, b: 2} - - # Generate rhs and ics for the Julia model - rhs_str, ics_str = model.generate_julia_diffeq() - self.assertIsInstance(rhs_str, str) - self.assertIn("ode_test_model", rhs_str) - self.assertIsInstance(ics_str, str) - self.assertIn("ode_test_model_u0", ics_str) - self.assertIn("(u0, p)", ics_str) - - # ODE model with input parameters - model = pybamm.BaseModel(name="ode test model 2") - a = pybamm.Variable("a") - b = pybamm.Variable("b") - p = pybamm.InputParameter("p") - q = pybamm.InputParameter("q") - model.rhs = {a: -a * p, b: a - b} - model.initial_conditions = {a: q, b: 2} - - # Generate rhs and ics for the Julia model - rhs_str, ics_str = model.generate_julia_diffeq(input_parameter_order=["p", "q"]) - self.assertIsInstance(rhs_str, str) - self.assertIn("ode_test_model_2", rhs_str) - self.assertIn("p, q = p", rhs_str) - - self.assertIsInstance(ics_str, str) - self.assertIn("ode_test_model_2_u0", ics_str) - self.assertIn("p, q = p", ics_str) - def test_set_initial_conditions(self): # Set up model model = pybamm.BaseModel() diff --git a/tests/unit/test_models/test_base_model_generate_julia_diffeq.py b/tests/unit/test_models/test_base_model_generate_julia_diffeq.py deleted file mode 100644 index b331b6b40a..0000000000 --- a/tests/unit/test_models/test_base_model_generate_julia_diffeq.py +++ /dev/null @@ -1,131 +0,0 @@ -# -# Tests for the base model class -# -import platform -import unittest -import pybamm - -have_julia = pybamm.have_julia() -if have_julia and platform.system() != "Windows": - from julia.api import Julia - - Julia(compiled_modules=False) - from julia import Main - - # load julia libraries required for evaluating the strings - Main.eval("using SparseArrays, LinearAlgebra") - - -@unittest.skipIf(not have_julia, "Julia not installed") -class TestBaseModelGenerateJuliaDiffEq(unittest.TestCase): - def test_generate_ode(self): - # ODE model with no input parameters - model = pybamm.BaseModel(name="ode test model") - a = pybamm.Variable("a") - b = pybamm.Variable("b") - model.rhs = {a: -a, b: a - b} - model.initial_conditions = {a: 1, b: 2} - - # Generate rhs and ics for the Julia model - rhs_str, ics_str = model.generate_julia_diffeq() - self.assertIsInstance(rhs_str, str) - self.assertIn("ode_test_model", rhs_str) - self.assertIn("(dy, y, p, t)", rhs_str) - self.assertIsInstance(ics_str, str) - self.assertIn("ode_test_model_u0", ics_str) - self.assertIn("(u0, p)", ics_str) - - # ODE model with input parameters - model = pybamm.BaseModel(name="ode test model 2") - a = pybamm.Variable("a") - b = pybamm.Variable("b") - p = pybamm.InputParameter("p") - q = pybamm.InputParameter("q") - model.rhs = {a: -a * p, b: a - b} - model.initial_conditions = {a: q, b: 2} - - # Generate rhs and ics for the Julia model - rhs_str, ics_str = model.generate_julia_diffeq(input_parameter_order=["p", "q"]) - self.assertIsInstance(rhs_str, str) - self.assertIn("ode_test_model_2", rhs_str) - self.assertIn("p, q = p", rhs_str) - - self.assertIsInstance(ics_str, str) - self.assertIn("ode_test_model_2_u0", ics_str) - self.assertIn("p, q = p", ics_str) - - def test_generate_dae(self): - # ODE model with no input parameters - model = pybamm.BaseModel(name="dae test model") - a = pybamm.Variable("a") - b = pybamm.Variable("b") - model.rhs = {a: -a} - model.algebraic = {b: a - b} - model.initial_conditions = {a: 1, b: 2} - - # Generate eqn and ics for the Julia model (semi-explicit) - eqn_str, ics_str = model.generate_julia_diffeq() - self.assertIsInstance(eqn_str, str) - self.assertIn("dae_test_model", eqn_str) - self.assertIn("(dy, y, p, t)", eqn_str) - self.assertIsInstance(ics_str, str) - self.assertIn("dae_test_model_u0", ics_str) - self.assertIn("(u0, p)", ics_str) - self.assertIn("[1.,2.]", ics_str) - - # Generate eqn and ics for the Julia model (implicit) - eqn_str, ics_str = model.generate_julia_diffeq(dae_type="implicit") - self.assertIsInstance(eqn_str, str) - self.assertIn("dae_test_model", eqn_str) - self.assertIn("(out, dy, y, p, t)", eqn_str) - self.assertIsInstance(ics_str, str) - self.assertIn("dae_test_model_u0", ics_str) - self.assertIn("(u0, p)", ics_str) - - # Calculate initial conditions in python - eqn_str, ics_str = model.generate_julia_diffeq( - get_consistent_ics_solver=pybamm.CasadiSolver() - ) - # Check that the initial conditions are consistent - self.assertIn("[1.,1.]", ics_str) - - def test_generate_pde(self): - # ODE model with no input parameters - model = pybamm.BaseModel(name="pde test model") - a = pybamm.Variable("a", domain="line") - b = pybamm.Variable("b", domain="line") - model.rhs = {a: pybamm.div(pybamm.grad(a)) + b, b: a - b} - model.boundary_conditions = { - a: {"left": (-1, "Dirichlet"), "right": (1, "Neumann")} - } - model.initial_conditions = {a: 1, b: 2} - - # Discretize - x = pybamm.SpatialVariable("x", domain=["line"]) - geometry = pybamm.Geometry( - {"line": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}} - ) - submesh_types = {"line": pybamm.Uniform1DSubMesh} - var_pts = {x: 10} - mesh = pybamm.Mesh(geometry, submesh_types, var_pts) - disc = pybamm.Discretisation(mesh, {"line": pybamm.FiniteVolume()}) - disc.process_model(model) - - # Generate rhs and ics for the Julia model - rhs_str, ics_str = model.generate_julia_diffeq() - self.assertIsInstance(rhs_str, str) - self.assertIn("pde_test_model", rhs_str) - self.assertIn("(dy, y, p, t)", rhs_str) - self.assertIsInstance(ics_str, str) - self.assertIn("pde_test_model_u0", ics_str) - self.assertIn("(u0, p)", ics_str) - - -if __name__ == "__main__": - print("Add -v for more debug output") - import sys - - if "-v" in sys.argv: - debug = True - pybamm.settings.debug_mode = True - unittest.main() diff --git a/tests/unit/test_solvers/test_julia_mtk.py b/tests/unit/test_solvers/test_julia_mtk.py deleted file mode 100644 index 0e8f0838f2..0000000000 --- a/tests/unit/test_solvers/test_julia_mtk.py +++ /dev/null @@ -1,240 +0,0 @@ -# -# Test for the evaluate-to-Julia functions -# -import pybamm - -import unittest -from platform import system - - -# julia imports -have_julia = pybamm.have_julia() - - -@unittest.skipIf(not have_julia, "Julia not installed") -@unittest.skipIf(system() == "Windows", "Julia not supported on windows") -class TestCreateSolveMTKModel(unittest.TestCase): - """ - These tests just make sure there are no errors when calling - pybamm.get_julia_mtk_model. TODO: add (comment out) tests that run and solve the - model. This needs (i) faster import of diffeqpy, (ii) working PDE discretisations - in Julia. - """ - - def test_exponential_decay_model(self): - model = pybamm.BaseModel() - v = pybamm.Variable("v") - model.rhs = {v: -2 * v} - model.initial_conditions = {v: 0.5} - - pybamm.get_julia_mtk_model(model) - - # Main.eval("using ModelingToolkit") - # Main.eval(mtk_str) - - # Main.tspan = (0.0, 10.0) - # # this definition of prob doesn't work, so we use Main.eval instead - # # prob = de.ODEProblem(Main.sys, Main.u0, Main.tspan) - - # Main.eval("prob = ODEProblem(sys, u0, tspan)") - # sol = de.solve(Main.prob, de.Tsit5()) - - # y_sol = np.concatenate(sol.u) - # y_exact = 0.5 * np.exp(-2 * sol.t) - # np.testing.assert_almost_equal(y_sol, y_exact, decimal=6) - - def test_lotka_volterra_model(self): - model = pybamm.BaseModel() - a = pybamm.InputParameter("a") - b = pybamm.InputParameter("b") - c = pybamm.InputParameter("c") - d = pybamm.InputParameter("d") - x = pybamm.Variable("x") - y = pybamm.Variable("y") - - model.rhs = {x: a * x - b * x * y, y: c * x * y - d * y} - model.initial_conditions = {x: 1.0, y: 1.0} - - pybamm.get_julia_mtk_model(model) - - # # Solve using julia - # Main.eval("using ModelingToolkit") - # Main.eval(mtk_str) - - # Main.tspan = (0.0, 10.0) - # Main.eval( - # """ - # begin - # p = [a => 1.5, b => 1.0, c => 3.0, d => 1.0] - # prob = ODEProblem(sys, u0, tspan, p) - # end - # """ - # ) - # sol_julia = de.solve(Main.prob, de.Tsit5(), reltol=1e-8, abstol=1e-8) - - # y_sol_julia = np.vstack(sol_julia.u).T - - # # Solve using pybamm - # sol_pybamm = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8).solve( - # model, sol_julia.t, inputs={"a": 1.5, "b": 1.0, "c": 3.0, "d": 1.0} - # ) - - # # Compare - # np.testing.assert_almost_equal(y_sol_julia, sol_pybamm.y, decimal=5) - - def test_dae_model(self): - model = pybamm.BaseModel() - x = pybamm.Variable("x") - y = pybamm.Variable("y") - - model.rhs = {x: -2 * x} - model.algebraic = {y: x - y} - model.initial_conditions = {x: 1.0, y: 1.0} - - pybamm.get_julia_mtk_model(model) - - # # Solve using julia - # Main.eval("using ModelingToolkit") - # Main.eval(mtk_str) - - # Main.tspan = (0.0, 10.0) - # Main.eval("prob = ODEProblem(sys, u0, tspan)") - # sol_julia = de.solve(Main.prob, de.Rodas5(), reltol=1e-8, abstol=1e-8) - - # y_sol_julia = np.vstack(sol_julia.u).T - - # # Solve using pybamm - # sol_pybamm = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8).solve(model, - # sol_julia.t) - - # # Compare - # np.testing.assert_almost_equal(y_sol_julia, sol_pybamm.y, decimal=5) - - def test_pde_model(self): - model = pybamm.BaseModel() - var = pybamm.Variable("var", domain="line") - - model.rhs = {var: pybamm.div(pybamm.grad(var))} - model.initial_conditions = {var: 1.0} - model.boundary_conditions = { - var: {"left": (1, "Dirichlet"), "right": (1, "Dirichlet")} - } - - geometry = {"line": {"x": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}} - - pybamm.get_julia_mtk_model(model, geometry=geometry, tspan=(0.0, 10.0)) - - # # Solve using julia - # Main.eval("using ModelingToolkit, DiffEqOperators") - # Main.eval(mtk_str) - - # Main.tspan = (0.0, 10.0) - # # Method of lines discretization - # Main.dx = 0.1 - # Main.order = 2 - # Main.eval("discretization = MOLFiniteDifference(dx,order)") - - # # Convert the PDE problem into an ODE problem - # Main.eval("prob = DiffEqOperators.discretize(pde_system,discretization)") - - # # Solve PDE problem - # sol_julia = de.solve(Main.prob, de.Tsit5(), reltol=1e-8, abstol=1e-8) - - # y_sol_julia = np.hstack(sol_julia.u) - - # # Check everything is equal to 1 - # # Just a simple test for now to get started - # np.testing.assert_equal(y_sol_julia, 1) - - def test_pde_model_spherical_polar(self): - model = pybamm.BaseModel() - var = pybamm.Variable("var", domain="particle") - - model.rhs = {var: pybamm.div(pybamm.grad(var))} - model.initial_conditions = {var: 1.0} - model.boundary_conditions = { - var: {"left": (1, "Dirichlet"), "right": (1, "Dirichlet")} - } - - geometry = { - "particle": {"r_n": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}} - } - - pybamm.get_julia_mtk_model(model, geometry=geometry, tspan=(0.0, 10.0)) - - # # Solve using julia - # Main.eval("using ModelingToolkit, DiffEqOperators") - # Main.eval(mtk_str) - - # Main.tspan = (0.0, 10.0) - # # Method of lines discretization - # Main.dx = 0.1 - # Main.order = 2 - # Main.eval("discretization = MOLFiniteDifference(dx,order)") - - # # Convert the PDE problem into an ODE problem - # Main.eval("prob = DiffEqOperators.discretize(pde_system,discretization)") - - # # Solve PDE problem - # sol_julia = de.solve(Main.prob, de.Tsit5(), reltol=1e-8, abstol=1e-8) - - # y_sol_julia = np.hstack(sol_julia.u) - - # # Check everything is equal to 1 - # # Just a simple test for now to get started - # np.testing.assert_equal(y_sol_julia, 1) - - def test_spm(self): - model = pybamm.lithium_ion.SPM() - parameter_values = model.default_parameter_values - parameter_values._replace_callable_function_parameters = False - sim = pybamm.Simulation(model, parameter_values=parameter_values) - sim.set_parameters() - pybamm.get_julia_mtk_model( - sim._model_with_set_params, geometry=sim.geometry, tspan=(0, 3600) - ) - - def test_spme(self): - model = pybamm.lithium_ion.SPMe() - parameter_values = model.default_parameter_values - parameter_values._replace_callable_function_parameters = False - sim = pybamm.Simulation(model, parameter_values=parameter_values) - sim.set_parameters() - pybamm.get_julia_mtk_model( - sim._model_with_set_params, geometry=sim.geometry, tspan=(0, 3600) - ) - - def test_dfn(self): - model = pybamm.lithium_ion.DFN() - parameter_values = model.default_parameter_values - parameter_values._replace_callable_function_parameters = False - sim = pybamm.Simulation(model, parameter_values=parameter_values) - sim.set_parameters() - pybamm.get_julia_mtk_model( - sim._model_with_set_params, geometry=sim.geometry, tspan=(0, 3600) - ) - - def test_exceptions(self): - model = pybamm.lithium_ion.SPM() - parameter_values = model.default_parameter_values - parameter_values._replace_callable_function_parameters = False - sim = pybamm.Simulation(model, parameter_values=parameter_values) - sim.set_parameters() - with self.assertRaisesRegex(ValueError, "must provide geometry"): - pybamm.get_julia_mtk_model( - sim._model_with_set_params, geometry=None, tspan=(0, 3600) - ) - with self.assertRaisesRegex(ValueError, "must provide tspan"): - pybamm.get_julia_mtk_model( - sim._model_with_set_params, geometry=sim.geometry, tspan=None - ) - - -if __name__ == "__main__": - print("Add -v for more debug output") - import sys - - if "-v" in sys.argv: - debug = True - pybamm.settings.debug_mode = True - unittest.main() diff --git a/tox.ini b/tox.ini index 57bd58fc6d..974efc5a14 100644 --- a/tox.ini +++ b/tox.ini @@ -1,5 +1,5 @@ [tox] -envlist = {windows}-{tests,unit,dev},tests,unit,dev,julia +envlist = {windows}-{tests,unit,dev},tests,unit,dev [testenv] skipsdist = true @@ -27,16 +27,6 @@ commands = dev-!windows-!mac: sh -c "echo export LD_LIBRARY_PATH={env:LD_LIBRARY_PATH} >> {envbindir}/activate" doctests: python run-tests.py --doctest -[testenv:julia] -platform = [linux, darwin] -skip_install = true -passenv = HOME -whitelist_externals = git -deps = - julia -commands = - python -c "import julia; julia.install()" - [testenv:pybamm-requires] platform = [linux, darwin] skip_install = true From 80e0e2c10dfc4e898a2716fc512440da263c05ec Mon Sep 17 00:00:00 2001 From: Alec Bills Date: Sun, 13 Nov 2022 15:57:41 -0500 Subject: [PATCH 02/10] symbol replacer --- pybamm/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pybamm/__init__.py b/pybamm/__init__.py index 699984b162..7013fb5876 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -95,6 +95,7 @@ from .expression_tree.operations.jacobian import Jacobian from .expression_tree.operations.convert_to_casadi import CasadiConverter from .expression_tree.operations.unpack_symbols import SymbolUnpacker + # # Model classes # From b5bebdc259e105eb5b32d26f6ac948d813b53722 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Sun, 13 Nov 2022 12:29:41 -0500 Subject: [PATCH 03/10] remove have julia --- pybamm/util.py | 25 ------------------------- tests/unit/test_util.py | 11 ----------- 2 files changed, 36 deletions(-) diff --git a/pybamm/util.py b/pybamm/util.py index 5821757026..d9cc806c06 100644 --- a/pybamm/util.py +++ b/pybamm/util.py @@ -288,31 +288,6 @@ def get_parameters_filepath(path): return os.path.join(pybamm.__path__[0], path) -def have_julia(): - """ - Checks whether the Julia programming language has been installed - """ - - # Try fetching info about julia - try: - info = JuliaInfo.load() - except (FileNotFoundError, subprocess.CalledProcessError): - return False - - # Compatibility: Checks - if not info.is_pycall_built(): # pragma: no cover - return False - if not info.is_compatible_python(): # pragma: no cover - return False - - # Confirm Julia() is callable - try: - Julia() - return True - except JuliaError: # pragma: no cover - return False - - def have_jax(): """Check if jax and jaxlib are installed with the correct versions""" return ( diff --git a/tests/unit/test_util.py b/tests/unit/test_util.py index 53b3ea8117..8c9c153f3e 100644 --- a/tests/unit/test_util.py +++ b/tests/unit/test_util.py @@ -92,17 +92,6 @@ def test_git_commit_info(self): self.assertIsInstance(git_commit_info, str) self.assertEqual(git_commit_info[:2], "v2") - @unittest.skipIf(not pybamm.have_julia(), "Julia not installed") - def test_have_julia(self): - # Remove julia from the path - with unittest.mock.patch.dict( - "os.environ", {"PATH": os.path.dirname(sys.executable)} - ): - self.assertFalse(pybamm.have_julia()) - - # Add it back - self.assertTrue(pybamm.have_julia()) - class TestSearch(unittest.TestCase): def test_url_gets_to_stdout(self): From 960d00a919c8b5ffd187f1a1a637df1c8bd4811c Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Sun, 13 Nov 2022 12:32:37 -0500 Subject: [PATCH 04/10] eliminate julia import --- pybamm/util.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pybamm/util.py b/pybamm/util.py index d9cc806c06..454ebe6f65 100644 --- a/pybamm/util.py +++ b/pybamm/util.py @@ -15,7 +15,6 @@ import timeit from platform import system import difflib -from julia.api import Julia, JuliaInfo, JuliaError import numpy as np import pkg_resources From d532df3ca6e8e07e4bcbb5fbc3df9b754c0f8487 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Sun, 13 Nov 2022 14:05:00 -0500 Subject: [PATCH 05/10] coverage --- pybamm/expression_tree/functions.py | 35 ------------------- .../test_expression_tree/test_functions.py | 4 --- 2 files changed, 39 deletions(-) diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py index f86a72e040..7fcc4bd7af 100644 --- a/pybamm/expression_tree/functions.py +++ b/pybamm/expression_tree/functions.py @@ -201,13 +201,6 @@ def _sympy_operator(self, child): """Apply appropriate SymPy operators.""" return child - @property - def julia_name(self): - "Return the name of the equivalent Julia function, for generating Julia code" - raise NotImplementedError( - "No julia name defined for function {}".format(self.function) - ) - def to_equation(self): """Convert the node and its subtree into a SymPy equation.""" if self.print_name is not None: @@ -256,14 +249,6 @@ def _function_new_copy(self, children): """See :meth:`pybamm.Function._function_new_copy()`""" return pybamm.simplify_if_constant(self.__class__(*children)) - @property - def julia_name(self): - """See :meth:`pybamm.Function.julia_name`""" - # By default, the julia name for a specific function is the class name - # in lowercase - # Some functions may overwrite this - return self.__class__.__name__.lower() - def _sympy_operator(self, child): """Apply appropriate SymPy operators.""" class_name = self.__class__.__name__.lower() @@ -281,11 +266,6 @@ def _function_diff(self, children, idx): """See :meth:`pybamm.Symbol._function_diff()`.""" return 1 / sqrt(children[0] ** 2 + 1) - @property - def julia_name(self): - """See :meth:`pybamm.Function.julia_name`""" - return "asinh" - def _sympy_operator(self, child): """Override :meth:`pybamm.Function._sympy_operator`""" return sympy.asinh(child) @@ -306,11 +286,6 @@ def _function_diff(self, children, idx): """See :meth:`pybamm.Function._function_diff()`.""" return 1 / (children[0] ** 2 + 1) - @property - def julia_name(self): - """See :meth:`pybamm.Function.julia_name`""" - return "atan" - def _sympy_operator(self, child): """Override :meth:`pybamm.Function._sympy_operator`""" return sympy.atan(child) @@ -426,11 +401,6 @@ class Max(SpecificFunction): def __init__(self, child): super().__init__(np.max, child) - @property - def julia_name(self): - """See :meth:`pybamm.Function.julia_name`""" - return "maximum" - def _evaluate_for_shape(self): """See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`""" # Max will always return a scalar @@ -451,11 +421,6 @@ class Min(SpecificFunction): def __init__(self, child): super().__init__(np.min, child) - @property - def julia_name(self): - """See :meth:`pybamm.Function.julia_name`""" - return "minimum" - def _evaluate_for_shape(self): """See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`""" # Min will always return a scalar diff --git a/tests/unit/test_expression_tree/test_functions.py b/tests/unit/test_expression_tree/test_functions.py index 67d5411d8e..fbc3625654 100644 --- a/tests/unit/test_expression_tree/test_functions.py +++ b/tests/unit/test_expression_tree/test_functions.py @@ -114,10 +114,6 @@ def test_exceptions(self): with self.assertRaises(pybamm.DomainError): pybamm.Function(test_multi_var_function, a, b) - fun = pybamm.Function(np.cos, pybamm.t) - with self.assertRaisesRegex(NotImplementedError, "No julia name"): - fun.julia_name - def test_function_unnamed(self): fun = pybamm.Function(np.cos, pybamm.t) self.assertEqual(fun.name, "function (cos)") From 0611f204275d0919c011e229f1efd484c8953679 Mon Sep 17 00:00:00 2001 From: Alec Bills <48105066+abillscmu@users.noreply.github.com> Date: Sun, 13 Nov 2022 16:53:24 -0500 Subject: [PATCH 06/10] Update CHANGELOG.md --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index c2ee1d24db..692b1714f8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,8 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) +# [v22.10.post1](https://github.com/pybamm-team/PyBaMM/tree/v22.10.post1) - 2022-10-31 +- Removed all julia generation code ([#2453](https://github.com/pybamm-team/PyBaMM/pull/2453)). Julia code will be hosted at [PyBaMM.jl]((https://github.com/tinosulzer/PyBaMM.jl) from now on. + # [v22.10](https://github.com/pybamm-team/PyBaMM/tree/v22.10) - 2022-10-31 ## Features From 641a4587cfdad0d6901bf9c82d8ec9676f782bfe Mon Sep 17 00:00:00 2001 From: Alec Bills <48105066+abillscmu@users.noreply.github.com> Date: Sun, 13 Nov 2022 16:53:55 -0500 Subject: [PATCH 07/10] Update CHANGELOG.md --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 692b1714f8..0425aff0a9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,7 +1,7 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) # [v22.10.post1](https://github.com/pybamm-team/PyBaMM/tree/v22.10.post1) - 2022-10-31 -- Removed all julia generation code ([#2453](https://github.com/pybamm-team/PyBaMM/pull/2453)). Julia code will be hosted at [PyBaMM.jl]((https://github.com/tinosulzer/PyBaMM.jl) from now on. +- Removed all julia generation code ([#2453](https://github.com/pybamm-team/PyBaMM/pull/2453)). Julia code will be hosted at [PyBaMM.jl](https://github.com/tinosulzer/PyBaMM.jl) from now on. # [v22.10](https://github.com/pybamm-team/PyBaMM/tree/v22.10) - 2022-10-31 From d9d45218a10598bd0ef35d93429195480ed772e2 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 13 Nov 2022 21:54:42 +0000 Subject: [PATCH 08/10] style: pre-commit fixes --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0425aff0a9..662cfbdd10 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,7 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) # [v22.10.post1](https://github.com/pybamm-team/PyBaMM/tree/v22.10.post1) - 2022-10-31 + - Removed all julia generation code ([#2453](https://github.com/pybamm-team/PyBaMM/pull/2453)). Julia code will be hosted at [PyBaMM.jl](https://github.com/tinosulzer/PyBaMM.jl) from now on. # [v22.10](https://github.com/pybamm-team/PyBaMM/tree/v22.10) - 2022-10-31 From f4b9cf195a8c1add60de772ac43634667e8507eb Mon Sep 17 00:00:00 2001 From: Alec Bills <48105066+abillscmu@users.noreply.github.com> Date: Sun, 13 Nov 2022 16:57:30 -0500 Subject: [PATCH 09/10] Update version.py --- pybamm/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pybamm/version.py b/pybamm/version.py index cbc21d2539..943ea9fb60 100644 --- a/pybamm/version.py +++ b/pybamm/version.py @@ -1 +1 @@ -__version__ = "22.10" +__version__ = "22.10.post1" From cbef36873bf88a515ef1863a5d31b1c28612e9f1 Mon Sep 17 00:00:00 2001 From: Ferran Brosa Planella Date: Fri, 2 Dec 2022 08:41:39 +0000 Subject: [PATCH 10/10] fixed version.py --- pybamm/version.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pybamm/version.py b/pybamm/version.py index bd31514bfd..969639bd2c 100644 --- a/pybamm/version.py +++ b/pybamm/version.py @@ -1,5 +1 @@ -<<<<<<< HEAD __version__ = "22.11" -======= -__version__ = "22.10.post1" ->>>>>>> main