Matrix compatibility and interdependent components

[tBuLi]

This PR makes two big new changes: support for sympy.MatrixExpr as variables and as an indirect consequence it also adds component interdependence.

Consider e.g.

a, b = parameters('a, b')
x, y, z = variables('x, y, z')
model_dict = {
    y: a**3 * x + b**2,
    z: y**2 + a * b
}

This could be a great asset when y is expensive to compute and used multiple times. In such a model, y is no longer a dependent variable, but rather an interdependent one.

The big challenge here was how to calculate the Jacobian and Hessian for such a model, since also there we would like D(z, a) to be a function of D(y, a):

{D(y, a): 3 * a**2 * x,
 D(y, b): 2 * b,
 D(z, a): b + 2 * y * D(y, a),
 D(z, b): a + 2 * y * D(y, b),
}

This PR so far addresses:

• [x] Introduce separate functions for building a model to represent the jacobian, hessian and least-squares.
• [x] Models now have an attribute self.jacobian_model and self.hessian_model, which are models to represent those quantities. self.eval_[...] calls the respective model and properly arranges the output. (The models just have a set of component representing the jacobian and hessian, so they need to be ordered properly.)
• [x] In addition to the previous point, evaluating models is now based purely on calling the eval_components, eval_jacobian, and eval_hessian. Unless changed by a subclass, eval_components calls numerical_components. The other two simply call self.jacobian_model or self.hessian_model.

• [x] All models now have an important new property: connectivity_mapping. This indicates for every component on which variables and parameters it depends. This is used in code generation when lambdifying symbolic expressions, such that every component is no longer assumed to depend on all independent variables and parameters, but is just provided with those it needs.

  For BaseNumericalModel subclasses such a connectivity_mapping should be provided, although the old syntax is also still supported.

• [x] Based on connectivity_mapping, a topological sort can be performed to find the right evaluation order of the components. This is stored in ordered_symbols. This added a dependency to the very lightweight package toposort.
• [x] In order to properly perform the chain rule we need to turn all variables into functions before calculating derivatives. To this end, there is now 'vars_as_functions' to provide a simple mapping between the two, and function_dict, which is equivalent to model_dict except that all variables have been turned into Function objects and the sorting is based on ordered_symbols, not alphabetical.
• [x] Removed a lot of commented out fossilized code which stemmed from a different era.
• [x] Constraints have been massively simplified, getting ever closer to the goal of eliminating them completely.
• [ ] Fixed #84, printing vector models in a more recognizable way. We could still consider printing the function_dict instead though.
• [ ] Move models to a separate file.
• [x] Rescale chi2 by 1/2 so we can always invert the Hessian immediately to get the covariance matrix, no worries about factors of 2

[pckroon]

Referencing #39

[tBuLi]

Don't worry, it's still very much a work in progress ;). The non-trivial comments I left open were things I wanted to think about carefully before implementing them.

For future reference: this mentiones how the hessian relates to the covariance matrix: https://www8.cs.umu.se/kurser/5DA001/HT07/lectures/lsq-handouts.pdf

My suggestion is to remove finite differences entirely, and to always use the Hessian of the objective to calculate covariances. This removes the discussion about using finite differences for optimization and will finally clean-up the mess that is covariance matrix calculation at this point.

[pckroon]

My suggestion is to remove finite differences entirely, and to always use the Hessian of the objective to calculate covariances. This removes the discussion about using finite differences for optimization and will finally clean-up the mess that is covariance matrix calculation at this point.

Bad plan, we still need it to get a covariance matrix for e.g. ODEModels and NumericalModels. What we can do, however, is deprecate the finite difference Jacobian, and change it for a finite difference Hessian method. Or keep both next to each other.