Issue 665 jac class #670

rtimms · 2019-10-18T14:26:31Z

Description

Adds a class Jacobian to re-use know jacobians of any expressions. Also changes the behavior of some binary simplifies which were returning dense matrices of zeros in some cases.

I was going to do Jacobian by equation as part of this (#651) but (I think) that needs to be done in discretisation so you can use _concatenate_in_order. I've added a method disc.create_jacobian which lets you do compute the jacobian, but by default I've left calculation of jac until the solve. The main reason being that it is much faster to compute the jacobian on the simplified model, but you dont want to simplify as part of discretisation as this prevents you from updating parameters. I've kept the create_jacobian method as there may be use cases where you want to compute the jacobian (once) of the unsimplified model, and then go ahead and solve multiple times but with different params etc.

Fixes #665

Type of change

Please add a line in the relevant section of CHANGELOG.md to document the change (include PR #) - note reverse order of PR #s.

New feature (non-breaking change which adds functionality)
Optimization (back-end change that speeds up the code)
Bug fix (non-breaking change which fixes an issue)

Key checklist:

No style issues: $ flake8
All tests pass: $ python run-tests.py --unit
The documentation builds: $ cd docs and then $ make clean; make html

You can run all three at once, using $ python run-tests.py --quick.

Further checks:

Code is commented, particularly in hard-to-understand areas
Tests added that prove fix is effective or that feature works

codecov · 2019-10-18T15:18:40Z

Codecov Report

Merging #670 into master will decrease coverage by 0.02%.
The diff coverage is 96.4%.

@@            Coverage Diff             @@
##           master     #670      +/-   ##
==========================================
- Coverage   98.59%   98.56%   -0.03%     
==========================================
  Files         169      170       +1     
  Lines        8450     8517      +67     
==========================================
+ Hits         8331     8395      +64     
- Misses        119      122       +3

Impacted Files	Coverage Δ
pybamm/solvers/scikits_dae_solver.py	`97.67% <ø> (ø)`	⬆️
pybamm/expression_tree/vector.py	`100% <ø> (ø)`	⬆️
pybamm/expression_tree/simplify.py	`94.98% <ø> (ø)`	⬆️
pybamm/expression_tree/scalar.py	`100% <100%> (ø)`	⬆️
pybamm/solvers/dae_solver.py	`99.23% <100%> (+0.01%)`	⬆️
pybamm/expression_tree/jacobian.py	`100% <100%> (ø)`
pybamm/expression_tree/independent_variable.py	`100% <100%> (ø)`	⬆️
pybamm/expression_tree/unary_operators.py	`96.18% <100%> (-0.02%)`	⬇️
pybamm/expression_tree/array.py	`100% <100%> (ø)`	⬆️
pybamm/expression_tree/functions.py	`100% <100%> (ø)`	⬆️
... and 9 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b21b738...ddc30ed. Read the comment docs.

valentinsulzer

Cheers Rob, looks like you found a few bugs that will make this more efficient. There are a few things to address, especially with the Jacobian class, see below

In order to reuse already-calculated jacobians, this function needs to keep calling its own jac method. For example, for a binary operator, call Jacobian.jac on left and right, and then call symbol._jac on processed_left and processed_right:

processed_left = self.jac(symbol.left)
processed_right = self.jac(symbol.right)
return symbol._jac(processed_left, processed_right)

Otherwise a new Jacobian object just gets created each time. See the Simplify class for example.
Apologies if you're doing this in a different way that I've missed.

valentinsulzer · 2019-10-18T23:04:09Z

pybamm/discretisations/discretisation.py

+        """Creates Jacobian of the discretised model.
+        Note that the model is assumed to be of the form M*y_dot = f(t,y), where
+        M is the (possibly singular) mass matrix. The Jacobian is df/dy.
+


Add a comment here referring to this PR to explain why this function isn't used. I could see us forgetting and wondering why we don't use this in future

also, if we make some kind of master simulation class we could use this method from the disretisation in the solver set up since we would be able to access the disc of the model when in the solver set up

valentinsulzer · 2019-10-18T23:06:39Z

pybamm/discretisations/discretisation.py

@@ -696,8 +754,15 @@ def _concatenate_in_order(self, var_eqn_dict, check_complete=False):
        ----------
        var_eqn_dict : dict
            Equations ({variable: equation} dict) to dicretise
+        check_complete : bool, optional
+            Whether to check keys in var_eqn_dict against self.y_slices. Defualt


got some defualt sneaking in again 👀

I have muscle memory for this typo...

haha I assumed it had got into your autocorrect

valentinsulzer · 2019-10-18T23:09:34Z

pybamm/expression_tree/binary_operators.py

+            if shape == ():
                return pybamm.Scalar(0)
+            else:
+                if len(shape) == 1 or shape[1] == 1:
+                    return pybamm.Vector(np.zeros(shape))
+                else:
+                    return pybamm.Matrix(csr_matrix(shape))


seeing them all together in review mode, these lines appear several times together and could be abstracted into a functon pybamm.zeros_of_shape(shape)

good suggestion

valentinsulzer · 2019-10-18T23:14:12Z

pybamm/expression_tree/jacobian.py

+        # jacobian removes the domain(s)
+        jac.domain = []
+        jac.auxiliary_domains = {}
+        return jac


In order to reuse already-calculated jacobians, this function needs to keep calling its own jac method. For example, for a binary operator, call Jacobian.jac on left and right, and then call symbol._jac on processed_left and processed_right:

processed_left = self.jac(symbol.left) processed_right = self.jac(symbol.right) return symbol._jac(processed_left, processed_right)

Otherwise a new Jacobian object just gets created each time. See the Simplify class for example.
Apologies if you're doing this in a different way that I've missed.

ah, thanks -- I had misunderstood how known_jacs was getting passed around.

valentinsulzer · 2019-10-18T23:15:48Z

pybamm/expression_tree/symbol.py

-            jac.domain = []
-            jac.auxiliary_domains = {}
-            return jac
+        return pybamm.Jacobian(known_jacs).jac(self, variable)


have just seen this known_jacs thing. I still think it's more efficient to only create the Jacobian class once

valentinsulzer · 2019-10-18T23:17:01Z

pybamm/models/base_model.py

@@ -52,7 +52,12 @@ class BaseModel(object):
        automatically
    jacobian : :class:`pybamm.Concatenation`
        Contains the Jacobian for the model. If model.use_jacobian is True, the
-        Jacobian is computed automatically during the set up in solve
+        Jacobian is computed automatically during discrisation


Suggested change

Jacobian is computed automatically during discrisation

Jacobian is computed automatically during discretisation

valentinsulzer · 2019-10-18T23:17:21Z

pybamm/models/base_model.py

+        Contains the Jacobian for the algebraic part of the model. This may be used
+        by the solver when calculating consistent initial conditions. If
+        model.use_jacobian is True, the Jacobian is computed automatically during
+        discrisation


Suggested change

discrisation

discretisation

valentinsulzer · 2019-10-18T23:18:05Z

pybamm/models/base_model.py

@@ -226,6 +232,14 @@ def jacobian(self):
    def jacobian(self, jacobian):
        self._jacobian = jacobian

+    @property
+    def jacobian_algebraic(self):


why is there a jacobian_algebraic but no jacobian_rhs?

only because we only ever use either the full jacobian (for the solver) or the algebraic part (for consistent ics). I'll add jacobian_rhs too for completeness.

oh right, I had forgotten about that, probably fine without jacobian_rhs then. Need to make sure we aren't calculating the jacobian for the algebraic part twice then (can't remember if already doing this)

valentinsulzer

Looks great now @rtimms !

valentinsulzer · 2019-10-24T12:50:47Z

pybamm/expression_tree/jacobian.py

+            left_jac = self.jac(left, variable)
+            right_jac = self.jac(right, variable)
+            # _outer_jac defined in pybamm.Outer
+            jac = symbol._outer_jac(left_jac, right_jac, variable)


not a big deal but you could have the isinstance(Outer) check on line 60 to avoid writing the process_children lines twice

valentinsulzer · 2019-10-24T12:51:37Z

tests/unit/test_expression_tree/test_jac.py

@@ -80,9 +80,9 @@ def test_linear(self):

        # test jac of outer if left evaluates to number
        func = pybamm.Outer(pybamm.Scalar(1), pybamm.Scalar(4))
-        jacobian = np.zeros((1, 4))
+        # jacobian = np.zeros((1, 4))


can be removed

rtimms added 6 commits October 16, 2019 14:24

#665 add jacobian class

d1acc43

#665 make binary simplifies sparse

8834f39

#665 move jac to discretise

c6d5717

#665 change default bahviour of jac

fbb18a4

#665 fix update params

d9b3cf0

#665 compute jac of simplified model before solve

6837960

valentinsulzer requested changes Oct 18, 2019

View reviewed changes

rtimms added 6 commits October 21, 2019 12:29

#665 tino comments

593370d

#665 merge master

719bf55

#665 jac of independent var

2559aef

#665 add tests

e4bfe1a

#665 debug jac of outer

9d1579c

#665 codacy

93d9862

valentinsulzer approved these changes Oct 24, 2019

View reviewed changes

rtimms added 2 commits October 24, 2019 14:48

#665 tino comments

d1b0b2d

Merge branch 'master' into issue-665-jac-class

ddc30ed

rtimms merged commit 46b47a1 into master Oct 24, 2019

rtimms deleted the issue-665-jac-class branch October 24, 2019 15:59

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Issue 665 jac class #670

Issue 665 jac class #670

rtimms commented Oct 18, 2019

codecov bot commented Oct 18, 2019 •

edited

valentinsulzer left a comment

valentinsulzer Oct 18, 2019

rtimms Oct 21, 2019

valentinsulzer Oct 18, 2019

rtimms Oct 21, 2019

valentinsulzer Oct 21, 2019

valentinsulzer Oct 18, 2019

rtimms Oct 21, 2019

valentinsulzer Oct 18, 2019

rtimms Oct 21, 2019

valentinsulzer Oct 18, 2019

valentinsulzer Oct 18, 2019

valentinsulzer Oct 18, 2019

valentinsulzer Oct 18, 2019

rtimms Oct 21, 2019

valentinsulzer Oct 21, 2019

valentinsulzer left a comment

valentinsulzer Oct 24, 2019

valentinsulzer Oct 24, 2019

	Jacobian is computed automatically during discrisation
	Jacobian is computed automatically during discretisation

Issue 665 jac class #670

Issue 665 jac class #670

Conversation

rtimms commented Oct 18, 2019

Description

Type of change

Key checklist:

Further checks:

codecov bot commented Oct 18, 2019 • edited

Codecov Report

valentinsulzer left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

valentinsulzer left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented Oct 18, 2019 •

edited