Need to flesh out docstrings for the solve functions (the one-off wrappers) and interpolation matrix.

Priors should be available through a shared library with the UQ code, since the framework is identical.

Ensure that the Dockerfile is correctly setup to use the released version on master.

The library should be able to handle time-dependent linear problems as a Linear FEM is solved at each time step in most treatments of time-dependent problems. At the moment this would have to be done manually, but perhaps some convenience functions could handle this for the user.

Related question, which I do not know the best answer to at the moment: How to do estimation for a time-dependent problem?

Convert talk into RST files for inclusion in documentation. I was unable to find a way to have either a notebook or a more version-controllable version of a notebook hold this and then convert to RST to include the notebook in the documentation (or at least not within the time I had available for putting the talk in the docs), so this will have to be the current solution.

Code needs some working examples.

Should there be a base class for the dressed up sparse matrix classes?

We never care about the forcing covariance matrix in isolation, only its action on a vector. Thus, we do not necessarily need to explicitly form it and would only need to implement some version of matrix-vector multiplication. (Looking at the code for the ForcingCovariance class should make this apparent that this is the only operation defined for the class.)

Allow the model to accept repeated measurements in some fashion and compute the appropriate posterior.

The Hessian of the log likelihood can be easily computed in closed form, add a method for doing this.

Revisit if memory deallocation and garbage collection is done carefully.

Add a changelog file to keep track of the evolution of the software from the start.

Building the documentation requires installing stat-fem to get the correct version number. However, the installation requires Firedrake, which has a complex build process. Work on fixing this. A few thoughts:

1. Could use the docker container the way that the tests currently do. Would be tricky to get the docs out of it, though, and the build scripts would need to be manually configured.

2. readthedocs builds in a virtual environment, but the Firedrake installer needs to create its own virtual environment. I imagine there is a workaround.

3. Or, write a manual install script to be run inside the virtualenv? Seems like the wrong solution.

The right answer is probably 2, but I need to test this out.

At the moment, the interpolation matrix does a projection to map the sensor data point to the FEM mesh points. This is correct for piecewise linear elements (1st order continuous Galerkin or Lagrange finite elements), but is not correct for the general case. To be truly correct, this should be an interpolation rather than a projection, but interpolation requires knowledge of the FEM basis functions and depends on more complex workings under the hood to calculate this correctly. Projection is easier to compute in the general case, and thus for the initial work on this I simply implemented projection to get the library working.

However, given that other approximations are made throughout the computation, it is unclear if this makes a significant difference and needs to be studied in more depth to determine how much error this introduces to the computation.

setuptools no longer allows specifyling license as a list of strings. Need to fix this.

What is the correct way to handle Dirichlet Boundary Conditions? BCs are disabled when solving the linear system as that was what was done in the other papers, but I don't completely understand why. Worth tracking this down and clearly explaining what is going on in the documentation.

Recent updates to interpolation in Firedrake should simplify the computation of the interpolation matrix and generalize it to other finite element spaces. This needs to be implemented and tested.

The Makefile is a bit opaque, as all commands are routed to Sphinx and not all of them show up in the help command (for instance clean). This is worth looking into.

Improve the docstrings and set up sphinx documentation.

At the moment, all processes have to collectively call the expensive _solve_forcing_covariance routine each time, even though some solves only need to be done on the root process of the ensemble. This presumably is best solved with some re-configuring of the MPI ensemble communicators to create a better master/worker division, where master does all "base" FEM solves, while the "sensor" solves are divided up among the remaining worker processors. Performance hit at the moment is probably minimal, though, as the number of sensor solves is presumably much larger in most cases.

Write routines for predicting values at new data points.

Firedrake gives more control over the PETSc options that are available, this should exend to the solves done by the LinearSolver class (can just instantiate a Firedrake version, which is reusable, thus we can just pass on all solves to that object).

Firedrake has updated their assemble options. To avoid the need to update these due to further upstream changes in the future, this capability should be deprecated.

Write some more details on ensemble parallelism (in particular give a full example).

Estimation routines currently do not cache things that are expensive to compute and are re-used. In particular, the forcing covariance matrix only needs to be computed once (it is currently computed twice at each step in the minimization routine), and the Cholesky factorization of the combined covariance matrix is computed twice (once in the log-likelihood, once in the derivatives). Caching these computations will dramatically improve the speed of estimation routines.

The paper that this software is based on has now published. The documentation (readme and docs) need to update the references based on this to link to the published paper.

I still have issues with the parallel tests hanging. It seems to occur reliably with 4 processes when running a sufficient fraction of the test suite and hangs without a clear reason why. (2 processes seems reliably fine.) When I just run one file at a time it appears to be fine, but once I add enough other tests it seems to cause issues. Changing the order doesn't seem to have a predictable pattern associated with it, but it is usually the LinearSolver tests or the ForcingCovariance tests that are most prone to this problem.

I don't think it is a memory usage issue -- the base test suite consistently uses ~100 MB of memory, and the parallel test suite uses roughly 4 times that as expected. This amount is nowhere near what would be needed to cause any issues.

For whatever reason, it never seems to have trouble on the Travis CI, but on my Mac or in my own version of the Docker container it is more problematic.

However, I have no idea at this point what else could be causing this other than some MPI issue that I have failed to uncover.

I can't seem to figure out why the coverage reports aren't uploading. A couple things to explore:

• Can I just install firedrake manually within the travis build? Coverage reports seem more straightforward if not in a docker container.
• I managed to get things working before on the travis branch, then they stopped. Look back at that commit and what changed after that to get some ideas.
• Keep tweaking? I don't know why I can never get the bash uploader to work...

Computing the forcing covariance is expensive as it is a large matrix (that is potentially dense, depending on the covariance parameters that are used). It cannot be stored for FEM problems of interest, so we use a sparse approximation to it in the computations. However, we do not know in advance what the structure of this sparse approximation of the matrix looks like, so we still have to compute all matrix elements in order to form the sparse version. A further complication is that memory allocation in PETSc is expensive, so the fact that we cannot know in advance how big this matrix is makes this a large computational expense.

At the moment, we deal with these issues by pre-computing all elements and only storing those that are needed in a Python dictionary. From the pre-computed version we can allocate memory and form the PETSc matrix. This has two issues: (1) this requires doing simple Python loops, which can be slow, and (2) the Python PETSc interface does not give as much flexibility when allocating memory for sparse matrices as the C interface. This means that significant improvements in speed and memory usage can be obtained by writing the matrix formation routines as a C extension. I believe that pre-computing the elements is still the best approach, which would mean that something like a linked list (as a replacement for the Python dictionary) is probably the best way to store the matrix elements prior to allocating memory and forming the PETSc matrix.

It appears that the new function signature for assemble has changed. Need to update the function wrapper based on this.

Finish implementation of estimation routines. Ideally, mimic behavior in the UQ platform where you can do either MLE or MCMC, and the MLE is used to initiate the MCMC and pick a good step size by inverting the Hessian. This will involve some (small) additional overhead to compute the necessary derivatives.

Docutils introduced a change that was not backwards compatible with Sphinx in a recent version:

https://blog.readthedocs.com/build-errors-docutils-0-18/

Fix this by bumping minimum version of sphinx in docs/requirements.txt

If multiple observations are made over time, provide a way to update the posterior by providing additional measurements to a linear solver object. Related to #11, and should be easier to implement as it can (I think) use existing machinery to do the update.

It is not clear if means computed in the data space are scaled by the discrepancy (I believe they are not, but I had to think about it which means that is is not clear). Getting the scaled mean is likely to be a common desired output, so the solver class should handle this for the user.

The log posterior derivative is not computed correctly for all cases (for non-unity values of rho, in particular). Missing a factor of rho (due to the fact that rho is on a logarithmic scale).

Parameter values are opaque, on a log scale, etc. We have a similar problem in the UQ library, so a shared solution in a 3rd party library may be the right answer.

More details should be provided to help users specify the PETSc options. At the very least, link to the Firedrake documentation on this.